Sensores de curvatura : optimización de su rendimiento

MarcosSotoBúa Tesisdoctoral SENSORESDECURVATURA: OPTIMIZACIÓNDESURENDIMIENTO DepartamentodeFísica Aplicada FacultaddeFísica Santiago de Compostela, 27 de noviembre de 2006 Eva Acosta Plaza, profesora titular de la Universidade de Santiago de Compostela, EXPONE: Que la presente memoria, titulada “Sensores de curvatura: optimizaci´on de su rendimiento”, ha sido realizada bajo mi direcci´on por el doctorando Marcos Soto B´ua en el ´ Area de ´ Optica del Departamento de F´ısica Aplicada de la Universidade de Santiago de Compostela y constituye la tesis doctoral que presenta para optar al t´ıtulo de Doctor en F´ısica por esta universidad. Fdo.: La directora, Fdo.: El doctorando, Eva Acosta Plaza Marcos Soto B´ua A todos pertenece lo que piensas. Solo es tuyo lo que sientes. Si quieres que sea tuyo lo que piensas, has de sentirlo. Friedrich von Schiller (1759-1805) A Josito A mis padres Agradecimientos Esta memoria es el resultado de un largo trabajo de investigaci´on realizado en el ´ Area de ´ Optica del Departamento de F´ısica Aplicada de la Universidade de Santiago de Compostela (USC). Quiero aqu´ı expresar mi agradecimiento a Eva Acosta Plaza, directora de la tesis doctoral, por su inestimable ayuda para la materializaci´on de este proyecto, y a Susana R´ıos Rodr´ıguez, por su implicaci´on y ´animo personal. A los miembros del despacho colectivo de la planta baja y sat´elites - al fin y al cabo, todo el mundo pasa por all´ı peri´odicamente -, por el inter´es que siempre han mostrado por mi investigaci´on. Pero, sobre todo, por estrechar v´ınculos con mi vida personal. Muy especialmente, a mi querido Jos´e Luis le agradezco el tremendo apoyo, en todos los ´ordenes, que me ha brindado para terminar este trabajo de investigaci´on. ´ El sabe que puede considerarlo como el suyo propio y a ´el se lo dedico. A mi familia al completo, su inefable ayuda y su confianza en m´ı a lo largo de todos estos a˜nos. A mis otros amigos, que siempre est´an ah´ı cuando los necesito para apoyarme y darme todo su cari˜no, les agradezco el compartir conmigo tantos buenos momentos en su tiempo libre, escaso y valioso. A los voluntarios y compa˜neros de Jos´e Luis en Interm´on-Oxfam (IO), su inter´es sincero en el progreso de la redacci´on de esta memoria. Para realizar este trabajo de investigaci´on he disfrutado de una beca predoctoral concedida por la Secretaria Xeral de Investigaci´on e Desenvolvemento de la Xunta de Galicia - hoy Direcci´on Xeral de Investigaci´on, Desenvolvemento e Innovaci´on (I+D+i) de la Conseller´ıa de Innovaci´on e Industria - y de otra perteneciente al Programa de Formaci´on de Personal Investigador (FPI) del Ministerio de Ciencia y Tecnolog´ıa - hoy Ministerio de Educaci´on y Ciencia -, organismos a los que les agradezco el apoyo econ´omico prestado. Por ´ultimo, indicar que el trabajo redactado en esta memoria se enmarca dentro de varios proyectos de investigaci´on financiados por la Xunta de Galicia (Ref. PGIDT00PXI20601PR), y por el Ministerio de Educaci´on y Ciencia (Ref. AYA20001565-C02-02 y AYA2004-07773-C02-02). 2Introducci´on curvatura en las configuraciones habituales del mismo. Segundo, c´omo mejorar la respuesta del sensor para que la fase recuperada tenga m´as resoluci´on. Tercero, debido a la mayor dificultad para realizar medidas precisas en la frontera, determinar el peso de la informaci´on recuperada en el borde de la pupila respecto del resto de informaci´on obtenida a partir de medidas en el interior. La memoria se ha estructurado en cinco cap´ıtulos. En el primero de ellos, se hace una revisi´on hist´orica de la ETI como t´ecnica de recuperaci´on de fase en diversas aplicaciones, se revisan los problemas principales que surgen al ponerla en pr´actica, las limitaciones de las soluciones utilizadas hasta la fecha y se proponen nuevas soluciones que constituyen los objetivos de este trabajo de investigaci´on. En el cap´ıtulo segundo, se trata la estimaci´on de la derivada axial de la irradiancia en el interior de la pupila suponiendo irradiancia uniforme en el plano de recuperaci´on, es decir, se estudia la estimaci´on del laplaciano de la fase. En concreto, se identifican las principales fuentes de error de las medidas de la curvatura del frente de onda y se modela su varianza en funci´on de estos factores y de la posici´on axial de los planos de irradiancia. Por ´ultimo, se procede al dise˜no de un m´etodo para optimizar el rendimiento del sensor en relaci´on con estas medidas, que se comprueba en un ejemplo ilustrativo. En el cap´ıtulo tercero, se introducen los m´etodos m´as habituales para obtener informaci´on de la fase en el borde de la pupila, se analizan sus principales limitaciones y se calcula una expresi´on para el error de la derivada normal de la fase en funci´on de estos factores limitantes y de la posici´on axial de los planos de irradiancia. Con esta informaci´on se procede al dise˜no de un procedimiento de optimizaci´on de este tipo de medidas, que, al igual que en el cap´ıtulo anterior, se ilustra con un ejemplo. Una vez identificadas las limitaciones de las implementaciones tradicionales de la ETI, en el cap´ıtulo cuarto se proponen otras configuraciones del sensor persiguiendo el objetivo de mejorar su rendimiento fijadas las condiciones de medida. Por ´ultimo, en el cap´ıtulo quinto, se aborda el problema de contorno que constituye la ETI aplicada a un haz de luz que ha atravesado una pupila, con la perpectiva de analizar la importancia en el resultado final de la informaci´on recuperada en el borde de la misma. La complejidad manifiesta de las medidas en el borde de la pupila justifica un estudio donde se pondere su peso en la recuperaci´on de la fase en comparaci´on con las medidas en el interior. Para concluir, se incluyen tres ap´endices no solo con herramientas de apoyo para la comprensi´on de los cap´ıtulos, sino tambi´en con informaci´on complementaria m´as orientada hacia aplicaciones concretas de la ETI, que no aparece en el cuerpo de la Introducci´on 3 memoria para dinamizar su lectura, pero no por ello menos importante. En el primero se deduce la ETI y se halla la versi´on de la misma particularizada a haces limitados a una abertura en el plano de recuperaci´on. En el segundo se obtienen las expresiones necesarias para aplicar el procedimiento que se desarrolla en esta memoria para optimizar la respuesta del sensor de curvatura, en aquellas aplicaciones donde haya una limitaci´on temporal que repercuta en la detecci´on de un n´umero finito de fotones. Por ´ultimo, el tercero concierne a la configuraci´on m´as habitual de esta t´ecnica de recuperaci´on de fase en el ´ambito de la ´ Optica adaptativa aplicada a la Astronom´ıa: la configuraci´on de Roddier. Cap´ıtulo 1 La Ecuaci´on de Transporte de Intensidad como t´ecnica de recuperaci´on de fase En este cap´ıtulo inicial, se aborda el papel de la Ecuaci´on de Transporte de Intensidad (ETI) en el problema de la recuperaci´on de fase a partir de medidas de irradiancia, para dar contexto al trabajo de investigaci´on que recoge esta memoria. Tras un repaso, en la primera secci´on, de las principales aplicaciones de la ETI, en la segunda secci´on se introduce esta ecuaci´on como relaci´on ´util para resolver un problema gen´erico de recuperaci´on de fase, as´ı como las necesidades y dificultades m´as importantes que surgen al aplicarla y algunas de las soluciones que se han propuesto hasta la fecha para reducirlas. Bas´andonos en estos antecedentes, en la tercera secci´on, se explican las cuestiones que han motivado este trabajo de investigaci´on y las propuestas generales para dar soluci´on a tales cuestiones, que constituyen los pilares del mismo. 1.1. Introducci´on Cuando una onda electromagn´etica se propaga en un medio determinado, su amplitud y su fase sufren modificaciones en mayor o menor grado, acumulando informaci´on acerca del medio que atraviesa. La variaci´on de amplitud provocada por la absorci´on de radiaci´on durante la propagaci´on proporciona datos acerca de la estructura del medio atravesado. No obstante, en ocasiones el medio es tan poco absorbente que la onda se transmite sin que su amplitud se modifique visiblemente. 6La Ecuaci´on de Transporte de Intensidad como t´ecnica de recuperaci´on de fase Es entonces cuando la obtenci´on de la fase de la onda a la salida del objeto o medio dispersivo es esencial para reunir informaci´on acerca del mismo, a costa de tener que emplear algoritmos y t´ecnicas complejas. Los dispositivos que se encargan de medir observables con los que obtener informaci´on acerca de la fase de la onda electromagn´etica se llaman sensores de frente de onda [1, 2]. Este tipo de aparatos se utilizan en muchas aplicaciones, como control de calidad de componentes ´opticos [3–5] o de haces l´aser de alta energ´ıa, pero su mayor desarrollo lleg´o dentro del campo de la ´ Optica Astron´omica de la mano de la ´ Optica activa y adaptativa, cuyo concepto fue descubierto en paralelo por dos investigadores en la d´ecada de los cincuenta (H. Babcock y V. P. Linnik) [6]. Sin embargo, los intentos serios de aplicar esta tecnolog´ıa en Astronom´ıa se hicieron esperar hasta mediados de los a˜nos setenta [7–9] y, hasta los noventa, la mayor parte del esfuerzo para desarrollarla se centr´o en las aplicaciones militares [6, 10, 11]. En Astronom´ıa la ´ Optica activa y adaptativa se utiliza para mejorar la resoluci´on de los telescopios terrestres compensando activamente sus aberraciones, inducidas t´ermicamente o causadas por defectos de fabricaci´on, as´ı como adapt´andose a las aberraciones provocadas por las turbulencias atmosf´ericas, empleando, para ello, un espejo deformable [6, 12–17]. Estos factores tienen m´as relevancia conforme aumenta el tama˜no de los telescopios [18], una necesidad fundamental a la hora de estudiar objetos celestes cada vez de menor tama˜no y tambi´en menos brillantes [19, 20]. La fase de una onda electromagn´etica puede obtenerse utilizando variedad de m´etodos. Por ejemplo, la holograf´ıa y la interferometr´ıa permiten caracterizar con precisi´on la forma del frente de onda a partir de patrones de interferencia. Cuando la aplicaci´on exige mayor robustez, estabilidad y mayor flexibilidad en los requerimientos de coherencia de la fuente, se emplean t´ecnicas no interferom´etricas, entre las que destacan los sensores Shack-Hartmann [1]. Alternativamente, pueden emplearse otros m´etodos no interferom´etricos que nos permitan deducir la amplitud compleja de la onda electromagn´etica a partir de medidas de irradiancia en eje [21– 38]. Este conjunto de t´ecnicas se conocen habitualmente como el problema de la recuperaci´on de fase. En este campo, hay m´etodos pr´acticos que extraen la fase como resultado de un procedimiento de optimizaci´on no lineal iterativo, basado, o bien en el conocimiento de la distribuci´on de intensidad en un segundo plano de difracci´on, o bien en la aplicaci´on de ciertas ligaduras [21, 22, 27]. Pero tambi´en existen m´etodos deterministas, que se fundamentan en una ecuaci´on matem´atica que describe una relaci´on m´as o menos precisa entre la amplitud y la fase de la onda [23–25, 31, 34, 37]. La recuperaci´on es determinista en el sentido de que la fase 1.1 Introducci´on 7 se obtiene directamente en funci´on de los datos de irradiancia y la cuesti´on de la unicidad de la soluci´on ya no surge aqu´ı [39–41], a diferencia de los m´etodos que emplean algoritmos iterativos [22, 27, 42]. En este trabajo nos centraremos en un subconjunto concreto de estos ´ultimos m´etodos, los que se basan en la Ecuaci´on de Transporte de Intensidad (ETI) [24, 25]. La ETI es una ecuaci´on diferencial de continuidad v´alida para describir, en aproximaci´on paraxial, la propagaci´on de luz en medios homog´eneos [43], tanto con iluminaci´on coherente como incoherente*[46, 47]. Si no hay ceros en la distribuci´on de irradiancia, esta ecuaci´on expresa una relaci´on ´unica entre la variaci´on transversal del frente de onda y la variaci´on axial de la intensidad [48, 49]. En particular, si la luz se distribuye uniformemente a la salida del medio atravesado, las variaciones axiales instant´aneas de la distribuci´on de irradiancia est´an provocadas por las curvaturas locales del frente de onda. Este comportamiento sienta las bases del funcionamiento de los conocidos sensores de curvatura, inventados por Roddier a finales de los a˜nos ochenta para su aplicaci´on en ´ Optica activa y adaptativa en Astronom´ıa, por sugerencia de Beckers [12, 50, 51]. Multitud de trabajos ilustran el rendimiento de estos sensores en aplicaciones como el control de calidad de telescopios terrestres y la compensaci´on de las aberraciones provocadas por las turbulencias en im´agenes de objetos celestes [10, 16, 17, 53, 54, 56–64]. Otros usos del sensor de curvatura en el campo de la ´ Optica Astron´omica son la alineaci´on en fase de espejos segmentados [65, 66], la obtenci´on del patr´on de centelleo de la luz de las estrellas [67] o la recuperaci´on de im´agenes solares contrastadas, que demuestra el ´exito de esta t´ecnica cuando opera con fuentes extensas [68–70]. En el aspecto de las comunicaciones por sat´elite, la ETI se ha aplicado en la detecci´on de la inclinaci´on del frente de onda para guiar con precisi´on un haz l´aser [71]. Tambi´en se ha demostrado su utilidad para el control de calidad de espejos [3] o lentes progresivas [72]. En el ´ambito de la ´ Optica Fisiol´ogica, el sensor de curvatura se ha empezado a aplicar recientemente para la medida de la topograf´ıa corneal y de la l´agrima en el ojo humano [73, 74], y de las aberraciones oculares [75]. Pero, junto con la ´ Optica Astron´omica, la Microscop´ıa es, sin duda, el otro campo donde las t´ecnicas de visualizaci´on de fase fundamentadas en la ETI han tenido una aceptaci´on notable [46, 49, 76–105]. Ello ha sido posible gracias a su simplicidad, as´ı como a la sencilla interpretaci´on de los resultados, que las distinguen de las *Recientemente, Paganin et al. han calculado la versi´on de esta ecuaci´on aplicable cuando el medio refractivo es no lineal [44]. Tambi´en se ha demostrado ´util la versi´on temporal de la ETI para el c´alculo de las propiedades no lineales de materiales ´opticos [45]. 8La Ecuaci´on de Transporte de Intensidad como t´ecnica de recuperaci´on de fase t´ecnicas interferom´etricas [80, 85]. En Microscop´ıa la recuperaci´on de la informaci´on acumulada en la fase es esencial en muestras poco absorbentes y delicadas, especialmente en el caso de espec´ımenes biol´ogicos, pues la adici´on de sustancias para incrementar el contraste por absorci´on puede da˜narlos irremediablemente [93]. Inicialmente, los objetos transparentes se visualizaban mediante el simple acto de la propagaci´on en el espacio libre, desenfocando la imagen del microscopio. El incremento de la defocalizaci´on aumenta el contraste, pero a costa de disminuir la resoluci´on. Para incrementar la resoluci´on, se idearon t´ecnicas que miden espec´ıficamente el desplazamiento de fase del haz que atraviesa la muestra, como el Contraste Diferencial por Interferencia (DIC) o la Microscop´ıa de fase de Zernike [106]. Persiguiendo este mismo objetivo e inspir´andose en la t´ecnica de visualizaci´on por desenfoque original, la ETI comienza a aplicarse plenamente en Microscop´ıa a finales de los a˜nos noventa para resolver objetos d´ebilmente absorbentes. En el ´ambito de la Microscop´ıa ´optica, la ETI se ha aplicado para caracterizar materiales y dispositivos ´opticos, como microlentes planas [77] y fibras ´opticas [80, 82, 85], as´ı como para visualizar estructuras de muestras biol´ogicas [80, 101]. Con el objetivo de resolver las micro o nanoestructuras que componen todo tipo de objetos, la luz visible se sustituye por otro tipo de radiaciones, como en el caso de la Microscop´ıa electr´onica [103], la Microscop´ıa por rayos X [76, 78, 79, 81, 87, 94] y la Microscop´ıa por haces de ´atomos o de neutrones [83, 107, 108]. Cuando se emplean haces de part´ıculas de este tipo, si la radiaci´on tiene poca energ´ıa, el objeto ha de someterse a dosis masivas para obtener im´agenes contrastadas por absorci´on, que pueden conducir a posibles da˜nos estructurales [84], o si, por el contrario, la energ´ıa de la radiaci´on es muy alta, el contraste por absorci´on puede ser despreciable. Al igual que en Microscop´ıa ´optica, la alternativa pasa por obtener im´agenes de fase contrastadas [76, 78, 109]. Existen sofisticadas t´ecnicas interferenciales que permiten obtener datos cuantitativos precisos, pero, cuando se utilizan radiaciones de este tipo, la relajaci´on de los requisitos de coherencia de la fuente es esencial con el fin de rebajar costes y simplificar t´ecnicamente el dispositivo [49, 99, 110]. La aplicaci´on de m´etodos de recuperaci´on de fase no interferenciales como los fundamentados en la ETI ha flexibilizado enormemente los requisitos t´ecnicos del sistema para muestras homog´eneas, alcanz´andose altas cotas de calidad en la recuperaci´on [92, 95]. En Microscop´ıa electr´onica fueron Lynch et al. quienes observaron que el contraste obtenido con una peque˜na defocalizaci´on es proporcional al laplaciano de la densidad de carga proyectada [93]. Este concepto tambi´en se ha aprovechado en Microscop´ıa electr´onica de Lorentz para obtener 1.2 Fundamentos 9 x y z z = 0 z -z Planode medida Planode medida Planode recuperación l Objetoo medio Onda plana Figura 1.1: Ilustraci´on de un t´ıpico experimento de recuperaci´on de fase basado en la ETI en un plano z= 0 a la salida de un objeto o medio difractivo. mapas de potenciales electrost´atico y magnetost´atico de materiales magn´eticos y superconductores delgados [86, 88, 89, 99, 102, 104, 105, 111–113]. Finalmente, la combinaci´on de las t´ecnicas de recuperaci´on de fase con la reconstrucci´on tomogr´afica ha mejorado mucho sus resultados, pues las im´agenes de fase contrastadas, menos afectadas por ruido que las im´agenes obtenidas por absorci´on, proporcionan una reconstrucci´on m´as limpia [98, 114, 115]. La recuperaci´on es cuantitativa y el resultado ha mostrado una calidad en ocasiones superior a otros m´etodos de medida empleados comercialmente [82, 85, 96, 97, 116]. 1.2. Fundamentos Consideremos la imagen de la Figura 1.1. La ilustraci´on describe un problema gen´erico de recuperaci´on de fase. Una onda electromagn´etica plana se propaga a lo largo de la direcci´on del eje ze incide en un medio u objeto d´ebilmente dispersivo que no absorbe significativamente la radiaci´on. El objeto o medio atravesado por la onda est´a localizado en el entorno del eje ´optico en el semiespacio z < 0 a una distancia ldel plano de recuperaci´on de fase, z= 0. Como el objeto o medio es d´ebilmente inhomog´eneo, la onda transmitida es paraxial, de la forma pI(r, z)eikW(r,z)eikz con 10 La Ecuaci´on de Transporte de Intensidad como t´ecnica de recuperaci´on de fase r= (x, y), y se propaga a lo largo del eje z. Aunque no se puede deducir informaci´on del objeto a partir de la absorci´on de radiaci´on, la fase de la onda, W(r, z), s´ı transporta con frecuencia importante informaci´on acerca del mismo. En esta memoria consideraremos el problema de la reconstrucci´on de la fase de la onda W(r) = W(r,0) en el plano z= 0, a partir de las medidas de la distribuci´on de irradiancia, I(r, z), en uno o m´as planos ortogonales al eje ´optico. De aqu´ı en adelante, obviaremos en este cap´ıtulo la dependencia en r de estas magnitudes. La ecuaci´on que describe c´omo cambia instant´aneamente la irradiancia de una onda electromagn´etica paraxial en funci´on de las variaciones transversales de la fase, es la Ecuaci´on de Transporte de Intensidad (ETI) [24]: −I′(z) = ∇·(I∇W) = ∇I·∇W+I∇2W(1.1) donde I′(z) = ∂I (z)/∂z es la derivada axial instant´anea en un plano zy∇= ∂xˆ x+∂yˆ yes el gradiente transversal. La Ec. (1.1) puede resolverse para recuperar la fase de la onda, W, en un plano zsi se conoce I′en dicho plano. Ichikawa et al. dan una interpretaci´on de los t´erminos de la derecha de la Ec. (1.1) que explica los fundamentos de los sensores de frente de onda basados en la ETI [117]. El primer t´ermino, ∇I· ∇W, representa la variaci´on de irradiancia causada por el desplazamiento transversal del haz no uniforme (∇I6= 0) debido a la inclinaci´on local de los frentes de onda, cuya normal o direcci´on del rayo viene dada por ∇W. El segundo t´ermino, I∇2W, puede interpretarse como las variaciones de irradiancia causadas por las convergencias o las divergencias locales del haz, que dependen de las curvaturas locales de los frentes de onda, ∇2W. Ambos sumandos determinan la variaci´on axial de la irradiancia que causan las deformaciones transversales de los frentes de onda a medida que el haz se propaga a lo largo del eje z. La Ec. (1.1) proporciona una relaci´on entre la irradiancia y la fase de la onda v´alida en el espacio libre. No obstante, en la pr´actica, la apertura finita del dispositivo que capta la luz se manifiesta con la presencia de una pupila, Po=P(r,0), en el plano de recuperaci´on de fase. Esta abertura influye en la evoluci´on axial de la irradiancia, lo que obliga a reformular la Ec. (1.1) de la forma (ver Ap´endice A): −I′(0) = [∇·(Io∇W)] Po−IoWrδc(1.2) donde Ioes la distribuci´on de irradiancia en el interior de la abertura en el plano de recuperaci´on o plano de la pupila, z= 0, y δces una distribuci´on lineal de Dirac que coincide con el borde de la pupila. Por tanto, en el plano z= 0, el primer t´ermino de 1.2 Fundamentos 11 Divergencia delhaz Convergencia delhaz I z I( )< o I z I( )> o Io Desplazamiento delborde Desplazamiento delborde Planodela pupila Planodemedida z= 0z z W Figura 1.2: Esquema para ilustrar el significado f´ısico de la ETI particularizada a intensidad uniforme en el plano de recuperaci´on de fase. Los rayos convergen y divergen debido a la influencia de las curvaturas locales del frente de onda, redistribuy´endose la irradiancia de los planos posteriores. En consecuencia, la intensidad aumenta en unas zonas y disminuye en otras, y se producen movimientos del borde que deforman la regi´on iluminada. la derecha describe la variaci´on axial instant´anea de la irradiancia en el interior de la abertura y el segundo en el borde o frontera. Si la distribuci´on de luz en el interior de la pupila es uniforme, entonces ∇Io· ∇Wse anula y la ETI toma la conocida forma [118]: −I′(0) = IoPo∇2W−IoWrδc(1.3) La Ec. (1.3) es una ecuaci´on de Poisson con condici´on de contorno de Neumann, que relaciona el incremento o decremento local de irradiancia que se produce durante la propagaci´on axial de la luz con la curvatura local de los frentes de onda en el interior de la pupila (ver Figura 1.2) y su derivada normal local en el borde. Constituye la base del funcionamiento de un tipo concreto de sensores que se derivan de la ETI, el llamado sensor de curvatura. Las Ecs. (1.2) y (1.3) tienen soluci´on anal´ıtica ´unica, salvo una constante o pist´on, siempre y cuando no haya ceros en la distribuci´on de irradiancia [48], ya que la presencia de puntos donde la intensidad se anula lleva a la aparici´on de dislocaciones de fase en forma de puntos rama o v´ortices [48, 90, 110, 119–122], que en este trabajo 18 La Ecuaci´on de Transporte de Intensidad como t´ecnica de recuperaci´on de fase Los errores de los datos de entrada y el comportamiento no lineal de la din´amica del borde en la propagaci´on de la luz reducen la precisi´on de esta medida. Disponer de una expresi´on donde se describa la dependencia del error de la estimaci´on de la derivada normal local con la separaci´on zfacilitar´ıa la determinaci´on objetiva de la posici´on de los planos que garantiza la mejor respuesta del sensor para unas condiciones de medida determinadas. En este trabajo nos proponemos calcularla. Paralelamente, estudiaremos hasta qu´e punto influye en la calidad de la medida el tama˜no de la regi´on de frontera. Tanto la ETI como las ecuaciones que describen la din´amica del borde de la regi´on iluminada pertenecen al ´ambito geom´etrico. En este trabajo no vamos a tener en cuenta el efecto de la difracci´on en la recuperaci´on, pues su influencia puede considerarse despreciable [52]. 1.5.3. Mejorar el rendimiento del sensor de curvatura bas´andose en tomar medidas de irradiancia en m´ultiples planos La respuesta de un sensor basado en la ETI que opere en la separaci´on ´optima podr´ıa ser insuficiente para garantizar el grado de resoluci´on en la fase recuperada que la aplicaci´on exige. Dentro del ´ambito de la Microscop´ıa, en los ´ultimos a˜nos se ha ganado resoluci´on gracias al dise˜no de algoritmos lineales alternativos, no basados en la ETI, que propagan menos error de no linealidad trabajando con datos de irradiancia tomados en planos m´as separados de z= 0, buscando reducir el efecto del ruido [36, 94, 139, 140]. Tambi´en se han dise˜nado algoritmos iterativos para refinar la fase inicialmente reconstruida por medio de la ETI a partir de las distribuciones de irradiancia en dos planos [16, 90, 98, 134, 141, 142]. Otra de las estrategias utilizadas consiste en emplear algoritmos que precisan de datos de irradiancia en m´ultiples planos [23, 121, 143–145]. Con estos m´etodos se incrementa la resoluci´on reduci´endose el efecto del ruido en la reconstrucci´on de la fase, pero a costa de incrementar las necesidades computacionales. Los algoritmos de este tipo basados en la ETI mejoran los resultados a base de promediar la derivada axial estimada a partir de m´ultiples pares de planos de irradiancia [144]. Siguiendo esta l´ınea de aumentar la resoluci´on espacial de la recuperaci´on a cambio de incrementar el n´umero de datos de irradiancia, proponemos en este trabajo aplicar f´ormulas en diferencias finitas para aproximar la derivada axial de la ETI 1.5 Objetivos 19 utilizando m´as de dos planos de irradiancia. Con ello se persigue el doble objetivo de minimizar la contribuci´on del ruido de medida y de reducir el error de no linealidad de la aproximaci´on, sin que ello implique necesariamente un incremento sustancial del gasto computacional al reconstruir la fase. 1.5.4. Analizar la contribuci´on de las medidas de frontera Se ha visto que para poder obtener una soluci´on ´unica del problema de recuperaci´on de fase representado por las Ecs. (1.2) ´o (1.3) es necesario estimar I′(0) en el interior de la regi´on iluminada y tambi´en en la frontera. Pero tambi´en se ha hecho constar la dificultad experimental de obtener resultados precisos para la derivada normal de la fase en el borde, pudiendo jugar la selecci´on del tama˜no de la regi´on de frontera un papel activo en ello. Las simetr´ıas del objeto o medio dispersivo, o bien el tipo de iluminaci´on, pueden fijar condiciones para evitar las medidas de frontera [24, 136, 146]. En la ´ultima parte de esta memoria, analizaremos con detalle el peso de la contribuci´on de las medidas en el contorno para estudiar c´omo obtener la m´axima informaci´on acerca de la fase con solo medidas en el interior del dominio. Cap´ıtulo 2 Optimizaci´on de la medida del interior En el cap´ıtulo anterior hemos introducido la Ecuaci´on de Transporte de Intensidad particularizada a irradiancia constante en el plano de recuperaci´on de fase, z= 0, como base del funcionamiento del sensor de curvatura. En este cap´ıtulo analizaremos el rendimiento de este tipo de dispositivos al recuperar la curvatura de los frentes de onda en el interior de la pupila, ∇2W, que depende de su capacidad para estimar la derivada axial de la irradiancia, I′(0). El ruido de las medidas de irradiancia en los dos planos de medida y el comportamiento no lineal intr´ınseco del sensor son las principales limitaciones que degradan su respuesta. Con esta informaci´on nos proponemos expresar matem´aticamente la varianza de la curvatura estimada en el interior en funci´on de estas dos contribuciones, que han de estar ligadas por medio de la distancia de separaci´on entre los planos. A partir de aqu´ı, dise˜naremos un m´etodo para deducir una posici´on de los planos de medida que garantice objetivamente la respuesta ´optima del sensor. Llevaremos a cabo el desarrollo de este procedimiento para la f´ormula de Stirling y para la f´ormula de Newton en dos planos. Dejaremos para m´as adelante el an´alisis de configuraciones alternativas del sensor que utilizan f´ormulas en diferencias finitas con m´as de dos datos de entrada. Finalmente, un ejemplo ilustrar´a la efectividad de este procedimiento de optimizaci´on, aplic´andolo a un sensor de curvatura acoplado a un dispositivo compensador de turbulencias atmosf´ericas. En concreto, utilizando los par´ametros caracter´ısticos del telescopio Gemini y de las condiciones atmosf´ericas promedio de su ubicaci´on, se determina la posici´on de los 22 Optimizaci´on de la medida del interior planos que proporciona la mejor medida de la curvatura de los frentes de onda*. 2.1. Factores que limitan la estimaci´on de la curvatura del frente de onda De acuerdo con la Ec. (1.3), el laplaciano de la fase en el interior de la pupila Po del plano z= 0 se relaciona con la derivada axial de la irradiancia en dicho plano, I′(0), de la forma: ∇2W=−I′(0) Io (2.1) donde Ioes la irradiancia uniforme en z= 0. Los dos agentes esenciales que deterioran la calidad del laplaciano de la fase recuperado son el ruido de los datos de entrada y el error del modelo en diferencias finitas utilizado para aproximar la derivada axial. En general, los detectores de irradiancia proporcionan se˜nales afectadas por ruido aditivo, que aumenta conforme la luminosidad se hace m´as escasa. El ruido se debe tanto a errores intr´ınsecos del detector como a errores provocados por la aleatoriedad espacio-temporal de los fotones al incidir en ´el [6, 14, 147]. El primer tipo se suele llamar ruido de detecci´on, que, en el caso de las c´amaras CCD, se conoce como ruido de lectura [147–149]. El segundo se llama ruido fot´onico y es fundamental en la detecci´on de luz y, por tanto, inevitable. Tampoco se puede evitar el efecto de la radiaci´on de fondo, que contribuye negativamente restando contraste a las im´agenes [149]. El ruido fot´onico depende del nivel de luz, a diferencia del ruido de detecci´on y del ruido de fondo, que son de la misma magnitud con independencia de la luminosidad. Por otro lado, como la derivada axial se construye aproxim´andola en diferencias finitas a primer orden, se presupone que el comportamiento del dispositivo es lineal en todo el rango axial donde se captan medidas de irradiancia [16, 53, 123, 131, 132]. Sin embargo, esto es solo cierto para valores de zpeque˜nos, donde el efecto del ruido es mayor. *Un c´alculo an´alogo se realiza en el Ap´endice C con la f´ormula de Roddier. Sus principios son semejantes, pero no completamente iguales, a los de la f´ormula central en diferencias finitas, requiriendo un estudio aparte, que no se incluye en el cuerpo principal de esta memoria para dinamizar su lectura. No obstante, debemos se˜nalar que este hecho no debe restar importancia a los resultados alcanzados, dada la trascendencia de la f´ormula de Roddier en la destacable expansi´on de los sensores de curvatura en el ´ambito de la ´ Optica Astron´omica. Es por ello que el ejemplo que se incluye en este cap´ıtulo ilustra los resultados obtenidos en el Ap´endice C para la f´ormula de Roddier. 2.2 Varianza de la estimaci´on de la derivada axial 23 2.2. Varianza de la estimaci´on de la derivada axial En las l´ıneas que siguen, calcularemos la varianza de la estimaci´on de la derivada axial en funci´on de la posici´on de los planos en el eje z. Modelaremos la presencia de ruido de medida: sin p´erdida de generalidad en los resultados, para calcular y analizar el comportamiento de la varianza de la curvatura, supondremos aqu´ı que las fluctuaciones de irradiancia en cada punto son de media cero e independientes del nivel de se˜nal. Esta situaci´on es equivalente a considerar como fluctuaciones por ruido de detecci´on las fluctuaciones provocadas por el ruido fot´onico, cometi´endose un error de estad´ıstica que no influye apreciablemente en el comportamiento cualitativo de la varianza. No obstante, en el Ap´endice B se incluye el c´alculo completo para condiciones de iluminaci´on limitadas, que es especialmente ´util en la aplicaci´on de ´ Optica Astron´omica. 2.2.1. F´ormula central de Stirling en diferencias finitas La f´ormula central en diferencias finitas permite estimar I′(0) = I′(r,0) en z= 0 a partir de las irradiancias en dos planos situados sim´etricamente a ambos lados de z= 0: ˆ I′(0) = i+−i− 2z(2.2) donde z > 0 y el circunflejo sobre ˆ I′denota una estimaci´on de I′.i+=i(r, z) e i−=i(r,−z) son las se˜nales proporcionadas por el dispositivo de detecci´on en los planos +zy−zrespectivamente, que en cada punto ry en cada instante de tiempo tson la suma de la irradiancia exacta sin ruido o determinista, I±=I(r,±z), y el ruido, n±=n±(r,±z), es decir: i±=I±+n±(2.3) Consideramos que los errores n±son de media cero, covarianza nula entre puntos distintos y varianza E[n±]2=σ2en todos los puntos de cada plano z, donde E{·} denota valor medio. Con el objetivo de obtener estimaciones precisas de la derivada axial de la intensidad en un punto cualquiera r, queremos encontrar la separaci´on entre los planos, z=z(r), tal que se verifica la siguiente condici´on: s2hˆ I′(0)i=Ehˆ I′(0) −I′(0)i2= m´ın (2.4) 24 Optimizaci´on de la medida del interior En definitiva, debemos calcular cu´anto se aleja la aproximaci´on experimental de la derivada axial, ˆ I′(0), del valor determinista de la derivada, I′(0), en el plano z= 0 y despu´es hallar el valor de zque minimiza la varianza s2hˆ I′(0)i. Suponiendo que la intensidad, I(r, z), es una funci´on univaluada, continua, cuya derivada primera es tambi´en continua y cuya derivada segunda existe en el rango (0, z) de operaci´on del sensor, entonces podemos aplicar el teorema de Taylor para la intensidad en cada punto ren torno a z= 0, que establece que [150]: I±=Io±I′(0)z+1 2I′′ (ζ±)z2(2.5) donde en el tercer sumando del miembro de la derecha, ζ+=ζ+(r) es un punto dentro del intervalo (0, z) (ζ−=ζ−(r) es otro tal que ζ−∈(−z, 0)). Este sumando representa el resto del desarrollo en serie de Taylor a primer orden y da cuenta del error de la aproximaci´on lineal de la derivada axial. Introduciendo el resultado anterior en la Ec. (2.2) y teniendo en cuenta la Ec. (2.3), la derivada axial estimada en cada punto rdel plano z= 0 es: ˆ I′(0) = I′(0) + n++n− 2z+1 4[I′′ (ζ+)−I′′ (ζ−)] z(2.6) A la vista de esta ecuaci´on, queda claro que la precisi´on de ˆ I′(0) depende del peso del segundo sumando, debido al error de las medidas, que limitar´a el valor de z inferiormente, y del tercer sumando, debido al error de la aproximaci´on en diferencias finitas, que lo limitar´a superiormente. Suponiendo que no hay correlaci´on entre el ruido en los puntos de ambos planos y que ambas fuentes de error contribuyen independientemente, se deduce que la varianza de la aproximaci´on de la derivada axial con la f´ormula de Stirling es: s2hˆ I′(0)i=σ2 2z2+1 16 [I′′ (ζ+)−I′′ (ζ−)]2z2(2.7) donde cabe recalcar que la varianza var´ıa de un punto a otro, ya que ζ±=ζ±(r). 2.2.2. F´ormula progresiva o regresiva de Newton en diferencias finitas La f´ormula progresiva o regresiva en diferencias finitas permite aproximar la derivada axial en un punto rdel plano de recuperaci´on de fase, I′(0) = I′(r,0), utilizando como datos de entrada la irradiancia en dicho plano, io=io(r,0), y la 2.3 Selecci´on del plano de medida 25 irradiancia en otro plano a una distancia za la derecha (+z, progresiva) o a la izquierda (−z, regresiva), i±=i±(r,±z): ˆ I′(0) = i±−io ±z(2.8) siendo z > 0. Las medidas de irradiancia se suponen degradadas por ruido de la misma estad´ıstica descrita en la secci´on anterior y obtendremos estimaciones precisas de I′(0) en ren la medida que la separaci´on z=z(r) verifique la condici´on (2.4). De forma an´aloga procedemos a aplicar el Teorema de Taylor a la intensidad en cada punto ren torno a z= 0, Ec. (2.5), y sustituimos en la Ec. (2.8) teniendo en cuenta la Ec. (2.3), obteni´endose una relaci´on similar a la Ec. (2.6) entre la aproximaci´on progresiva o regresiva y el valor determinista de la derivada axial. Suponiendo que la contribuci´on del ruido y del error de la aproximaci´on de la derivada en diferencias finitas son independientes y teniendo en cuenta la hip´otesis de no correlaci´on del ruido en los puntos rde ambos planos se deduce que la varianza de la derivada axial estimada en cada rcon la f´ormula de Newton es: s2hˆ I′(0)i=2σ2 z2+1 4[I′′ (ζ±)]2z2(2.9) donde ζ+=ζ+(r) es un punto tal que ζ+∈(0, z) (ζ−∈(−z, 0)). 2.3. Selecci´on del plano de medida Las varianzas (2.7) y (2.9) demuestran que en id´enticas condiciones de medida, la precisi´on que se alcanza al estimar la curvatura del frente de onda depende de la f´ormula concreta que se escoja para aproximar la derivada axial. En estas expresiones generales de las varianzas, el primer sumando est´a relacionado con el ruido de las medidas de irradiancia y aumenta conforme la separaci´on zdisminuye. El segundo sumando depende exclusivamente de la forma del frente de onda a recuperar y describe, en funci´on de z, cu´anto se desv´ıa el modelo en diferencias finitas respecto de la curvatura exacta del frente de onda. Las Ecs. (2.7) y (2.9) re´unen caracter´ısticas que impiden emplearlas directamente para hallar una ´unica separaci´on entre los planos que optimice la respuesta del sensor en la configuraci´on correspondiente. Por un lado, informan de que, en general, el error de la aproximaci´on de la derivada axial depende del punto rconcreto. Por otro lado, el segundo t´ermino de estas varianzas es desconocido en cada punto r, ya que ζ=ζ(r). El m´etodo de selecci´on de la separaci´on entre los planos deber´a eludir 26 Optimizaci´on de la medida del interior estos problemas para proporcionar una ´unica posici´on de los mismos que optimice globalmente la respuesta del sensor en el plano de recuperaci´on de fase, z= 0. Antes de proceder a dise˜narlo, es necesario realizar un an´alisis previo de la propagaci´on axial de la irradiancia en funci´on de las caracter´ısticas del frente de onda. 2.3.1. An´alisis de la evoluci´on axial de la intensidad El problema de la propagaci´on de un haz difractado por una abertura en z= 0 da lugar a integrales dobles de la forma [151, 152]: ZZ g(x, y, z)eikf(x,y,z)dx dy (2.10) donde, en la aproximaci´on de Fresnel: f(x, y, z) = 1 2z(x−xc)2+ (y−yc)2+W(x, y) (2.11) y donde ges la amplitud de la onda en el plano de la abertura. Tanto gcomo fson independientes de ky el dominio de integraci´on est´a determinado por la abertura. En el rango de la aproximaci´on de Fresnel, la distribuci´on de irradiancia en cada plano axial depende, en esencia, de las caracter´ısticas de la fase a segundo orden en los llamados puntos cr´ıticos de primera clase, rc= (xc, yc), del plano z= 0 [152, 153], es decir: I(r, z) = Io(rc) 1 + z∇2W(rc) + z2H[W(rc)] (2.12) donde ∇2W(rc) es el operador laplaciano y H[W(rc)] es el operador hessiano de la fase evaluados en dichos puntos. Los puntos cr´ıticos son puntos del plano z= 0 donde el argumento de la exponencial del integrando de la Ec. (2.10), f=f(r), verifica: ∂f ∂x =∂f ∂y = 0 (2.13) La Ec. (2.12) es v´alida siempre y cuando r= (x, y) no se encuentre en las inmediaciones de las c´austicas [153]. En realidad, su denominador deber´ıa estar en valor absoluto, pero aqu´ı se ha eliminado y se ha tomado el signo positivo porque trabajamos con planos de irradiancia situados antes de las c´austicas. Esta ecuaci´on nos muestra que la suposici´on de que la irradiancia evoluciona linealmente en la variable zes una simplificaci´on, que solo es cierta para valores de zsuficientemente peque˜nos, pues, en realidad, la evoluci´on axial de la intensidad en el rango de Fresnel depende de las caracter´ısticas de la fase de forma m´as compleja. 2.3 Selecci´on del plano de medida 27 En general, a cada punto rle corresponde un punto cr´ıtico distinto del plano z= 0 y la ecuaci´on que relaciona rcon rces complicada de resolver. Sin embargo, la Ec. (2.12) no presenta una dependencia expl´ıcita en los puntos cr´ıticos en el caso particular de que la fase Wen el plano z= 0 sea cuadr´atica, esto es, que tenga la forma: Wo=a+bx +cy +dx2+ey2+fxy (2.14) donde arepresenta el pist´on, bycson inclinaciones y d,eyfest´an relacionados con el operador laplaciano y con el operador hessiano de la fase cuadr´atica Wo. Dado que estos dos operadores son invariantes, podemos simplificar Worotando los ejes de coordenadas de tal modo que la Ec. (2.14) pueda reescribirse en su forma can´onica [154]: Wo=A+Bx +Cy +Dx2+Ey2(2.15) De las Ecs. (2.14) y (2.15), se deduce que ∇2Wo= 2(d+e) = 2(D+E) y H(Wo) = 4de −f2= 4DE. De este modo, la intensidad evoluciona en cada punto r= (x, y) de acuerdo con: I(z) = Io 1 + 2(D+E)z+ 4DEz2(2.16) donde no hay dependencia en los puntos cr´ıticos, pues cada uno de los puntos de la pupila tiene el mismo valor del laplaciano y del operador hessiano. Por tanto, dado un plano z, la irradiancia es la misma para cada punto rde dicho plano. La aplicaci´on de la Ec. (2.16) en las f´ormulas (2.7) y (2.9) de la §2.2 proporciona, en este caso, la misma varianza de ˆ I′(0) para todos los puntos rde un mismo plano axial, cuya minimizaci´on permite calcular un ´unico valor para la separaci´on ´optima de los planos de medida. Aunque esta no es la situaci´on que normalmente nos encontraremos globalmente en el plano de la pupila, como iremos viendo, el estudio detallado de estas varianzas en funci´on de los par´ametros DyEnos permitir´a dise˜nar una estrategia para resolver problemas m´as complejos. 2.3.2. Varianza de la estimaci´on del laplaciano de una fase cuadr´atica En este apartado analizaremos la evoluci´on del error de la estimaci´on de la curvatura de un frente de onda cuadr´atico Wocon f´ormulas en diferencias finitas, en funci´on de la posici´on axial de los planos de irradiancia, z. En concreto, haremos un an´alisis detallado para distintos valores de las curvaturas principales en el eje x,D, y 34 Optimizaci´on de la medida del interior 0.05 0.09 0.13 0.17 0.21 0.25 0.29 0.33 0.37 0.41 0.45 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 (a) z ǫo↓ zopt =0.17 ↓ zmin =0.21 ←zmin =0.21 ↑ zmin =0.29 −0.34 −0.3 −0.26−0.22−0.18−0.14 −0.1 −0.06 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 (b) −z ǫo ↓ zopt =−0.11 ↓ zmin =−0.13 ↑ zmin =−0.13 ↑ zmin =−0.19 0.06 0.1 0.14 0.18 0.22 0.26 0.3 0.34 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 (c) z ǫo ↓ zmin =0.13 ↓ zmin =0.15 ↑ zmin =0.15 ↑ zmin =0.22 D= 1/2, E = 0 2 (D+E) = 1 D= 1/6, E = 1/3 2 (D+E) = 1 D= 1/4, E = 1/4 2 (D+E) = 1 D= 1/3, E =−1/6 2 (|D|+|E|) = 1 Figura 2.5: Evoluci´on del error de las curvaturas de los frentes de onda cuadr´aticos que verifican 2(|D|+|E|) = 1 en funci´on de la posici´on de los planos, z, suponiendo que se estiman con la f´ormula central (a), con la f´ormula regresiva (b) y con la f´ormula progresiva (c), bajo la influencia de ruido de detecci´on con desviaci´on relativa σ/Io= 0,01. y 2.5(c) para la f´ormula progresiva de Newton, se muestra la desviaci´on ǫosolo para aquellas combinaciones de DyEtales que 2(|D|+|E|) iguala la cota superior, V= 1. En estas figuras tambi´en se indica la posici´on del plano donde se minimiza ǫo, en las mismas unidades que la posici´on de la l´ınea focal. Es decir, supongamos que el frente es un cilindro parab´olico con D= 1/(4f) = 1/2 y E= 0. Si V= 1 m−1, entonces la l´ınea focal se halla en f = 0,5 m y el m´ınimo en zopt = 0,17 m, para la f´ormula central, y en zopt = 0,11 m, para las f´ormulas de Newton. A partir de ellas se deduce que de todas las combinaciones de DyEtales que 2(|D|+|E|) = V, la curva que est´a por encima de todas las dem´as en la zona del m´ınimo corresponde al cilindro parab´olico, seguida, en orden decreciente, por la curva del paraboloide el´ıptico, la del paraboloide parab´olico y, por ´ultimo, la del hiperb´olico. 2.3 Selecci´on del plano de medida 35 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 z I(z) D=1/2;E=0 D=1/6;E=1/3 D=1/4;E=1/4 Figura 2.6: Evoluci´on axial de la intensidad para las fases cuadr´aticas de laplaciano unidad analizadas. En l´ınea de puntos se representa la aproximaci´on lineal proporcional al laplaciano. La raz´on para este comportamiento es que, aunque consideremos frentes cuadr´aticos con el mismo valor del laplaciano, la contribuci´on del error de no linealidad var´ıa entre unos y otros. En concreto, de entre los cuatro est´andares de cu´adricas estudiados, la irradiancia asociada al cilindro parab´olico es la que evoluciona axialmente m´as alejada del comportamiento lineal (ver Figura 2.6). Adicionalmente, cabe destacar que la curva asociada al frente de onda plano se encuentra siempre por debajo de cualquier otra curva, sea cual sea el valor de Vy el tipo de fase cuadr´atica con que comparemos. Por ´ultimo, la Figura 2.5 nos indica tambi´en que en las mismas condiciones de medida, la f´ormula central proporciona mejores resultados que las f´ormulas de Newton. 2.3.3. Selecci´on de la separaci´on ´optima y cotas superior e inferior de error En la secci´on anterior, hemos calculado y representado frente a zla desviaci´on est´andar de la estimaci´on experimental de la curvatura de un frente de onda 36 Optimizaci´on de la medida del interior cuadr´atico en funci´on del nivel de ruido, σ/Io, y de las caracter´ısticas del frente expresadas en sus curvaturas principales DyE. El an´alisis de estas curvas de error nos permitir´a deducir ahora la separaci´on entre los planos que proporciona la mejor respuesta del sensor ante frentes de onda cuadr´aticos. En concreto, la Figura 2.5 nos informa de que si lo ´unico que conocemos de la fase cuadr´atica en el plano z= 0 es que su laplaciano alcanza como m´aximo un valor V, entonces la peor respuesta del sensor de curvatura se da ante los frentes con forma de cilindro parab´olico de curvatura V. Por eso, la mejor respuesta del sensor ante frentes de onda cuadr´aticos se garantizar´ıa separando los planos de irradiancia la distancia z=zopt que viene determinada por el m´ınimo de la curva de error correspondiente al cilindro parab´olico de curvatura V. Bas´andonos en este resultado alcanzado con este caso sencillo, encontraremos en esta secci´on la separaci´on entre planos que optimiza la respuesta global del sensor de curvatura en el interior ante frentes de onda m´as complicados. Si suponemos que el frente de onda, W, en los puntos del plano z= 0 es suave, continuo y diferenciable, entonces existen entornos Ω donde Wpuede aproximarse localmente por superficies cuadr´aticas, Wo, de la forma (2.14), donde ahora a representa el pist´on, b=∂W/∂x yc=∂W/∂y son las inclinaciones locales de Wen Ω, y d= (1/2)(∂2W/∂x2), e= (1/2)(∂2W/∂y2) y f=∂2W/∂x∂y est´an relacionados con el operador laplaciano y con el operador hessiano locales de la fase Wen el entorno Ω. La uni´on de todos los entornos Ω recubre la superficie del frente de onda completamente. Como los operadores laplaciano y hessiano son invariantes, podemos simplificar la cu´adrica (2.14) rotando los ejes de coordenadas y reescribiendo Woen su forma can´onica, es decir, en la forma descrita por la Ec. (2.15). Los par´ametros DyEde la aproximaci´on cuadr´atica local Wo, representan ahora las curvaturas principales locales del frente de onda Wen el interior del entorno Ω. Por tanto, el operador laplaciano y el operador hessiano locales de Wasociados a cada entorno Ω son constantes en cada plano axial para todos los puntos de Ω y la evoluci´on de la irradiancia en el eje zse puede expresar localmente como en la Ec. (2.12), en funci´on de los par´ametros DyEdel entorno. Estas curvaturas de cada entorno Ω provocan convergencias o divergencias de los haces, que se manifiestan con cambios de la distribuci´on de luz dentro de la regi´on iluminada en funci´on de la posici´on axial z(ver Figura 1.2). En definitiva, de acuerdo con esta hip´otesis, podemos afirmar que las Ecs. (2.18) y (2.20) describen el error de la recuperaci´on de la curvatura local del frente de onda Wen un entorno Ω. De este modo, concentr´emonos ahora en la estimaci´on de las curvaturas locales de 2.3 Selecci´on del plano de medida 37 un frente de onda, suponiendo que existen entornos Ω donde Wpuede aproximarse por superfices cuadr´aticas Wotales que 2 (|D|+|E|) es menor o igual que un cierto valor normalizado a la unidad, que denotamos por V. En la Figura 2.5, zopt denota la posici´on de los planos que minimiza la desviaci´on est´andar de la estimaci´on de la curvatura local del frente de onda en un entorno con forma de cilindro parab´olico con V= 1. Si en Ω la cu´adrica fuese de otro tipo, entonces, de acuerdo con esa figura, situando los planos en zopt est´a garantizado que su laplaciano se estimar´ıa con menos desviaci´on que si fuese un cilindro parab´olico. Veamos qu´e implicaciones tiene este comportamiento a nivel local en la recuperaci´on del laplaciano a nivel global, es decir, en la uni´on de todos los entornos del plano z= 0. Si los planos de medida se sit´uan en las posiciones z < zopt, entonces las curvaturas locales estimadas en cualquier parte del interior de la pupila tendr´an siempre m´as error que si los planos se colocan en las posiciones z=zopt. Por otro lado, si las irradiancias se detectan en z > zopt, aunque mejora la estimaci´on en los entornos donde la aproximaci´on de Wes un paraboloide o un plano, empeora la recuperaci´on de las curvaturas locales en entornos donde Wse aproxima por cilindros parab´olicos. De este modo, si hubiese una densidad significativa de entornos dentro del dominio definido por la pupila donde el frente de onda se aproximase mejor por medio de superficies cil´ındricas de curvatura m´axima (V= 1), si z6=zopt, globalmente empeorar´ıa la estimaci´on del laplaciano de la fase en z= 0. En consecuencia, si la ´unica informaci´on a priori disponible acerca de los frentes de onda es que su curvatura m´axima es V, minimizando ǫ2 opara un cilindro parab´olico de curvatura V, es decir, D=V/2 y E= 0, obtendr´ıamos zopt, la mejor posici´on de los planos que es compatible con esta informaci´on. Todo esto independientemente de la f´ormula en diferencias finitas que se utilice para configurar el sensor. En t´erminos generales, con el objetivo de escoger adecuadamente sus componentes principales, cuando se procede al dise˜no de cualquier dispositivo sensorial de fase debe disponerse de un cierto grado de conocimiento a priori acerca de los frentes de onda que se quieren medir. Por ejemplo, en el caso de los sensores de frente de onda ShackHartmann, es necesaria cierta informaci´on previa de la fase para dise˜nar las matrices de microlentes con el espaciado y la apertura num´erica adecuados. En el caso de los sensores de curvatura, una vez que se ha determinado el nivel de ruido de los planos de irradiancia, el paso siguiente es fijar su posici´on, para lo que es necesario hacer uso de la informaci´on de la fase disponible. Si conocemos una estimaci´on de las m´aximas curvaturas locales de la fase en la abertura, tenemos ya la capacidad para hallar objetivamente una separaci´on ´optima siguiendo el m´etodo descrito hasta 38 Optimizaci´on de la medida del interior ahora. No obstante, normalmente estos datos no siempre se conocen. En su lugar podr´ıamos disponer de cierta informaci´on estad´ıstica de los frentes de onda que nos oriente a la hora de escoger un valor para V. Ilustr´emoslo de la siguiente forma. Por ejemplo, supongamos que la fase en el plano z= 0 var´ıa en el tiempo y se conoce lo suficiente acerca de su comportamiento estad´ıstico como para que podamos: estimar L=|h∇2W(r)i|, donde h·i indica valor medio espacio-temporal en el total de la pupila. estimar ∆L=qh[∇2W(r)−h∇2W(r)i]2i, la desviaci´on cuadr´atica media espacio-temporal calculada en la totalidad de la pupila. Si no hay informaci´on adicional, proponemos tomar los planos de irradiancia en la posici´on zopt que minimiza ǫo(z) tomando D=V/2 y E= 0, siendo V=L+α∆L, donde αes un par´ametro que cuantifica la probabilidad de que V∈[L−α∆L, L +α∆L]*. De este modo, nos aseguramos de que la desviaci´on de la curvatura local de Wen un entorno Ω cualquiera, ǫo(z), oscila, para todo z, entre las siguientes cotas: ǫinf(z)≤ǫo(z)≤ǫsup(z) (2.21) donde ǫinf(z) es la mayor cota inferior, que est´a fijada por la curva de un frente de onda plano, y ǫsup(z) es la menor cota superior, que viene indicada por la curva del cilindro parab´olico de curvatura V. La magnitud de estas cotas ser´a distinta seg´un se emplee una u otra f´ormula en diferencias finitas. A continuaci´on, expresaremos las ecuaciones de las cotas en funci´on de la variable normalizada h=V z, que permite generalizar el error y la posici´on ´optima de los planos a cualquier valor de V: a) F´ormula central de Stirling ǫo(z) de la Ec. (2.18) para D=V/2 y E= 0 representa para cada zla menor cota superior de error de la estimaci´on de la curvatura, ǫsup, que, con la normalizaci´on h=V z, es: ǫsup V2=1 2h2σ Io2 +h2 1−h22 (2.22) *Si α= 1, la probabilidad es del 68,3 %, si α= 2, la probabilidad es del 95,4 % y si α= 3, la probabilidad llega al 99,7 % [156]. 2.3 Selecci´on del plano de medida 39 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 h ǫsup,inf V (a) ↓ hopt =0.079 ↓ hopt =0.17 ↓ hopt =0.35 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 h ǫsup,inf V (b) ↓ hopt =0.037 ↓ hopt =0.11 ↓ hopt =0.29 σ /Io=0.001 σ /Io=0.01 σ /Io=0.1 Figura 2.7: Cota superior de error, ǫsup, y cota inferior de error, ǫinf, de la estimaci´on de la curvatura local de la fase en Ω con la f´ormula central (a) y con las f´ormulas de Newton (b), para varios niveles del ruido de detecci´on, σ/Io. mientras que ǫo(z) para D=E= 0 representa la mayor cota inferior de error de la estimaci´on de la curvatura, ǫinf: ǫinf V2=1 2h2σ Io2 (2.23) En la Figura 2.7(a) representamos estas cotas para algunos ejemplos de σ/Io, apreci´andose la disminuci´on del rendimiento del sensor conforme aumenta el 40 Optimizaci´on de la medida del interior nivel de ruido en las medidas de irradiancia. Para hallar una relaci´on entre la separaci´on ´optima, hopt, y el nivel de ruido de las medidas, σ/Io, derivamos la Ec. (2.22) respecto de he igualamos a cero, obteniendo: −σ Io2 1−h2 opt3+ 4h6 opt = 0 (2.24) b) F´ormulas progresiva y regresiva de Newton La Ec. (2.20) particularizada a D=V/2 y E= 0 proporciona la menor cota superior de error para cada h=V z: ǫsup V2=2 h2σ Io2 +h 1−h2 (2.25) mientras que si se particulariza a D=E= 0, proporciona la mayor cota inferior de error: ǫinf V2=2 h2σ Io2 (2.26) En la Figura 2.7(b) representamos estas cotas para los mismos ejemplos de σ/Io que en el apartado anterior. Comparando esta figura con la Figura 2.7(a), se ve con claridad el aumento de la incertidumbre de las curvaturas recuperadas cuando se usan las f´ormulas de Newton en comparaci´on con la f´ormula central. De la minimizaci´on de la Ec. (2.25) respecto de hobtenemos la siguiente ecuaci´on, que relaciona la separaci´on ´optima, hopt, con el nivel de ruido en las medidas de irradiancia, σ/Io: −2σ Io2 (1 −hopt)3+h4 opt = 0 (2.27) Por ´ultimo, la Figura 2.8(a) muestra la dependencia de la separaci´on ´optima, hopt, con el nivel de ruido, σ/Io, suponiendo que se implementa la f´ormula central. En los apartados (b) y (c), se representan las cotas superior, ǫsup/V , e inferior, ǫinf/V , de error de la curvatura en la posici´on ´optima. Asimismo, se indican los datos de los coeficientes β,γyδdel ajuste de las cotas y de las posiciones ´optimas a la siguiente ecuaci´on en funci´on del nivel de ruido: fun = βrσ Io+γrσ Io2 +δrσ Io3 (2.28) Sustituyendo en la Ec. (2.28) los coeficientes del apartado (a) de la Figura 2.8, podemos calcular f´acilmente la separaci´on entre los planos, hopt, que optimiza la 2.3 Selecci´on del plano de medida 41 0 0.12 0.24 0.36 0.48 0.6 hopt (a) β= 1,620 ±0,002 γ=−1,942 ±0,007 δ= 1,000 ±0,006 0 0.2 0.4 0.6 0.8 1 ǫsup V (b) β= 0,4581 ±0,0005 γ= 1,077 ±0,002 δ=−0,250 ±0,001 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.18 0.36 0.54 0.72 0.9 ǫinf V σ Io (c) β= 0,3692 ±0,0004 γ= 0,924 ±0,001 δ=−0,168 ±0,001 Figura 2.8: Posici´on ´optima de los planos y cotas superior e inferior de error de la curvatura local de la fase en Ω, en funci´on del nivel de ruido en el plano de la abertura, σ/Io, para la f´ormula central. Con el s´ımbolo “o” se indican algunos de los puntos utilizados para realizar el ajuste, que se obtienen a partir de las Ecs. (2.22), (2.23) o (2.24) para cada valor de σ/Io. La l´ınea continua representa el ajuste de los puntos. respuesta del sensor de curvatura para unas condiciones de medida determinadas. El mismo procedimiento para los coeficientes de los apartados (b) y (c) nos proporciona las cotas entre las que oscila el error de las curvaturas locales estimadas con el sensor. En la Figura 2.9 se hace la representaci´on an´aloga para las f´ormulas de Newton. Las Figuras 2.8 y 2.9 constituyen uno de los resultados m´as importantes de este trabajo de investigaci´on, pues son una referencia objetiva para determinar no solo la posici´on de los planos que optimiza la respuesta del sensor para unas condiciones de operaci´on dadas, sino tambi´en las cotas superior e inferior entre las que oscila el error de las curvaturas locales estimadas, que ning´un otro criterio previo proporcionaba. Para poder aplicar este procedimiento solo resulta necesario tener informaci´on sobre el nivel de ruido de medida, σ/Io, y un m´ınimo conocimiento a priori acerca de la fase, es decir, una cota superior de las curvaturas locales, V. 42 Optimizaci´on de la medida del interior 0 0.12 0.24 0.36 0.48 0.6 hopt (a) β= 1,1600 ±0,0003 γ=−0.8563 ±0,0009 δ= 0,3040 ±0,0007 0 0.46 0.92 1.38 1.84 2.3 ǫsup V (b) β= 1,67750 ±0,00005 γ= 1,0281 ±0,0001 δ= 0,0897 ±0,0001 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.4 0.8 1.2 1.6 2 ǫinf V σ Io (c) β= 1,18410 ±0,00005 γ= 1,0947 ±0,0002 δ= 0,0833 ±0,0001 Figura 2.9: Posici´on ´optima de los planos y cotas superior e inferior de error de la curvatura local de la fase en Ω, en funci´on del nivel de ruido en el plano de la abertura, σ/Io, para las f´ormulas de Newton. Con el s´ımbolo “o” se indican algunos de los puntos utilizados para realizar el ajuste, que se obtienen a partir de la Ec. (2.25), (2.26) o (2.27) para cada valor de σ/Io. La l´ınea continua representa el ajuste de los puntos. 2.4. Ejemplo: Selecci´on de la separaci´on ´optima para el telescopio Gemini En las l´ıneas siguientes, se ilustra con un ejemplo el m´etodo de selecci´on de la posici´on de los planos descrito en la §2.3. El problema escogido consiste en calcular la separaci´on entre los planos que optimiza la respuesta del sensor de curvatura acoplado al sistema de ´optica adaptativa del telescopio Gemini. Este sensor de curvatura, al igual que la mayor´ıa de los que son utilizados en ´ Optica Astron´omica, se ha dise˜nado de acuerdo a la configuraci´on de Roddier (ver §1.3). El sensor se emplea para medir las aberraciones de frentes de onda distorsionados por turbulencias atmosf´ericas y trabaja con un n´umero de fotones limitado, siendo significativa la influencia del ruido fot´onico. A lo largo de este cap´ıtulo, hemos 2.4 Ejemplo: Selecci´on de la separaci´on ´optima para el telescopio Gemini 43 supuesto que las medidas de irradiancia solo est´an afectadas por ruido de detecci´on independiente del nivel de se˜nal. Los resultados que se han mostrado para las f´ormulas t´ıpicas en diferencias finitas son cualitativamente an´alogos a los que se obtienen para un n´umero de fotones disponible limitado, seg´un puede comprobarse en el Ap´endice B. No obstante, la diferente estad´ıstica del ruido de medida en un caso y en el otro se traduce en cambios cuantitativos en los resultados, de ah´ı la utilidad de los datos presentados en ese ap´endice. Por su parte, los c´alculos espec´ıficos para la f´ormula de Roddier, adecuados para su aplicaci´on en este ejemplo, se muestran en el Ap´endice C. Para seleccionar los planos de irradiancia, nos basaremos solamente en los datos t´ecnicos publicados en el trabajo de Rigaut et al. de 1997, relativos a los componentes del telescopio Gemini y a las condiciones de turbulencia medias de su ubicaci´on en Mauna Kea (Hawaii) [56]. 2.4.1. Caracter´ısticas del sistema y de la turbulencia En la Tabla 2.1 se resumen los par´ametros principales que caracterizan el telescopio Gemini y las turbulencias atmosf´ericas medias de su emplazamiento [56]. Los frentes turbulentos que se quieren compensar provienen de una estrella natural de magnitud mR= 14,7 en la banda R*. El ritmo de emisi´on de fotones indicado est´a ya calculado para la totalidad de la abertura del telescopio Gemini. De los datos de la tabla se desprende que el sistema de ´ Optica adaptativa es real, con 56 subaberturas repartidas en una pupila de 7,9 m de di´ametro con una obstrucci´on central de 1,2 m de di´ametro. Cada subabertura tiene asociado un detector y un actuador que modifica la forma de un espejo deformable. Se emplean fotodiodos de avalancha (APD), limitados solo por ruido fot´onico, que operan en una longitud de onda λ= 0,7µm. El 50 % de los fotones que emite la estrella no se detectan por efecto de la atm´osfera, de la transmisi´on instrumental del telescopio y de la eficiencia cu´antica del detector. Por ´ultimo, el comportamiento estad´ıstico de las turbulencias atmosf´ericas en Mauna Kea est´a caracterizado por el par´ametro de Fried, ro= 25 cm, para una longitud de onda λ= 0,55 µm. Con los par´ametros de la Tabla 2.1 calcularemos un valor para Vy para el n´umero de fotones promedio disponibles en cada subabertura, N, necesarios para hallar zopt de acuerdo con los resultados del Ap´endice C. *La banda astron´omica Rest´a centrada en una longitud de onda λ= 0,7µm con un ancho de banda efectivo de 0,22 µm [149]. 50 Optimizaci´on de la medida de frontera planos de medida. En estas situaciones, todo apunta a que las medidas de frontera se realizan empleando las mismas posiciones de los planos que las medidas del interior. No obstante, a priori no existe ninguna raz´on que garantice que la separaci´on elegida para optimizar la medida del laplaciano permita obtener tambi´en medidas de frontera de buena calidad. Debido a los problemas que surgen al poner en pr´actica el procedimiento anterior, se han propuesto otros m´etodos para estimar Wr, que utilizan como datos de entrada el radio de la regi´on iluminada en dos planos adyacentes ortogonales al eje z[131]. Mediante este tipo de m´etodos, se elimina la contribuci´on del ruido de las medidas de irradiancia, pero se sustituye por la inevitable incertidumbre de las medidas de la posici´on del borde. En este cap´ıtulo desarrollaremos un procedimiento, an´alogo al realizado para el interior, para optimizar la se˜nal de frontera proporcionada por el sensor de curvatura. Comenzaremos describiendo con m´as detalle los fundamentos de estos m´etodos y proseguiremos analizando las principales limitaciones que presentan y que determinan la precisi´on del resultado final. Por ´ultimo, modelaremos el error de Wren funci´on de la posici´on de los planos de medida en el eje z, que nos permitir´a encontrar la separaci´on entre ellos que minimiza las contribuciones del error en conjunto. 3.1. M´etodos de medida de la derivada normal de la fase en la frontera Cualquiera de los dos m´etodos que se emplean para obtener Wren el borde de la pupila tiene como finalidad cuantificar el cambio de dimensi´on radial de la zona geom´etricamente iluminada. El primero, consiste en medir directamente la deformaci´on radial de las regiones iluminadas en dos planos paralelos perpendiculares al eje ´optico del sistema [131]. Lo llamaremos m´etodo directo. En el segundo, el desplazamiento radial del borde se mide indirectamente, como resultado de implementar la ETI aproximando la derivada axial con diferencias finitas a primer orden de medidas de irradiancia en dos planos promediadas radialmente en el contorno [136]. Nos referiremos a ´el como m´etodo indirecto*. *En este cap´ıtulo estudiaremos exclusivamente las aproximaciones cl´asicas en diferencias finitas de la derivada axial. La f´ormula de Roddier, semejante en algunos aspectos a la f´ormula central, no se analiza en este cap´ıtulo para facilitar su comprensi´on. No obstante, dado su destacable protagonismo en el campo de la ´ Optica Astron´omica, merece un estudio aparte en este trabajo, que se incluye en 3.1 M´etodos de medida de la derivada normal de la fase en la frontera 51 En esta secci´on deduciremos las ecuaciones que relacionan las medidas que lleva a cabo el sensor con la derivada normal de la fase en la frontera. Partiendo de estas ecuaciones, describiremos las simplificaciones que se realizan para poderlas implementar experimentalmente en el sensor de curvatura. Este paso es muy importante para identificar, m´as adelante, las principales fuentes de error que deterioran la respuesta del sensor en relaci´on con las medidas de frontera. 3.1.1. M´etodo directo de medida de la deformaci´on de la regi´on iluminada La hip´otesis de partida consiste en que la deformaci´on radial de las regiones iluminadas en dos planos zy−z, es proporcional a la derivada de la fase normal al borde de la pupila, Wr. De este modo, en cada ´angulo θse toman, experimentalmente, las seguientes relaciones para estimar Wren la frontera [131]: R+−R−= 2zˆ Wr R±−R=±zˆ Wr (3.1) donde R±=R±(θ) es el radio de la regi´on iluminada en las posiciones ±zyˆ Wr= ˆ Wr(θ) es la estimaci´on de la derivada normal en el punto (R, θ) del borde. Dado que en la primera ecuaci´on se usan datos en dos planos situados sim´etricamente a ambos lados de z= 0, la llamaremos f´ormula central. En cuanto a la segunda ecuaci´on, la llamaremos f´ormula progresiva o bien f´ormula regresiva seg´un se utilice el radio en el plano +zo bien en en el plano −z, respectivamente. En lo que sigue, demostraremos que las Ecs. (3.1) son, en realidad, una aproximaci´on a primer orden de la din´amica del borde de la regi´on iluminada en planos perpendiculares al eje z, v´alida solo muy cerca de z= 0. Partimos de las ecuaciones de los rayos asociadas a una onda paraxial de fase W que se propaga a lo largo del eje z[152, 153]: x=xc±zWxc y=yc±zWyc(3.2) donde Wxc= (∂W/∂x)(xc,yc)yWyc= (∂W/∂y)(xc,yc). Estas ecuaciones describen la trayectoria del rayo que pasa por el punto (xc, yc) del plano de la pupila, z= 0, en funci´on de las coordenadas y de las inclinaciones de la fase, WxcyWyc, en dicho punto. el Ap´endice C. 52 Optimizaci´on de la medida de frontera R+R q qox y dz= 0 z Figura 3.1: Ejemplo de evoluci´on de la zona geom´etricamente iluminada desde el plano de la abertura, z= 0, a un plano desplazado, z. En coordenadas polares se expresan de la forma:    rcos(θ−θc) = rc±zWrc rsen(θ−θc) = ±zWθc rc (3.3) donde Wrc= (∂W/∂r)(rc,θc)yWθc= (∂W/∂θ)(rc,θc), siendo (rc, θc) las coordenadas polares del punto del plano z= 0 por donde pasa el rayo. Si la distancia de separaci´on escogida, z, es lo suficientemente peque˜na como para que el desplazamiento angular de los rayos, θ−θc, no sea muy grande, podemos truncar a primer orden los desarrollos en serie de potencias de las funciones trigonom´etricas de las ecuaciones anteriores, obteni´endose:    r=rc±zWrc θ=θc±zWθc rrc (3.4) Pero para llegar a las Ecs. (3.1) es necesario hacer una aproximaci´on m´as fuerte, es decir, hay que suponer que la deformaci´on de la regi´on iluminada es estrictamente radial, esto es, θ=θc. En definitiva, se supone que el efecto de la derivada angular de la fase en los puntos de la frontera no es significativo desde el plano de recuperaci´on, z= 0, al plano z. En la Figura 3.1 se muestra un ejemplo para ilustrar la evoluci´on axial de la zona iluminada, donde el vector dindica el desplazamiento del borde, que, en general, 3.1 M´etodos de medida de la derivada normal de la fase en la frontera 53 tiene componente radial y tangencial. Si se trabaja con distancias zlejos de las c´austicas, podemos asegurar que los puntos del borde en el plano z= 0, (R, θc), dan lugar a puntos del borde en un plano ±z, (R±, θ) [153]. Por tanto, a partir de las Ecs. (3.3), se deduce que el radio de la regi´on iluminada evoluciona en cada ´angulo θde acuerdo con la siguiente ecuaci´on: R±=s(R±zWrc)2+zWθc R2 (3.5) Si la distancia zes muy peque˜na, entonces la zona iluminada se deforma linealmente con zy proporcionalmente a Wrc, lo que posibilita el m´etodo descrito por las Ecs. (3.1), es decir: R±−R=s(R±zWrc)2+zWθc R2 −R∼ =Rr1±2zWrc R−R ∼ =R1±zWrc R−R=±zWrc(3.6) De acuerdo con este modelo, en cada ´angulo θ, la derivada normal ser´ıa la constante de proporcionalidad que liga el desplazamiento radial del borde con el desplazamiento longitudinal del plano de medida en el eje z. Pero, en realidad, la Ec. (3.6) informa de que la linealidad del proceso de deformaci´on de la zona iluminada solo se puede suponer en un rango limitado del eje z. El motivo, ya adelantado en la Figura 3.1, es que, en realidad, la derivada angular de la fase, Wθc, interviene tambi´en en la traslaci´on del borde de la regi´on iluminada. Como las medidas de los radios de las dos zonas iluminadas se llevan a cabo para un ´angulo θconstante, la derivada radial calculada como diferencia entre ambos radios se asigna al punto (R, θ) de la frontera de la pupila cuando, en realidad, corresponde a otro punto (R, θc) de las inmediaciones (ver Ecs. (3.4)). Esta asignaci´on incorrecta induce un error que crece con la separaci´on z, debido a la contribuci´on no siempre despreciable del t´ermino zWθc/R de la Ec. (3.5). 3.1.2. Metodo indirecto de medida de la deformaci´on de la regi´on iluminada Con este m´etodo se trata de obtener la derivada de la fase normal al borde de la pupila a partir de medidas de la derivada axial de la intensidad, de acuerdo con la relaci´on entre ambas magnitudes descrita por la ETI (Ec. (1.3)): I′(0) = IoWrδc(3.7) 54 Optimizaci´on de la medida de frontera En esta ecuaci´on, Wrva modulada por una delta de Dirac en el borde de la pupila, δc=δ(r−R). De este modo, una forma sencilla de extraer la derivada normal es promediando la ETI radialmente en un intervalo que incluya la delta de Dirac del borde [136], es decir: I′(0) = 1 LZR+L/2 R−L/2 I′(0) dr =1 LZR+L/2 R−L/2 IoWrδ(r−R)dr =Io Wr L(3.8) donde I′(0) representa, en un ´angulo θ, el promedio radial de la derivada axial dentro de una regi´on alrededor del borde de anchura radial L: la regi´on de frontera. Experimentalmente, se implementa alguna de las aproximaciones en diferencias finitas de la derivada axial, cuyos datos de entrada son las distribuciones de irradiancia de dos planos axiales ±zy/o z= 0. En la Figura 3.2(a) se representan ejemplos de los perfiles de irradiancia en un ´angulo θdel plano z,I(z), del plano z= 0, I(0), y del plano −z,I(−z). Matem´aticamente pueden expresarse de la forma:        I(0) = Io[1 −H(r−R)] I(z) = I+[1 −H(r−R+)] I(−z) = I−[1 −H(r−R−)] (3.9) siendo R±el radio de la regi´on iluminada asociado al ´angulo θen el plano ±z(Ec. (3.5)), I±la irradiancia en el punto (r, θ) de dicho plano en el supuesto de que no hubiese una pupila en el plano z= 0 (en aproximaci´on de Fresnel, Ec. (2.12)) y H es la funci´on de salto unitario, es decir [158]: H(r−R) =    1 si r > R 0 si r≤R ,(3.10) que representa el efecto que causa la pupila del plano z= 0. Por tanto, las aproximaciones en diferencias finitas de la derivada axial que t´ıpicamente se implementan para estimar la derivada normal de la fase en la frontera, ˆ Wr, pueden expresarse de la forma siguiente: a) F´ormula Central o de Stirling En un punto (R, θ) del borde de la abertura, la derivada axial se estima como: ˆ I′(0) = I+−I− 2z(3.11) 3.1 M´etodos de medida de la derivada normal de la fase en la frontera 55 r r r r I(0) I z(- ) Io Io Io r I z( ) R+ RR (a) rr I z I z( )+ (- ) Io R+ RL r I z I z( )- (- ) Io R+ RL (b) Figura 3.2: (a) Ejemplo de perfiles de irradiancia en planos ±zyz= 0 perpendiculares al eje del sistema. (b) Suma y resta de las irradiancias de los planos +zy−z. con z > 0, donde el circunflejo sobre ˆ I′indica que se trata de una estimaci´on de la derivada axial y donde: I±=1 LZR+L/2 R−L/2 I±[1 −H(r−R±)] dr =1 LZR± R−L/2 I±dr (3.12) 56 Optimizaci´on de la medida de frontera El algoritmo es m´as preciso a medida que se reduce la distancia entre los planos, es decir: ˆ I′(0) = 1 2z1 LZR+ R−L/2 I+dr −1 LZR− R−L/2 I−dr∼ =Io LR+−R−(R−−R) 2z ∼ =Io L zWr−(−zWr) 2z=Io Wr L=I′(0) (3.13) donde tambi´en se comprueba que este m´etodo y el descrito por las Ecs. (3.1) son equivalentes cuando los planos se sit´uan muy cerca de z= 0. b) F´ormula progresiva y regresiva de Newton En este caso, la derivada axial se aproxima en cada punto (R, θ) del borde de la forma: ˆ I′(0) = I±−Io ±z(3.14) con I±de la Ec. (3.12) e Io: Io=1 LZR+L/2 R−L/2 Io[1 −H(r−R)] dr =Io LR−R+L 2=Io 2(3.15) Llevando a cabo un proceso an´alogo al del apartado anterior se deduce f´acilmente que: ˆ I′(0) ∼ =Io Wr L=I′(0) (3.16) 3.2. Factores que limitan la medida de frontera En la secci´on anterior, hemos introducido los dos m´etodos que se utilizan para obtener la derivada normal de la fase en el borde de la pupila. Las f´ormulas empleadas experimentalmente son aproximaciones a primer orden de las ecuaciones exactas y la estimaci´on de Wrpor estos m´etodos puede llegar a ser muy imprecisa. Concretamente, como ya vimos, los algoritmos funcionan correctamente siempre y cuando pueda suponerse que un punto (R, θ) del borde de la pupila en z= 0 se mueve hasta (R±zWr, θ) en el plano ±z. No obstante, la deformaci´on real del borde difiere, en general, de este comportamiento lineal. Acabamos de comprobar que la precisi´on de cualquiera de los dos m´etodos se incrementa conforme los planos se acercan [16, 131, 136]. Pero tambi´en hay que tener en cuenta la presencia de errores en los datos de entrada, cuyos efectos se hacen notar sobre todo en las posiciones cercanas al plano z= 0. En el caso de 3.3 Varianza de la estimaci´on de la derivada normal en la frontera 57 los m´etodos directos, los errores de los datos son las incertidumbres de las medidas de la posici´on del borde de la zona iluminada. En los m´etodos indirectos, son las fluctuaciones de las irradiancias debido al ruido. En definitiva, los errores de las medidas y el efecto de la derivada angular de la fase en la din´amica del borde de la regi´on iluminada son dos problemas fundamentales, ´ıntimamente ligados por la posici´on de los planos, z, que determina la precisi´on del resultado final. Por otro lado, para que el m´etodo basado en la ETI sea preciso, la derivada que debemos medir en la frontera debe ser infinitamente estrecha, seg´un indica la delta de Dirac, δc(Ecs. (3.7) y (3.8)). No obstante, en la pr´actica, el salto de ˆ I′(0) = [I(z)−I(−z)] /(2z) en la regi´on de frontera tiene un perfil aproximadamente rectangular, que solo converge a la delta de Dirac en el l´ımite en que zes peque˜no [158] (ver Figura 3.2(b)). Adem´as, puesto que Wrse estima promediando ˆ I′(0) en la frontera, la precisi´on del m´etodo indirecto tambi´en depende de la anchura radial seleccionada para la regi´on de frontera, L. 3.3. Varianza de la estimaci´on de la derivada normal en la frontera A continuaci´on calcularemos la varianza de la derivada normal estimada siguiendo los m´etodos directo e indirecto en funci´on de la distancia de separaci´on z. 3.3.1. M´etodo directo En el m´etodo directo solo distinguimos dos fuentes de error. La primera es el error del modelo, es decir, cu´anto se aleja de la ecuaci´on exacta la ecuaci´on simplificada que emplea el algoritmo de medida. La segunda es la propagaci´on de los errores de los datos de entrada a la derivada normal. Vamos a suponer que estas dos contribuciones son independientes y la varianza total de ˆ Wrser´a, por tanto, la suma de ambas. Supongamos que el m´etodo de detecci´on de la posici´on del borde de las regiones iluminadas proporciona medidas de los radios R±yRen cada ´angulo θcon una incertidumbre s(R±) = s(R) = σ. Podemos estimar la derivada normal en un punto (R, θ) utilizando las Ecs. (3.1), cuyas varianzas son, para cada f´ormula: a) F´ormula central Las medidas son R+yR−en el plano +zy el plano −z, respectivamente. 58 Optimizaci´on de la medida de frontera Entonces, la varianza de ˆ Wres: s2(ˆ Wr) = ∂ˆ Wr ∂R+!2 s2(R+) + ∂ˆ Wr ∂R−!2 s2(R−) + ˆ Wr−Wr2 =σ2 2z2+R+−R− 2z−Wr2 (3.17) con R±de la Ec. (3.5). b) F´ormula progresiva y regresiva Las medidas son R±en el plano ±zyRen z= 0. La varianza de la derivada normal de la fase en el borde de la pupila es, entonces: s2(ˆ Wr) = ∂ˆ Wr ∂R±!2 s2(R±) + ∂ˆ Wr ∂R !2 s2(R) + ˆ Wr−Wr2 =2σ2 z2+R±−R ±z−Wr2 (3.18) con R±de la Ec. (3.5). 3.3.2. M´etodo indirecto En este caso la derivada normal de la fase en la frontera se estima a partir de medidas de irradiancia. La luz se capta utilizando detectores situados en el borde de la pupila en una regi´on que abarca desde R−L/2 a R+L/2, como se muestra en la Figura 3.3 para una pupila circular [16]. Cada detector, de anchura radial L, cubre un ´angulo ∆θdonde realiza medidas promedidas de irradiancia. El resultado para la derivada normal, ˆ Wr, es una medida promediada en el ´angulo ∆θ, con centro en θ, que se asigna al punto (R, θ) [136]. No obstante, en este c´alculo te´orico vamos a suponer que en cada ´angulo θdel contorno tenemos un detector que proporciona medidas de irradiancia, es decir, supondremos que podemos estimar Wren cada ´angulo θ. Llamemos I+eI−a las irradiancias promedio en las posiciones zy−z, respectivamente, captadas por el detector del ´angulo θ. La cantidad de luz que capta este detector var´ıa dependiendo de la distancia zporque las aberraciones de los frentes de onda act´uan deformando la regi´on iluminada, como se muestra en la Figura 3.3. En el m´etodo indirecto las dos fuentes principales de error son el error del modelo utilizado para aproximar la derivada axial y el ruido de las medidas de irradiancia, 3.3 Varianza de la estimaci´on de la derivada normal en la frontera 59 q Dq z=0 z L R Figura 3.3: Ejemplo de matriz de detectores en la regi´on de frontera. En l´ınea de puntos se dibujan las zonas iluminadas en el plano de la pupila, z= 0, y en un plano desplazado, z. que consideramos contribuyen independientemente. Sin p´erdida de generalidad en los resultados, para calcular y analizar el comportamiento de la varianza de ˆ Wr, s2(ˆ Wr), supondremos que las fluctuaciones de irradiancia en cada punto son de media cero e independientes del nivel de se˜nal. Esta es una simplificaci´on, ya que as´ı las fluctuaciones por ruido fot´onico se consideran como fluctuaciones por ruido de detecci´on, cometi´endose un error de estad´ıstica que no influye apreciablemente en el comportamiento cualitativo de s2(ˆ Wr). No obstante, en el Ap´endice B se incluye este c´alculo en funci´on del n´umero de fotones que llegan al detector, de utilidad cuando se precisen datos cuantitativos en aplicaciones donde la disponibilidad de luz es limitada, como, por ejemplo, en ´ Optica Astron´omica. Si n±es el ruido de las medidas de irradiancia en cada punto r(ver §2.2.1), el ruido del promedio radial de la irradiancia en la frontera, η±, es tambi´en de media cero y de varianza σ2: s2(I±) = E[η±]2=σ2(3.19) donde E{·} denota valor medio. Teniendo en cuenta estas hip´otesis, calcularemos s2(ˆ Wr) partiendo de la varianza de la derivada axial de la intensidad estimada en diferencias finitas. An´alogamente al procedimiento de optimizaci´on de las medidas del interior, con el objetivo de obtener estimaciones precisas de I′(0) en cada ´angulo θde la frontera, encontraremos la 66 Optimizaci´on de la medida de frontera A continuaci´on analizaremos el comportamiento de ǫoen funci´on de zy de los par´ametros de Wo(Ec. (3.25)). Clasificaremos los frentes de onda cuadr´aticos utilizando la misma notaci´on que en la P´agina 28 del Cap´ıtulo 2, en funci´on de los valores de sus curvaturas principales DyE. Es decir, distinguiremos entre frentes de onda con forma de cilindro parab´olico, paraboloide el´ıptico, paraboloide parab´olico y paraboloide hiperb´olico. Como antes, consideraremos solo combinaciones de los par´ametros tales que 2 (|D|+|E|)≤V, siendo Vuna cota superior, para garantizar que los frentes de onda cuadr´aticos son suaves con peque˜nas desviaciones en torno a un frente de onda plano. No obstante, como Wordepende tambi´en de los coeficientes ByC, debemos acotar el valor m´aximo de estas inclinaciones para mantener la condici´on de paraxialidad de los frentes de onda. En concreto, analizaremos frentes de onda cuadr´aticos tales que √B2+C2≤U, donde Ues otra cota superior. En las Figuras 3.4(a), 3.4(b) y 3.4(c) representamos ǫofrente a zen la configuraci´on central, la configuraci´on progresiva y la configuraci´on regresiva, respectivamente, del m´etodo directo. Para cada plano z, el valor de ǫo(z) representado en las curvas es el resultado de calcular el m´aximo de la varianza de ˆ Wor de todos los ´angulos θen dicho plano z. Suponemos que la medida del radio en un ´angulo θ,R± o, tiene una incertidumbre σ= 0,01, siendo el radio de la pupila R= 1, y que las cotas de las curvaturas principales y de las inclinaciones son V=U= 1*. No obstante, tan solo hemos representado los resultados para aquellos valores de DyEtales que 2 (|D|+|E|) = 1 y aquellos de ByCtales que √B2+C2= 1. Esto implica que las curvas representadas corresponden a los frentes de onda menos suaves, para los que, en vista de los resultados obtenidos para la medida del interior, se da la respuesta m´as deteriorada del sensor de curvatura. Tomando estos mismos valores de los par´ametros de Wo, de las cotas UyVy del radio de la pupila, R, en la Figura 3.5 se muestra la varianza de la derivada normal estimada utilizando el m´etodo indirecto en sus tres configuraciones, en funci´on de la posici´on axial de los planos, z. El error se ha calculado suponiendo σ/Io= 0,01 para la fluctuaci´on de la irradiancia promedio detectada en un ´angulo θdel plano z= 0. Hemos tomado L= 0,3Rpara la dimensi´on radial del anillo de detectores, pues solo vamos a estudiar el comportamiento de las curvas de error en funci´on de zy del tipo de cu´adrica. La influencia de la elecci´on del tama˜no de la regi´on de frontera *En la pr´actica, normalmente se verificar´a 1/V > R. Sin embargo, se han tomado estos valores sencillos para U,VyRporque nuestro objetivo aqu´ı es aclarar el comportamiento cualitativo de las curvas de error y no tanto dar resultados cuantitativos. Por otro lado, en la §3.6 se ver´a que para frentes de onda turbulentos, U≈V. 3.4 Selecci´on del plano de medida 67 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.1 0.2 0.3 0.4 0.5 ǫo (a) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.1 0.2 0.3 0.4 0.5 ǫo (b) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.1 0.2 0.3 0.4 0.5 z ǫo (c) D= 1/2, E = 0 2 (D+E) = 1 D= 1/3, E = 1/6 2 (D+E) = 1 D= 1/4, E = 1/4 2 (D+E) = 1 D= 1/3, E =−1/6 2 (|D|+|E|) = 1 D= 0, E = 0 2 (D+E) = 0 Figura 3.4: Curvas de error de ˆ Worestimada con la configuraci´on central (a), con la configuraci´on progresiva (b) y con la regresiva (c), utilizando el m´etodo directo para σ/Io= 0,01, L= 0,3R, R= 1 y fases cuadr´aticas con U=V= 1. Tambi´en se representa el resultado para el frente de onda plano. se analizar´a m´as adelante. Por ´ultimo, al igual que en las curvas anteriores para el m´etodo directo, en todas las figuras se ha calculado la desviaci´on est´andar en cada uno de los ´angulos y se ha representado el m´aximo de todos los errores obtenidos en cada plano z. En definitiva, podemos asegurar que cada curva muestra la menor 68 Optimizaci´on de la medida de frontera 0 0.015 0.03 0.045 0.06 0.075 0.09 0.105 0.12 0.135 0.15 0 0.1 0.2 0.3 0.4 0.5 ǫo (a) 0 0.015 0.03 0.045 0.06 0.075 0.09 0.105 0.12 0.135 0.15 0 0.1 0.2 0.3 0.4 0.5 ǫo (b) 0 0.015 0.03 0.045 0.06 0.075 0.09 0.105 0.12 0.135 0.15 0 0.1 0.2 0.3 0.4 0.5 z ǫo (c) D= 1/2, E = 0 2 (D+E) = 1 D= 1/3, E = 1/6 2 (D+E) = 1 D= 1/4, E = 1/4 2 (D+E) = 1 D= 1/3, E =−1/6 2 (|D|+|E|) = 1 D= 0, E = 0 2 (D+E) = 0 Figura 3.5: Curvas de error de ˆ Worestimada con la configuraci´on central (a), con la configuraci´on progresiva (b) y con la regresiva (c), utilizando el m´etodo indirecto para σ/Io= 0,01, L= 0,3R, R= 1 y fases cuadr´aticas con U=V= 1. Tambi´en se representa el resultado para el frente de onda plano. cota superior del error de ˆ Wor. En ambos m´etodos e independientemente del tipo de f´ormula empleada, el comportamiento cualitativo de estas curvas es similar al observado para las curvas de error del laplaciano recuperado en el interior. Sus principales caracter´ısticas son: 3.4 Selecci´on del plano de medida 69 A medida que disminuye la distancia entre los planos el error tiende a infinito y a medida que uno cualquiera de los dos planos se aproxima a la posici´on de la l´ınea focal en xo en ym´as peque˜na, el error se incrementa tambi´en dr´asticamente. Manteniendo constantes las condiciones de detecci´on, σ, siempre hay una desviaci´on m´ınima absoluta, ǫo(zm´ın), que decrece a medida que disminuyen √B2+C2o 2 (|D|+|E|). La posici´on correspondiente de los planos, zm´ın, se aleja de z= 0. Es importante indicar que si las caracter´ısticas del frente de onda se mantienen constantes, la posici´on del m´ınimo absoluto, zm´ın, se aleja de z= 0 a medida que aumenta la incertidumbre de los datos de entrada, σ. Entre las curvas asociadas a una configuraci´on del sensor y las asociadas a otra, tambi´en se aprecian diferencias. Si hablamos del m´etodo indirecto, a igualdad de valores de los par´ametros de la fase, del ruido y de la dimensi´on radial de los detectores, la f´ormula m´as precisa para estimar la derivada normal es la f´ormula central, seguida de la f´ormula regresiva y, por ´ultimo, la f´ormula progresiva (ver Figura 3.6). En el caso del m´etodo directo, la configuraci´on central proporciona el sensor con la mejor calidad de respuesta. No obstante, hay que indicar que si todos los par´ametros tienen signos opuestos, entonces los resultados representados para la f´ormula progresiva pasan a ser los de la regresiva y a la inversa. Cuando los planos se sit´uan sim´etricamente a ambos lados del plano z= 0, el resultado no cambia. Teniendo todo esto en cuenta, desde el punto de vista cuantitativo e independientemente del tipo de configuraci´on del sensor, se observa que la curva que est´a por encima de todas las dem´as corresponde al cilindro parab´olico, seguida de la correspondiente al paraboloide el´ıptico, al paraboloide parab´olico y, por ´ultimo, al hiperb´olico. Es decir, an´alogamente a lo visto para las medidas del interior, la peor respuesta del sensor en la frontera se da ante los frentes de onda cil´ındricos porque la contribuci´on del error del modelo es en este caso superior al de los otros tipos de frentes cuadr´aticos. Por ´ultimo, es necesario recordar que para representar ǫohemos calculado el m´aximo error en cada ´angulo θ. En realidad, las curvas dibujadas coinciden precisamente con las curvas de error de ˆ Woren los ´angulos θdonde |Wor| ≪ |Woθ|, es decir, para obtener la mejor estimaci´on de Wordebemos minimizar el error de 70 Optimizaci´on de la medida de frontera 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.13 0.15 0.17 0.19 0.21 0.23 0.25 0.27 0.29 0.31 0.33 z ǫo Fórmula central Fórmula progresiva Fórmula regresiva Figura 3.6: Representaci´on frente a zde la desviaci´on est´andar m´axima de la derivada normal de la fase cuadr´atica con √B2+C2= 1, D= 1/2, E= 0, estimada con el m´etodo indirecto empleando las f´ormulas en diferencias finitas indicadas, siendo R= 1 y L= 0,3Ren cada caso. ˆ Woren dichos ´angulos. En la Figura 3.7 se representan WoryWoθpara un cilindro parab´olico con B= 1, C= 0, D= 1/2 y E= 0. Todas las curvas dibujadas para este tipo de cu´adrica corresponden a ´angulos θdonde Wor≈0 y Woθes m´axima en valor absoluto (en este ejemplo, en torno a 5π/4 ´o 7π/4). 3.5. Selecci´on de la separaci´on ´optima El an´alisis realizado en la secci´on anterior para fases cuadr´aticas en el plano z= 0 nos permitir´a a continuaci´on dise˜nar un m´etodo para optimizar la estimaci´on de la derivada normal de una fase cualquiera en la frontera. Este m´etodo complementa el procedimiento de optimizaci´on de la recuperaci´on del laplaciano del interior desarrollado en el Cap´ıtulo 2. En concreto, las figuras representadas nos informan de que si lo ´unico que conocemos acerca de la fase en el plano de la pupila, z= 0, es que se trata de una fase cuadr´atica (Ec. (3.25)) y que sus par´ametros verifican 3.5 Selecci´on de la separaci´on ´optima 71 0 pi/4 pi/2 3pi/4 pi 5pi/4 3pi/2 7pi/4 2pi −1.5 −1 −0.5 0 0.5 1 1.5 θ W or W oθ Figura 3.7: Derivada normal y angular del frente cil´ındrico parab´olico con B= 1, C= 0, D= 1/2 yE= 0 en la frontera. que 2 (|D|+|E|) alcanza como m´aximo un valor Vy que √B2+C2est´a tambi´en acotada por un valor U, entonces la peor respuesta del sensor de curvatura se da ante los frentes con forma de cilindro parab´olico de curvatura Vy de inclinaci´on U. Por este motivo, separando los planos de irradiancia la distancia z=zmin, que viene determinada por el m´ınimo de la curva de error correspondiente al cilindro parab´olico de curvatura Ve inclinaci´on U, se garantiza la mejor respuesta del sensor ante cualquier frente de onda cuadr´atico como los considerados. Llegados a este punto, nos basaremos en el procedimiento que se utiliz´o para optimizar las medidas del laplaciano del interior, para dise˜nar un m´etodo an´alogo para seleccionar la posici´on de los planos que optimiza las medidas de frontera. As´ı, supongamos que el frente de onda W=W(x, y) en los puntos del plano z= 0 es suave, continuo y diferenciable. Entonces, tambi´en puede aproximarse localmente en entornos, Ω, situados en el borde de la regi´on iluminada, por superficies cuadr´aticas 72 Optimizaci´on de la medida de frontera de la forma siguiente: Wo=a+bx +cy +dx2+ey2+fxy (3.36) donde arepresenta el pist´on, b=∂W/∂x yc=∂W/∂y son las inclinaciones en xy en y, respectivamente, y d= (1/2) (∂2W/∂x2), e= (1/2) (∂2W/∂y2) y f=∂2W/∂x∂y est´an relacionados con el laplaciano y el hessiano de Wen Ω. A diferencia de la forma can´onica de Wo(Ec. (3.25)), aqu´ı aparece el par´ametro fporque, en general, los ejes del entorno Ω pueden tener cualquier orientaci´on. Para cada punto (R, θ) del borde suponemos que existe un entorno Ω = Ω (θ) donde el frente de onda re´une estas caracter´ısticas. Esto quiere decir que en cada ´angulo θ, suponemos que los movimientos del borde de la regi´on iluminada est´an regidos por las deformaciones que provocan las inclinaciones y las curvaturas locales del frente de onda en Ω. Es decir, en cada ´angulo θ, el borde de la regi´on iluminada evoluciona de la misma forma que el borde de la elipse asociada a la aproximaci´on cuadr´atica local del frente de onda en Ω. Para ilustrar esta hip´otesis emplearemos un ejemplo. Supongamos que en el plano de la pupila, z= 0, la fase es: W(x, y) = Ax3−3xy2(3.37) Queremos analizar el movimiento del borde de la zona iluminada en un punto (xc, yc) del contorno, desde el plano z= 0 al plano z. Las Ecs. (3.2) describen la evoluci´on de la posici´on del punto (xc, yc) de z= 0, a un punto (x, y) del plano zque, para este caso particular, son las siguientes:    x=xc+ 3Az(x2 c−y2 c) y=yc−6Azxcyc (3.38) Calculemos la fase local cuadr´atica Wode la Ec. (3.36) evaluada en puntos (xo, yo) de z= 0 dentro del entorno Ω asociado al punto (xc, yc). Para ello debemos obtener los valores de los par´ametros b,c,d,eyfde Woen el entorno Ω. En concreto: b= 3A(y2 o−x2 o) c= 6Axoyo d= 3Axo e=−3Axo f=−6Ayo (3.39) 3.5 Selecci´on de la separaci´on ´optima 73 −1.5 −1 −0.5 0 0.5 1 1.5 −1 −0.5 0 0.5 1 x y z=0.25 θ=π/ 3 −1.5 −1 −0.5 0 0.5 1 1.5 −1 −0.5 0 0.5 1 x y z=0.5 θ=π/ 3 −1.5 −1 −0.5 0 0.5 1 1.5 −1 −0.5 0 0.5 1 x y z=0.75 θ=π/ 3 −1.5 −1 −0.5 0 0.5 1 1.5 −1 −0.5 0 0.5 1 x y z=1 θ=π/ 3 Figura 3.8: Secuencia que muestra la evoluci´on del borde de la regi´on iluminada desde el plano z= 0 al plano z. En l´ınea delgada continua se representa la pupila en el plano z= 0. En l´ınea gruesa continua se muestra el borde de la regi´on iluminada en cada plano z, que corresponde a la fase W. Por ´ultimo, en l´ınea discontinua se muestra la evoluci´on del borde de la regi´on iluminada que predice la aproximaci´on cuadr´atica Woen la direcci´on θ=π/3. y las coordenadas de los puntos (xo, yo) del interior del entorno Ω, incluidos los del borde de la zona iluminada, evolucionar´an axialmente de acuerdo con:    x=xo+z(b+ 2dxo+fyo) y=yo+z(c+ 2exo+fxo) (3.40) En la Figura 3.8 se muestra la pupila circular en el plano z= 0 de radio R= 1, junto con la evoluci´on, en cuatro planos z, del borde de la regi´on iluminada correspondiente a Wy tambi´en a la aproximaci´on cuadr´atica local Woasociada al entorno Ω del punto con R= 1 y θ=π/3, tomando A= 1/(24R). En ella se observa la t´ıpica zona iluminada correspondiente al trifoil junto con la elipse correspondiente a la fase cuadr´atica, que evolucionan en contacto en la direcci´on θ=π/3. Esta secuencia confirma que el radio de la regi´on iluminada evoluciona en torno a cada ´angulo θen funci´on de las caracter´ısticas locales de la fase a primer y segundo orden. En la secci´on anterior, hemos representado curvas que describen el rendimiento del 74 Optimizaci´on de la medida de frontera sensor a la hora de recuperar la derivada normal de una fase cuadr´atica en la frontera, Wor, en funci´on de la posici´on de los planos en el eje z. Pero, en vista del razonamiento que acabamos de hacer, tambi´en nos informan de la precisi´on con la que el sensor estima la derivada normal local de un frente de onda cualquiera en un entorno Ω de un punto del borde. Ahora bien, los errores representados en la secci´on anterior para Woson la menor cota superior del error en cualquier ´angulo θy, por tanto, son independientes de la orientaci´on de la elipse del plano ±z. Entonces, podremos aplicar las varianzas de la secci´on anterior para todo entorno Ω independientemente de la orientaci´on de sus ejes, en funci´on de los par´ametros B,C,DyEde la forma can´onica de la fase cuadr´atica Wo. Otra consecuencia de que las curvas representadas sean cotas superiores de error, donde se verifica 2 (|D|+|E|) = Vy√B2+C2=U, es que aquellas curvas donde 2 (|D|+|E|)< V y√B2+C2< U se situar´ıan por debajo de las mostradas. Fij´emonos, entonces, en la curva asociada a un frente de onda con forma de cilindro parab´olico. Si llamamos zm´ın a la posici´on del m´ınimo de esta curva, si los planos de medida se sit´uan en las posiciones z < zm´ın, entonces las derivadas normales locales se estimar´an con menos precisi´on que si se mide en las posiciones z=zm´ın, cualquiera que sea el frente de onda cuadr´atico. Asimismo, si los planos se colocan en posiciones z > zm´ın, aunque mejora la estimaci´on en los entornos donde la aproximaci´on cuadr´atica Woes el´ıptica, parab´olica o hiperb´olica, empeora la recuperaci´on de la derivada normal en entornos donde Woes cil´ındrica. Siguiendo el mismo razonamiento que en las medidas del interior, si existe una densidad significativa de entornos donde el frente de onda tiene forma fundamentalmente cil´ındrica, entonces la recuperaci´on de la derivada normal empeorar´ıa tomando z > zm´ın. En definitiva, si de la fase solo conocemos las cotas UyV, minimizando la m´axima varianza ǫ2 opara un frente de onda con forma de cilindro parab´olico de curvatura Ve inclinaci´on U, obtendr´ıamos la mejor posici´on de los planos compatible con esta informaci´on. En la Figura 3.6 se representa la curva asociada a D= 1/2 y E= 0, con √B2+C2= 1, R= 1 y L= 0,3R, para las f´ormulas en diferencias finitas del m´etodo indirecto, de donde se deduce que la f´ormula central proporciona la mejor respuesta del sensor. Es importante indicar que el error se mantiene en niveles bajos en un rango de zm´as amplio para la f´ormula central que para las f´ormulas de Newton. Esta es una caracter´ıstica muy a tener en cuenta, pues la elecci´on de la posici´on de los planos se muestra m´as flexible para la f´ormula central que para las f´ormulas de Newton. 3.5 Selecci´on de la separaci´on ´optima 75 A diferencia de las medidas del interior, para optimizar la recuperaci´on de Wr con el m´etodo indirecto, hay que escoger adicionalmente una dimensi´on radial de la regi´on de frontera, L, que conjuntamente con la selecci´on apropiada de zminimice la varianza de ˆ Wr. Es decir, para cada valor de Lexiste uno de z,zm´ın, que minimiza la varianza de ˆ Wr. Pero para seleccionar la separaci´on ´optima, zopt, debemos hallar el par (L, zm´ın) que garantiza la mejor respuesta del sensor. En concreto, se trata de ajustar La la deformaci´on m´axima de la regi´on iluminada en el plano de medida, evitando la p´erdida de informaci´on que se produce cuando el borde de la regi´on iluminada sobrepasa los l´ımites interior o exterior del anillo de detectores de la regi´on de frontera. Solo as´ı podremos decir que zopt =zm´ın. En la Figura 3.9 se representa ǫo(z) en la configuraci´on central y en la de Newton, para R= 1, √B2+C2= 1, D= 1/2 y E= 0, y para varias desviaciones de las medidas de irradiancia, σ/Io, cuando Lse ha ajustado a la deformaci´on radial m´axima de la regi´on iluminada en zopt. El proceso que hemos seguido es el siguiente: Dado un valor de σ/Io, simplemente hemos variado Ly, en paralelo, zm´ın, hasta que Lcoincide con el valor m´aximo de R+ o−R− oevaluado en zm´ın, en un proceso iterativo. Por ejemplo, supongamos que σ/Io= 0,01. Partiendo de L= 1 se obtiene zm´ın = 0,056. Pero situando los planos en esa posici´on, se obtiene que R+ o−R− o= 0,22 6= 1, lo que indica que la dimensi´on radial de la regi´on de frontera no est´a ajustada correctamente. Si ahora tomamos L= 0,22, entonces zm´ın = 0,041 y R+ o−R− o= 0,16 6= 0,22. Como L sigue siendo todav´ıa distinto a R+ o−R− o, tomamos L= 0,16 y calculamos de nuevo zm´ın. El proceso contin´ua iterativamente hasta que el valor de Lpreseleccionado coincide con R+ o−R− o. En este ejemplo esta situaci´on se produce cuando L= 0,14 y podremos decir que la posici´on del m´ınimo error, zm´ın, nos proporciona el plano ´optimo, es decir, zopt =zm´ın = 0,036. Comparando con la Figura 3.5, se aprecia que la magnitud del error en el plano ´optimo disminuye apreciablemente cuando Lse ajusta a la deformaci´on radial m´axima de la regi´on iluminada. No obstante, habitualmente no dispondremos de la informaci´on manejada hasta aqu´ı acerca de los par´ametros B,C,DyEen Ω, sino que en su lugar tendremos acceso a cierta informaci´on a priori acerca de la estad´ıstica de los frentes de onda que llegan al sensor, de utilidad para proceder a seleccionar los planos. An´alogamente al procedimiento desarrollado para las medidas del interior, supongamos que nuestro conocimiento previo de la estad´ıstica de los frentes de onda en los puntos del plano z= 0, nos permite: estimar V=|h∇2W(r)i|+αq[∇2W(r)−h∇2W(r)i]2, 82 Sensor de curvatura multiplano costa de incrementar enormemente las necesidades computacionales [23, 121, 143– 145]. Tambi´en se han aplicado algorimos no lineales iterativos para refinar la fase inicialmente recuperada por medio de la ETI a partir de medidas de irradiancia en dos planos [16, 90, 98, 134, 141]. Como una aportaci´on m´as a la t´ecnica de recuperaci´on de fase basada en la ETI, proponemos en este cap´ıtulo aplicar f´ormulas en diferencias finitas m´as precisas para aproximar la derivada axial. Las configuraciones del sensor que se derivan de estas f´ormulas utilizan indefectiblemente medidas en m´as de dos planos de irradiancia. Escogiendo adecuadamente la posici´on y el peso de los datos de entrada de la f´ormula, se minimiza la contribuci´on del ruido de medida y se reduce el error de no linealidad de la aproximaci´on, sin que ello implique necesariamente un incremento sustancial del gasto computacional al reconstruir la fase. 4.2. Varianza de la estimaci´on de la derivada axial con f´ormulas en m´ultiples planos En el Cap´ıtulo 2, calculamos la varianza de la respuesta del sensor para las dos aproximaciones cl´asicas de la derivada en diferencias finitas con dos planos de irradiancia: la f´ormula central o de Stirling y las f´ormulas de Newton progresiva y regresiva. El error de no linealidad de la aproximaci´on en diferencias finitas se dedujo a partir del resto de la f´ormula de Taylor de la irradiancia en torno al plano de medida, z= 0, truncado a primer orden. Por otro lado, el efecto del ruido se model´o suponiendo que las fluctuaciones de irradiancia en cada punto son de media cero e independientes del nivel de se˜nal. Siguiendo un procedimiento an´alogo, en esta secci´on calcularemos la varianza de estas dos aproximaciones de la derivada axial en diferencias finitas generalizadas a M≥2 planos de irradiancia en funci´on de la separaci´on interplanos, z. 4.2.1. F´ormula central de Stirling En un instante de tiempo t, la intensidad en cada punto r= (x, y) de un plano jz,ij=i(r, jz), se expresa como la suma de la intensidad exacta sin ruido o determin´ıstica, Ij=I(r, jz), y el ruido, nj=nj(r, jz), con j∈Z+: ij=Ij+nj(4.1) 4.2 Varianza de la estimaci´on de la derivada axial con f´ormulas en m´ultiples planos 83 donde njson errores aditivos con media cero, varianza σ2en todos los puntos de cada plano jz y correlaci´on nula entre puntos distintos. Se quiere estimar I′(0) = I′(r,0) y para ello utilizaremos la aproximaci´on central en diferencias finitas usando M= 2m datos de entrada. Se define: ˆ I′(0) = 1 z m X j=−m ajij(4.2) donde m∈Z+,z > 0 y el circunflejo sobre ˆ I′denota una estimaci´on de I′. Para cada punto rse quieren encontrar los pesos ajque verifican las dos condiciones siguientes: l´ım z→0Enˆ I′(0)o=I′(0), s2hˆ I′(0)i=Ehˆ I′(0) −I′(0)i2= m´ın (4.3) donde E {·} denota valor medio. Para cada punto ry en torno al plano z= 0, se verifica: ˆ I′(0) = 1 z m X j=−m ajnj+1 z m X j=−m ajIo+z z m X j=−m ajjI′(0) + z2 z m X j=−m aj j2 2I′′(ζj) (4.4) con ζj∈(0, jz)∀j, que se obtiene sustituyendo el desarrollo de Taylor de Ijhasta primer orden (Ec. (2.5)) en la aproximaci´on (4.2). Entonces, en cada punto rse verifica: ˆ I′(0) − m X j=−m ajjI′(0) = 1 z m X j=−m ajnj+1 z m X j=−m ajIo+z m X j=−m aj j2 2I′′(ζj) (4.5) Haciendo: m P j=−m ajj= 1 m P j=−m aj= 0 (4.6) la diferencia entre la derivada exacta y la estimada es: ˆ I′(0) −I′(0) = 1 z m X j=−m ajnj+z m X j=−m aj j2 2I′′(ζj) (4.7) Es decir, la aproximaci´on en diferencias finitas se aleja del valor real de la derivada axial una cantidad que depende del ruido de los datos de entrada y del resto del desarrollo en serie de Taylor de la irradiancia truncado a primer orden, estando ambas contribuciones pesadas por el correspondiente coeficiente ajasociado a cada plano jz y satisfaci´endose la primera condici´on de la Ec. (4.3). La desviaci´on cuadr´atica media es: s2hˆ I′(0)i=σ2 z2 m X j=−m a2 j!+z2"m X j=−m aj j2 2I′′(ζj)#2 (4.8) 84 Sensor de curvatura multiplano M m a±1a±2a±3a±4a±5a±6 2 1 ±1/2 4 2 ±1/5±1/10 6 3 ±3/28 ±1/14 ±1/28 8 4 ±1/15 ±1/20 ±1/30 ±1/60 10 5 ±1/22 ±2/55 ±3/110 ±1/55 ±1/110 12 6 ±3/91 ±5/182 ±2/91 ±3/182 ±1/91 ±1/182 Tabla 4.1: Coeficientes ajcorrespondientes a la f´ormula central en funci´on del n´umero de planos. donde se han tenido en cuenta las hip´otesis asumidas para el ruido, junto con la independencia entre la contribuci´on de los errores de los datos y la del error de truncamiento. Para que la varianza de ˆ I′(0) satisfaga la segunda condici´on de la Ec. (4.3), necesariamente ha de escogerse el peso de los datos de entrada de forma que se minimice la contribuci´on del ruido de medida, es decir, han de escogerse los coeficientes ajtales que: m X j=−m a2 j= m´ın (4.9) Aplicando el m´etodo de los multiplicadores de Lagrange para resolver la Ec. (4.9) con las condiciones (4.6) [159], se obtiene: aj=3j m(m+ 1)(2m+ 1) (4.10) En la Tabla 4.1, se indican algunos de estos coeficientes para el n´umero de planos Mindicado. El c´alculo de la suma de la Ec. (4.9) es: m X j=−m a2 j=3 m(m+ 1)(2m+ 1) = ˜am(4.11) Por tanto, la expresi´on general para la varianza de la aproximaci´on central de la derivada (4.2) es: s2hˆ I′(0)i=σ2 z2˜am+z2"m X j=−m aj j2 2I′′(ζj)#2 (4.12) 4.2 Varianza de la estimaci´on de la derivada axial con f´ormulas en m´ultiples planos 85 4.2.2. F´ormula progresiva y regresiva de Newton En este caso la derivada axial se aproxima en r= (x, y) como diferencias finitas de medidas de la irradiancia en el plano z= 0 y en M−1 = mplanos equiespaciados situados, o bien a la derecha (diferenciaci´on progresiva), o bien a la izquierda de z= 0 (diferenciaci´on regresiva). De este modo, definimos una estimaci´on de I′(0) = I′(r,0) como: ˆ I′(0) = 1 ±z m X j=1 bj(ij−io) (4.13) donde m∈Z+,z > 0 y el denominador es +zpara la aproximaci´on progresiva y −z en el caso regresivo. El ruido que afecta a los datos de irradiancia (Ec. (4.1)), re´une las mismas caracter´ısticas que en la secci´on anterior. Para cada punto se quieren hallar los coeficientes bjtales que la aproximaci´on (4.13) verifica las Ecs. (4.3). Como antes, sustituimos en la Ec. (4.13) la f´ormula de Taylor para Ijen torno az= 0 y obtenemos la diferencia entre la derivada exacta y la estimada. En esta ocasi´on, para que se verifiquen las condiciones (4.3), los coeficientes bjsolo han de cumplir: m P j=1 bjj= 1 m P j=1 b2 j= m´ın (4.14) As´ı, despejamos la diferencia entre la derivada exacta y la estimada en una ecuaci´on an´aloga a la Ec. (4.7). Considerando la hip´otesis de que la contribuci´on del ruido de los datos de irradiancia y la del error de truncamiento son independientes, llegamos a una ecuaci´on para la varianza de la aproximaci´on (4.13), an´aloga a la Ec. (4.8): s2hˆ I′(0)i=2σ2 z2 m X j=1 b2 j!+z2"m X j=1 bj j2 2I′′(ζj)#2 (4.15) con ζj∈(0, jz)∀jpara la aproximaci´on progresiva y con ζj∈(−jz, 0) ∀jpara la regresiva. Para calcular los coeficientes bj, imponemos las condiciones (4.14) para bj y hallamos el correspondiente multiplicador de Lagrange, obteni´endose: bj=6j m(m+ 1)(2m+ 1) (4.16) En la Tabla 4.2, se indican algunos de los coeficientes bjasociados al n´umero de planos Mindicado. Con la Ec. (4.16) ya podemos calcular: m X j=1 b2 j=6 m(m+ 1)(2m+ 1) =˜ bm(4.17) 86 Sensor de curvatura multiplano M m b1b2b3b4b5b6 2 1 1 3 2 1/5 2/5 4 3 1/14 1/7 3/14 5 4 1/30 1/15 1/10 2/15 6 5 1/55 2/55 3/55 4/55 1/11 7 6 1/91 2/91 3/91 4/91 5/91 6/91 Tabla 4.2: Coeficientes bjcorrespondientes a la f´ormula progresiva y regresiva en funci´on del n´umero de planos. es decir, ˜ bm= 2˜am. Sustituyendo la Ec. (4.17) en la Ec. (4.15), se tiene la expresi´on general para la varianza de las aproximaciones progresiva y regresiva de la derivada en un punto r: s2hˆ I′(0)i=2σ2 z2˜ bm+z2"m X j=1 bj j2 2I′′(ζj)#2 (4.18) 4.3. Selecci´on del plano de medida Las Ecs. (4.12) y (4.18) describen, para cada punto r, las varianzas de ambas f´ormulas en diferencias finitas con m´ultiples datos de entrada en funci´on de la distancia interplanos, z. Las varianzas calculadas en el Cap´ıtulo 2 para f´ormulas en dos planos son casos particulares de estas ecuaciones. Esencialmente, las f´ormulas aqu´ı calculadas presentan caracter´ısticas muy similares a aquellas. El primer sumando describe la contribuci´on de los errores de los datos de entrada y es inversamente proporcional a z2. El segundo sumando es directamente proporcional az2y da cuenta del error de no linealidad de las f´ormulas en funci´on del resto del desarrollo en serie de Taylor de la irradiancia en torno al plano z= 0. La precisi´on de la f´ormula en diferencias finitas y, por tanto, la calidad de respuesta del sensor de curvatura depende de la contribuci´on de estos dos sumandos. Los coeficientes del primero, ˜amy˜ bm, se han escogido para minimizar la propagaci´on del ruido de los datos de irradiancia. En la Tabla 4.3 se indican algunos de ellos para distintos valores de m. An´alogamente a lo visto en las f´ormulas en dos planos, el segundo sumando es m´as complejo. El punto ζj∈(0, jz) es desconocido y, en general, distinto para cada punto rde cada plano jz. Esto nos impide la aplicaci´on directa 4.3 Selecci´on del plano de medida 87 m˜am 1 1/2 2 1/10 3 1/28 4 1/60 5 1/110 25 1/11050 50 1/676700 Tabla 4.3: Coeficientes ˜ampara distintos valores de m. de las ecuaciones halladas para calcular una separaci´on interplanos que optimice la respuesta del sensor globalmente en el plano z= 0. Siguiendo un procedimiento an´alogo al de los Cap´ıtulos 2 y 3, hallaremos en esta secci´on expresiones obtenidas a partir de las Ecs. (4.12) y (4.18) para optimizar el funcionamiento de un sensor que utilice m´ultiples planos. 4.3.1. Varianza de la estimaci´on de la curvatura local y separaci´on ´optima Seg´un lo visto en la §2.3.1, en la integral de difracci´on de Fresnel contribuyen fundamentalmente regiones en la vecindad de ciertos puntos de la pupila de entrada, rcj. De este modo, la irradiancia Ij=I(r, jz) en un punto revoluciona en cada plano jz esencialmente en funci´on de las caracter´ısticas de la fase en torno a dichos puntos, rcj, del plano z= 0 (Ec. (2.12)). Pero si trabajamos en entornos Ω del frente de onda donde este se comporta fundamentalmente como su recubrimiento cuadr´atico Wo(Ec. (2.14) o su equivalente en forma can´onica, Ec. (2.15)), toda dependencia en rcjimpl´ıcita en la Ec. (2.12) desaparece. As´ı, en todos los puntos rdel interior de Ω, la irradiancia evoluciona en cada plano jz de la forma: Ij=Io 1 + 2 (D+E)jz + 4DEj2z2(4.19) donde DyEson las curvaturas principales del frente en el entorno Ω y donde los planos jz est´an situados antes de las c´austicas. Haciendo uso de la ecuaci´on anterior, podemos calcular varianzas de las f´ormulas en diferencias finitas con m´ultiples planos v´alidas en puntos rdel interior de Ω, de 88 Sensor de curvatura multiplano forma an´aloga al c´alculo realizado para las f´ormulas en dos planos del Cap´ıtulo 2, y sustituirlas en la siguiente ecuaci´on para obtener la varianza de la curvatura local en Ω, ǫ2 o: ǫ2 o=1 I2 on4 (D+E)2σ2+s2hˆ I′(0)io (4.20) Hay que indicar aqu´ı que, mientras que con dos planos de irradiancia se consider´o que la varianza de Io,σ2, es despreciable comparado con s2hˆ I′(0)i, cuando se utilizan m´as de dos planos podr´ıa darse el caso de que la varianza de la derivada axial disminuya hasta llegar a alcanzar valores del orden de σ2. De ah´ı el no despreciarla en este caso. Entonces, las varianzas asociadas a cada f´ormula son: a) F´ormula central A partir de la Ec. (4.7), calculamos la contribuci´on del error de truncamiento de la aproximaci´on central en diferencias finitas: z"m X j=−m aj j2 2I′′(ζj)#=1 z m X j=−m ajIj−I′(0) (4.21) De este modo, la Ec. (4.12) en funci´on de los par´ametros DyEdel entorno Ω es: s2hˆ I′(0)i=σ2˜am z2+I2 o"2 (D+E) + 1 z m X j=−m aj 1 + 2 (D+E)jz + 4DEj2z2#2 (4.22) Teniendo en cuenta la Ec. (4.10), a−j=−aj, y podemos simplificar el segundo t´ermino de la ecuaci´on anterior y sustituirla en la Ec. (4.20), obteni´endose la siguiente expresi´on para la varianza de la curvatura local estimada con la f´ormula central en el entorno Ω: ǫ2 o=4 (D+E)2+˜am z2σ Io2 + 4 (D+E)2 ×"1−˜am m X j=−m j2 (1 + 4DEj2z2)2−4 (D+E)2j2z2#2 (4.23) b) F´ormula progresiva y regresiva De modo an´alogo al apartado anterior, calculamos el segundo t´ermino de la Ec. (4.18) y obtenemos la varianza de esta f´ormula en funci´on de los par´ametros 4.3 Selecci´on del plano de medida 89 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0 0.012 0.024 0.036 0.048 0.06 0.072 0.084 0.096 0.108 0.12 z ǫo (a) D= 1/2, E = 0 2 (D+E) = 1 D= 1/3, E = 0 2 (D+E) = 2/3 D= 1/6, E = 0 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0 0.012 0.024 0.036 0.048 0.06 0.072 0.084 0.096 0.108 0.12 z ǫo (b) D= 1/3, E = 1/6 2 (D+E) = 1 D= 2/9, E = 1/9 2 (D+E) = 2/3 D= 1/9, E = 1/18 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0 0.012 0.024 0.036 0.048 0.06 0.072 0.084 0.096 0.108 0.12 z ǫo (c) D= 1/4, E = 1/4 2 (D+E) = 1 D= 1/6, E = 1/6 2 (D+E) = 2/3 D= 1/12, E = 1/12 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0 0.012 0.024 0.036 0.048 0.06 0.072 0.084 0.096 0.108 0.12 z ǫo (d) D= 1/3, E =−1/6 2 (|D|+|E|) = 1 D= 2/9, E =−1/9 2 (|D|+|E|) = 2/3 D= 1/9, E =−1/18 2 (|D|+|E|) = 1/3 D= 0, E = 0 2 (|D|+|E|) = 0 Figura 4.1: Desviaci´on est´andar de la curvatura local del frente de onda en un entorno Ω en funci´on de la posici´on de los planos, z, suponiendo que se estima con la f´ormula central de Stirling particularizada a M= 8 planos de medida, siendo σ/Io= 0,01. El apartado (a) corresponde a aproximaciones cuadr´aticas con forma de cilindro parab´olico, el (b) a aproximaciones con forma de paraboloide el´ıptico, el (c) a aproximaciones con forma de paraboloide parab´olico y el (d) a aproximaciones de tipo paraboloide hiperb´olico. DyEcaracter´ısticos de la aproximaci´on cuadr´atica del frente de onda en el entorno Ω. Por ´ultimo, sustituyendo en la Ec. (4.20), calculamos la varianza de la curvatura local de Westimada en el entorno Ω con la f´ormula progresiva y regresiva: ǫ2 o="4 (D+E)2+ 2˜ bm z2#σ Io2 +"2 (D+E)±1 z m X j=1 bj∓2 (D+E)jz −4DEj2z2 1±2 (D+E)jz + 4DEj2z2#2 (4.24) En las Figuras 4.1, 4.2 y 4.3, se representan las curvas de error de la curvatura local estimada en el entorno Ω con las f´ormulas central, progresiva y regresiva, 90 Sensor de curvatura multiplano 0 0.010.020.03 0.040.050.060.070.08 0.09 0.1 0 0.024 0.048 0.072 0.096 0.12 0.144 0.168 0.192 0.216 0.24 z ǫo (a) D= 1/2, E = 0 2 (D+E) = 1 D= 1/3, E = 0 2 (D+E) = 2/3 D= 1/6, E = 0 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0 0.010.020.03 0.040.05 0.060.07 0.080.09 0.1 0 0.024 0.048 0.072 0.096 0.12 0.144 0.168 0.192 0.216 0.24 z ǫo (b) D= 1/3, E = 1/6 2 (D+E) = 1 D= 2/9, E = 1/9 2 (D+E) = 2/3 D= 1/9, E = 1/18 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0 0.010.020.03 0.040.050.060.070.08 0.09 0.1 0 0.024 0.048 0.072 0.096 0.12 0.144 0.168 0.192 0.216 0.24 z ǫo (c) D= 1/4, E = 1/4 2 (D+E) = 1 D= 1/6, E = 1/6 2 (D+E) = 2/3 D= 1/12, E = 1/12 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0 0.010.020.03 0.040.05 0.060.07 0.080.09 0.1 0 0.024 0.048 0.072 0.096 0.12 0.144 0.168 0.192 0.216 0.24 z ǫo (d) D= 1/3, E =−1/6 2 (|D|+|E|) = 1 D= 2/9, E =−1/9 2 (|D|+|E|) = 2/3 D= 1/9, E =−1/18 2 (|D|+|E|) = 1/3 D= 0, E = 0 2 (|D|+|E|) = 0 Figura 4.2: Desviaci´on est´andar de la curvatura local del frente de onda en un entorno Ω en funci´on de la posici´on de los planos, z, suponiendo que se estima con la f´ormula progresiva de Newton particularizada a M= 8 planos de medida, siendo σ/Io= 0,01. El apartado (a) corresponde a aproximaciones cuadr´aticas con forma de cilindro parab´olico, el (b) a aproximaciones con forma de paraboloide el´ıptico, el (c) a aproximaciones con forma de paraboloide parab´olico y el (d) a aproximaciones de tipo paraboloide hiperb´olico. respectivamente, con M= 8 planos de irradiancia*. Suponemos que las irradiancias proporcionadas por el detector fluct´uan con desviaci´on relativa σ/Io= 0,01. Los valores de los par´ametros DyEse han seleccionado manteniendo la condici´on 2(|D|+|E|)≤V= 1, al igual que en las figuras de la §2.3.2 asociadas a f´ormulas con datos de irradiancia en solo dos planos. El apartado (a) de cada figura muestra el comportamiento de ǫopara frentes de onda con forma de cilindro parab´olico en Ω, el (b) para frentes con forma de paraboloide el´ıptico, el (c) para frentes con forma de paraboloide esf´erico y el (d) para aquellos con forma de paraboloide hiperb´olico. De nuevo, tal y como observamos en las curvas para dos planos, de todas las *Podr´ıa utilizarse otro n´umero de planos y se llegar´ıa a las mismas conclusiones, pues el comportamiento cualitativo general de las curvas en funci´on de zno cambia por escoger un n´umero de planos distinto en la f´ormula en diferencias finitas. 4.3 Selecci´on del plano de medida 91 0 0.010.020.03 0.040.050.060.070.08 0.09 0.1 0 0.024 0.048 0.072 0.096 0.12 0.144 0.168 0.192 0.216 0.24 z ǫo (a) D= 1/2, E = 0 2 (D+E) = 1 D= 1/3, E = 0 2 (D+E) = 2/3 D= 1/6, E = 0 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0 0.010.020.03 0.040.05 0.060.07 0.080.09 0.1 0 0.024 0.048 0.072 0.096 0.12 0.144 0.168 0.192 0.216 0.24 z ǫo (b) D= 1/3, E = 1/6 2 (D+E) = 1 D= 2/9, E = 1/9 2 (D+E) = 2/3 D= 1/9, E = 1/18 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0 0.010.020.03 0.040.050.060.070.08 0.09 0.1 0 0.024 0.048 0.072 0.096 0.12 0.144 0.168 0.192 0.216 0.24 z ǫo (c) D= 1/4, E = 1/4 2 (D+E) = 1 D= 1/6, E = 1/6 2 (D+E) = 2/3 D= 1/12, E = 1/12 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0 0.010.020.03 0.040.05 0.060.07 0.080.09 0.1 0 0.024 0.048 0.072 0.096 0.12 0.144 0.168 0.192 0.216 0.24 z ǫo (d) D= 1/3, E =−1/6 2 (|D|+|E|) = 1 D= 2/9, E =−1/9 2 (|D|+|E|) = 2/3 D= 1/9, E =−1/18 2 (|D|+|E|) = 1/3 D= 0, E = 0 2 (|D|+|E|) = 0 Figura 4.3: Desviaci´on est´andar de la curvatura local del frente de onda en un entorno Ω en funci´on de la posici´on de los planos, z, suponiendo que se estima con la f´ormula regresiva de Newton particularizada a M= 8 planos de medida, siendo σ/Io= 0,01. El apartado (a) corresponde a aproximaciones cuadr´aticas con forma de cilindro parab´olico, el (b) a aproximaciones con forma de paraboloide el´ıptico, el (c) a aproximaciones con forma de paraboloide parab´olico y el (d) a aproximaciones de tipo paraboloide hiperb´olico. combinaciones de DyEtales que 2(|D|+|E|) = V= 1 (curvas en l´ınea continua de las figuras) la curva que est´a por encima de todas las dem´as corresponde al cilindro parab´olico, seguida, en orden decreciente, por la curva del paraboloide parab´olico, la del paraboloide el´ıptico y, por ´ultimo, la del hiperb´olico. Si llamamos zopt a la separaci´on interplanos que coincide con el m´ınimo de la curva con V= 1 del cilindro parab´olico, si los planos de medida se posicionan separados entre s´ı una distancia z < zopt, entonces las curvaturas locales estimadas en cualquier entorno Ω del dominio donde haya otro tipo de cu´adrica tendr´an siempre m´as error que si se separan z= zopt. Por otro lado, si las irradiancias se detectan en planos separados z > zopt, aunque mejora la estimaci´on de las curvaturas locales en los entornos donde la aproximaci´on cuadr´atica es un paraboloide o un plano, empeora en entornos donde el frente de onda puede aproximarse por cilindros parab´olicos. Consecuentemente, la distancia 98 Sensor de curvatura multiplano b) F´ormula progresiva y regresiva Calculamos la menor cota superior tomando D=V/2 y E= 0 en la Ec. (4.24) y evaluando en h=V z: ǫsup V2= 1 + 2˜ bm h2!σ Io2 +"1− m X j=1 bjj 1−jh#2 (4.37) La mayor cota inferior se obtiene sustituyendo D=E= 0: ǫinf V2=2˜ bm h2σ Io2 (4.38) y minimizando la Ec. (4.37) respecto de he igualando a cero, obtenemos una relaci´on entre la posici´on ´optima de los planos, hopt, y el nivel de ruido de medida, σ/Io: 2˜ bm h3 opt σ Io2 −"1− m X j=1 bjj 1−jhopt #"m X j=1 bjj 1−jhopt 2#= 0 (4.39) En la Figura 4.5(b) se muestran estas cotas para los mismos valores de σ/Io que en la Figura 4.5(a). Al igual que ocurr´ıa con las f´ormulas con dos planos, el error en la separaci´on ´optima es mayor en las curvas correspondientes a las f´ormulas de Newton, como se deduce comparando aquellos pares de curvas de ambas figuras asociados al mismo nivel de ruido. En definitiva, independientemente del n´umero de planos, los errores de los datos de entrada y el error de modelizaci´on se propagan m´as a la medida de curvatura si se usan las f´ormulas de Newton en lugar de la f´ormula central. Por otro lado, la Figura 4.5, para f´ormulas en M= 8 planos, junto con la Figura 2.7, para f´ormulas en M= 2 planos, nos permiten comparar el rendimiento del sensor en relaci´on al n´umero de planos de irradiancia que utiliza para estimar la curvatura local. Por ejemplo, en las curvas con σ/Io= 0,1 para la f´ormula central, se aprecia un descenso de la cota superior de error en la separaci´on ´optima desde ǫsup/V = 0,24 con dos planos hasta ǫsup/V = 0,20 con ocho planos. En las curvas con σ/Io= 0,01, el error disminuye desde ǫsup/V = 0,051 con dos planos hasta ǫsup/V = 0,039 con ocho planos. Para las f´ormulas de Newton, si σ/Io= 0,1, el error se reduce desde ǫsup/V = 0,64 con dos planos hasta ǫsup/V = 0,43 con ocho planos. Si σ/Io= 0,01, la reducci´on va desde ǫsup/V = 0,18 con dos planos hasta ǫsup/V = 0,12 con ocho planos. 4.4 Rendimiento del sensor de curvatura multiplano 99 0 4 8 12 16 20 24 28 32 36 40 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055 (a) M ǫo(hopt) V 0 4 8 12 16 20 24 28 32 36 40 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 (b) M ǫo(hopt) V ǫsup ǫinf ǫsup ǫinf Figura 4.6: Evoluci´on de la cota superior de error de la curvatura local en funci´on del n´umero de planos de medida, suponiendo que se estima con la f´ormula central en diferencias finitas, (a), y con la f´ormula progresiva y regresiva de Newton, (b). En ambos casos σ/Io= 0,01. La evoluci´on de las cotas superior e inferior de error en hopt en funci´on del n´umero de planos puede verse m´as claramente en la Figura 4.6. En esta representaci´on se muestra c´omo a consecuencia del uso de m´as planos de irradiancia el error de la estimaci´on de la curvatura local podr´ıa hacerse tan peque˜no como para llegar a ser comparable a σ. El empleo de f´ormulas con m´as planos de irradiancia tambi´en trae consigo la reducci´on de la distancia interplanos a medida que aumenta M(ver Figura 4.7). Esta situaci´on pone un l´ımite superior en el n´umero de medidas de irradiancia, pues el valor de zpodr´ıa reducirse tanto que la imprecisi´on de la posici´on de los planos, que en este trabajo no se est´a teniendo en cuenta, podr´ıa llegar a ser significativa. 100 Sensor de curvatura multiplano 0 4 8 12 16 20 24 28 32 36 40 0 0.017 0.034 0.051 0.068 0.085 0.102 0.119 0.136 0.153 0.17 (a) M hopt 0 4 8 12 16 20 24 28 32 36 40 0 0.011 0.022 0.033 0.044 0.055 0.066 0.077 0.088 0.099 0.11 (b) M hopt Figura 4.7: Evoluci´on de la separaci´on interplanos conforme crece el n´umero de planos de medida, para la f´ormula central en diferencias finitas, (a), y para la f´ormula progresiva y regresiva de Newton, (b), suponiendo σ/Io= 0,01. Adicionalmente, si zllega a ser suficientemente peque˜no, la hip´otesis de no correlaci´on entre las medidas de irradiancia en un plano y los adyacentes podr´ıa dejar de ser consistente. La Figura 4.8 es una representaci´on que explica c´omo el plano m´as alejado de la pupila est´a cada vez m´as cerca de z= 0 conforme aumenta el n´umero de planos. Es decir, la disminuci´on de hopt a consecuencia del aumento del n´umero de planos, va pareja a una reducci´on del rango axial donde se toman medidas de irradiancia. A ra´ız de estos resultados, podemos concluir que en experimentos donde la irradiancia solo est´a afectada por ruido de media cero independiente del nivel de se˜nal, la estimaci´on de la curvatura es m´as precisa si se emplean f´ormulas con m´as 4.4 Rendimiento del sensor de curvatura multiplano 101 0 4 8 12 16 20 24 28 32 36 40 0.14 0.144 0.148 0.152 0.156 0.16 0.164 0.168 0.172 0.176 0.18 (a) M M 2hopt 0 4 8 12 16 20 24 28 32 36 40 0.065 0.07 0.075 0.08 0.085 0.09 0.095 0.1 0.105 0.11 0.115 (b) M (M−1)hopt Figura 4.8: Evoluci´on de la posici´on de los planos m´as alejados de z= 0 conforme crece el n´umero de planos de medida, para la f´ormula central en diferencias finitas, (a), y para la f´ormula progresiva y regresiva de Newton, (b), suponiendo σ/Io= 0,01. de dos planos. A cambio de tener que sustituir el detector por otro m´as avanzado, mejorar el rendimiento del sensor exige llevar a cabo medidas en m´as de dos planos para ganar resoluci´on. Podemos verlo en la Figura 4.9. Esta es una representaci´on de la calidad del detector que deber´ıamos utilizar en un sensor de curvatura que trabaje con solo dos planos, para que en la posici´on ´optima estimemos las curvaturas locales con la misma precisi´on que un sensor multiplano optimizado que utilice un detector que proporcione medidas de irradiancia con desviaci´on est´andar σ/Io= 0,01. 102 Sensor de curvatura multiplano 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 3 4 5 6 7 8 9 10 x 10−3 M σ Io Figura 4.9: Nivel de ruido de medida, σ/Io, asumible en un sensor que funciona con dos planos de medida, para garantizar la recuperaci´on de una curvatura local de magnitud Ven el esquema central con la misma precisi´on que un sensor de curvatura multiplano que utiliza un detector afectado por ruido de detecci´on con desviaci´on relativa σ/Io= 0,01. 4.4.2. Calidad de la derivada normal recuperada en funci´on del n´umero de planos Si la aproximaci´on cuadr´atica del frente de onda en cada entorno Ω verifica que √B2+C2≤Uy 2 (|D|+|E|)≤V, en la §4.3.2 se ha visto que la separaci´on interplanos est´a fijada por el m´ınimo de la varianza asociada a un frente cil´ındrico de curvatura m´axima, D=V/2 y E= 0, y con la m´axima inclinaci´on, √B2+C2=U. Pero para conseguir la mejor respuesta del sensor, adicionalmente debemos adaptar el tama˜no de la regi´on de frontera, L, para que coincida con la deformaci´on radial m´axima de la regi´on iluminada. As´ı se minimiza la contribuci´on del ruido de las medidas de irradiancia evitando, al mismo tiempo, la p´erdida de informaci´on que implica que el borde de la regi´on iluminada sobrepase los l´ımites de los detectores que constituyen la regi´on de frontera. 4.4 Rendimiento del sensor de curvatura multiplano 103 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z ǫo (a) 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z ǫo (b) σ / Io= 0.1 ; L= 0.47R;zopt = 0.033 σ / Io= 0.03 ; L= 0.23R;zopt = 0.017 σ / Io= 0.01 ; L= 0.11R;zopt = 0.0077 σ / Io= 0.1 ; L= 0.25R;zopt = 0.0089 σ / Io= 0.03 ; L= 0.074R;zopt = 0.0026 σ / Io= 0.01 ; L= 0.023R;zopt = 0.00082 Figura 4.10: Curvas de error de la derivada normal de la fase en la frontera estimada con la configuraci´on central (a) y la de Newton (b) particularizadas a M= 8 planos de medida, cuando el tama˜no de la zona de detecci´on, L, se ajusta a la deformaci´on radial m´axima de la regi´on iluminada, suponiendo el nivel de ruido en el plano de recuperaci´on, σ/Io, indicado. Adicionalmente, se incluye la separaci´on ´optima, zopt. Tomemos entonces un frente de onda cil´ındrico con U= 1 y V= 1. En las Figuras 4.10(a), para la f´ormula central, y 4.10(b), para la f´ormula de Newton, se representa ǫo(z) cuando Lse ha ajustado a la deformaci´on radial m´axima de la regi´on iluminada en zopt para los valores de σ/Ioindicados. Comparando con la figura correspondiente para dos planos, Figura 3.9, se aprecia una disminuci´on del error de la estimaci´on de Woren zopt. En el caso de la f´ormula central, la disminuci´on va desde ǫo= 0,36 con dos planos a ǫo= 0,27 con ocho planos, para σ/Io= 0,1, y desde ǫo= 0,16 con 104 Sensor de curvatura multiplano dos planos a 0,12 con ocho planos, para σ/Io= 0,01. En el caso de la f´ormula de Newton, el error disminuye desde 0,49 con dos planos a ǫo= 0,3 con ocho planos, para σ/Io= 0,1, y desde ǫo= 0,15 con dos planos a ǫo= 0,09 con ocho planos, para σ/Io= 0,01. Al igual que con el laplaciano, el empleo de m´as datos de entrada tambi´en reduce el error de la recuperaci´on de Wor. En definitiva, en experimentos donde las medidas de irradiancia est´an afectadas solo por ruido de media cero e independiente del nivel de se˜nal, la estimaci´on de la derivada normal es m´as precisa si aumenta el n´umero de datos de irradiancia. Por ´ultimo, an´alogamente a lo observado con las f´ormulas en dos planos, ǫo permanece en valores cercanos al m´ınimo en un rango axial mayor en la f´ormula central que en las f´ormulas de Newton. Por tanto, la elecci´on de zopt con la f´ormula central se muestra m´as flexible que con las f´ormulas de Newton, independientemente del n´umero de planos utilizado para aproximar la derivada axial. Cap´ıtulo 5 El problema de contorno 5.1. Introducci´on Desde el punto de vista matem´atico, la reconstrucci´on de fase en la t´ecnica basada en la ETI se reduce a resolver una ecuaci´on diferencial en derivadas parciales con condici´on de contorno de Neumann. Sin embargo, existen dificultades pr´acticas manifiestas a la hora de obtener medidas de frontera de calidad, debido, fundamentalmente, a las limitaciones espec´ıficas que afectan a este tipo de medidas respecto de las del interior. Este hecho hace que nos planteemos en este cap´ıtulo el objetivo de analizar el peso de las medidas de la derivada normal de la fase en la frontera en este tipo de m´etodos de recuperaci´on de fase. Asimismo, se estudiar´an soluciones para minimizar el error a˜nadido por las medidas de contorno, viendo cuando estas son estrictamente necesarias y cuando no lo son para poder recuperar la fase de forma ´unica. 5.2. Estimaci´on del frente de onda a partir de medidas de irradiancia Los algoritmos que se usan para reconstruir la fase a partir de la informaci´on proporcionada por un sensor de frente de onda se suelen clasificar en algoritmos zonales y algoritmos modales [160]. En el primer caso, el frente de onda se recupera en un conjunto discreto de puntos de medida, mientras que en el segundo se reconstruye en funci´on de los coeficientes de un desarrollo en serie funcional. Los principales algoritmos de reconstrucci´on de tipo zonal basados en la ETI han 106 El problema de contorno sido desarrollados para el caso particular de distribuciones de intensidad uniformes en el plano de recuperaci´on de fase. Estos m´etodos resuelven la Ec. (1.3) empleando funciones de Green [24], diferencias finitas [3, 52, 126] o transformadas de Fourier [161, 162]. El principal problema que caracteriza los m´etodos de estimaci´on zonal, para el caso que nos ocupa, es la dificultad de introducir en los algoritmos las condiciones de contorno, lo que generalmente conlleva errores que se a˜naden a los de las medidas de frontera. En los m´etodos que utilizan transformadas de Fourier se resuelve el problema imponiendo condiciones peri´odicas ficticias, como en el caso de G. Vdovin [3], o bien definiendo un dominio mayor que el real con el objetivo de tener en cuenta la contribuci´on de la frontera, como en el caso de Roddier et al. [161]. En el m´etodo utilizado por M. R. Teague, se calcula una soluci´on para la fase que depende de una funci´on de Green y su derivada normal [24]. No obstante, ambas presentan singularidades integrables que hacen perder precisi´on en la recuperaci´on de la fase. Las estimaciones modales de magnitudes f´ısicas son bastante habituales y aparecen en varias ´areas de la F´ısica (por ejemplo, el desarrollo multipolar del potencial electrost´atico en funci´on de polinomios de Legendre [163] o la soluci´on de la ecuaci´on de Schr¨odinger para el oscilador arm´onico, donde la funci´on de onda se expresa como combinaci´on lineal de funciones propias constituidas por el producto de un polinomio de Hermite y una funci´on exponencial de Gauss [164]). La principal ventaja de los m´etodos de estimaci´on modal respecto de los de estimaci´on zonal es que proporcionan una expresi´on de la magnitud como desarrollo en un conjunto de funciones. La recuperaci´on modal del frente de onda tiene por objetivo estimar los coeficientes del desarrollo de la fase en un conjunto de funciones base, generalmente polinomios, suponiendo as´ı un comportamiento anal´ıtico de la fase que permite extrapolarla a cualquier punto de la abertura. Aun cuando es posible utilizar como funciones base cualquier conjunto de polinomios, con frecuencia se escogen en ´ Optica los polinomios de Zernike, debido fundamentalmente a que son f´aciles de relacionar con las aberraciones cl´asicas [2, 165]. Se han desarrollado algunos m´etodos de estimaci´on modal del frente de onda espec´ıficos para su aplicaci´on combinada con la ETI [135–137, 146]. Sin embargo, estos m´etodos no son ortogonales, es decir, los coeficientes del desarrollo de la fase en la base escogida de funciones ortogonales no pueden obtenerse independientemente porque se expresan en funci´on de derivadas de dichas funciones. Aunque la fase se exprese de forma modal mediante un desarrollo en serie de funciones ortogonales, la 5.3 Estimaci´on modal ortogonal 107 actuaci´on de los operadores lineales laplaciano y derivada normal sobre las funciones base del desarrollo no garantiza que se conserve la ortogonalidad entre las funciones obtenidas, produci´endose un acoplamiento entre los coeficientes. Se han propuesto varias estrategias que permiten resolver el problema del acoplamiento entre los coeficientes modales para geometr´ıas concretas de los planos de medida [166, 167]. No obstante, el m´etodo de los proyectores ortogonales da una soluci´on m´as general y m´as adecuada a los objetivos de este cap´ıtulo [168]. Para ello se introduce un conjunto de funciones, los proyectores ortogonales o funciones de peso, que son ortogonales a las funciones que resultan de aplicar un operador lineal (gradiente, Laplaciano, ...) a las funciones base. Esto permite estimar cada coeficiente modal de forma independiente como la proyecci´on de las medidas del sensor sobre la funci´on base elegida. Este m´etodo no es necesariamente el mejor m´etodo de recuperaci´on de fase. Sin embargo, permite separar la informaci´on de la fase relativa al interior de la correspondiente a la frontera y la necesidad de imponer condiciones de contorno para un problema particular es mucho m´as intuitiva que en otros algoritmos, pues las condiciones de frontera aparecen en el m´etodo de forma natural. Por este motivo, este procedimiento resulta muy adecuado para estudiar la influencia de las condiciones de contorno en el problema de recuperaci´on de fase a partir de la Ecuaci´on de Transporte de Intensidad. En lo que sigue, describiremos brevemente los fundamentos en los que se basa el m´etodo. 5.3. Estimaci´on modal ortogonal Sea f(r) un observable que est´a relacionado con una determinada magnitud W(r) por medio de la ecuaci´on o conjunto de ecuaciones: L[r] = f(r) (5.1) donde Lrepresenta un operador lineal actuando sobre Wyrel vector de posici´on en un espacio N-dimensional. Sea una funci´on de aproximaci´on de W(r) construida como combinaci´on lineal de un conjunto de funciones base ψi(r) con coeficientes ai(i= 1,2,...M): QM(r) = M X i=1 aiψi(r) (5.2) 114 El problema de contorno A continuaci´on se ver´a qu´e implica el hecho de imponer la condici´on de contorno adicional ξk(ρ= 1, θ) = 0 en el c´alculo de las funciones de peso. (En lo que sigue se utilizar´a la notaci´on de Noll para los polinomios de Zernike [157]). (a) Modos rotacionalmente sim´etricos, Z0 n, con npar Los polinomios de Zernike que describen los modos rotacionalmente sim´etricos se pueden obtener a partir de la siguiente expresi´on: Z0 n=√n+ 1 n/2 X s=0 (−1)s(n−s)! s!n 2−s!2(5.19) Una soluci´on de la forma: ξ0 n(ρ) = Zρ 1 1 ρ′"Zρ′ 1 Z0 n(ρ′′)ρ′′dρ′′#dρ′(5.20) satisface la Ec. (5.18) y las condiciones de contorno de Neumann y de Dirichlet, [dξ0 n(ρ)/dρ]ρ=1 =ξ0 n(ρ= 1) = 0, a excepci´on de la funci´on de peso asociada al pist´on, con n=m= 0. Por tanto, los modos rotacionalmente sim´etricos, excepto el pist´on, pueden calcularse midiendo tan solo el Laplaciano de la fase. (b) Modos no sim´etricos, Zm n,m6= 0 Los polinomios de Zernike que describen los modos no sim´etricos se obtienen a partir de la siguiente f´ormula: Zm n=p2 (n+ 1) (n−m)/2 X s=0 (−1)s(n−s)! s! [(n+m)/2−s]! [(n−m)/2−s]!ρn−2s ×   sen (mθ) si m > 0 cos (mθ) si m < 0 (5.21) La funci´on de peso se puede calcular en este caso mediante el m´etodo de separaci´on de variables, obteni´endose: ξm n(ρ, θ) = Aρ|m|+Bρ−|m|+ρ|m| 2|m|Zρ 0 ρ′(1−|m|)Z|m| n(ρ′)dρ′ −ρ−|m| 2|m|Zρ 0 ρ′(1+|m|)Z|m| n(ρ′)dρ′   sen (mθ) si m > 0 cos (mθ) si m < 0 (5.22) donde AyBson constantes arbitrarias (correspondientes a la soluci´on de la ecuaci´on homog´enea). Con el fin de obtener soluciones anal´ıticas en todo el dominio, se impone B= 0. 5.4 An´alisis de la necesidad de las medidas en el contorno 115 (b.1) Modos cuya frecuencia azimutal, |m|, es diferente a la radial, n Si en la expresi´on (5.22) se toma: A=−1 2|m|Z1 0 ρ(1−|m|)Z|m| n(ρ)dρ (5.23) se puede comprobar f´acilmente que las funciones de peso verifican ambas condiciones homog´eneas, [∂ξm n(ρ, θ)/∂ρ]ρ=1 = 0 y ξm n(ρ= 1, θ) = 0, y, por tanto, estos modos tambi´en pueden calcularse a partir de medidas de la curvatura del frente solamente. (b.2) Modos cuya frecuencia azimutal, |m|, es igual a la radial, n:Zm |m| Este tipo de modos corresponden a funciones arm´onicas, soluci´on de la ecuaci´on de Poisson, separables en parte radial y parte angular. La expresi´on de las funciones de peso se reduce a: ξm |m|(ρ, θ) = Aρ|m|+ρ|m|+2 4 (|m|+ 1)   sen (mθ) si m > 0 cos (mθ) si m < 0 (5.24) siendo: A=−(|m|+ 2) 4|m|(|m|+ 1) (5.25) para obtener h∂ξm |m|(ρ, θ)/∂ρiρ=1 = 0. Sin embargo, en este caso ξm |m|(ρ= 1) = −1/[2 |m|(|m|+ 1)] 6= 0, implicando que para estimar estos modos se necesitan medidas en el contorno, hecho esperable debido a su car´acter arm´onico. Por tanto, es posible recuperar todos los coeficientes modales del desarrollo de la fase en polinomios de Zernike a partir de solamente datos de curvatura del frente de onda, a excepci´on de aquellos correspondientes a polinomios de frecuencia azimutal y radial iguales, es decir, los modos arm´onicos. En la Tabla 5.1 se incluyen algunas funciones de peso correspondientes al desarrollo de la fase en una base de polinomios de Zernike e iluminaci´on uniforme en z= 0. Condiciones de contorno para las funciones de peso asociadas a las funciones pseudoanal´ıticas de Karhunen-Lo`eve, Ki(ρ, θ) En aplicaciones relacionadas con turbulencias atmosf´ericas, los polinomios de Zernike no son un conjunto ´optimo de funciones base porque sus coeficientes 116 El problema de contorno n m Zm n(ρ, θ)ξm n(ρ, θ) 1 1 2ρsen (θ)−1 4Io(3ρ−ρ3) sen (θ) 2 0 √3 (2ρ2−1) −√3 8Io(−ρ4+ 2ρ2−1) 2 2 √6ρ2sen (2θ)−√6 12Io(2ρ2−ρ4) sen (2θ) 3 1 √8 (3ρ3−ρ) sen (θ)−√8 8Io(ρ5−2ρ3+ρ) sen (θ) 3 3 √8ρ3sen (3θ)−√8 48Io(−3ρ5+ 5ρ3) sen (3θ) 4 0 √5 (6ρ4−6ρ2+ 1) −√5 4Io(−ρ6+ 9ρ4−6ρ2+ 1) 4 2 √10 (4ρ4−3ρ2) sen (2θ)−√10 8Io(−ρ6+ 2ρ4−ρ2) sen (2θ) Tabla 5.1: Polinomios de Zernike y funciones de peso asociadas para iluminaci´on uniforme. asociados no son independientes estad´ısticamente. La necesidad de realizar un estudio adecuado de las turbulencias atmosf´ericas llev´o a la b´usqueda de unas funciones base cuyos coeficientes asociados no estuvieran correlacionados. Estas funciones son las denominadas funciones de Karhunen-Lo`eve. Al igual que los polinomios de Zernike, las funciones de Karhunen-Lo`eve pueden factorizarse en productos de funciones radiales y angulares [172]: Kl(ρ, θ) = Λq p(ρ) Θq(θ) (5.26) donde Θq(θ) = (sen (qθ) si q > 0 cos (qθ) si q < 0y donde las partes radiales se obtienen resolviendo la ecuaci´on integral: Z1 0 Yq(ρ, ρ′) Λq p(ρ)dρ =α2 kΛq p(ρ′) (5.27) siendo hα2 kila varianza normalizada de los coeficientes de Karhunen-Lo`eve e Yq(ρ, ρ′) la funci´on n´ucleo q-´esima dependiente de la funci´on de estructura normalizada. Para cada Yqdada existe un conjunto de autofunciones ortonormales, Λq p(ρ), cuya condici´on de ortogonalidad es: Z1 0 Λq p(ρ) Λq p′(ρ)ρdρ =1 2δpp′(5.28) La soluci´on de la Ec. (5.27), aun no siendo anal´ıtica, admite una buena aproximaci´on en serie de potencias de la forma: Λq p(ρ) = ρq∞ X s=0 gsρ2s(5.29) 5.4 An´alisis de la necesidad de las medidas en el contorno 117 p q gog1g2g3g4g5g6g7 1 0 2,1558 6,4545 4,1230 1,3654 0,2192 0,0320 1 1 1,5545 0,2444 0,0438 1 2 3,2273 2,7488 1,0245 0,1758 0,0269 1 3 6,0555 8,1895 4,5994 1,3254 0,2139 1 4 10,729 19,053 14,215 5,7788 1,3699 0,1279 2 0 3,1888 32,765 85,091 97,412 62,759 24,974 6,1831 0,7158 2 1 6,5287 19,114 17,922 8,6614 2,3235 0,3526 2 2 16,382 56,727 73,054 50,190 20,788 5,2854 0,6234 Tabla 5.2: Coeficientes de ajuste de la parte radial de las funciones pseudoanal´ıticas de KarhunenLo`eve [172]. donde la r´apida convergencia a cero de los coeficientes gshace posible obtener una aproximaci´on polin´omica precisa de las funciones Λq p(ρ). En la Tabla 5.2 aparecen tabulados los coeficientes asociados a funciones de Karhunen-Lo`eve, tomando una soluci´on aproximada de la Ec. (5.27) truncando (5.29) a 8 t´erminos. Teniendo en cuenta entonces que la parte radial de las funciones de KarhunenLo`eve admite un desarrollo anal´ıtico, esta se puede tambi´en desarrollar como combinaci´on lineal de polinomios de Zernike con frecuencias azimutales id´enticas de la forma: Λq p(ρ) = ∞ X k=2 UklRm n(ρ) (5.30) donde mrepresenta la frecuencia azimutal, el ´ındice kordena los polinomios de Zernike y se relaciona con los ´ındices (n, m), y el ´ındice lordena las funciones de Karhunen-Lo`eve y se relaciona con los ´ındices (p, q). Adem´as, en ciertas aplicaciones esta funci´on se puede aproximar por medio de un n´umero finito de polinomios de Zernike, con lo que la serie anterior aparece truncada. Las funciones de KarhunenLo`eve aproximadas de esta forma se llaman funciones pseudoanal´ıticas de KarhunenLo`eve. Como consecuencia, se pueden calcular las funciones de peso para este caso como combinaci´on lineal de las funciones de peso previamente calculadas en el apartado anterior. Sin embargo, este tipo de aproximaci´on para toda funci´on de KarhunenLo`eve, a excepci´on de las rotacionalmente sim´etricas, contiene siempre un polinomio de Zernike con car´acter arm´onico, de forma que no pueden evaluarse sin la se˜nal de 118 El problema de contorno k m Ukm k m Ukm 0 0 1∗2 1 y3+ 3x2y−y 1 0 x∗1 2 x3+ 3xy2−x 0 1 y∗2 2 5x4+ 5y4−6x2−6y2+ 18x2y2+ 1 1 1 xy∗3 0 5x3+ 3xy2−3x 2 0 3x2+y2−1 0 3 5y3+ 3x2y−3y 0 2 x2+ 3y2−1 3 1 5x3y+ 3xy3−3xy Tabla 5.3: Polinomios de Apell y Kamp´e de F´eriet. frontera. Por tanto, si se escogen las funciones de Karhunen-Lo`eve para desarrollar la fase, a partir ´unicamente de medidas del Laplaciano, solo se pueden recuperar los modos con m= 0. Condiciones de contorno para las funciones de peso asociadas a polinomios de Appell y Kamp´e de F´eriet, Ui(x, y) Estos polinomios nunca han sido utilizados en ´ Optica, pero aqu´ı nos servir´an para ilustrar otra faceta m´as del problema de contorno que representa la recuperaci´on de fase por medio de la ETI. Constituyen una base ortogonal en el c´ırculo de radio unidad, pero no todos son separables en parte radial y parte angular. Obedecen la expresi´on general [173]: Ukm =∂k+m ∂xk∂ymx2+y2−1k+m(5.31) En la Tabla 5.3 figura una lista de los 12 primeros polinomios de Apell y Kamp´e de F´eriet, se˜nal´andose los que tienen car´acter arm´onico por medio del s´ımbolo *. Se puede demostrar f´acilmente que no existen funciones de peso asociadas a estos polinomios tales que se anulen en el contorno porque el contenido arm´onico y el no arm´onico no est´a separado en esta base de polinomios. Adem´as, para algunos de ellos tampoco existe soluci´on al problema de Neumann o Dirichlet y, por tanto, se necesitan ambas condiciones de contorno para resolver el problema de acuerdo con el m´etodo de los proyectores ortogonales. De este modo, se concluye que dependiendo de la base ortogonal elegida para 5.4 An´alisis de la necesidad de las medidas en el contorno 119 desarrollar la fase, existe la posibilidad de recuperar algunos de los modos sin necesidad de medidas de la derivada del frente de onda en la frontera. Se ha demostrado que cuando la base est´a constituida por polinomios de Apell y Kamp´e de F´eriet se necesitan condiciones de contorno de Neumann y de Dirichlet para recuperar la fase sin ambig¨uedades. En cambio, cuando la fase se desarrolla bas´andose en funciones pseudoanal´ıticas de Karhunen-Lo`eve, solo se puede calcular la contribuci´on de los modos rotacionalmente sim´etricos sin medidas en la frontera. Por ´ultimo, para una base de polinomios de Zernike, solo los modos arm´onicos necesitan condiciones de frontera para poder ser recuperados. La m´axima informaci´on acerca del frente de onda se puede obtener sin condiciones de contorno mediante la elecci´on de una base ortogonal para desarrollar la fase que distinga modos arm´onicos y no arm´onicos, pues la informaci´on de la parte arm´onica de la fase se encuentra exclusivamente en la frontera. Los polinomios de Zernike son una base que re´une esta caracter´ıstica, adem´as de la de su ortogonalidad en el c´ırculo de radio unidad. En realidad, sin haber llevado a cabo un estudio te´orico riguroso, sino por simple inspecci´on de los primeros modos, I. Han ya hab´ıa observado que desarrollando la fase en una base de polinomios de Zernike pod´ıa obtener el contenido arm´onico a partir de la medida de Wren la frontera y del contenido no arm´onico obtenido midiendo ∇2Wen el interior [136]. Ya que existen bases para las que algunos modos pueden obtenerse sin necesidad de medidas en el contorno, se concluye que algunas de las aberraciones presentes en el frente de onda, relacionadas con la parte no arm´onica, pueden estimarse usando solo informaci´on del Laplaciano de la fase. En la secci´on siguiente, se ver´a si la elecci´on de una distribuci´on adecuada de irradiancia en el plano de medida, permite estimar la fase sin necesidad de medidas en el contorno, para lo cual se analizar´an dos ejemplos de distribuciones de irradiancia. 5.4.2. Distribuci´on no uniforme de irradiancia En esta secci´on se resolver´a la ETI para distribuciones de irradiancia no uniformes en el plano de medida. Nuestro objetivo es estudiar la posibilidad de obtener la fase de forma ´unica sin necesidad de medidas en la frontera dependiendo del comportamiento local de la distribuci´on de irradiancia en el contorno. Con el fin de simplificar los c´alculos, pero sin p´erdida de generalidad en las conclusiones principales, se analizar´an dos distribuciones rotacionalmente sim´etricas que corresponden a los dos perfiles de intensidad siguientes: 120 El problema de contorno (a) Perfil de intensidad gaussiano: I(ρ) =    I0e−βρ2 0 ρ≤1 ρ > 1 (5.32) (b) Perfil de intensidad parab´olico: I(ρ) =    I0(1 −ρ2) 0 ρ≤1 ρ > 1 (5.33) Sea cual sea la distribuci´on de irradiancia en el plano de recuperaci´on, las fases que no ser´ıan detectables por medio de medidas de irradiancia en el interior del dominio son aquellas que sean soluci´on de la ecuaci´on homog´enea, es decir, aquellas que verifiquen: ∇·(I∇W) = 0 (5.34) pues I′(0) = 0, con lo que los coeficientes modales solo podr´ıan calcularse en funci´on de medidas de Wo de Wren la frontera. Perfil gaussiano de intensidad En la Ref. [23], Teague propone la elecci´on de un perfil de intensidad que caiga pr´acticamente a cero en la frontera, con el fin de minimizar la influencia del contorno en el proceso de recuperaci´on de fase a partir de medidas de irradiancia en ausencia de medidas en el borde. En concreto, se propone un perfil de la forma Ioe−βρ2N. Para realizar en este trabajo un estudio cualitativo riguroso de la necesidad de las medidas en el contorno cuando la distribuci´on de irradiancia tiene esta forma, ser´a suficiente con escoger un perfil gaussiano (N= 1). Entonces, sup´ongase que la distribuci´on de irradiancia en el plano de recuperaci´on viene determinada por el perfil gaussiano de la Ec. (5.32). Para cada polinomio de Zernike del desarrollo de la fase, resolveremos la ecuaci´on que proporciona las funciones de peso, es decir: ∇·hIoe−βρ2∇ξm n(ρ, θ)i=Zm n(ρ, θ) (5.35) Se distinguen dos casos: 5.4 An´alisis de la necesidad de las medidas en el contorno 121 (a) Modos rotacionalmente sim´etricos Z0 n, con npar La funci´on de peso asociada tiene la forma: ξ0 n(0) = Zρ 1 eβρ′2 Ioρ′ Zρ′ 1 ρ′′Z0 n(ρ′′)dρ′′!dρ′(5.36) y satisface la Ec. (5.35) con condiciones de contorno de Neumann y de Dirichlet: [dξ0 n(ρ)/dρ]ρ=1 =ξ0 n(ρ= 1) = 0 . De esta forma, para recuperar los modos rotacionalmente sim´etricos, salvo el pist´on, solo son necesarias medidas del gradiente axial de intensidad en el dominio, de acuerdo con la ecuaci´on: ak=−ZΣ ξkI′(0) dΣ (5.37) (b) Modos no sim´etricos Zm n En este caso, la soluci´on de la Ec. (5.35) no es tan sencilla como la anterior y, para calcularla, es necesario simplificar la ecuaci´on realizando diversos cambios de variable. Primero resolveremos la ecuaci´on homog´enea y despu´es encontraremos una soluci´on particular con el fin de hallar la soluci´on general de la ecuaci´on de las funciones de peso. Suponiendo que las funciones de peso son factorizables en parte radial y parte angular, ξm n(ρ, θ) = Rm n(ρ) Θ (mθ) con Θ (mθ) = (sen (mθ) cos (mθ) si si m > 0 m < 0, el problema se reduce a encontrar la soluci´on de la ecuaci´on homog´enea para la parte radial: ρ2d2R dρ2+ρ1−2ρ2α2dR dρ −m2R= 0 (5.38) es decir: R(ρ) = ρ|m|φ|m| 2,|m|+ 1; βρ2(5.39) donde φes la funci´on hipergeom´etrica confluente [159]. Por otra parte, se puede demostrar que: ξm n(ρ, θ) = eβρ2Pm n(ρ) Θ (mθ) (5.40) es una soluci´on particular de la Ec. (5.35), donde Pm n(ρ) es un polinomio de grado nque se puede calcular por el m´etodo de variaci´on de par´ametros [150] 122 El problema de contorno y, por tanto, una funci´on de la forma: ξm n(ρ, θ) = Aρ−|m|φ|m| 2,|m|+ 1; βρ2+Bρ|m|φ|m| 2,|m|+ 1; βρ2 +eβρ2Pm n(ρ)iΘ (mθ) (5.41) es la soluci´on de la Ec. (5.35). Las constantes AyBse eligen dependiendo de las condiciones de contorno necesarias para resolver el problema y de las propiedades de las funciones de peso. Por ejemplo, las funciones de peso han de poderse integrar f´acilmente y, por eso, las soluciones que se buscan son anal´ıticas. As´ı, se toma A= 0, mientras que Bse elige de forma que ξm nverifique condiciones de contorno de Neumann, es decir, [∂ξm n(ρ, θ)/∂ρ]ρ=1 = 0. Una vez escogida Bde esta forma, se puede demostrar que ξm nno verifica simult´aneamente la condici´on de Dirichlet, es decir, ξm n(ρ= 1, θ)6= 0. Por tanto, se necesitan condiciones de contorno de Neumann para la fase, con el fin de poder evaluar la contribuci´on de este tipo de modos de acuerdo con la expresi´on: ak=−ZΣ ξkI′(0) dΣ−IC WrIξkdθ (5.42) De estos resultados se concluye que cuando la distribuci´on de irradiancia en el plano de recuperaci´on se puede describir por un perfil gaussiano, se necesitan medidas de la pendiente radial de la fase en la frontera con el fin de encontrar soluci´on ´unica anal´ıtica al problema. En la Tabla 5.4 se incluye un listado de las funciones de peso para los 12 primeros polinomios de Zernike, correspondientes a esta distribuci´on de irradiancia. No obstante, en la Ec. (5.42) se observa que el comportamiento del producto I(ρ= 1, θ)ξk(ρ= 1, θ) juega un papel muy importante en la contribuci´on de la frontera en la fase recuperada. Es decir, la integral de frontera HCWrIξkdθ puede llegar a ser despreciable respecto de la integral de superficie en el dominio, para valores lo suficientemente peque˜nos del producto Iξk. En la Figura 5.1 se muestra el comportamiento de la parte radial de dicho producto en la frontera para los modos arm´onicos ξm |m|en funci´on del par´ametro β, que define la anchura de la intensidad gaussiana. A medida que el par´ametro βde la gaussiana aumenta, el producto I(ρ= 1) Rm |m|(ρ= 1) disminuye, con lo que para valores de βlo suficientemente grandes se puede despreciar la integral de contorno y estimar la fase de forma 5.4 An´alisis de la necesidad de las medidas en el contorno 123 n m Zm n(ρ, θ)ξm n(ρ, θ) 1 1 2ρsen (θ)−ρ 3Ioβ(1 + 2β)φ(1 2,2;βρ2) φ(3 2,2;β)eβ−eβρ2sen (θ) 2 0 √3 (2ρ2−1) √3 4Ioβ2heβ−(1 + β−βρ2)eβρ2i 2 2 √6ρ2sen (2θ)−√6 8Ioβρ2(1 + β)φ(1,3;βρ2) φ(2,3;β)eβ−eβρ2sen (2θ) 3 1 √8 (3ρ3−ρ) sen (θ)−√8 30Ioβ2(−2β2+ 5β−6) φ(1 2,2;βρ2) φ(3 2,2;β)ρeβ−[9βρ3−(10β+ 6) ρ]eβρ2sen (θ) 3 3 √8ρ3sen (3θ)−√8 30Ioβρ3(3 + 2β)φ(3 2,4;βρ2) φ(5 2,4;β)eβ−3eβρ2sen (3θ) 4 0 √5 (6ρ4−6ρ2+ 1) √5 4Ioβ3n(β−4) eβ−[(4 + 3β)βρ2−2β2ρ4−(β2+ 3β+ 4)] eβρ2o 4 2 √10 (4ρ4−3ρ2) sen (2θ)−√10 24Ioβ2(4 + 3β−β2)φ(1,3;βρ2) φ(2,3;β)ρ2eβ−(8βρ4−13ρ2)eβρ2sen (2θ) Tabla 5.4: Polinomios de Zernike y funciones de peso asociadas para un perfil gaussiano de intensidad. Ap´endice A La Ecuaci´on de Transporte de Intensidad (ETI) A continuaci´on deduciremos la Ecuaci´on de Transporte de Intensidad (ETI) a partir de la Ecuaci´on Paraxial de Ondas y calcularemos la versi´on para haces que atraviesan una pupila en el plano z= 0. Consideremos una onda electromagn´etica monocrom´atica que se propaga en la direcci´on z, cuya amplitud compleja viene dada por: u(r, z) = A(r, z)eikW (r,z)eikz (A.1) donde r= (x, y), W(r, z) es la fase en el plano z,kes el n´umero de ondas y A(r, z) = pI(r, z), siendo I(r, z) la intensidad de dicha onda en un plano z. La amplitud compleja, u(r, z), es una soluci´on de la Ecuaci´on Paraxial de Ondas: 2ik ∂ ∂z +∇2u(r, z) = 0 (A.2) siendo ∇2=∂2 ∂x2+∂2 ∂y2el laplaciano bidimensional, es decir: 2ik ∂ ∂z +∇2A(r, z)eikW(r,z)eikz = 0 (A.3) Operando se obtiene la ecuaci´on siguiente: −2k2A∂W ∂z +∇2A−kA2|∇W|2eikW +i2ik∂A ∂z +kA∇2W+ 2∇A·∇WeikW = 0 (A.4) 132 La Ecuaci´on de Transporte de Intensidad (ETI) Igualando a cero la parte imaginaria y teniendo en cuenta que ∇A=∇I/ 2√I, se deduce que: I∇2W+∇I·∇W=−∂I ∂z (A.5) que no es m´as que la Ecuaci´on de Transporte de Intensidad. Normalmente se ver´a, sin embargo, de la forma: ∇·(I∇W) = −∂I ∂z (A.6) La ETI relaciona la variaci´on axial instant´anea de la irradiancia en un plano z con la variaci´on transversal de la fase en dicho plano. La versi´on de la Ec. (A.6) es adecuada para estudiar la propagaci´on de ondas que se extienden en todo el espacio. Calculemos ahora una relaci´on similar para un haz que atraviesa una pupila en el plano z= 0. Denotaremos los ejes transversales en cada plano zcomo xy y los del plano z= 0 como xoyo. Si trabajamos en el plano z= 0, la ETI es: ∂I(z) ∂z z=0 =−∇·[I(0)∇W] (A.7) donde I(z) = I(x, y, z) es la irradiancia en cada punto (x, y) de un plano z,I(0) = I(xo, yo,0) es la irradiancia en z= 0 y W=W(xo, yo) es la fase en cada punto (xo, yo) del plano z= 0. La onda (A.1) se propaga a lo largo del eje zhasta que en el plano z= 0 atraviesa la pupila circular Po=P(xo, yo). En concreto, la pupila modula la intensidad del plano de la abertura y ejerce un efecto en la propagaci´on de la luz, que ha de describir la ecuaci´on que proporciona la variaci´on axial de irradiancia. De este modo, la distribuci´on de irradiancia en el plano z= 0, I(0), se expresa en funci´on de la pupila Pode la forma: I(0) = I(xo, yo,0) = I(xo, yo)P(xo, yo) = IoPo(A.8) donde Io=I(xo, yo) es la irradiancia en los puntos del interior de Po. Si se sustituye en el miembro de la derecha de la ETI se tiene: −∇·[I(0)∇W] = −∇·[IoPo(∇W)] =−(∇Io)(∇W) + Io(∇2W)Po−Io(∇W)(∇Po) (A.9) Para calcular (∇W) (∇Po) empleamos coordenadas polares planas (ro, θo), de modo que la pupila circular se puede modelar como Po=P(ro) = 1 −H(ro−R), donde 133 H(ro−R) es la funci´on de salto unitario en el punto ro=R, siendo Rel radio de la pupila, es decir [158]: H(ro−R) =    1 si ro> R 0 si ro≤R (A.10) De este modo: (∇W) (∇Po) = Wroˆ ro+Wθo ro ˆ θo(Proˆ ro) = WroPro(A.11) donde Wro=∂W/∂ro,Wθo=∂W/∂θoyPro=∂Po/∂ro. Como Po= 1 −H(ro−R) entonces: Pro=∂ ∂ro [1 −H(ro−R)] = −δ(ro−R) = −δc(A.12) donde el sub´ındice cindica que la delta de Dirac coincide con el borde de la pupila, C. Si sustituimos en la Ec. (A.9): −∇·[IoPo(∇W)] = −∇·(Io∇W)Po+IoWroδc(A.13) De este modo, la ETI pasa a tener la forma usual: ∂I (z) ∂z z=0 =−∇·(Io∇W)Po+IoWroδc(A.14) Ap´endice B Rendimiento del sensor de curvatura en condiciones de iluminaci´on limitadas A lo largo de esta memoria, se han dado los pasos necesarios para dise˜nar un procedimiento que ayude a optimizar la respuesta de un sensor de curvatura de acuerdo a criterios objetivos. Para ello, hemos analizado cualitativamente el rendimiento del sensor, suponiendo que todas las fluctuaciones de las medidas de irradiancia que realiza se deb´ıan a ruido de detecci´on. Sin embargo, cuando el sensor se emplea para recuperar aberraciones din´amicas como, por ejemplo, frentes de onda turbulentos, las medidas de irradiancia pueden constar de un n´umero bajo de fotoelectrones y, consecuentemente, el ruido fot´onico adquiere m´as protagonismo. La diferente estad´ıstica de un tipo de fluctuaciones y otro implica algunos cambios en el c´alculo de las varianzas de las curvaturas y de las derivadas normales locales, aunque, como veremos, su comportamiento cualitativo en funci´on de la posici´on de los planos de irradiancia en el eje zno se modifica apreciablemente. No obstante, las diferencias pueden traducirse en cambios cuantitativos significativos del error y de la posici´on ´optima de los planos. En este ap´endice nos proponemos calcular estas varianzas en funci´on del nivel de ruido fot´onico y del ruido de detecci´on, que nos permitir´an hallar la posici´on ´optima de los planos conocido el n´umero de fotones promedio detectados durante el tiempo de exposici´on. 136 Rendimiento del sensor de curvatura en condiciones de iluminaci´on limitadas B.1. Varianza de la estimaci´on de la curvatura local, cotas de error y posici´on ´optima Supongamos que durante un cierto tiempo de exposici´on, τ, un detector Σ capta en un plano jz un n´umero de fotones, ij, tal que: ij=Ij+nrms j+nph j(B.1) donde Ijes el n´umero de fotones que se detectan en promedio en el plano jz durante el tiempo de exposici´on, τ,nrms json errores de detecci´on, aditivos, con media cero, varianza σ2y covarianza nula entre detectores y planos distintos, mientras que nph j es ruido fot´onico, aditivo, de media cero, con una varianza σ2 j=Ijy covarianza nula entre detectores y planos distintos. De acuerdo con el Cap´ıtulo 3, a partir de estas medidas de irradiancia se construye una f´ormula en diferencias finitas para aproximar la derivada axial, ˆ I′(0), en el plano z= 0. La relaci´on entre ˆ I′(0) y la curvatura del frente de onda que proporciona la ETI se puede emplear para calcular la varianza de la curvatura del frente a partir de s2hˆ I′(0)i. Si en lugar de la Ec. (4.1) se introduce la Ec. (B.1) en los c´alculos de la §4.2, entonces en la varianza de la derivada axial aproximada en diferencias finitas interviene un t´ermino m´as que informa de la contribuci´on del ruido fot´onico, que all´ı no tuvimos en cuenta. En concreto, obtenemos la siguiente ecuaci´on para la varianza de ˆ I′(0) estimada en Σ con la f´ormula central con M= 2mplanos (Ec. (4.2)): s2hˆ I′(0)i=1 z2 σ2˜am+ m X j=−m a2 jIj!+z2"m X j=−m aj j2 2I′′(ξj)#2 (B.2) La ecuaci´on correspondiente a la f´ormula progresiva y regresiva con M=m+ 1 planos (Ec. (4.13)) es: s2hˆ I′(0)i=1 z2"2σ2+Io˜ bm+ m X j=1 b2 jIj#+z2"m X j=1 bj j2 2I′′(ξj)#2 (B.3) Los coeficientes aj,bj, ˜amy˜ bmson los mismos que los del Cap´ıtulo 4 (ver §4.2). Supongamos que durante el tiempo de exposici´on un detector Σ capta en promedio un n´umero finito de fotones, N, que se distribuyen entre los distintos planos de medida, es decir: N= m P j=−m Ij, N =Io+ m P j=1 Ij(B.4) B.1 Varianza de la estimaci´on de la curvatura local, cotas de error y posici´on ´optima 137 para el esquema central y para las f´ormulas de Newton, respectivamente. Ahora supongamos que los frentes de onda son lo suficientemente suaves como para que en el ´area de cada detector Σ se pueda considerar que la fase se comporta esencialmente como una fase cuadr´atica, Wo, (Ec. (2.15)). Introduciendo s2hˆ I′(0)i en la Ec. (4.20) y teniendo en cuenta la Ec. (2.15), calculamos la varianza de la curvatura local estimada en Σ en funci´on de Npara ambos tipos de f´ormulas: a) F´ormula central ǫ2 o=σ N2˜am z2+ 4(D+E)2"m X j=−m 1 1 + 2(D+E)jz + 4DEj2z2#2 +1 N"m X j=−m 1 1 + 2(D+E)jz + 4DEj2z2# ×"4(D+E)2+1 z2 m X j=−m a2 j 1 + 2(D+E)jz + 4DEj2z2# +4(D+E)2"1−˜am m X j=−m j2 [1 + 4DEj2z2]2−[2(D+E)jz]2#2 (B.5) b) F´ormula progresiva y regresiva ǫ2 o=σ N2"2˜ bm z2+ 4(D+E)2#"1 + m X j=1 1 1±2(D+E)jz + 4DEj2z2#2 +1 N"1 + m X j=1 1 1±2(D+E)jz + 4DEj2z2# ×"˜ bm z2+ 4(D+E)2+1 z2 p X j=1 b2 j 1±2(D+E)jz + 4DEj2z2# +"2(D+E)±1 z m X j=1 dj∓2(D+E)jz −4DEj2z2 1±2(D+E)jz + 4DEj2z2#2 (B.6) En la Figura B.1, para dos planos de medida, y en la Figura B.2, para ocho planos, se muestran ejemplos de las curvas de error que se obtienen a partir de la Ec. (B.5) para los mismos valores de los par´ametros de la fase cuadr´atica local en Σ, Wo, que en las curvas an´alogas de las Figuras 2.2 y 4.1. Se supone que el n´umero de fotones detectados es N= 104y que el detector es cuasiideal, esto es, σ= 0. 138 Rendimiento del sensor de curvatura en condiciones de iluminaci´on limitadas 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 0.3 z ǫo (a) D= 1/2, E = 0 2 (D+E) = 1 D= 1/3, E = 0 2 (D+E) = 2/3 D= 1/6, E = 0 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 0.3 z ǫo (b) D= 1/3, E = 1/6 2 (D+E) = 1 D= 2/9, E = 1/9 2 (D+E) = 2/3 D= 1/9, E = 1/18 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 0.3 z ǫo (c) D= 1/4, E = 1/4 2 (D+E) = 1 D= 1/6, E = 1/6 2 (D+E) = 2/3 D= 1/12, E = 1/12 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 0.3 z ǫo (d) D= 1/3, E =−1/6 2 (|D|+|E|) = 1 D= 2/9, E =−1/9 2 (|D|+|E|) = 2/3 D= 1/9, E =−1/18 2 (|D|+|E|) = 1/3 D= 0, E = 0 2 (|D|+|E|) = 0 Figura B.1: Desviaci´on est´andar de la curvatura local del frente de onda en Σ en funci´on de la posici´on de los planos, z, suponiendo que se estima con la f´ormula central particularizada a M= 2 planos de medida, con N= 104yσ= 0. El apartado (a) corresponde a aproximaciones cuadr´aticas con forma de cilindro parab´olico, el (b) a aproximaciones con forma de paraboloide el´ıptico, el (c) a aproximaciones con forma de paraboloide parab´olico y el (d) a aproximaciones con forma de paraboloide hiperb´olico. Esto quiere decir que las fluctuaciones relativas por ruido fot´onico, con desviaci´on 1/√N= 0,01, son de la misma magnitud que las consideradas en las Figuras 2.2 y 4.1, causadas por ruido de detecci´on con desviaci´on σ/Io= 0,01, y los cambios entre unas curvas y otras se deben a la diferente estad´ıstica de las fluctuaciones de los datos de irradiancia. A pesar de ello, el comportamiento cualitativo de unas curvas y otras es id´entico y podremos aplicar el procedimiento dise˜nado en el Cap´ıtulo 2 para optimizar la estimaci´on de las curvaturas locales en Σ. A partir de las Figuras B.3, para dos planos de medida, y B.4, para ocho planos, se puede llegar a las mismas conclusiones para las f´ormulas de Newton. Continuando con el esquema seguido en esta memoria, calculamos ahora las cotas superior e inferior de error obtenidas particularizando las Ecs. (B.5) y (B.6) a D= B.1 Varianza de la estimaci´on de la curvatura local, cotas de error y posici´on ´optima 139 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 0.3 z ǫo (a) D= 1/2, E = 0 2 (D+E) = 1 D= 1/3, E = 0 2 (D+E) = 2/3 D= 1/6, E = 0 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 0.3 z ǫo (b) D= 1/3, E = 1/6 2 (D+E) = 1 D= 2/9, E = 1/9 2 (D+E) = 2/3 D= 1/9, E = 1/18 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 0.3 z ǫo (c) D= 1/4, E = 1/4 2 (D+E) = 1 D= 1/6, E = 1/6 2 (D+E) = 2/3 D= 1/12, E = 1/12 2 (D+E) = 1/3 D= 0, E = 0 2 (D+E) = 0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 0.3 z ǫo (d) D= 1/3, E =−1/6 2 (|D|+|E|) = 1 D= 2/9, E =−1/9 2 (|D|+|E|) = 2/3 D= 1/9, E =−1/18 2 (|D|+|E|) = 1/3 D= 0, E = 0 2 (|D|+|E|) = 0 Figura B.2: Desviaci´on est´andar de la curvatura local del frente de onda en Σ en funci´on de la posici´on de los planos, z, suponiendo que se estima con la f´ormula central particularizada a M= 8 planos de medida, con N= 104yσ= 0. El apartado (a) corresponde a aproximaciones cuadr´aticas con forma de cilindro parab´olico, el (b) a aproximaciones con forma de paraboloide el´ıptico, el (c) a aproximaciones con forma de paraboloide parab´olico y el (d) a aproximaciones con forma de paraboloide hiperb´olico. V/2 y E= 0, y D=E= 0, respectivamente: a) F´ormula central a.1) Cota superior de error ǫsup V2=σ N2˜am h2+ 1 m X j=−m 1 1 + jh!2 +1 N m X j=−m 1 1 + jh! × 1 + 1 h2 m X j=−m a2 j 1 + jh!+"1−˜am m X j=−m j2 1−(jh)2#2 (B.7) a.2) Cota inferior de error ǫinf V2= 2m˜am h22mσ N2+1 N (B.8)