Repositorio Institucional de Documentos

2014 13 Elvis Javier Lacruz Calderón Problema de Lambert para órbitas perturbadas : aplicación a la búsqueda de órbitas cuasiestacionarias Departamento Director/es Matemática Aplicada Abad Medina, Alberto José Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Elvis Javier Lacruz Calderón PROBLEMA DE LAMBERT PARA ÓRBITAS PERTURBADAS : APLICACIÓN A LA BÚSQUEDA DE ÓRBITAS CUASIESTACIONARIAS Director/es Matemática Aplicada Abad Medina, Alberto José Tesis Doctoral Autor 2014 Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Problema de Lambert para ´orbitas perturbadas. Aplicaci´on a la b´usqueda de ´orbitas cuasiestacionarias Elvis Javier Lacruz Calder´on Problema de Lambert para ´orbitas perturbadas. Aplicaci´on a la b´usqueda de ´orbitas cuasiestacionarias Elvis Javier Lacruz Calder´on Memoria presentada para optar al grado de Doctor en Ciencias Matem´aticas. Dirigida por Dr. Alberto J. Abad Medina. ´ Indice general Introducci´on V 1. Movimiento orbital 1 1.1. Movimiento kepleriano y orbital . . . . . . . . . . . . . . . . . . . . . 1 1.1.1. Movimiento kepleriano . . . . . . . . . . . . . . . . . . . . . . 1 1.1.2. Movimiento orbital . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.3. Orbitas keplerianas . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.4. ´ Orbitas .............................. 7 1.2. Sistemas de referencia . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2.1. Sistemas de referencia espacial y planetoc´entrico . . . . . . . . 8 1.2.2. ´ Orbitarelativa .......................... 10 1.2.3. Posici´on, velocidad y aceleraci´on en un sistema no inercial . . 13 1.3. Ecuaciones fundamentales . . . . . . . . . . . . . . . . . . . . . . . . 16 1.3.1. Ecuaciones del movimiento de un sat´elite artificial . . . . . . 16 1.3.2. Ecuaciones variacionales . . . . . . . . . . . . . . . . . . . . . 19 2. Diferenciaci´on autom´atica y c´alculo de perturbaciones 23 ii ´ Indice general 2.1. Diferenciaci´on autom´atica . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2. Potencial gravitacional de un planeta . . . . . . . . . . . . . . . . . . 28 2.2.1. Modelos gravitacionales . . . . . . . . . . . . . . . . . . . . . 29 2.2.2. Esquema de c´alculo del potencial gravitatorio de un planeta . 33 2.2.3. C´alculo de derivadas parciales del potencial gravitatorio . . . . 41 2.2.4. Test num´ericos . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.3. Otrasfuerzas ............................... 51 2.3.1. Fuerzas conservativas . . . . . . . . . . . . . . . . . . . . . . . 51 2.3.2. Fuerzas no conservativas . . . . . . . . . . . . . . . . . . . . . 54 3. Arcos orbitales y problema de Lambert 59 3.1. Arcos keplerianos. Problema cl´asico de Lambert . . . . . . . . . . . . 59 3.1.1. ´ Orbitas keplerianas que pasan por dos puntos . . . . . . . . . 59 3.1.2. Transferencias orbitales . . . . . . . . . . . . . . . . . . . . . . 60 3.1.3. Problema y teorema de Lambert . . . . . . . . . . . . . . . . . 64 3.1.4. M´etodo de Battin para la resoluci´on del problema de Lambert 66 3.2. Arcos orbitales. Problema de Lambert generalizado . . . . . . . . . . 72 3.2.1. Problema de Lambert para un modelo orbital perturbado . . . 72 3.2.2. M´etodo de correcci´on de ´orbitas peri´odicas . . . . . . . . . . . 72 3.2.3. Extensi´on del m´etodo de correcci´on de ´orbitas . . . . . . . . . 76 3.2.4. M´etodo de Lambert generalizado . . . . . . . . . . . . . . . . 78 ´ Indice general iii 3.2.5. Aplicaci´on: b´usqueda de arcos orbitales alrededor de la Tierra ydelaLuna............................ 79 4. Arcos orbitales cerrados 85 4.1. Arcos keplerianos cerrados . . . . . . . . . . . . . . . . . . . . . . . . 86 4.1.1. Coordenadas planetoc´entricas de un punto de la ´orbita relativa 88 4.1.2. Algunas propiedades de los arcos keplerianos cerrados . . . . . 90 4.1.3. ´ Orbitas que repiten la traza con un arco kepleriano cerrado: ´orbitas tipo Molniya . . . . . . . . . . . . . . . . . . . . . . . 93 4.2. ´ Orbitas cuasiestacionarias . . . . . . . . . . . . . . . . . . . . . . . . 96 4.2.1. ´ Orbitas cuasiestacionarias en la Luna . . . . . . . . . . . . . . 100 4.3. Arcos orbitales cerrados: mantenimiento de ´orbitas cuasiestacionarias 107 Conclusiones y trabajo futuro 113 Introducci´on Kepler, a comienzos del siglo XVII, enunci´o tres leyes que describen el movimiento de los planetas. A finales del mismo siglo, Newton da un paso m´as all´a y formula la ley de gravitaci´on universal, que re´une los principios de la Mec´anica y del C´alculo Diferencial para dar una s´olida teor´ıa sobre el movimiento de los cuerpos en el espacio. La Mec´anica Celeste recoge todos los esfuerzos de los principales matem´aticos de los siglos XVIII al XX en el estudio de las caracter´ısticas de este movimiento. El a˜no 1957 supone un punto de inflexi´on en este tema. El 4 de octubre de 1957, la antigua Uni´on Sovi´etica lanz´o al espacio, con ´exito, el sat´elite Sputnik1-1, desde la base de Kazajist´an. Este hecho marca el inicio de la era espacial y con ´el nace la Astrodin´amica, que extiende los conocimientos de la Mec´anica Celeste a movimientos en el espacio de cuerpos no naturales. La Astrodin´amica no solo extiende el rango y complejidad de los problemas de la Mec´anica Celeste a˜nadiendo nuevas perturbaciones al modelo propuesto por Kepler, sino que tambi´en introduce otros problemas nuevos como es, por ejemplo, el dise˜no de complejas trayectorias que permitan a una nave viajar por todo el sistema solar o incluso, en el futuro, modificar la ´orbita de objetos naturales como los asteroides potencialmente peligrosos para la Tierra. Desde entonces, y a lo largo de estas ´ultimas cinco d´ecadas, se han realizado una gran variedad de misiones espaciales, algunas con objetivos muy diferentes. Una de las partes fundamentales de cada proyecto es el an´alisis de misi´on o la determinaci´on previa de la ´orbita que seguir´a el orbitador para cumplir las especificaciones de la misi´on (Wertz and Larson, 2010). El poder imponer para una ´orbita el cumplimiento de determinadas cualidades, particulares o generales, ha permitido alcanzar con ´exito muchos de los objetivos de las misiones ejecutadas hasta hoy. Con ello ha habido un 1Cuatro d´ıas antes del lanzamiento, en el Comit´e Sp´ecial de l’Ann´e G´eophysique Internationale organizado en Whashington, el Profesor Segei M. Poloskov en su presentaci´on titulada “Sputnik Zemli”, di´o a entender que eran capaces de lanzar un sat´elite al espacio en las siguientes semanas. vi Introducci´on incremento progresivo del espectro de posibilidades para misiones espaciales. De entre toda la variedad de ´orbitas de sat´elites artificiales terrestres las m´as conocidas son las llamadas ´orbitas geoestacionarias. Estas ´orbitas pertenecen a un conjunto m´as amplio, que llamaremos ´orbitas s´ıncronas2, caracterizado por su sincron´ıa con la rotaci´on del cuerpo central, esto es, porque el valor de su periodo orbital coincide con el de el periodo de rotaci´on del planeta (o la Luna en su caso). Si el cuerpo central es la Tierra las ´orbitas son llamadas geos´ıncronas y su periodo orbital debe ser igual a un d´ıa sid´ereo, o lo que es igual a 23h56m4. s09, lo que equivale a decir que su semieje mayor debe medir 42164 km. Las ´orbitas estacionarias3son aquellas ´orbitas s´ıncronas cuya excentricidad e inclinaci´on son nulas. Las ´orbitas estacionarias aparecen siempre, para un observador situado en la superficie del cuerpo central, como un punto fijo en el ecuador celeste. Las primeras ideas sobre la ´orbita geoestacionaria fueron publicadas por el austroh´ungaro Herman Potoˇcnik en 1928 en su trabajo titulado “El problema del viaje espacial-El motor cohete”. A˜nos despu´es, en 1945, el brit´anico Arthur C. Clarke hace menci´on a las ventajas que puede tener el uso de esta ´orbita para las comunicaciones, pues con solo tres sat´elites es suficiente establecer comunicaci´on con cualquier punto del planeta, exceptuando regiones cerca de los polos. En el a˜no 1963 es lanzado, y puesto en ´orbita, el Syncom-2 que es el primer sat´elite de comunicaciones en ´orbita geoestacionaria. Esto confirma las ideas de Arthur C. Clarke y constituye un avance fundamental en distintas disciplinas o ´areas del conocimiento como son: las telecomunicaciones, climatolog´ıa, oceanograf´ıa, geodesia, militar, entre muchas otras m´as. Desde el ´exito de la misi´on del Syncom-2 se han puesto centenares de sat´elites en esta ´orbita. Debido a las peculiaridades de las ´orbitas geoestacionarias, podemos pensar en ´estas como restringidas a una zona del espacio que forma un aro o anillo ecuatorial alrededor de la Tierra, delimitado por unos pocos kil´ometros desde la propia ´orbita geoestacionaria (Flohrer et al., 2011), por ello, existen grandes limitaciones al n´umero m´aximo de sat´elites geostacionarios posibles, lo que obliga a retirar los que quedan inactivos a una ´orbita cementerio para dejar un hueco a nuevos sat´elites. Cincuenta a˜nos despu´es del primer sat´elite geoestacionario, la franja ecuatorial, donde se encuentran, comienza a estar congestionada (a d´ıa de hoy existen 416 sat´elites geoestacionarios activos). Estudios realizados por Milani et al. 2Geos´ıncronas, areos´ıncronas y selenos´ıncronas para la Tierra, Marte y la Luna respectivamente. 3Geoestacionarias, areoestacionaria y selenoestacionarias para la Tierra, Marte y la Luna respectivamente. vii (2011) describen c´omo es la distribuci´on de los orbitadores en esta zona del espacio, en base a datos obtenidos por medios ´opticos y telem´etricos, los cuales han permitido catalogarlos como sat´elites activos, no activos y los que no est´an controlados con tama˜nos superiores a 1 cm. A fines pr´acticos la zona geoestacionaria, o el anillo que forma, se ha dividido a su vez en peque˜nas ventanas continuas con espacio disponible para albergar uno o varios sat´elites. Existen razones que justifican este hecho como son: evitar el cruce de las se˜nales de emisi´on y recepci´on entre una estaci´on terrena y el orbitador (Maral and Bousquet, 2009), tambi´en para evitar colisiones entre s´ı, realizando las maniobras necesarias para corregir su posici´on nominal, producto de la deriva ocasionada por las perturbaciones orbitales. Estos desplazamientos han sido estudiados con modelos te´oricos o pueden ser calculados con m´etodos observacionales (Montojo et al., 2011; Montojo et al.; Lacruz and Abad, 2008). Un estudio detallado de las caracter´ısticas fundamentales que determinan la ´orbita geoestacionaria lo podemos encontrar en Soop (1994). El tama˜no habitual de las ventanas suele ser del orden de 1◦para la longitud (direcci´on Este-Oeste) y de [−0·1◦,+0·1◦] para latitud (direcci´on NorteSur) (Capderou, 2005), sin embargo, algunas ventanas llegan a ser de tama˜nos m´as reducidos, del orden de 0·05◦en longitud y latitud, equivalente a unos (35×35) km., siendo ´estos, casos muy especiales como lo indica Evans (1999). En virtud del escenario que se presenta actualmente en relaci´on al uso de la zona geoestacionaria, y en vista del crecimiento de la demanda requerida, principalmente en el uso de las telecomunicaciones, que adem´as se incrementa si consideramos la mala cobertura de estos sat´elites para zonas de latitud alta, han surgido planteamientos que conducen a indagar si existe alguna otra ´area espacial que permita extender la zona geoestacionaria, y si existe, bajo qu´e condiciones est´a sujeta o condicionada. La alternativa m´as usada para comunicaciones en lugares de latitud alta han sido las ´orbitas de tipo Molniya. Estas ´orbitas fueron ideadas en la Uni´on Sovi´etica en el a˜no 1960, aunque la primera misi´on con ´exito fue lanzada en el a˜no 1963. Sus caracter´ısticas fundamentales son: un periodo de medio d´ıa sid´ereo, excentricidad muy alta e inclinaci´on cr´ıtica. El periodo consigue una traza que se repite cada dos vueltas, sobrevolando siempre las mismas zonas de la superficie terrestre. La alta excentricidad indica que el sat´elite se encuentra, durante la mayor parte de su periodo, cerca del apoastro. Finalmente, debido a la inclinaci´on cr´ıtica, la posici´on del apoastro se mantiene estable. Si situamos tres sat´elites en esta ´orbita nos aseguramos de que en todo momento uno de los sat´elites se encuentra muy pr´oximo al apoastro, que supondremos situado por encima del lugar de la Tierra donde estable- viii Introducci´on cer las comunicaciones. Con este m´etodo nos aseguramos una misi´on con un coste de lanzamiento mucho menor que la ´orbita geoestacionaria, aunque se precisan tres sat´elites para una completa cobertura, sin embargo, para su seguimiento se necesita una antena m´ovil, pues la ventana que ocupa el sat´elite, desde el punto de vista del observador, es muy grande, de hecho, un sat´elite determinado termina saliendo de esta ventana aunque siempre habr´a alguno de los tres visible. Este tipo de sat´elite permite avanzar en la idea de que las ´orbitas que repiten la traza pueden constituir una buena soluci´on para las comunicaciones, siempre que se consiga una ventana de visibilidad que sea peque˜na o en la que siempre aparezca alg´un sat´elite de la misi´on. Este tipo de misiones ha sido tambi´en propuesta (La´ınez and Romay, 2009; La´ınez et al., 2009) como una alternativa a las grandes constelaciones, tipo GPS, para dar cobertura precisa a una regi´on peque˜na con un peque˜no n´umero de sat´elites. Otra alternativa, que recientemente ha comenzado a estudiarse, se basa en el uso de fuerzas externas (motores de bajo impulso, velas solares, etc.) para construir ´orbitas ex´oticas, no-keplerianas, que McInnes (1999) llama ´orbitas desplazadas. Estas ´orbitas tienen la peculiaridad de situarse en un plano orbital paralelo al ecuatorial con una rotaci´on sincronizada con el planeta para aparecer, desde la superficie, como un punto estacionario fuera del ecuador. En los trabajos realizados por McKay et al. (2009); Anderson and Macdonald (2010), se estudia la estabilidad de familias de ´orbitas desplazadas y se demuestra que es posible dise˜narlas aplicando una peque˜na fuerza producida por un sistema de propulsi´on de bajo empuje, cuya magnitud sea constante y su direcci´on sea siempre la misma. La posibilidad de utilizar impulsos discretos para conseguir ´orbitas desplazadas aparece en un art´ıculo de Nock (1984). Con objeto de observar los anillos de Saturno desde las sondas Voyager4, Nock propuso el uso de peque˜nos impulsos discretos en intervalos de tiempos cortos que desplazaban la ´orbita, en las proximidades del anillo, pero impidi´endole cruzarlo. Posteriormente McInnes (2011) considera el uso de impulsos discretos, en lugar de usar motores de bajo empuje de forma continua, para generar desplazamientos m´as grandes y poder formar familias de ´orbitas desplazadas no ecuatoriales y circulares. La propuesta de C.R. McInnes consiste en generar una ´orbita desplazada aplicando una serie de maniobras orbitales realizadas a intervalos de tiempo constantes, cada vez que el sat´elite pase por un punto fijo del sistema de referencia planetoc´entri4Las sondas espaciales Voyager 1 y 2 fueron lanzadas en 1977 desde Cabo Ca˜naveral. ix co (que rota con el planeta). De esta forma, la ´orbita est´a formada por arcos cerrados en el sistema planetoc´entrico que se convierten en arcos abiertos en el sistema espacial (inercial). La trayectoria final en el sistema espacial no es exactamente paralela al ecuador, como en el caso de empuje continuo, sino que est´a formada por una especie de “corona” representada por arcos que tienen su v´ertice en un paralelo. Visto desde la superficie, el sat´elite no ocupa un punto fijo, sino que se desplaza describiendo un arco cerrado cuya magnitud puede variar. La obtenci´on de cada uno de los arcos, y la maniobra necesaria para recorrerlos, viene expresada en t´erminos del movimiento relativo a un punto de la trayectoria. Para ello, si llamamos u0a un punto cualquiera de la trayectoria, expresado en el sistema planetoc´entrico en un cierto instante inicial t0, y ua la posici´on del orbitador en otro instante cualquiera, podremos llamar ρ=u−u0, a la posici´on relativa del orbitador respecto de su posici´on inicial. La ecuaci´on diferencial que rige la evoluci´on del vector ρ, expresada en el sistema rotante, ser´a: ¨ ρ+ (2ω×˙ ρ)+(ω×(ω×ρ)) = −∂∇VK ∂uu=u0 ρ−∇VK(u0)−ω×(ω×u0), donde VKes el potencial kepleriano, −∂∇VK ∂uu=u0 es la matriz hessiana evaluada en u0y el vector ω= (0,0, ω), donde ωes la velocidad angular que define el movimiento del planeta en torno a su eje de rotaci´on. Como se ha dicho antes, la ´orbita propuesta por McInnes representa, en el sistema planetoc´entrico, un arco cerrado, por lo que si partimos de un punto u0, el valor inicial de ρser´a cero, mientras que al final del arco ρdebe ser de nuevo igual a cero. Para buscar las ´orbitas, McInnes desarrolla un m´etodo num´erico que parte de la linealizaci´on de la ecuaci´on diferencial anterior y obtiene la expresi´on final en t´erminos de la matriz de transici´on. Este tipo de ´orbitas vistas desde el sistema planetoc´entrico no aparecen como un punto del espacio sino como un arco que se cierra. Este arco depende del punto donde se cierra, al que llamaremos v´ertice del arco y del tiempo que se tarde en recorrerlo. A este tipo de ´orbita le llamaremos ´orbita cuasiestacionaria si desde un lugar de la Tierra el orbitador siempre es visible. En cierto modo, es similar a una ´orbita estacionaria situada en cualquier punto de la Tierra, no solo en el ecuador, pero cuya ventana es en general mucho m´as amplia que la ventana de una ´orbita geoestacionaria. McInnes presenta un m´etodo num´erico para el c´alculo de este tipo de ´orbitas, xIntroducci´on v´alido para el modelo kepleriano, pero no realiza un estudio de c´omo son dichas ´orbitas, ni responde a la pregunta de cu´antas hay. En esta memoria hemos pretendido partir de la idea de McInnes y realizar un estudio profundo de dichas ´orbitas que permita saber cu´antas existen y clasificarlas en funci´on de sus propiedades, as´ı como relacionar estos conceptos y propiedades con los de otro tipo de ´orbitas como las ´orbitas que repiten la traza, etc. Para ello, en lugar de un m´etodo num´erico adaptado para este problema, hemos utilizado un m´etodo cl´asico muy estudiado y contrastado y del que se conocen sus propiedades de forma exhaustiva, el problema de Lambert. Adem´as, hemos extendido las ´orbitas, no solo a modelos orbitales keplerianos, sino a cualquier modelo orbital que considere todo tipo de perturbaciones. El proceso de estudio de dichas ´orbitas nos ha obligado a considerar otros problemas en los que hemos realizado nuevas aportaciones, que complementan las herramientas necesarias para completar nuestro estudio. Una de las nuevas cuestiones abordadas es la extensi´on del problema de Lambert a modelos orbitales perturbados. Por otro lado, con vistas a la evaluaci´on del modelo de fuerzas, hemos desarrollado un nuevo m´etodo de evaluaci´on de las derivadas, basado en el m´etodo de diferenciaci´on autom´atica. Este m´etodo, que ha sido publicado recientemente (Abad and Lacruz, 2013), permite evaluar las derivadas de cualquier orden del potencial, de una forma muy r´apida, lo que nos da las componentes de la fuerza de un modelo completo de potencial planetario, as´ı como la matriz hessiana, ´util para la integraci´on de las ecuaciones variacionales, sin necesidad de programar largas y complejas expresiones anal´ıticas. Una gran parte del trabajo realizado para completar esta memoria ha sido la creaci´on del software necesario para la implementaci´on de los m´etodos desarrollados, as´ı como su aplicaci´on sistem´atica para reproducir los ejemplos y casos mostrados. Algunos de los algoritmos programados dentro de este software est´an contenidos en muchos otras aplicaciones o librer´ıas, sin embargo, debido en unas ocasiones a su alto coste por ser programas profesionales, en otras a su “desconocida” fiabilidad para librer´ıas de dominio p´ublico y al lenguaje de programaci´on que, en ocasiones, resulta inadecuado para nuestros prop´ositos, hemos preferido, salvo en caso del m´etodo de integraci´on num´erica, desarrollar nuestro propio software, que ha sido escrito en lenguaje C, y al que hemos llamado OrbitsC. OrbitsC utiliza el integrador dopri85, desarrollado por Hairer et al. (1993), para integrar las ecuaciones del movimiento de un sat´elite artificial, tanto en un sistema inercial como en un sistema rotante. Adem´as, integra las ecuaciones variacionales, 5http://www.unige.ch/ hairer/software.html Movimiento kepleriano y orbital 3 de la ´orbita, con mOymPlas masas respectivas del cuerpo central y el orbitador. Si introducimos el vector velocidad Xpodemos expresar estas ecuaciones como un conjunto de ecuaciones diferenciales de orden uno en la forma: ˙x=X,˙ X=Fk.(1.3) F´acilmente puede demostrarse que el movimiento del orbitador, resultante de la integraci´on de estas ecuaciones, cumple las leyes de Kepler, (Abad, 2012). 1.1.2. Movimiento orbital En la realidad, adem´as de la fuerza kepleriana, aparecen una serie de perturbaciones que modifican las ecuaciones y la trayectoria de la ´orbita, entre las que podemos destacar: Potencial gravitatorio del cuerpo central debido a la forma no esf´erica de los cuerpos celestes. Atracci´on gravitacional de otros cuerpos del sistema distintos del central: la Luna y el Sol en el caso de ´orbitas terrestre, la Tierra para ´orbitas lunares, etc. Rozamiento o resistencia de la atm´osfera en aquellos cuerpos que la poseen. Presi´on de radiaci´on solar. Efectos de mareas, relatividad y otros. Las anteriores perturbaciones pueden formularse por medio de una fuerza adicional Pque modifica las ecuaciones del movimiento (1.3), en la forma: ˙x=X,˙ X=F=Fk+P.(1.4) El movimiento que se deduce de dichas ecuaciones ser´a llamado movimiento orbital y coincide con el movimiento kepleriano cuando el vector perturbaci´on vale cero, P= 0. En general, se verifica la relaci´on kPk<< kFkk, por lo que la descripci´on y propiedades del movimiento orbital son muy pr´oximas a las del movimiento kepleriano, sirviendo ´este como primera aproximaci´on para comprender y estudiar el movimiento de los cuerpos en el espacio. Por ello, en las secciones que siguen ser´a estudiado en detalle el concepto de ´orbita kepleriana. 4Movimiento orbital 1.1.3. Orbitas keplerianas Llamaremos ´orbita kepleriana y la denotaremos con el s´ımbolo O, a la soluci´on de las ecuaciones del problema kepleriano (1.3) para unas condiciones iniciales dadas. Entenderemos por ´orbita, no solo la trayectoria del orbitador, sino todos sus par´ametros, tanto est´aticos o constantes, como din´amicos o variables. Las ecuaciones del problema kepleriano (1.3) constituyen un sistema de seis ecuaciones diferenciales de orden uno. De acuerdo con la teor´ıa de ecuaciones diferenciales ordinarias una soluci´on de dicho sistema vendr´a dada como x=x(t, C), donde C= (C1, C2, C3, C4, C5, C6) representa un vector de seis constantes independientes que llamaremos variables de estado porque permiten determinar cualquier par´ametro de la ´orbita en cualquier instante, es decir, caracterizan la ´orbita. Los seis elementos que componen las variables de estado son constantes de la ´orbita o variables din´amicas particularizadas para un instante dado. En este ´ultimo caso hay que dar el valor de ´estas as´ı como el instante t0en que han sido calculadas. Una vez determinado el conjunto de variables de estado, la ´orbita quedar´a caracterizada por ´este y, pondremos O(C) si los elementos del vector de estado son constantes de la ´orbita y O(t0,C) si son variables particularizadas en t0. Las variables de estado pueden ser elegidas de diversas maneras. La m´as natural, desde el punto de vista de las ecuaciones diferenciales, es a trav´es de los valores del vector de posici´on, x0, y velocidad, X0, para un instante dado. Al vector de dimensi´on seis compuesto por las componentes de los vectores x0yX0se le llama vector de estado. De esta forma una ´orbita kepleriana podr´a ser representada como O(t0,x0,X0). Cada aspecto y propiedad de una ´orbita kepleriana Opuede ser representado por un par´ametro orbital ovariable din´amica. Estos par´ametros pueden ser constantes, como la excentricidad eo la norma del momento angular Go variables como el vector de posici´on xo la anomal´ıa verdadera f. Una vez conocidos los seis o siete elementos que caracterizan la ´orbita, ´esta queda completamente determinada junto con todos sus par´ametros. En el caso de que un par´ametro, que de forma gen´erica llamaremos σ, sea constante, utilizaremos la notaci´on σ(O), pues este par´ametro solo depende de la ´orbita, sin embargo, cuando el par´ametro sea variable depender´a a su vez del instante ten que sea calculado, por lo que pondremos σ(t, O). Movimiento kepleriano y orbital 5 Para representar una ´orbita de forma mucho m´as descriptiva, tanto desde un punto de vista geom´etrico, como astron´omico y din´amico, se utilizan unas variables de estado que se adaptan completamente a la descripci´on del movimiento dada por las leyes de Kepler: los elementos orbitales. En primer lugar tomaremos los dos elementos que caracterizan la forma de la c´onica, esto es, el semieje mayor a(o el semilado recto4p) y la excentricidad e, que caracterizan la forma y dimensiones de la c´onica. Para completar la informaci´on sobre la trayectoria necesitaremos situarla en el espacio, para lo cual basta observar la figura 1.1 y recordar que la ´orbita est´a contenida en un plano perpendicular al vector momento angular G=x×Xo lo que es igual a su direcci´on n. Supondremos, por ahora, que la ´orbita no coincide con el plano Oxy del sistema espacial, esto es, que n×e36= 0. O P x A l ua n e1 e2 e3 Ω ω f i i Figura 1.1: ´ Orbita kepleriana en el espacio. Puesto que el plano de la ´orbita y el plano fundamental del sistema espacial Oxy no son paralelos, necesariamente se cortar´an en una recta que pasa por Oy pertenece a ambos planos y que llamaremos l´ınea de los nodos. Tomaremos como direcci´on positiva de dicha recta la que contiene el nodo ascendente, o punto de la ´orbita en el que el orbitador pasa de coordenadas znegativas a positivas. El vector unitario ldefine la l´ınea de los nodos y forma un ´angulo Ω, ´angulo del nodo, con e1. El ´angulo Ω puede tomar cualquier valor entre 0 y 2π. El ´angulo que forman el vector ncon e3ser´a llamado inclinaci´on, y denotado 4El semieje mayor no est´a definido para la par´abola. 6Movimiento orbital por i, y representa tambi´en el ´angulo entre el plano Oxy y el de la ´orbita. El ´angulo ipuede tomar un valor cualquiera entre 0 y π. El vector nrepresenta tambi´en el sentido de la rotaci´on de la part´ıcula alrededor del eje definido por n, pues debido a su definici´on ´esta tiene siempre lugar en sentido contrario a las agujas del reloj si se observa desde el extremo de n. As´ı pues, el ´angulo que forma ncon e3indica tambi´en el sentido de giro observado desde un punto cualquiera de la parte positiva del eje Oz. Un ´angulo ientre 0 y π/2 indicar´a una ´orbita directa (sentido de giro contrario a las agujas del reloj), mientras que una inclinaci´on entre π/2 y πindicar´a una ´orbita retr´ograda (sentido de giro igual al de las agujas del reloj). Los dos ´angulos, Ω e irepresentan la posici´on del plano de la ´orbita en el espacio, pero para poder representar con exactitud la forma de la c´onica, hay que situar la direcci´on del eje de la misma dentro de su plano. El eje de la c´onica lleva la direcci´on de la l´ınea de los ´apsides, a, que forma un ´angulo ωcon la l´ınea de los nodos. Dicho ´angulo ser´a llamado argumento del periastro, representa la posici´on relativa de la c´onica en su plano y es la tercera variable angular de la ´orbita. El argumento del periastro toma un valor cualquiera entre 0 y 2π. Se han completado as´ı los cinco elementos que caracterizan la geometr´ıa de la ´orbita, esto es, la forma, dimensiones y situaci´on de la curva que recorre el orbitador. Para completar la caracterizaci´on de la ´orbita bastar´a un elemento que describa la din´amica, o lo que es igual, que nos informe en qu´e punto de la curva o trayectoria se encuentra el orbitador en cada instante t. O P0 P fE r Figura 1.2: Anomal´ıas en el movimiento kepleriano. Para ello recordemos que si se elige un sistema de coordenadas polares, (r, f), en el plano de la ´orbita (ver figura 1.2), con centro en el cuerpo central O, y eje de coordenadas polares la direcci´on del periastro, entonces la ecuaci´on de una c´onica Movimiento kepleriano y orbital 7 nos da la relaci´on r=p 1 + ecos f,(1.5) donde fes llamada anomal´ıa verdadera. Para encontrar la realci´on de fcon tse introducen dos nuevas variables angulares E(anomal´ıa exc´entrica) y `(anomal´ıa media) a trav´es de las relaciones `=E−esen E, y tan f 2=r1 + e 1−etan E 2,(1.6) que permiten obtener fen funci´on de `. Finalmente, para encontrar la posici´on del orbitador en cada instante basta tener en cuenta la relaci´on `=n(t−T), siendo n=pµ/a3el movimiento medio y Tla ´epoca o instante de paso del orbitador por el periastro. Para caracterizar su din´amica basta considerar finalmente la constante5T. Llamaremos elementos orbitales al conjunto de seis constantes (a, e, i, Ω, ω, T). La obtenci´on de los elementos orbitales a partir de las condiciones iniciales (t0,x0,X0) y viceversa, cuya demostraci´on puede verse en Abad (2012), demuestra la equivalencia entre ambos, y en consecuencia los elementos orbitales constituyen un conjunto de variables de estado. Por tanto la ´orbita kepleriana se puede caracterizar como O(a, e, i, Ω, ω, T) o bien como O(t0,x0,X0). 1.1.4. ´ Orbitas Cuando aparezca una peque˜na perturbaci´on en el modelo, en la forma formulada en las ecuaciones (1.4), dejaremos de usar la palabra kepleriana y pasaremos a llamar ´orbita perturbada o simplemente ´orbita a la soluci´on. El peque˜no valor de la fuerza perturbadora frente a la fuerza kepleriana hace que una ´orbita pueda ser considerada como una ´orbita osculatriz instant´aneamente kepleriana, esto es, que sus propiedades se pueden estudiar, para un instante dado, como los de una ´orbita kepleriana, pero de forma que los par´ametros que caracterizan esta ´orbita kepleriana (los elementos orbitales) dejan de ser constantes y var´ıan con el tiempo: a(t), e(t), i(t),Ω(t), ω(t) y T(t). 5Aunque el elemento Tes constante hay que tener en cuenta que, para ´orbitas el´ıpticas, ´este var´ıa de una vuelta a otra aumentando en una cantidad igual al periodo orbital P. 8Movimiento orbital Siguiendo con la notaci´on introducida en el apartado anterior podremos caracterizar una ´orbita bien a partir de sus condiciones iniciales O(t0,x0,X0) (al igual que para ´orbitas keplerianas) o bien a partir del valor de los elementos orbitales para el instante inicial O(t0, a0, e0, i0,Ω0, ω0, T0). 1.2. Sistemas de referencia 1.2.1. Sistemas de referencia espacial y planetoc´entrico Para representar el movimiento de una part´ıcula Pen el espacio es necesario determinar, en cada instante t, el vector OP ∈R3que une un punto fijo O, que se toma como origen, con la posici´on de la part´ıcula en dicho instante. Si elegimos una base ortonormal (s1,s2,s3) de R3, existe un conjunto de tres componentes: y1, y2, y3, de forma que se puede expresar como combinaci´on lineal de los vectores directores de la base, como sigue: OP = 3 X i=0 yisi.(1.7) De aqu´ı en adelante usaremos indistintamente la notaci´on anterior y la expresi´on de un vector dada por la matriz columna 3 ×1 dada por OP = (y1, y2, y3)T. Llamaremos sistema de referencia al conjunto S={O, s1,s2,s3}formado por el origen y la base. En lo que sigue tomaremos siempre el centro de masas del cuerpo central como origen, por lo que no volveremos a considerarlo e identificaremos el sistema de referencia ´unicamente con la base elegida. La misi´on para la que est´a dise˜nado un sat´elite artificial depende, en buena medida, de su posici´on sobre la superficie del cuerpo que orbita, por lo que es imprescindible considerar la forma y movimientos de este cuerpo a la hora de analizar la ´orbita. En general todos los cuerpos del Sistema Solar sobre los que, ahora o en un futuro, se pretende situar sat´elites artificiales tienen aproximadamente la misma forma y movimiento. Podemos considerar un planeta6como un elipsoide de revoluci´on (figura 1.3) que gira con velocidad angular constante ω, alrededor de un eje fijo, perpendicular a un plano fijo que, de forma gen´erica llamaremos ecuador. 6Hablaremos de planeta en forma gen´erica aunque esto se puede extender de forma gen´erica a cualquiera de las grandes lunas de los planetas. Sistemas de referencia 9 p1 p2 p3 λ ψ φ a b P S polo del planeta ecuador del planeta Figura 1.3: Sistema de referencia planetoc´entrico. Para determinar la posici´on de un sat´elite con respecto a la superficie del planeta es necesario establecer un sistema de coordenadas basado en un sistema de referencia que sea fijo con respecto al planeta. Para ello estableceremos el llamado sistema planetoc´entrico,{p1,p2,p3}, donde: p3representa la direcci´on del eje de rotaci´on fijo; el plano formado por p1yp2coincide con el plano del ecuador; la direcci´on de p1representa la direcci´on de un meridiano de referencia llamado primer meridiano omeridiano cero; y finalmente p2=p3×p1. Por otro lado, las ecuaciones (1.3) o (1.4), est´an referidas a un sistema de referencia inercial7,{e1,e2,e3}, que llamaremos sistema espacial8, tal que el plano fundamental, formado por e1ye2, representa el ecuador del planeta, con e1se˜nalando una direcci´on fija en el espacio9. El vector e3representa la direcci´on del polo del planeta, por lo que coincide con p3. De esta manera el vector x(t), que se˜nala la posici´on del orbitador, OP, en cada instante, tiene sus componentes referidas al sistema espacial como: x= 3 X i=0 xiei.(1.8) La relaci´on entre los dos sistemas de referencia viene dada a trav´es de una rotaci´on elemental, respecto al eje Oz, y de ´angulo θ, que en forma vectorial puede ponerse como p1= cos θe1+ sen θe2,p2=−sen θe1+ cos θe2,p3=e3,(1.9) donde el ´angulo θentre los vectores e1yp1, viene dado por la expresi´on θ=ωt+θ0,(1.10) 7Ver el siguiente apartado. 8Si consideramos ´orbitas alrededor del Sol el sistema espacial coincide con el sistema ecl´ıptico 9La direcci´on de e1es la intersecci´on del ecuador celeste con el ecuador del planeta. Para la Tierra, es el equinoccio γ. 10 Movimiento orbital siendo ωla velocidad angular constante del planeta y θ0el ´angulo entre e1yp1en el instante inicial. El ´angulo θ, que var´ıa de 0 a 2πen una rotaci´on o d´ıa sid´ereo del planeta, representa el reloj de tiempo sid´ereo del planeta10. En lo que sigue llamaremos u= 3 X i=0 uipi,(1.11) al vector de posici´on del orbitador, OP, referido al sistema planetoc´entrico. 1.2.2. ´ Orbita relativa Las ´orbitas, vistas como la curva x(t) soluci´on de las ecuaciones (1.3) o (1.4), representan una c´onica, que es una curva cerrada. Puesto que esta curva se observa desde el sistema espacial, la llamaremos, de aqu´ı en adelante ´orbita espacial. La parte superior izquierda de la figura 1.4, muestra una de estas ´orbitas donde se ve que al final del periodo orbital la posici´on del orbitador coincide con el punto inicial. Cuando la misma curva se representa en el sistema planetoc´entrico, u(t), la rotaci´on del sistema produce que el punto final de la ´orbita no coincida con el inicial, por lo que la curva no se cierra al final del periodo, como puede verse en la parte derecha de la figura 1.4. A la ´orbita referida al sistema planetoc´entrico le llamaremos ´orbita planetoc´entrica. La proyecci´on de la ´orbita relativa sobre un mapa plano que representa la superficie del planeta, en este caso la Tierra, se le llama traza de la ´orbita. La parte inferior de la misma figura 1.4, muestra la traza de la misma ´orbita espacial donde se ve m´as claramente el hecho de la p´erdida de periodicidad de la ´orbita planetoc´entrica. Debido a lo dicho, es conveniente clarificar el concepto de periodo que ser´a usado a lo largo de esta memoria. En general, cuando hablemos de periodo orbital, nos referimos al periodo de una ´orbita espacial kepleriana, que siempre existe. Para una ´orbita (no kepleriana) el concepto de periodicidad puede generalizarse igual que el resto de par´ametros constantes de la ´orbita que se transforman en par´ametros variables, as´ı podremos hablar del una funci´on P(t) que representa en cada instante el periodo de la ´orbita kepleriana osculatriz. 10Para la Tierra θrepresenta el tiempo sid´ereo medio en Greenwich. Sistemas de referencia 11 Figura 1.4: Izquierda: ´orbita referida al sistema espacial. Derecha: ´orbita relativa o referida al sistema planetoc´entrico. Abajo: traza o proyecci´on de la misma ´orbita sobre un mapa plano de la superficie terrestre. Una ´orbita planetoc´entica puede ser peri´odica o no serlo, pero dicho periodo, si existe, no ser´a igual al periodo orbital correspondiente. La ´unica excepci´on a la afirmaci´on anterior la constituyen las ´orbitas de periodo orbital igual al periodo de rotaci´on del planeta u ´orbitas planetos´ıncronas. Para que una ´orbita planetoc´entrica sea peri´odica tienen que existir dos n´umeros enteros nym, tales que nP −mR = 0, siendo Pel periodo orbital y Rel periodo de rotaci´on del planeta. 12 Movimiento orbital La figura 1.5 representa dos ´orbitas planetoc´entricas diferentes. La de la parte superior izquierda no es peri´odica y puede verse, tanto en la figura tridimensional como en la traza, que al cabo de seis periodos orbitales no se ha cerrado y, adem´as va formando una curva que va llenando todos los puntos de la Tierra entre dos latitudes l´ımite. Por ello, a una ´orbita de este tipo se le llama ´orbita densa. En la parte superior derecha de la figura 1.5 aparece una ´orbita planetoc´entrica peri´odica (en este caso geos´ıncrona), que aunque representa tambi´en seis periodos orbitales siempre se repite, su proyecci´on sobre un mapa plano sobre la superficie de la Tierra se aprecia en la parte inferior derecha de la misma figura. Figura 1.5: Izquierda: la imagen superior corresponde a seis vueltas de una ´orbita densa y su respectiva proyecci´on sobre la superficie en la parte inferior. Derecha: la imagen superior corresponde a seis vueltas de una ´orbita relativa peri´odica y su respectiva proyecci´on equivale a la imagen inferior. Ecuaciones fundamentales 19 Para sistemas conservativos se utiliza habitualmente una formulaci´on hamiltoniana, en la que el movimiento en el sistema rotante viene caracterizado por un hamiltoniano de la forma H(u,U) = 1 2U2−ω·(u×U) + V(u),(1.41) donde las componentes de urepresentan las coordenadas, que corresponden a la posici´on expresada en el sistema rotante, mientras que las componentes del vector Uson sus momentos asociados. Las ecuaciones de Hamilton se expresar´an en la forma u0=U−ω×u, U0=−ω×U−∇uV,(1.42) y, como puede observarse, no coinciden con las ecuaciones (1.35) que describen el movimiento en el sistema relativo. Para comprender esta discrepancia basta observar la primera de las expresiones (1.42) y compararla con (1.31), donde U=u0, lo que lleva a concluir que el vector U, de momentos, coincide con RPE X, esto es la velocidad absoluta expresada en el sistema rotante y no con la velocidad relativa U. 1.3.2. Ecuaciones variacionales La matriz de transici´on de estado, esto es, la variaci´on de las componentes de la soluci´on del sistema con respecto a sus condiciones iniciales resulta de gran inter´es en la correcci´on y determinaci´on de la estabilidad de las ´orbitas peri´odicas. El m´etodo de Lambert generalizado, que desarrollaremos en el cap´ıtulo 3, requiere tambi´en del uso de esta matriz, que ser´a obtenida por integraci´on de las ecuaciones variacionales del problema. En este apartado plantearemos las ecuaciones variacionales del movimiento de un sat´elite. Para ello partiremos de las ecuaciones del problema como ecuaciones (aut´onomas o no) de orden uno: ˙ y=f(t, y),y0=˜ y(t0),y∈R6,(1.43) donde ypuede ser (x,X) o (u,U) seg´un formulemos el problema en el sistema espacial o planetoc´entrico respectivamente. La expresi´on de fser´a distinta en cada caso. 20 Movimiento orbital Para encontrar la variaci´on de la soluci´on de (1.43) con respecto a las condiciones iniciales, derivaremos, con respecto a ellas, la propia ecuaci´on (1.43) obteni´endose la ecuaci´on ˙ yy0 =fy·yy0 ,(1.44) que, tras llamar Φ=yy0 , se transforma en ˙ Φ=fy·Φ,(1.45) y representa la ecuaci´on variacional del sistema (1.43). Para obtener las ecuaciones variacionales en el sistema espacial tendremos en cuenta en primer lugar que Φ=  xx0xX0 Xx0XX0  ,(1.46) mientras que, atendiendo a la expresi´on (1.33), tendremos fy= O I FxFX ,(1.47) siendo OeIlas matrices nulas e identidad de orden 3, FxyFX, respectivamente, las matrices de derivadas del vector fuerza respecto de la posici´on y de la velocidad expresadas en el sistema espacial. Si el sistema es conservativo la fuerza deriva de un potencial V(x) y la matriz fytomar´a la forma fy= O I Vxx O ,(1.48) siendo Vxx la matriz hessiana de V. En el sistema planetoc´entrico tendremos que Φ=  uu0uU0 Uu0UU0  ,(1.49) mientras que, atendiendo a la expresi´on (1.35), tendremos fy= O I Fu−W2FU−2W ,(1.50) Ecuaciones fundamentales 21 siendo, OeIlas matrices nula e identidad de orden 3, FuyFU, respectivamente, las matrices de derivadas del vector fuerza respecto de la posici´on y de la velocidad expresadas en el sistema planetoc´entrico, y Wla matriz (1.22), asociada a la rotaci´on del sistema. Si el sistema es conservativo la fuerza deriva de un potencial V(u) y la matriz fytomar´a la forma fy= O I Vuu −W2−2W ,(1.51) siendo, Vuu la matriz hessiana de V. Las expresiones (1.36) y (1.38) permiten expresar la fuerza y el potencial en el sistema de referencia adecuado cuando la expresi´on de la fuerza o el potencial es calculada en un sistema distinto al de la formulaci´on de las ecuaciones diferenciales o cuando se mezclan fuerzas o potenciales expresados en ambos sistemas. Con las ecuaciones variacionales podemos tener el mismo problema, por lo que debemos poder expresar las matrices FxyFX, en t´erminos de FuyFUy viceversa, as´ı como Vxyx, en t´erminos de Vuyu. Para ello, si atendemos a la primera expresi´on dada en (1.36) obtendremos: Fx=REP Fu·ux+FU·Ux, FX=REP Fu·ux+FU·Ux,(1.52) donde debemos sustituir ux,uX,UxyUXpor sus valores obtenidos por derivaci´on de las relaciones (1.27), (1.30) y (1.31). As´ı, tendremos las siguientes relaciones: ux=RPE ,uX= 0, Ux=−RPE ·W,UX=RPE ,(1.53) por lo que finalmente, podremos poner Fx=REP ·Fu−FU·W·RPE , FX=REP ·FU·RPE ,(1.54) as´ı como, Fu=RPE ·Fx+FX·W·RPE , FU=RPE ·FX·REP .(1.55) 22 Movimiento orbital Cuando las fuerzas derivan de un potencial tendremos las relaciones Vxx =REP ·Vuu ·RPE , Vuu =RPE ·Vxx ·RPE ,(1.56) que se obtienen aplicando las propiedades del hessiano de una funci´on. La expresi´on de fy, dada por las ecuaciones (1.47), (1.48), (1.50) y (1.51), se obtiene de forma anal´ıtica para cada problema. Esto conduce a reformular las ecuaciones variacionales del problema del sat´elite en cada modelo de perturbaci´on. El m´etodo de diferenciaci´on autom´atica desarrollado en el siguiente cap´ıtulo (cap´ıtulo 2), permite el c´alculo num´erico de las derivadas, de cualquier orden, de una funci´on cualquiera sin necesidad de sus expresiones anal´ıticas ni de m´etodos aproximados. De esta forma, la formulaci´on dada en este apartado nos permite obtener autom´aticamente las ecuaciones variacionales para cualquier modelo de perturbaciones. Cap´ıtulo 2 Aplicaci´on de la diferenciaci´on autom´atica a la evaluaci´on de las perturbaciones orbitales 2.1. Diferenciaci´on autom´atica El c´alculo simb´olico de la derivada de una funci´on resulta, en general, una tarea relativamente sencilla, cuya mayor complicaci´on radica en la complejidad de la expresi´on expl´ıcita de la funci´on a derivar y de la propia derivada. Por otro lado, la obtenci´on num´erica de la derivada en un punto, basada en el m´etodo de las diferencias finitas, resulta mucho menos precisa que otros m´etodos num´ericos. La t´ecnica conocida como diferenciaci´on autom´atica permite simplificar notablemente la obtenci´on del valor num´erico de la derivada sin utilizar la expresi´on simb´olica expl´ıcita de ´esta, ni aplicar ninguna t´ecnica num´erica de aproximaci´on. La diferenciaci´on autom´atica, tambi´en conocida como diferenciaci´on computacional o algor´ıtmica es una herramienta mediante la cual, utilizando t´ecnicas propias del c´alculo simb´olico, podemos construir, sin necesidad de aplicar ninguna t´ecnica de aproximaci´on, el valor num´erico de la derivada de una funci´on en un punto. Para ello la funci´on a derivar debe descomponerse en una serie de operaciones y funciones elementales (e.g. suma, producto, sin,cos,log,exp), a las que poder aplicar sistem´aticamente, y de forma num´erica, la regla de la cadena, (Neidinger, 1992; Rall and Corliss, 1996; Tsukanov and Hall, 2003; Griewank and Walther, 2008). 24 Diferenciaci´on autom´atica y c´alculo de perturbaciones Esta t´ecnica puede tambi´en aplicarse al c´alculo de derivadas parciales y de esta forma, aplic´andola junto con propiedades sencillas de c´alculo vectorial, extenderla para la obtenci´on del gradiente de una funci´on f(x) de varias variables, con x∈ Rn. Para ello, llamando ∇xfal vector de derivadas parciales de frespecto de x, bastar´a extender las reglas de derivaci´on basadas en la regla de la cadena al gradiente de una funci´on, en la forma ∇x(u v) = u∇xv+v∇xu, ∇x(u+v) = ∇xu+∇xv, ∇x(αu) = α∇xu, ∇x(uα) = α uα−1∇xu, ∇x(sin u) = cos u∇xsin u, ∇x(cos u) = −sin u∇xcos u, (2.1) donde, αes un n´umero real constante, uyvfunciones de x. Para poder evaluar el gradiente es necesario, al igual que para el c´alculo de una derivada, la descomposici´on de la funci´on en funciones y operaciones simples, unarias y binarias, para las que tengamos alguna relaci´on como las dadas en (2.1). A fin de ilustrar como opera en m´etodo de diferenciaci´on autom´atica para la evaluaci´on de una funci´on y su gradiente, consideremos la funci´on f(x, y, z), con x, y, z ∈R, una funci´on f:R3−→ R, dada por la expresi´on: f(x, y, z) = x2+z√y+ 5z2 sin x2+ 3 .(2.2) En la parte izquierda del esquema (2.3) vemos una sucesi´on de expresiones si,con i= 1, . . . ,14 que sucesivamente evaluadas, partiendo de los valores de s1=x, s2=y, s3=z, permiten llegar finalmente a la expresi´on de s14, que representa la funci´on f. Observemos que dicho esquema de evaluaci´on es ´optimo en el sentido de que no se repite ninguna operaci´on. Por ejemplo, en la funci´on aparece dos veces el t´ermino x2, una en el numerador y otra en el denominador, en el esquema de evaluaci´on este t´ermino se calcula una ´unica vez en s4y, su valor se usa posteriormente en el c´alculo del denominador, s10, y el numerador, s12. La obtenci´on de un esquema de c´alculo eficiente, que efect´ue el n´umero m´ınimo de operaciones, constituye una parte muy importante del esquema de evaluaci´on. El siguiente esquema muestra como puede realizarse la evaluaci´on de la funci´on, Diferenciaci´on autom´atica 25 as´ı como la evaluaci´on de su gradiente. s1←x, ∇xs1←(1,0,0), s2←y, ∇xs2←(0,1,0), s3←z, ∇xs3←(0,0,1), s4=x2←s1s1,∇xs4← ∇x(s1s1), s5=√y←s(1/2) 2,∇xs5← ∇x(s(1/2) 2), s6=z2←s3s3,∇xs6← ∇x(s3s3), s7= 5z2←5s6,∇xs7← ∇x(5s6), s8=z√y←s3s5,∇xs8← ∇x(s3s5), s9= sin x2←sin s4,∇xs9← ∇x(sin s4), s10 = sin x2+ 3 ←s9+ 3,∇xs10 ← ∇x(s9+ 3), s11 = 1/(sin x2+ 3) ←s(−1) 10 ,∇xs11 ← ∇x(s(−1) 10 ), s12 =x2+z√y←s4+s8,∇xs12 ← ∇x(s4+s8), s13 =x2+z√y+ 5z2←s8+s7,∇xs13 ← ∇x(s8+s7), s14 = (x2+z√y+ 5z2)/(sin x2+ 3) ←s13s11,∇xs14 ← ∇x(s13s11). (2.3) Si en lugar de comenzar con s1=x, s2=y, s3=z, comenzamos con valores num´ericos en lugar de (x, y, z) el resultado final ser´a el valor de la funci´on fevaluada para esos valores num´ericos. Este esquema, que permite evaluar esta funci´on de forma eficiente, puede extenderse a la evaluaci´on de cualquier operador del que conozcamos las reglas de evaluaci´on para las funciones elementales en que se descompone f. Las expresiones (2.1) extienden a este caso el operador gradiente, lo que permite calcular el valor del gradiente de la funci´on sin m´as que aplicar el mismo esquema con este operador. Para ello, como se ve en la parte derecha de (2.3) partiremos del valor del gradiente de las variables, esto es ∇xx= (1,0,0),∇xy= (0,1,0),∇xz= (0,0,1), y continuaremos el mismo esquema de evaluaci´on hasta llegar al gradiente de la funci´on. En este caso las relaciones (2.1) dependen tanto del gradiente como del valor de las funciones por lo que para obtener el resultado final debe calcularse, en primer lugar, la parte izquierda de (2.3), y con estos valores evaluar, posteriormente la parte derecha. En lo que sigue generalizaremos este proceso para extenderlo al c´alculo de las derivadas parciales de cualquier orden y, al contrario de lo que se ha explicado hasta aqu´ı, realizar en cada paso, para cada si, el c´alculo simult´aneo de la funci´on siy todas sus derivadas hasta el orden deseado. Para ello en primer lugar introduciremos la notaci´on que nos ayudar´a en este proceso. 26 Diferenciaci´on autom´atica y c´alculo de perturbaciones Sea f(x) una funci´on diferenciable, donde x∈Rn, entonces usaremos el vector i= (i1, i2, . . . , in)∈Nn 0, para denotar el ´ındice que representa la derivada parcial : fi=∂O(i)f ∂xi1 1∂xi2 2. . . ∂xin n ,(2.4) donde, O(i) = i1+i2+. . .+in, corresponde al orden de la derivada. De esta forma, el conjunto de los ´ındices de todas las derivadas hasta el orden o∈N, quedar´a definido por I(o) = {i= (i1, i2, . . . , in)|O(i)≤o, 0≤ij≤o}.(2.5) Podemos definir en este conjunto por una relaci´on de orden total, tal que: i≺j⇐⇒ (O(i)<O(j), O(i) = O(j), ik=jk, k = 0, . . . , m, im+1 < jm+1. Esta relaci´on de orden nos permite identificar cualquier vector icon un n´umero entero entre 0 y card(I(o))1. Dados dos elementos, i,j, o ´ındices de I(o), diremos, de aqu´ı en adelante que: •i−j, es el ´ındice correspondiente al vector (i1−j1, i2−j2, . . . , in−jn). •i≤j⇐⇒ ik≤jk, k = 1,2, . . . , n. Esto representa una nueva relaci´on de orden, usada en la proposici´on 2.1, que da el conjunto de las derivadas parciales necesarias para evaluar una en particular. •i j=i1 j1i2 j2. . . in jn. •i∗=i−1k, donde 1k, es un vector con todos sus elementos igual a 0, excepto el k-´esimo elemento que es igual a 1. kpuede ser cualquiera de las componentes distintas de cero del vector i. Tsukanov and Hall (2003) muestran que la mejor elecci´on para kes la posici´on del m´ınimo de las componentes idistintas de cero. Con la notaci´on introducida para representar el conjunto de los ´ındices de las derivadas, podemos agrupar todas las derivadas parciales de la funci´on f(x), hasta un orden de derivada o∈N, en un conjunto que las represente. Para ello, definimos el conjunto de vectores de derivadas parciales de una funci´on como: Do(f) = {f0, f1, . . . , fi, . . . , f`}, 1Cardinalidad o n´umero de elementos del conjunto I(o). Diferenciaci´on autom´atica 27 donde, f0corresponde a la funci´on f, y con 1,iy`representando el primer, i-´esimo y ´ultimo elemento del conjunto ordenado I(o). Podemos notar que el sub´ındice `es igual a card(I(o)) −1. As´ı, una vez que tengamos la funci´on descompuesta en sus operaciones elementales podemos aplicar las reglas de diferenciaci´on a cualquier orden de derivadas. De esta manera, tenemos las siguientes proposiciones cuyas demostraciones pueden encontrarse en (Neidinger, 1992; Tsukanov and Hall, 2003; Griewank and Walther, 2008; Abad et al., 2012). Proposici´on 2.1 Supongamos tres funciones U=U(x),V=V(x),W=W(x),x∈ Rn, y α, β ∈R, y Do(U) = {U0,U1, . . . , Ui, . . . , U`}, Do(V) = {V0,V1, . . . , Vi, . . . , V`}, Do(W) = {W0,W1, . . . , Wi, . . . , W`}, sus vectores de derivadas. Entonces tendremos: •Si W=αU+βV, el i-´esimo elemento del vector Do(W)vendr´a dado por Wi=αUi+βVi. •Si W=U V, el i-´esimo elemento del vector Do(W)vendr´a dado por Wi=X v≤ii vUvVi−v. •Si W=Uα, α 6= 0, el i-´esimo elemento del vector Do(W)vendr´a dado por W0=Uα 0, Wi=1 U0X v≤i∗i∗ v(αWvUi−v−UvWi−v). Considerando las mismas condiciones de la proposici´on 2.1, se puede extender las reglas de la diferenciaci´on autom´atica para funciones elementales como el sin y el cos. Para estas funciones tenemos la siguiente proposici´on. 28 Diferenciaci´on autom´atica y c´alculo de perturbaciones Proposici´on 2.2 Considerando las mismas condiciones de la proposici´on 2.1, si V= sin UyW= cos U, entonces el i-´esimo elemento de los vectores Do(V)y Do(W)est´a dado por las siguientes expresiones: V0= sin U0,W0= cos U0, Vi=X v≤i∗i∗ vWvUi−v,Wi=−X v≤i∗i∗ vVvUi−v. En esta proposici´on se observa que el c´alculo de las derivadas de la funciones seno y coseno debe ser simult´aneo, pues no se puede obtener uno sin el otro. En lo que sigue aplicaremos las anteriores proposiciones a funciones de tres o seis variables. En la tabla 2.1 se muestran los conjuntos de´ındices necesarios para realizar todo el proceso de c´alculo de derivadas hasta orden tres para funciones de tres variables. La primera columna presenta los ´ındices de las derivadas correspondientes hasta el tercer orden de la funci´on. La segunda muestra el n´umero entero, asignado por la relaci´on de orden total al vector i. En la tercera, cuarta y quinta aparecen, respectivamente, los ´ındices i∗, los coeficientes {i v|v≤i},y el conjunto de los ´ındices {v|v≤i}. 2.2. Potencial gravitacional de un planeta El potencial gravitacional de un planeta2de masa Mque act´ua sobre un punto material externo (orbitador) definido en un sistema de coordenadas esf´ericas, puede ser representado por arm´onicos esf´ericos (Heiskanen and Moritz, 1967; Torge, 2001; Abad, 2012) a trav´es de la expresi´on V(r, λ, ψ) = −µ rX n≥0rp rnn X m=0 (Cnm cos mλ +Snm sin mλ)Pnm(sin ψ),(2.6) siendo, r, λ yψlas coordenadas planetoc´entricas, o coordenadas polares esf´ericas relativas al sistema planetoc´entrico, µ=GMel par´ametro gravitacional del planeta, definido a partir de la constante de gravitaci´on universal, G, y la masa del planeta, M, rpel radio medio ecuatorial del planeta y finalmente Pnm los polinomios asociados de Legendre. Los ´ındices n, m ∈Nson llamados, respectivamente, grado y orden. 2Podemos extender los resultados de este apartado a cualquier cuerpo s´olido que act´ue como cuerpo central de una ´orbita, p. e. la Luna. Potencial gravitacional de un planeta 35 La evaluaci´on de cada elemento de la expresi´on (2.16) puede efectuarse de forma iterativa. Para ello, tendremos en cuenta, por un lado, la relaci´on ρ0=rp r, ρn=ρn−1ρ0,(2.18) para los t´erminos ρn, mientras que para calcular um,yvmharemos uso de las propiedades de las funciones circulares, que conducen a las siguientes relaciones: u0= 1, v0= 0, u1=u/r, v1=v/r, um=um−1u1−vm−1v1, vm=vm−1u1+um−1v1, m > 1. (2.19) Para completar la evaluaci´on del potencial Vdebemos sumar los t´erminos Vnm, lo que, teniendo en cuenta los valores entre los que var´ıa el doble ´ındice (n, m) puede hacerse de varias formas distintas. Para comprender esto mejor representaremos los t´erminos Vnm en una estructura bidimensional. Con esta representaci´on podemos ver, en primer lugar, que la suma de t´erminos puede realizarse fila a fila (horizontalmente), o lo que es igual sumando cada vez todos los t´erminos del mismo grado. Esta suma, que puede verse en la figura 2.1, se corresponde con la expresi´on dada por el doble sumatorio (2.10) o lo que es igual con la expresi´on (2.20), donde hemos llamado Vnal t´ermino que representa la suma de todos los elementos de grado n. Vn0 . . . V30 V20 V21 V22 V31 V32 V33 . . .. . .. . . . . .... Vn1Vn2Vn3. . . Vnn   - - - N X n=2 Vn, Figura 2.1: Suma por grados (horizontal). Vn= m´ın(n,M) X m=0 Vnm,(2.20) V3=V30 +V31 +. . . La suma puede realizarse tambi´en por columnas (verticalmente), o lo que es igual por ´ordenes, de esta forma tendremos el esquema mostrado en la figura 2.2. La suma por ´ordenes se representa por (2.21), donde hemos llamado Vma la suma de todos los elementos del orden m. Finalmente, podemos sumar cada una de las diagonales en la forma en que se muestra en la figura 2.3. La suma por diagonales se representa por el sumatorio (2.22), donde hemos llamado Vda la suma de la diagonal d. 36 Diferenciaci´on autom´atica y c´alculo de perturbaciones Vn0 . . . V30 V20 V21 V22 V31 V32 V33 . . .. . .. . . . . .... Vn1Vn2Vn3. . . Vnn ? ? ? 6 6 6 M X m=0 Vm, Figura 2.2: Suma por ´ordenes (vertical). Vm= N X n=m´ax(2,m)Vnm,(2.21) V1=V21 +V31 +. . . Vn0 . . . V30 V20 V21 V22 V31 V32 V33 . . .. . .. . . . . .... Vn1Vn2Vn3. . . Vnn N X d=0 Vd, @ @ @ @ @ @ @ @ @R @R @R I I I Figura 2.3: Suma por diagonales. Vd= m´ın(N,M+d) X i≥m´ax(2,d)Vi,i−d,(2.22) V2=V20 +V31 +. . . Cualquiera de estas elecciones permitir´a la evaluaci´on del potencial planetario. Adem´as si somos capaces de implementar cualquiera de ellas de manera que calcular un elemento (fila, columna o diagonal) sea independiente del c´alculo de elementos del mismo tipo, entonces podremos construir un c´odigo que pueda ser paralelizado, con la consiguiente mejora en tiempo de CPU si lo ejecutamos en cualquier ordenador con m´ultiples procesadores y/o n´ucleos. La elecci´on de una u otra forma de sumar los t´erminos Vnm est´a fuertemente relacionada con el m´etodo de iteraci´on elegido para evaluar las funciones derivadas de Legendre Qnm. Seg´un muestran Lundberg and Schutzf (1988) hay muchas f´ormulas de recursi´on distintas para evaluar las estas funciones, sin embargo, imponiendo una serie de restricciones algebraicas, el n´umero recursiones diferentes se reduce a siete, que pueden verse en la tabla 2.3. Los coeficientes constantes αi nm yβi nm, con i={1,2,3,4,5,6,7}que acompa˜nan a los dos t´erminos en la combinaci´on lineal de cada recurrencia dependen exclusivamente de los valores que tomen el grado y del orden. Cada αi nm yβi nm, se muestra en la tabla 2.4 correspondiente a cada formulaci´on de la tabla 2.3. Potencial gravitacional de un planeta 37 Tabla 2.3: F´ormulas de recurrencia para la evaluaci´on de los Qnm, para m<n−1. I: Qnm =α1 nm t Q(n−1)m-β1 nm Q(n−2)m, II: Qnm =α2 nm t Qn(m+1) +β2 nm (t2−1) Qn(m+2), III: Qnm =α3 nm t Q(n−1)m+β3 nm (t2−1) Q(n−1)(m+1), IV: Qnm =α4 nm t Qn(m+1) -β4 nm Q(n−1)(m+1), V: Qnm =α5 nm Q(n−2)m+β5 nm (t2−1) Q(n−1)(m+1), VI: Qnm =α6 nm Q(n−1)(m+1) +β6 nm (t2−1) Qn(m+2), VII: Qnm =α7 nm Qn(m+1)/t +β7 nm (t2−1) Qn(m+1). Tabla 2.4: Coeficientes αi nm yβi nm, correspondientes a los siete esquemas de recurrencia de las funciones Qnm con m<n−1 en la tabla 2.3. I: α1 nm =(2n−1) n−m,β1 nm =(n+m−1) n−m, II: α2 nm =2(m+1) (n−m)(n+m+1) ,β2 nm =1 (n−m)(n+m+1), III: α3 nm =(n+m) (n−m),β3 nm =1 (n−m), IV: α4 nm =1 (n−m),β4 nm =1 (n−m), V: α5 nm =(n+m)(n+m−1) (n−m−1)(n−m),β5 nm =(2n−1) (n−m−1)(n−m), VI: α6 nm =2(m+1) (n−m)(n−m−1),β6 nm =1 (n−m)(n+m+1), VII: α7 nm =(n+m) (n−m),β7 nm =1 (n−m). En todos los casos las recurrencias deben ser inicializadas por medio de las expresiones Qnn = (2n−1) Qn−1,n−1= (2n−1)!!,(2.23) para los t´erminos de la diagonal Qnn, y con las expresiones Qn(n−1) =t Qnn = (2n−1) t Qn−1,n−1,(2.24) para los t´erminos de la subdiagonal Qn(n−1). Observemos que los t´erminos de la diagonal son constantes num´ericas que pueden obtenerse por iteraci´on o bien direc- 38 Diferenciaci´on autom´atica y c´alculo de perturbaciones tamente a partir de la definici´on del doble factorial6, mientras que la subdiagonal se obtiene a partir del t´ermino anterior en la fila o en la columna. De los siete esquemas de recurrencia de la tabla 2.3 eliminamos el esquema VII porque el t´ermino taparece en el denominador de la expresi´on, lo que conducir´a a una singularidad cuando intentemos evaluar la funci´on Qnm en el polo, esto es, para sin ψ= 0. Podemos ver en un gr´afico bidimensional, figuras 2.4 y 2.5, como funcionan los esquemas V y VI. En ellos se observa que para calcular un elemento necesitamos otro un grado y orden menor y otro del mismo grado y dos ´ordenes anteriores (esquema V) y el mismo orden y dos grados anteriores (esquema VII). Qn,m Qn−1,m Qn−2,m Qn,m+1 Qn−1,m+1 Qn−2,m+1 Qn,m+2 Qn−1,m+2 Qn−2,m+2     - Figura 2.4: Esquema V. Qn,m Qn−1,m Qn−2,m Qn,m+1 Qn−1,m+1 Qn−2,m+1 Qn,m+2 Qn−1,m+2 Qn−2,m+2     6 Figura 2.5: Esquema VI. Para los esquemas III y IV las expresiones permiten calcular de forma m´as simple cada elemento disminuyendo en una unidad, en cada caso, el grado o el orden, como puede verse en las figuras 2.6 y 2.7. Qn,m Qn−1,m Qn,m+1 Qn−1,m+1   ? Figura 2.6: Esquema III. Qn,m Qn−1,m Qn,m+1 Qn−1,m+1    Figura 2.7: Esquema IV. En cualquiera de estos casos para calcular un elemento necesitamos otros dos de alguna fila, columna o diagonal anterior (o posterior), de esta forma ,estos esquemas de c´alculo hacen depender el c´alculo de cada fila (columna o diagonal) de otras filas (columnas o diagonales), por lo que el algoritmo no podr´a ser paralelizado. Los anteriores motivos nos permite excluir los esquemas III, IV, V, VI y VII, por lo que nos centraremos ´unicamente en los esquemas I y II. En las figuras 2.8 y 2.9, puede verse que en el esquema I se calcula cada elemento a partir de los anteriores 6(2n−1)!! = 1 3 . . . (2n−1) = (2n−1)(n−1)!!. Potencial gravitacional de un planeta 39 de su misma columna (mismo orden), mientras que en el esquema II cada elemento es calculado a partir de dos posteriores de su misma fila (mismo grado). Qn,m Qn−1,m Qn−2,m Qn,m+1 Qn−1,m+1 Qn−2,m+1 Qn,m+2 Qn−1,m+2 Qn−2,m+2 ? ? Figura 2.8: Esquema I. Qn,m Qn−1,m Qn−2,m Qn,m+1 Qn−1,m+1 Qn−2,m+1 Qn,m+2 Qn−1,m+2 Qn−2,m+2   Figura 2.9: Esquema II. Podemos extender la definici´on de las funciones derivadas de Legendre haciendo Qnm = 0,con n < m. De esta forma, la relaci´on (2.24) puede expresarse en la forma Qn(n−1) = (2n−1) t Qn−1,n−1=α1 n(n−1) t Qn−1,n−1+β1 n(n−1)t Qn−2,n−1, =t Qnn =α2 n,n−1t Qn,n +β2 n(n−1) (t2−1) Qn,n+1, lo que equivale a extender las formulas de recursi´on I y II a todos los pares de ´ındices (n, m) con m<n, e inicializar la recursi´on con los valores Qn−1,n = 0, Qnn = (2n−1)!!,para n≥0. Con esta extensi´on, si combinamos el esquema I, con el m´etodo de suma por orden o el esquema II con el m´etodo de suma por grado podremos conseguir un c´odigo paralelo. No tendremos ninguna posible paralelizaci´on basada en la suma diagonal. Por otro lado, se han realizado numerosos estudios sobre la estabilidad num´erica de las formulaciones I y II. Por un lado, Wiggins and Saito (1971) y posteriormente por Olvers and Smith (1983) realizan dichos estudios aplicados a las recursiones de los polinomios asociados de Legendre. En el art´ıculo de Lundberg and Schutzf (1988) y posteriormente Holmes and Featherstone (2002) demuestran que los esquemas I y II son tambi´en los m´as eficientes desde el punto de vista de su estabilidad num´erica en la evaluaci´on para ´ordenes altos, en particular el primero, por lo que el esquema I de recurrencia ser´a el elegido para aplicar el algoritmo de diferenciaci´on autom´atica y calcular las derivadas del potencial. El factor (n−m)!/(n+m)!, que est´a presente, de forma impl´ıcita, en las expresiones de Qnm, hace que ´este alcance magnitudes muy grandes para valores altos del orden n. Al mismo tiempo los valores de los arm´onicos Cnm ySnm se hacen muy peque˜nos. La necesidad de operar simult´aneamente con valores muy peque˜nos y muy grades en la evaluaci´on de los t´erminos Vnm produce una gran inestabilidad num´erica en el proceso de c´alculo del potencial y sus derivadas. 40 Diferenciaci´on autom´atica y c´alculo de perturbaciones Podemos evitar este fen´omeno introduciendo el factor de normalizaci´on Nnm (Kaula, 1966; Lorrell, 1969) definido por la siguiente expresi´on: Nnm =s(2 −δ0m)(2n+ 1)(n−m)! (n+m)! ,(2.25) donde δ0mrepresenta la delta de Kronecker. Mediante el factor Nnm podemos sustituir los elementos Cnm, Snm y, Qnm por sus expresiones normalizadas definidas a partir de las relaciones (¯ Cnm ¯ Snm )=1 Nnm (Cnm Snm ),¯ Qnm(t) = NnmQnm(t),(2.26) mediante las cuales la expresi´on de Vnm se pondr´a como Vnm =ρn(¯ Cnmum+¯ Snmvm)¯ Qnm(w1).(2.27) Los valores ¯ Cy,¯ Sde los arm´onicos normalizados son los que normalmente presentan los modelos de potencial planetario, por lo que no deben ser transformados para la evaluaci´on de V. El esquema I para el c´alculo de las funciones derivadas normalizadas de Legendre se resume en las expresiones ¯ Qnm(w) =                0, n < m, 1, n =m= 0, ¯γm¯ Qm−1,m−1, n =m > 0, ¯αnm w¯ Qn−1,m(w) + ¯ βnm ¯ Qn−2,m(w), n > m, (2.28) donde, ¯γmesta dado por ¯γm=δmr2m+ 1 2m, δm=(√2, m = 1, 1, m 6= 1,(2.29) y los coeficientes ¯αnm y¯ βnm, que antes llam´abamos ¯α1 nm y¯ β1 nm, pueden expresarse como ¯αnm =s(2n+ 1)(2n−1) (n−m)(n+m), ¯ βnm =−s(2n+ 1)(n+m−1)(n−m−1) (2n−3)(n−m)(n+m). (2.30) Potencial gravitacional de un planeta 41 De esta forma, el esquema de evaluaci´on puede representarse gr´aficamente a trav´es de la figura 2.10. En ´esta partimos de los dos primeros elementos de cada columna (orden) que forman las dos primeras diagonales y que se han inicializado a partir de los valores dados por las tres primeras l´ıneas de la expresi´on (2.28). El resto de valores puede obtenerse columna a columna usando la expresi´on de la cuarta l´ınea de (2.28) que precisa ´unicamente de los dos elementos anteriores de dicha columna. ¯ Qn0 . . . ¯ Q30 ¯ Q20 ¯ Q10 ¯ Q00 ¯ Q−10 ¯ Q11 ¯ Q01 ¯ Q21 ¯ Q22 ¯ Q12 ¯ Q31 ¯ Q32 ¯ Q33 ¯ Q23 . . .. . .. . . . . .... ¯ QN1¯ QN2¯ QN3. . . ¯ QNM ¯ QN−1,M ? ? ? ? ? ? ? ? ? ? ? ? ? ? Figura 2.10: Evaluaci´on de las funciones derivadas normalizadas de Legendre. 2.2.3. C´alculo de derivadas parciales del potencial gravitatorio El esquema de c´alculo de potencial visto en el apartado anterior, junto con las reglas de derivaci´on descritas en el apartado de diferenciaci´on autom´atica, permiten el c´alculo de las derivadas de cualquier orden del potencial planetario para cualquier modelo de potencial. En el presente apartado describiremos las distintas partes del algoritmo que realiza dicho c´alculo, al que hemos llamado GPDC (Gravity Potential Derivatives Calculator), y cuyo c´odigo, escrito en lenguaje C, puede descargarse libremente desde la p´agina web http://gme.unizar.es/software/gpdc. El objetivo de este algoritmo es calcular el vector Do(V) de derivadas del potencial Vtanto respecto a las coordenadas cartesianas como respecto a las coordenadas esf´ericas planetoc´entricas. 42 Diferenciaci´on autom´atica y c´alculo de perturbaciones En nuestro caso tenemos que Do(V) = {V0,V1,...,Vi,...,V`}, siendo Vi=∂O(i)V ∂ui1∂vi2∂wi3,o bien Vi=∂O(i)V ∂ri1∂λi2∂ψi3,(2.31) donde los ´ındices irepresentan los elementos del conjunto ordenado I(o) = {i= (i1, i2, i3)|O(i)≤o, 0≤ij≤o}. Por el momento nos centraremos en el c´alculo de las derivadas de Vcon respecto a las coordenadas cartesianas (u, v, w), expresi´on a la izquierda de (2.31). Para el c´alculo de Do(V) deberemos calcular por separado el vector Do(VK) y el vector Do(VP). La expresi´on del potencial kepleriano, dada por VK=−µ/r, con r=√u2+v2+w2, permite un sencillo esquema para el c´alculo de Do(VK) que se muestra en la parte 1 del algoritmo GPDC. Algoritmo GPDC parte 1: C´alculo de Do(VK), con respecto a las coordenadas (u, v, w). Data:µ, u, v, w. Result:Do(VK), con respecto a (u, v, w) Do(u)←(u, 0,0,1,0,...,0) Do(v)←(v, 0,1,0,0,...,0) Do(w)←(w, 1,0,0,0,...,0) s1← Do(u u) s2← Do(v v) s3← Do(w w) s4← Do(s1+s2+s3) Do(1/r)← Do(s−1/2 4) Do(VK)← Do(−µ(1/r)) Adem´as del elemento Do(VK), con esta parte del algoritmo se han calculado, y almacenado, los elementos Do(u),Do(v),Do(w) y Do(1/r). Partiendo de la constante rpy de Do(1/r) podremos calcular los elementos Do(ρi), con i= 0, . . . , N por medio de la parte 2 del algoritmo. Potencial gravitacional de un planeta 43 Algoritmo GPDC parte 2: C´alculo de Do(ρi). Data:rp,Do(1/r). Result:Do(ρi), i = 0, . . . , N Do(ρ0)← Do(rp(1/r)) for i←1to Ndo Do(ρi)← Do(ρi−1ρ0) end Partiendo de los elementos Do(u),Do(v),Do(w) y Do(1/r) podremos calcular Do(w1), Do(ui) con Do(vi),con i= 0, . . . , M por medio de la parte 3 del algoritmo. Puede observarse en este algoritmo que las constantes iniciales u0= 1 y v0= 0, dadas en (2.19), se representan como dos vectores con todos sus elementos iguales a cero (derivada de una constante igual a cero) excepto el primero que toma el valor de la constante. Algoritmo GPDC parte 3: C´alculo de los t´erminos Do(um),Do(vm). Data:Do(u),Do(v),Do(w),Do(1/r). Result:Do(w1),Do(ui),Do(vi), i = 0, . . . , M Do(u0)←(1,0,0,0,0. . . , 0) Do(v0)←(0,0,0,0,0. . . , 0) Do(u1)← Do(u(1/r)) Do(v1)← Do(w(1/r)) Do(w1)← Do(w(1/r)) for m←2to Mdo Do(s1)← Do(um−1u1) Do(s2)← Do(vm−1v1) Do(um)← Do(s1−s2) Do(s1)← Do(vm−1u1) Do(s2)← Do(um−1v1) Do(vm)← Do(s1+s2) end La parte m´as importante del algoritmo es el proceso de c´alculo de cada sumando Vmde la expresi´on (2.21), para ello partimos de los escalares µ, rp,¯ Cij,¯ Sij,¯ Qmm, αij, βij y los vectores de derivadas Do(ρi),Do(uj),Do(vj) y Do(w) previamente calculados. Luego, a partir de los vectores anteriormente obtenidos, construiremos el procedimiento de c´alculo del vector Do(VP) como se muestra en la parte 4 del algoritmo 44 Diferenciaci´on autom´atica y c´alculo de perturbaciones GPDC. Algoritmo GPDC parte 4: C´alculo de las derivadas del potencial Do(VP). Data:µ, rp,¯ Cnm,¯ Snm,¯ Qmm, αnm, βnm,Do(ρn),Do(um),Do(vm),Do(w). Result:Do(VP) Do(VP)←(0,0,0,0,...,0) for m←0to Mdo Do(Vm)←(0,0,0,0,...,0) for n←m´ax(2, m)to Ndo Do(s1)← Do(w¯ Q(n−1)m) Do(¯ Qnm)← Do(αnm s1) if n6=m+ 1 then Do(s1)← Do(βnm ¯ Q(n−2)m) Do(¯ Qnm)← Do(¯ Qnm +s1) end Do(s1)← Do(¯ Cnm um) Do(s2)← Do(¯ Snm vm) Do(s1)← Do(s1+s2) Do(s1)← Do(s1¯ Qnm) Do(s1)← Do(s1ρn) Do(Vm)← Do(Vm+s1) end Do(VP)← Do(VP+Vm) Do(VP)← Do(−(µ/rp)VP) end Finalmente, una vez calculados Do(VK) y Do(VK) podemos calcular Do(VK) a partir de la parte 5 del algoritmo. Algoritmo GPDC parte 5: C´alculo de Do(V). Data:Do(VK),Do(VP). Result:Do(V) Do(V)← Do(VK+VP) El algoritmo GPDC de c´alculo de las derivadas respecto de las coordenadas cartesianas est´a completo con la uni´on de las partes 1, 2, 3, 4 y 5. Si queremos calcular el conjunto de derivadas Do(V) respecto de las coordenadas polares esf´ericas (r, λ, ψ), parte izquierda de (2.31), bastar´a recordar las relaciones (2.15) y sustituir la parte 1 del algoritmo por la parte 6, que calcula las derivadas del Otras fuerzas 51 Las tres tablas muestran los beneficios de la paralelizaci´on a medida que aumenta el grado de los modelos gravitacionales y el orden de las derivadas. Sin embargo, cuando el n´umero de hilos aumenta, E(p) disminuye. Esto se debe al gran n´umero de variables compartidas necesarias en la paralelizaci´on, las cuales generan esta p´erdida de eficiencia. 2.3. Otras fuerzas 2.3.1. Fuerzas conservativas Fuerza radial Una fuerza de tipo radial puede describirse por medio de la expresi´on αx/r, con r=kxk, y un par´ametro α∈Rque, al igual que en las fuerzas tangencial y normal, puede ser constante o funci´on de t. Una fuerza de este tipo es una fuerza conservativa cuyo potencial ser´a VR(r) = VR(r) = α r, (2.33) donde hemos resaltado el hecho de que dicha fuerza, que depende ´unicamente de la distancia r, puede ser calculada indistintamente en el sistema espacial y en el planetoc´entrico. Para calcular este potencial, y sus derivadas, aplicaremos un algoritmo similar al del c´alculo del potencial kepleriano, visto en apartados anteriores. El esquema para la fuerza radial es el que se muestra en el algoritmo 7. Para calcular Do(VR) usaremos el mismo algoritmo sustituyendo (x, y, z) por (u, v, w). Si el par´ametro αdepende de tser´a necesario a˜nadir una llamada a una funci´on que calcule el valor del par´ametro para un determinado instante. Perturbaci´on de un tercer cuerpo El efecto de un tercer cuerpo sobre la ´orbita de un sat´elite viene expresado, en un sistema espacial, por medio de la funci´on potencial V3b=−µP1 ||xP−x|| −x·xP ||xP||3,(2.34) 52 Diferenciaci´on autom´atica y c´alculo de perturbaciones Algoritmo 7: C´alculo de Do(VR), en funci´on de las coordenadas (x, y, z). Data:α, x, y, z. Result:Do(VR), e.f.c. (x, y, z) Do(x)←(x, 0,0,1,0,...,0) Do(y)←(y, 0,1,0,0,...,0) Do(z)←(z, 1,0,0,0,...,0) s1← Do(x x) s2← Do(y y) s3← Do(z z) s4← Do(s1+s2+s3) Do(r)← Do(s1/2 4) Do(VR)← Do(α(1/r)) donde µP=GmP, siendo mPyxPla masa y posici´on del tercer cuerpo P. Si observamos dicha expresi´on podemos concluir que el primer sumando depende de la distancia entre el sat´elite y el tercer cuerpo, mientras que el segundo es el cociente entre el coseno del ´angulo entre la direcci´on del sat´elite y el tercer cuerpo y la distancia del tercer cuerpo al principal. Teniendo esto en cuenta, podemos concluir que el potencial puede ser calculado indistintamente en el sistema espacial y planetoc´entrico y se obtendr´a el mismo valor V3b=V3b. Para el c´alculo del potencial y sus derivadas usaremos el algoritmo 8, donde hemos llamado (p1, p2, p3) a las componentes del vector xP, mientras que pnrepresenta el valor 1/||xP||. En el caso general las componentes de xPy su norma son funciones de t, por lo que para usar el algoritmo 8 debe calcularse previamente el valor de p1, p2, p3ypn para un instante dado con una funci´on externa al algoritmo. Un caso particular muy ´util e interesante se presenta en aquellos casos en los que se investiga la ´orbita de un sat´elite artificial en torno a la luna de un planeta cuando la rotaci´on y el movimiento orbital de la luna respecto al planeta est´en sincronizados. Esto ocurre en la mayor parte de las lunas de planetas del Sistema Solar, incluida la Luna terrestre que, como es bien sabido, presenta siempre la misma cara a la Tierra debido a que el periodo orbital de la Luna respecto a la Tierra y su periodo de rotaci´on son iguales. Simplificaremos el problema suponiendo, en lo que sigue, que la luna orbita en Otras fuerzas 53 Algoritmo 8: C´alculo de Do(V3b), en funci´on de las coordenadas (x, y, z). Data:µk, p1, p2, p3, pn, x, y, z. Result:Do(V3b), e.f.c. (x, y, z) Do(x)←(x, 0,0,1,0,...,0) Do(y)←(y, 0,1,0,0,...,0) Do(z)←(z, 1,0,0,0,...,0) s1← Do(x p1) s2← Do(y p2) s3← Do(z p3) s4← Do(s1+s2+s3) s5← Do(s4pn) s6← Do(p1−x) s7← Do(p2−y) s8← Do(p3−z) s9← Do(s6s6) s10 ← Do(s7s7) s11 ← Do(s8s8) s12 ← Do(s9+s10 +s11) s13 ← Do(s−1/2 12 ) s14 ← Do(s5−s13) Do(V3b)← Do(µks14) torno al planeta en una ´orbita circular, de radio d, y ecuatorial y elegimos el meridiano cero de la luna suponiendo que el planeta ocupa la posici´on xP=−dp3. Si el vector de posici´on del sat´elite viene dado por la expresi´on x=x1p1+x2p2+ x3p3, entonces la perturbaci´on producida por el planeta sobre el sat´elite en ´orbita lunar podr´a expresarse (en el sistema planetoc´entrico) como V=µ"z d2 −1 px2+y2+ (z−d)2#,(2.35) siendo µ=Gm, con mla masa del planeta. En los cap´ıtulos siguientes presentaremos ejemplos donde se aplica esta expresi´on para el caso de ´orbitas lunares perturbadas por la Tierra para las que denotaremos por V⊕al potencial perturbador. 54 Diferenciaci´on autom´atica y c´alculo de perturbaciones 2.3.2. Fuerzas no conservativas Las fuerzas que se presentan a continuaci´on son fuerzas no conservativas que dependen no solo de la posici´on, sino tambi´en de la velocidad, esto es, dependen de seis variables (x,X) o (u,U). En el caso del potencial gravitacional hemos construido un esquema de c´alculo que permite calcular las derivadas del potencial hasta cualquier orden, en este caso, puesto que las necesidades de c´alculo se extienden ´unicamente a problemas de tipo orbital, nos limitaremos al c´alculo de derivadas hasta orden uno, puesto que para la propagaci´on necesitamos ´unicamente la expresi´on de la fuerza, sin derivar, mientras que para las ecuaciones variacionales son necesarias las derivadas de orden uno de la fuerza. Para ello calcularemos el vector Do(F) de derivadas de una funci´on F, tanto respecto a la posici´on como a la velocidad, esto es Do(F) = {F0,F1,F2,...,F6}, siendo Fi=∂O(i)F ∂xi1∂yi2∂zi3∂Xi4∂Y i5∂Zi6,oFi=∂O(i)F ∂ui1∂vi2∂wi3∂Ui4∂V i5∂Wi6, donde los ´ındices irepresentan los elementos del conjunto ordenado I(o) = {i= (i1, i2, i3, i4, i5, i6)|O(i)≤o, 0≤ij≤o}. En este caso no podemos seguir el esquema de ´ındices dado en la tabla 2.1, sino que tendremos que usar el mostrado en la tabla 2.8. Tabla 2.8: ´ Indices de las derivadas parciales (hasta orden uno) y sus valores asociados para una funci´on de seis variables. (i1, i2, . . . , i6)i i∗{i v|v≤i} {v|v≤i} (1,0,0,0,0,0) 1 0 {1,1} {0,1} (0,1,0,0,0,0) 2 0 {1,1} {0,2} (0,0,1,0,0,0) 3 0 {1,1} {0,3} (0,0,0,1,0,0) 4 0 {1,1} {0,4} (0,0,0,0,1,0) 5 0 {1,1} {0,5} (0,0,0,0,0,1) 6 0 {1,1} {0,6} Otras fuerzas 55 Fuerza tangencial Una fuerza tangencial lleva la direcci´on del vector velocidad. Como se ha visto en el cap´ıtulo anterior existen dos vectores velocidad diferentes: velocidad absoluta X o velocidad relativa U, que representan, respectivamente, la velocidad medida en el sistema espacial o en el sistema planetoc´entrico. De esta forma podremos hablar de la fuerza tangencial medida en el sistema espacial FTE = (F1 TE, F2 TE, F3 TE) = αX vE , vE=||X||, α ∈R,(2.36) y la fuerza tangencial medida en el sistema planetoc´entrico, FTP = (F1 TP ,F2 TP ,F3 TP ) = αU vP , vP=||U||, α ∈R.(2.37) Algoritmo 9: C´alculo de Do(FTE), en funci´on del vector velocidad X. Data:α, X, Y, Z. Result:Do(F1 TE),Do(F2 TE),Do(F3 TE), e.f.c. (X, Y, Z) Do(X)←(X, 0,0,0,1,0,0) Do(Y)←(Y, 0,0,0,0,1,0) Do(Z)←(Z, 0,0,0,0,0,1) s0← Do(X X) s1← Do(Y Y ) s2← Do(Z Z) s3← Do(s0+s1+s2) Do(1/v)← Do(s−1/2 3) s4← Do(α(1/v)) Do(F1 TE)← Do(s4X) Do(F2 TE)← Do(s4Y) Do(F3 TE)← Do(s4Z) Esta fuerza no es conservativa por lo que debemos calcular por separado sus tres componentes. El algoritmo 9 muestra como calcular la fuerza FTE, y sus derivadas hasta cualquier orden por medio de la diferenciaci´on autom´atica. Para calcular la fuerza FTP puede usarse el mismo algoritmo 9 sustituyendo el vector Xpor el vector U. 56 Diferenciaci´on autom´atica y c´alculo de perturbaciones Fuerza normal Una fuerza normal es aquella proporcional a la direcci´on del vector normal x×X, esto es FN=αx×X ||x×X||, α ∈R.(2.38) El vector direcci´on normal representa un ´unico vector, por lo que la relaci´on entre la fuerza normal, calculada en el sistema espacial y la calculada en el sistema planetoc´entrico siguen la misma regla que cualquier otro vector, esto es FN(x,X) = REP ·FN(u,U).(2.39) Para calcular la fuerza FN, que al igual que la fuerza tangencial no es conservativa, utilizaremos el algoritmo 10, que nos da las componentes de la fuerza y todas sus derivadas, hasta el orden deseado. El c´alculo de FNser´a id´entico sustituyendo x,Xpor u,U. Otras fuerzas 57 Algoritmo 10: C´alculo de Do(FN), en funci´on de los vectores de posici´on y velocidad x,X. Data:α, x, y, z, X, Y, Z. Result:Do(F1 N),Do(F2 N),Do(F3 N), e.f.c. (x, y, z, X, Y, Z) Do(x)←(x, 1,0,0,0,0,0) Do(y)←(y, 0,1,0,0,0,0) Do(z)←(z, 0,0,1,0,0,0) Do(X)←(X, 0,0,0,1,0,0) Do(Y)←(Y, 0,0,0,0,1,0) Do(Z)←(Z, 0,0,0,0,0,1) s0← Do(x Y ) s1← Do(X y) s2← Do(x Z) s3← Do(z X) s4← Do(y Z) s5← Do(z Y ) s6← Do(s4−s5) s7← Do(s3−s2) s8← Do(s0−s1) s9← Do(s6s6) s10 ← Do(s7s7) s11 ← Do(s8s8) s12 ← Do(s9+s10 +s11) s13 ← Do(s−1/2 12 ) s14 ← Do(α s13) Do(F1 N)← Do(s14s6) Do(F2 N)← Do(s14s7) Do(F3 N)← Do(s14s8) Cap´ıtulo 3 Arcos orbitales y problema de Lambert 3.1. Arcos keplerianos. Problema cl´asico de Lambert 3.1.1. ´ Orbitas keplerianas que pasan por dos puntos Llamaremos arco kepleriano al segmento de curva x(t), comprendido entre dos instantes de tiempo t∈[t0, t0+T]1, donde x(t) representa la trayectoria de la ´orbita soluci´on del sistema kepleriano (1.3). De acuerdo con lo visto en el cap´ıtulo primero de esta memoria, el problema de la determinaci´on de un arco kepleriano consiste en la integraci´on, num´erica o anal´ıtica, del problema de valores iniciales (posici´on y velocidad) de la ecuaci´on diferencial ordinaria dada por (1.3). Si atendemos a la definici´on de arco kepleriano podemos formular un nuevo problema, distinto del problema de valor inicial anterior, y en el que no conozcamos la posici´on y velocidad sino, ´unicamente, las dos posiciones que representan los extremos del arco. Este problema de contorno puede reformularse a trav´es de la siguiente pregunta: ¿ Cu´antos, y cu´ales son, los arcos keplerianos que unen dos puntos 1De aqu´ı en adelante tomaremos como origen de tiempos t0= 0, esto es, el instante del primer punto del arco. 60 Arcos orbitales y problema de Lambert x1,x2∈R3? Este problema es muy ´util desde el punto de vista astron´omico, pues permite desarrollar m´etodos de determinaci´on de ´orbitas en el Sistema Solar para cuerpos menores como asteroides, cometas, etc. En la actualidad, debido a la necesidad de dise˜no y planificaci´on de misiones espaciales se ha extendido a otro m´as amplio como es el problema de las transferencias orbitales que agrupa problemas como el de llevar una nave de un punto a otro del espacio, los encuentros orbitales entre naves o rendevouz, optimizaci´on de trayectorias, etc. El problema de las transferencias orbitales consiste en conectar dos puntos pero sin fijar el tiempo de tr´ansito. 3.1.2. Transferencias orbitales El problema de las transferencias orbitales tiene dos partes: B´usqueda de todas las ´orbitas keplerianas que conectan dos puntos en el espacio por medio de un arco kepleriano. Elecci´on, entre todas las ´orbitas obtenidas en el punto anterior, de aquellas a las que podemos acceder por medio de una maniobra orbital (o cambio de velocidad por aplicaci´on de un impulso o ∆v) de coste m´ınimo. El segundo problema, que no ser´a considerado aqu´ı, se traduce en un problema de optimizaci´on cuando formulamos la expresi´on matem´atica que traduce el coste de la maniobra a una funci´on objetivo a minimizar. Para buscar las ´orbitas keplerianas que conectan dos posiciones x1yx2, consideraremos que ´estas representan los puntos P1yP2, de manera que se tiene la relaci´on x1=OP1,yx2=OP2, siendo Oel cuerpo central. Los puntos O, P1yP2forman un plano, por lo que cualquier ´orbita kepleriana de transferencia, que pase de P1aP2, debe estar contenida en ese mismo plano. Como puede verse en la figura 3.1, el recorrido para ir de P1aP2puede realizarse, bien recorriendo el trayecto m´as corto (´angulo agudo θ1= cos−1(x1·x2/kx1kkx2k)∈ [0, π], trayecto en negrita en la figura) o bien el m´as largo, θ2= (2π−θ1)∈[π, 2π] (trayecto de trazo discontinuo). En el primer caso el vector nque representa la norma del momento angular y la direcci´on del plano ser´a n1= (x1×x2)/kx1×x2k, mientras que en el segundo caso se tendr´a n=n2= (x2×x1)/kx1×x2k. Arcos keplerianos. Problema cl´asico de Lambert 67 De acuerdo con el teorema de Lambert, la expresi´on (3.6) contiene la informaci´on necesaria para construir la soluci´on al problema de Lambert. Gran parte de los m´etodos cl´asicos de resoluci´on del problema de Lambert se basan en una adecuada formulaci´on de la funci´on Φ de dicha ecuaci´on. As´ı, por ejemplo, Lagrange la expresa en la forma √µT=a3/2[(α−sin α)−(β−sin β)],(3.9) donde introduce los par´ametros α= (φ+ψ) y β= (φ−ψ). En Prussing (1979) encontramos un estudio geom´etrico detallado de αyβ. Por otro lado Gauss logra eliminar el t´ermino cos φde las expresiones dadas en (3.7) y reducir ´estas a la ecuaci´on siguiente: √µT=a3/2(2ψ−sin 2ψ)+2λ∆a1/2sin ψ. (3.10) Los detalles del c´alculo de esta expresi´on se describen en su cl´asica memoria de funciones hipergeom´etricas y sus expansiones en fracciones continuas, publicadas tres a˜nos despu´es de la Theoria Motus. Al ser introducidos los par´ametros constantes lym, que dependen exclusivamente de la geometr´ıa del tri´angulo de transferencia, el intervalo de tiempo Ty la constante de gravitaci´on µ, a trav´es de las expresiones l=(1 −λ)2 4λ, m =µT2 (2λs)3,(3.11) y una nueva variable, y. Gauss transforma la ecuaci´on del tiempo dada en (3.10) en las dos ecuaciones siguientes: y2=m l+ sin2(ψ/2),(3.12) y3−y2=m2ψ−sin 2ψ sin3ψ.(3.13) Gauss introduce, posteriormente otra expresi´on para el par´ametro l, dado por: l=sin2(θ/2) + tan2(2w) cos(θ/2) ,con tan2(2w) = r1rr2 r1 +r2rr2 r1−2 4r2 ,(3.14) el cual mejora los resultados cuando el ´angulo de transferencia es muy peque˜no. Para buscar la soluci´on a las ecuaciones (3.12) y (3.13), basta con evaluar de forma simult´anea las variables yyψ. Para ello, el t´ermino de la derecha de la ecuaci´on c´ubica en (3.13), es expresado (Gauss, 1809), como 2ψ−sin 2ψ sin3ψ=4 3F(3,1; 5/2; x),(3.15) 68 Arcos orbitales y problema de Lambert siendo F, una funci´on hipergeom´etrica con la variable x= sin2(ψ/2). As´ı, las ecuaciones (3.12) y (3.13) son reformuladas como sigue: x=m y2−l, y3−y2−hy =h 9,(3.16) donde h=m/(5/6 + l+ξ), y ξuna funci´on de la variable xdeterminada a partir de la siguiente fracci´on continua : ξ(x) = 2x2 35 1 + 2x 35 − 40x 63 1− 4x 99 1−··· .(3.17) Esto permite, posteriormente, el c´alculo de xde manera sistem´atica hasta que deje de cambiar dentro de los l´ımites de una cierta tolerancia establecida. De aqu´ı, la ecuaci´on c´ubica admite exactamente una ra´ız real que puede ser determinada por un proceso iterativo de sustituciones sucesivas. Las ecuaciones dadas por Gauss pueden generalizarse a cualquier otro tipo de movimiento no el´ıptico si se extiende la definici´on de la variable xen la forma x=                sin21 4(E0−ET),elipse, 0,par´abola, sinh21 4(H0−HT),hip´erbola, (3.18) siendo, H0yHTanomal´ıas hiperb´olicas. Battin et al. (1978) parten del m´etodo desarrollado por el propio Gauss con el objetivo de eliminar la singularidad del m´etodo (cuando θ=π) y acelerar su convergencia. Para ello, combinan las formulaciones de Lagrange y la de Gauss, y alteran la geometr´ıa del problema manteniendo fijos los puntos P1,P2y el tiempo Tde tr´ansito, pero moviendo el foco de atracci´on hacia un punto en la recta perpendicular a la l´ınea que une los dos puntos. El m´etodo de resoluci´on se construye a partir de las expresiones del movimiento kepleriano del problema transformado y su relaci´on con las expresiones del problema original. La figura 3.4 muestra la geometr´ıa de la ´orbita transformada. En ella pueden observarse las siguientes propiedades: el semieje mayor es perpendicular a P1P2; Arcos keplerianos. Problema cl´asico de Lambert 69 el tiempo de transferencia desde el pericentro xop hasta el punto P2es la mitad del intervalo de tiempo T; la distancia del nuevo foco O0hasta P2es (r1+r2)/2; y, finalmente, la anomal´ıa verdadera fqueda definida en t´erminos θ, a trav´es del cos f= (2 λ∆)/(r1+r2). O0 f rorop xop c 2 P2 P1 (r1+r2) 2 Figura 3.4: Geometr´ıa de la ´orbita kepleriana transformada. Con esto, la ecuaci´on del tiempo de transferencia se puede expresar a partir de la ecuaci´on de Kepler aplicada al problema transformado, tal que 1 2rµ a3T=E−eosin(E),(3.19) siendo, Ela anomal´ıa exc´entrica del punto P2en la ´orbita transformada y, eosu excentricidad. De acuerdo con la figura 3.4, el punto medio del arco kepleriano que conecta a los puntos P1yP2est´a indicado por xop, de norma rop =kxopky relacionado con ro, a trav´es de las propiedades geom´etricas del tri´angulo O0P1P2, (Battin and Vaughan, 1984), de forma que rop =1 4r1+r2+ 2√r1r2cos θ 2,(3.20) por lo que, la ecuaci´on (3.19) se puede expresar como 1 2rµ a3T=E−sin E+ (1 −eo) sin E, (3.21) donde, E= (1/2)(ET−E0) = ψ= (1/2)(α−β). Al ser considerada, por un lado, la relaci´on de la anomal´ıa verdadera y la exc´entrica seg´un la ecuaci´on (1.6), y por otro que (1 −eo)=(rop/a) sec2(E/2), e introdu- 70 Arcos orbitales y problema de Lambert ciendo la variable x, conjuntamente con los par´ametros lym, en la forma x= tan2E 2, l = tan2f 2, m =µT2 (2rop)3,(3.22) se llega al factor 2rop a=4x (1 + x)(l+x). Si sustituimos ´este y efectuamos una reordenaci´on de la expresi´on (3.20), ´esta se transforma en las dos ecuaciones siguientes: y2=m (l+x)(1 + x),(3.23) y3−y2=m2E−sin(E) 4 tan3(E/2) ,(3.24) similares a las desarrolladas por Gauss. Battin and Vaughan (1984) reformulan el par´ametro lcomo, l=sin2(θ/2) + tan2(2w) sin2(θ/4) + tan2(2w) + cos(θ/2),(3.25) lo que permite remover la singularidad en θ=π, y mejorar los resultados cuando este ´angulo es peque˜no. Por otro lado, el t´ermino que acompa˜na al par´ametro mde la ecuaci´on c´ubica es expresado en t´erminos de una funci´on hipergeom´etrica equivalente a la dada para el algoritmo de Gauss, tal que 2E−sin E 4 tan3(E/2) =1 2 tan2(E/2) E/2 tan(E/2) −1 1 + tan2(E/2),(3.26) se convierte en 2E−sin E 4 tan3(E/2) =−d dxF(1/2,1; 3/2; −x) = 2x+ξ (1 + x)[4x+ξ(3 + x)],(3.27) al ser introducida la variable xdefinida en la expresi´on (3.22). Con esto, las ecuaciones (3.23) y (3.24) son funciones exclusivamente de xey, justamente como el algoritmo de Gauss, y pueden ser aplicadas para cualquier c´onica extendiendo la definici´on de xen la forma x=                tan21 4(E0−ET),elipse, 0,par´abola, −tanh21 4(H0−HT),hip´erbola. (3.28) Arcos keplerianos. Problema cl´asico de Lambert 71 El valor de xas´ı definido tiene un rango entre -1 y +∞. En un proceso iterativo de sustituciones sucesivas se obtiene la soluci´on para xe y. Para ello, como valor de inicio es suficiente con tomar x=lpara el caso el´ıptico yx= 0, para cualquier otro. Para mejorar la convergencia del m´etodo puede introducirse la variable η, como η=x (√1 + x)2,−1< η < 1,(3.29) y la funci´on ξ(x) y las variables auxiliares h1yh2como ξ(x) = 8√1 + x+ 1 3 + 1 η+ξ(η) ,(3.30) h1=(x+l)2(1 + 3x+ξ(x)) (1 + 2x+l)(4x+ξ(3 + x)),(3.31) h2=m(x−l+ξ(x)) (1 + 2x+l)(4x+ξ(3 + x)),(3.32) y de esta forma obtener unas ecuaciones an´alogas a la ecuaci´on c´ubica obtenida por Gauss, esto es x+1 + l 2=s1−l 22 +m y2,(3.33) y3−y2(1 + h1) = h2.(3.34) Finalmente, se puede obtener el valor del semieje mayor y el semiper´ımetro por a=ms(1 + λ)2 8xy2, p =4r1r2sin2(θ/2) c2po,(3.35) con po=c2(1 + x)2 16 a x . En Vallado (2001), puede encontrarse un seudoc´odigo del algoritmo desarrollado por Battin en el que incluye el c´alculo de la velocidad a partir de los elementos orbitales ayp. 72 Arcos orbitales y problema de Lambert 3.2. Arcos orbitales. Problema de Lambert generalizado 3.2.1. Problema de Lambert para un modelo orbital perturbado Podemos extender el concepto de arco kepleriano si sustituimos el modelo orbital kepleriano, dado por las ecuaciones diferenciales (1.3), por el dado por las ecuaciones perturbadas (1.4). De esta forma llamaremos arco orbital al segmento de curva x(t), comprendido entre dos instantes de tiempo t∈[0,T], donde x(t) representa la trayectoria de la ´orbita soluci´on del sistema (1.4). El problema cl´asico de Lambert, analizado en el apartado anterior, as´ı como cualquiera de sus m´etodos de resoluci´on, buscan el arco kepleriano que conecta dos puntos del espacio, pero para un modelo orbital perturbado este arco no es sino una aproximaci´on del arco orbital que los une. En el siguiente apartado introduciremos un m´etodo que llamaremos m´etodo de Lambert generalizado que mejora la soluci´on cl´asica kepleriana del problema de Lambert transformando el arco kepleriano en el arco orbital. Para ello nos basaremos en una extensi´on del m´etodo de mejora de ´orbitas peri´odicas desarrollado en Abad et al. (2011). 3.2.2. M´etodo de correcci´on de ´orbitas peri´odicas En Abad et al. (2011) se desarrolla un algoritmo corrector de mejora de ´orbitas peri´odicas que usa el m´etodo cl´asico de Newton-Raphson con extensiones basadas en el uso del software TIDES4(Abad et al., 2012), que permite obtener ´orbitas peri´odicas con cualquier precisi´on, y en el m´etodo de resoluci´on de sistemas lineales SVD (descomposici´on en valores singulares), que permite que el sistema lineal que debe resolverse no sea necesariamente cuadrado. En nuestro caso prescindiremos del uso de TIDES porque no necesitamos trabajar en m´ultiple precisi´on y porque TIDES no puede ser adaptado a la integraci´on de alguno de los modelos de perturbaciones que tratamos en esta memoria. 4Taylor series Integrator for Differential EquationS. Arcos orbitales. Problema de Lambert generalizado 73 Supongamos un sistema din´amico aut´onomo dado por el sistema de ecuaciones diferenciales de orden uno ˙ y=f(y),y0=y(0),y∈Rn,(3.36) donde y0representa el vector de condiciones iniciales5. La soluci´on del sistema vendr´a dada por el vector y=y(t;y0), t ∈R,y∈Rn,(3.37) Una ´orbita peri´odica del sistema dado en (3.36) est´a caracterizada por un vector y0de condiciones iniciales dadas y un periodo de tiempo T ∈ R, tal que se cumple la siguiente condici´on: y(T;y0)−y0= 0.(3.38) El problema de buscar ´orbitas peri´odicas consiste en la b´usqueda de un conjunto de valores (T;y0) que verifiquen la condici´on anterior y, por tanto, se convierten en las condiciones iniciales de la ´orbita peri´odica. En Astrodin´amica existen diversos m´etodos que nos dan ´orbitas que aunque no sean peri´odicas se encuentran pr´oximas a una peri´odica, de forma que en lugar de la condici´on (3.38) se cumple y(T;y0)−y0≈0. El m´etodo que aqu´ı se presenta se basa en la b´usqueda de una correcci´on (∆T,∆y0) a la ´orbita preliminar (T,y0), de forma que y(T+ ∆T;y0+ ∆y0)−(y0+ ∆y0),(3.39) sea cero, 0, o un valor aproximado a cero. Este m´etodo ser´a aplicado de forma iterativa hasta el momento en que el valor anterior sea lo suficientemente pr´oximo a cero para considerar la ´orbita peri´odica. Para calcular la correcci´on (∆T,∆y0) desarrollaremos la expresi´on (3.39) en una serie de Taylor multivariable, hasta el primer orden y la igualaremos a cero, esto es y(T;y0) + yt∆T+yy0·∆y0−(y0+ ∆y0) = 0.(3.40) La matriz Φ=yy0 representa matriz de transici´on del sistema, esto es, la matriz de derivadas parciales de la soluci´on respecto de las condiciones iniciales y es la soluci´on de las ecuaciones variacionales del sistema din´amico. Esta matriz, evaluada en (T;y0), es la llamada matriz de monodrom´ıa,M=Φ(T;y0) y juega un papel muy importante en el estudio de las ´orbitas peri´odicas, no solo por su correcci´on sino tambi´en para el estudio de estabilidad. 5Usaremos, sin p´erdida de generalidad, el valor t0= 0. 74 Arcos orbitales y problema de Lambert Por otro lado, el t´ermino yt, equivale al vector columna que determina la derivada de la soluci´on respecto al tiempo, es decir, ˙ y=f(y), y corresponde a la expresi´on de fevaluada en yT=y(T;y), cuando t=T, esto es, f(yT). Por lo tanto, la ecuaci´on (3.40), queda en forma matricial como: f(yT)M−I ∆T ∆y0!= (y0−yT),(3.41) donde, Ies una matriz identidad de orden n, y representa un sistema lineal de n ecuaciones con n+ 1 inc´ognitas. Usualmente, los m´etodos de mejora de ´orbitas a˜naden al sistema de ecuaciones anterior una nueva condici´on (f(y0))T∆y0= 0,(3.42) que, por un lado evita posibles desplazamientos tangentes que impidan la mejora de la ´orbita ya que siguen la direcci´on de la misma, y por otro permiten expresar el sistema en la forma  f(yT)M−I 0 (f(y0))T  ∆T ∆y0 = y0−yT 0 ,(3.43) con lo que el n´umero de ecuaciones es igual al n´umero de inc´ognitas, esto es un sistema lineal de (n+ 1) ×(n+ 1). En ocasiones se desea que la ´orbita peri´odica encontrada mantenga constante un par´ametro6, o vector de par´ametros, definido por la condici´on Q(t;y) = q, por lo que debemos a˜nadir esta condici´on a la condici´on de periodicidad dada en (3.38). As´ı, para mantener la restricci´on desarrollamos hasta el primer orden e igualamos a cero la expresi´on Q(T+ ∆T;y0+ ∆y0)−q,(3.44) lo que nos llevar´a a imponer la condici´on Qt∆T+Qy·∆y0= (q−Q(T;y0)),(3.45) donde, QyyQtest´an evaluadas en (T,y). Entonces, a˜nadiendo las ecuaciones (3.45) al sistema (3.43), obtendremos final6Por ejemplo, la energ´ıa o la constante de Jacobi. Arcos orbitales. Problema de Lambert generalizado 75 mente el sistema:      f(yT)M−I 0 (f(y0))T QtQy      ∆T ∆y0 =     y−yT 0 q−Q(T;y0)     ,(3.46) el cual es un sistema lineal de n+1+decuaciones con n+ 1 inc´ognitas, siendo del n´umero de par´ametros (dimensi´on de Q) a conservar. El sistema lineal que verifica la correcci´on de la ´orbita peri´odica puede expresarse como A α =b, siendo Auna matriz m×n, donde mno tiene por qu´e ser igual a n. De esta forma no podremos asegurar la existencia de una soluci´on del sistema, sin embargo, puesto que el objetivo es corregir la ´orbita iterativamente, en varios pasos, nos conformaremos con encontrar la soluci´on de m´ınima norma, esto es, buscaremos el vector αque hace m´ınima la diferencia d=kAα −bk. El m´etodo de descomposici´on en valores singulares permite factorizar una matriz A, de dimensiones (m×n) como, A=MΣNT,(3.47) donde la matriz M=AATes una matriz ortogonal (m×m), N=ATAuna matriz ortogonal (n×n), y Σuna matriz (m×n) cuyos elementos son todos iguales a cero excepto los primeros elementos de la diagonal Σii =σi≥0, con i≤m´ın(n, m). Los elementos σisatisfacen una ordenaci´on tal que: σ1≥σ2≥ ··· ≥ σn, y son los valores singulares de la matriz A. El rango de la matriz Aes igual al n´umero de valores singulares distintos de cero determinado por el n´umero r≤m´ın(n, m). Demostraci´on y propiedades de la ecuaci´on (3.47) la podemos encontrar en (Stoer and Bulirsch, 1993). Por otro lado, tenemos que el sistema matricial es equivalente a, Aα −b=MΣNTα−b,(3.48) que al ser multiplicado por la transpuesta de M, se transforma en el sistema Σβ−b∗, donde β=NTαyb∗=MTb. Como las matrices ortogonales preservan la norma (Demmel, 1997) tenemos d=kAα −bk=kΣβ−b∗k. Finalmente, puesto que Σes una matriz diagonal, podemos encontrar f´acilmente los elementos del vector βque minimizan a d. La soluci´on vendr´a dada por las expresiones βi=b∗ i/σi, si σi6= 0 y βicualquier valor en otro caso. Por lo tanto, la soluci´on de menor norma vendr´a dada por α=Nβy, la distancia a la soluci´on vendr´a dada por d= (Pi>r(b∗ i)2)1/2), (Demmel, 1997). 76 Arcos orbitales y problema de Lambert 3.2.3. Extensi´on del m´etodo de correcci´on de ´orbitas El m´etodo descrito en la anterior subsecci´on corresponde a la mejora de ´orbitas peri´odicas, sin embargo, puede ser extendido para modificar las condiciones a cumplir por la ´orbita mejorada en dos sentidos diferentes: 1. Fijar alguno de los elementos de las condiciones iniciales. Por ejemplo, que el per´ıodo T, o alguna de las componentes del vector de condiciones iniciales no var´ıe. 2. Cambiar la condici´on de periodicidad por otra diferente. Por ejemplo si no queremos que una ´orbita sea peri´odica sino que su trayectoria se cruce en un punto, basta restringir la condici´on (3.38) a la posici´on y no a la velocidad. El m´etodo de correcci´on busca un conjunto de par´ametros {T ,y0= (y0 1, y0 2,···, y0 n)} que representen el periodo y las condiciones iniciales de una ´orbita peri´odica. Con objeto de fijar indistintamente cualquiera de estos elementos introduciremos el vector extendido ˜ y0∈Rn+1 tal que ˜ y0= (y0 0, y0 1, y0 2,···, y0 n), con y0 0≡ T ∈ R. Llamaremos J ⊆ {0,1, . . . , n}al conjunto de ´ındices de los elementos de ˜ yque pueden variar. Por ejemplo, si queremos que el periodo quede fijo, entonces tendremos J={1, . . . , n}. La cardinalidad7de Jser´a nJ≤(n+ 1), donde el signo igual significa que todos los elementos de ˜ ypueden variar. Para reformular la condici´on (3.38) llamaremos, en primer lugar, I ⊆ {1,2, . . . , n} al conjunto de ´ındices de las componentes a las que afecta la condici´on (3.38). La cardinalidad de Iser´a nI≤n, donde el signo igual significa que la restricci´on afecta a todas las componentes del vector. En el caso de buscar ´orbitas que se crucen en R3tendremos que n= 6, siendo las tres primeras componentes la posici´on y las tres ´ultimas la velocidad, por tanto I={1,2,3}. Finalmente, la condici´on (3.38) puede ser reformulada en la forma yI(˜ y0) = ψI(˜ y0,z),con ψI,z∈RnI,(3.49) siendo, zun vector de par´ametros constantes8. En este caso, puesto que no se busca una ´orbita peri´odica, el elemento y0 0=Tdeja de ser el periodo y se transforma en el 7N´umero de elementos del conjunto. 8La necesidad de introducir este vector de par´ametros se ver´a cuando apliquemos este m´etodo para construir una versi´on del m´etodo de Lambert aplicado a ´orbitas no keplerianas. Arcos orbitales. Problema de Lambert generalizado 83 importante y se han calculado los siguientes par´ametros: Velocidad XK iinicial que debe darse al orbitador en el instante inicial para llegar al punto xfa trav´es de un arco kepleriano, esto es en una ´orbita kepleriana. Esta velocidad se calcula con el m´etodo de Lambert cl´asico. Velocidad XT iinicial que debe darse al orbitador en el instante inicial para llegar al punto xfa trav´es de un arco orbital que sigue la ´orbita de un sat´elite que se mueve de acuerdo con el modelo gravitacional LP165P completo (165 ×165). Esta velocidad se calcula con el m´etodo de Lambert generalizado. Velocidad XK+V⊕ iinicial que debe darse al orbitador en el instante inicial para llegar al punto xfa trav´es de un arco orbital para un modelo orbital que considera la fuerza gravitacional kepleriana de la Luna y la perturbaci´on producida por la atracci´on gravitacional terrestre. Esta velocidad, as´ı como las dos siguientes, se calcula con el m´etodo de Lambert generalizado. El potencial perturbador de la Tierra ser´a tomado a partir de la expresi´on de V⊕dada por la ecuaci´on (2.35), donde la distancia entre la Luna y la Tierra ha sido tomada aproximadamente como 221.37 radios lunares. Velocidad XJ+V⊕ iinicial que debe darse al orbitador en el instante inicial para llegar al punto xfa trav´es de un arco orbital obtenido para un modelo que considera el problema principal de un sat´elite lunar perturbado por la atracci´on gravitacional terrestre. Velocidad XT+V⊕ iinicial que debe darse al orbitador en el instante inicial para llegar al punto xfa trav´es de un arco orbital que sigue la ´orbita de un sat´elite que se mueve de acuerdo con el modelo gravitacional lunar LP165P completo (165 ×165) al que le a˜nadimos la atracci´on gravitacional terrestre. Siguiendo la analog´ıa con los casos anteriores se presenta en la tabla 3.7 las velocidades iniciales obtenidas en los cinco casos, mientras que en la tabla 3.8, podemos apreciar el error en metros obtenido entre el vector final y el vector obtenido de la propagaci´on combinando las distintas condiciones iniciales con los distintos modelos orbitales. De los datos mostrados en la tabla 3.8, vemos que si utilizamos condiciones iniciales obtenidas sin la inclusi´on de la perturbaci´on del tercer cuerpo, al ser propagadas 84 Arcos orbitales y problema de Lambert xi(8.715027353750699, 4.683953615172548, 1.445295040567537) xf(4.200843302825212, 8.370251693456526, 3.503658536291552) XK i(-0.0088885095462878493, 0.014290897509099339, 0.0072824480070416611) XT i( -0.0088885035106045941, 0.014290911545181533, 0.0072824584317568389) XK+V⊕ i(-0.0089524274990118100, 0.014315675017697833, 0.0072920030886850637) XJ+V⊕ i(-0.0089524182594439197, 0.014315682236216251, 0.0072920137460271029) XT+V⊕ i(-0.0089524214748922335, 0.014315689059034071, 0.0072920135224772992) Tabla 3.7: Datos del arco correspondiente a una ´orbita lunar media. Tiempo de tr´ansito T= 342. (Unidades: radios lunares y minutos). C.I. K LP165P K +V⊕J2+V⊕LP165P +V⊕ (xi,XK i) 1.12×10−8121.9 39778.74 39775.09 39777.80 (xi,XT i) 1.33×10−839779.67 39776.02 39778.73 (xi,XK+V⊕ i) 8.35×10−910.31 12.06 (xi,XJ+V⊕ i) 1.72×10−94.32 (xi,XT+V⊕ i) 1.12×10−8 Tabla 3.8: Datos del arco correspondiente a la ´orbita media considerando el potencial gravitacional de la Luna y la Tierra como perturbaci´on de un tercer cuerpo. con los distintos modelos que incorporan esta perturbaci´on, el error llega a ser de kil´ometros, de un orden aproximado a los 39 km. Por el contrario, cuando se obtienen condiciones iniciales considerando el potencial del tercer cuerpo, al propagar con cada modelo que incluyan esta perturbaci´on el error est´a por debajo de los 12 m, en el caso m´as completo y menor para otros modelos incompletos. El mejor resultado lo aportan las propias condiciones iniciales al ser propagadas con su respectivo modelo, al igual que los casos anteriormente expuestos, el error llega a ser de 10−8m. Tras analizar los resultados presentados en este apartado podemos concluir que si aplicamos los resultados obtenidos por medio del m´etodo de Lambert para un modelo de perturbaciones, incompleto la propagaci´on con el modelo completo nos aleja del punto final en una distancia que oscila de metros a kil´ometros dependiendo del modelo. Sin embargo, aplicando el m´etodo generalizado de Lambert conseguimos una aproximaci´on que en el peor de los casos es menos que 10−8m, esto es, suficientemente buena en cualquier caso. Cap´ıtulo 4 Arcos orbitales cerrados Las definiciones de arco kepleriano y arco orbital, estudiados en el cap´ıtulo anterior, est´an asociadas a la soluci´on de las ecuaciones (1.3) y (1.4) y como consecuencia al sistema de referencia espacial. La definici´on de estos conceptos puede ser f´acilmente extendida al sistema de referencia planetoc´entrico si se sustituyen esas ecuaciones diferenciales por las ecuaciones (1.34) o (1.35). A estos arcos, definidos en el sistema planetoc´entrico, les llamaremos arcos relativos. Toda misi´on espacial de un sat´elite artificial est´a dise˜nada para un determinado tipo de observaci´on asociado siempre a un observador u observable situado sobre la superficie del planeta. As´ı, el concepto de periodicidad no se asocia al concepto de periodicidad orbital, sino a la combinaci´on de ´este y de la rotaci´on del planeta, lo que conduce a t´erminos como los de ´orbitas que repiten la traza o, lo que es igual, a´orbitas peri´odicas en el sistema planetoc´entrico. Ejemplos de la importancia de este tipo de ´orbitas nos los dan las ´orbitas planetoc´entricas o planetoestacionarias, ´orbitas de tipo Molniya, ´orbitas congeladas, etc. Si nos restringimos a ´orbitas keplerianas podemos decir que una ´orbita kepleriana que repite la traza no es sino un arco kepleriano relativo en el que al cabo de un tiempo el punto inicial del arco coincide con el punto final del mismo, lo que de ahora en adelante llamaremos arco kepleriano cerrado. La periodicidad de estos arcos en el sistema relativo exige que en el punto final no solo coincida la posici´on sino tambi´en la velocidad del punto. El concepto de arco kepleriano cerrado puede extenderse si atendemos ´unicamente a la posici´on y no a la velocidad. Tal como se ha dicho en la introducci´on de 86 Arcos orbitales cerrados esta memoria, McInnes (2011), introduce este concepto en su b´usqueda de ´orbitas no keplerianas desplazadas. En este cap´ıtulo se redefine el concepto de arco kepleriano cerrado y se propone un nuevo m´etodo de b´usqueda, basado en el m´etodo de Lambert. Adem´as se efect´ua un exhaustivo estudio de sus propiedades. Finalmente, aplicando el m´etodo de Lambert generalizado podemos extender el m´etodo de b´usqueda de arcos keplerianos cerrados a los arcos orbitales cerrados, en los que se considera las perturbaciones del movimiento orbital. 4.1. Arcos keplerianos cerrados Comenzaremos el estudio de los arcos orbitales cerrados restringi´endonos a un modelo orbital kepleriano. Llamaremos arco kepleriano cerrado de periodo Ty v´ertice u, a un segmento de ´orbita relativa kepleriana que verifica la relaci´on uT=u0=u,(4.1) siendo u0yuTlas posiciones del orbitador en los instantes t= 0,yt=T, referidas al sistema planetoc´entrico, respectivamente. La parte izquierda de la figura 4.1 muestra un arco cerrado. Figura 4.1: Izquierda: arco kepleriano cerrado. Derecha: el mismo arco mostrado en el sistema espacial. Un arco kepleriano cerrado est´a caracterizado, al igual que el resto de arcos keplerianos, por dos puntos y el tiempo empleado en recorrerlos. A diferencia de los arcos keplerianos, los dos puntos que lo caracterizan est´an referidos al sistema planetoc´entrico, y por ello no podemos aplicar directamente el m´etodo de Lambert1 1Cualquiera de los m´etodos que resuelven el problema de Lambert est´an formulados en el sistema espacial. Arcos keplerianos cerrados 87 para buscar la ´orbita que lo caracteriza. Sin embargo, puesto que conocemos los instantes de tiempo de ambos puntos podemos obtener las posiciones de los dos puntos en el sistema espacial sin m´as que aplicar la rotaci´on correspondiente, esto es, x0=R3(θ0)u0,xT=R3(θ0+ωT)u0,(4.2) donde, ωyθ0representan, respectivamente, la velocidad angular de rotaci´on del planeta y el ´angulo de rotaci´on del planeta en el instante t= 0. De esta forma, el arco cerrado de la izquierda de la figura 4.1 se transforma en el arco orbital de la derecha de la misma figura, y la b´usqueda de un arco kepleriano cerrado se transforma en la b´usqueda de un arco kepleriano, por lo que el m´etodo de Lambert puede ser aplicado, lo que demuestra, simult´aneamente, que dado un punto del espacio uy un tiempo Texisten ´unicamente dos arcos cerrados (uno con una ´orbita directa y otro con una ´orbita retr´ograda) que tienen como v´ertice dicho punto y se recorren en ese tiempo. En la figura 4.2 se pueden apreciar las dos soluciones, directa y retr´ograda, correspondientes a un arco cerrado terrestre con v´ertice en un punto a una distancia 3.5r⊕, longitud λ= 0◦y latitud geoc´entrica ψ= 45◦y periodo T= 8h. Se muestra en rojo la soluci´on directa y en blanco la soluci´on retr´ograda, tanto en la parte izquierda superior como en la inferior que representa la proyecci´on de ambos arcos sobre la superficie terrestre. Puede observarse en la figura 4.2 que la proyecci´on del arco directo sobre la superficie terrestre est´a confinada en una regi´on peque˜na de la Tierra, mientras que la de la ´orbita retr´ograda recorre todas las longitudes2. Por otro lado, la parte superior derecha de la figura muestra la ´orbita correspondiente al arco directo extendida a periodo de tiempo de varias veces el periodo del arco. Podemos observar que la ´orbita no es peri´odica, esto es, no repite su traza. A lo largo de la ´orbita el arco cerrado se repite m´ultiples veces, sin embargo, los v´ertices de estos arcos no coincidan. 2Para otros arcos se invierte esta propiedad. 88 Arcos orbitales cerrados Figura 4.2: Parte superior izquierda e inferior: arcos cerrados directo (rojo) y retr´ogrado (blanco). Parte superior derecha: extensi´on de la ´orbita del arco cerrado directo. 4.1.1. Coordenadas planetoc´entricas de un punto de la ´orbita relativa Para comenzar a estudiar las propiedades de los arcos keplerianos cerrados, trabajaremos con las coordenadas polares esf´ericas (r, λ, ψ), esto es, con la distancia rentre el orbitador y el centro de masas del cuerpo central, la longitud, λ, y lati- Arcos keplerianos cerrados 89 tud3,ψ, planetoc´entricas del orbitador en lugar de las coordenadas cartesianas. Por ello, en primer lugar buscaremos la relaci´on entre estas coordenadas y los elementos orbitales. Para ello, observemos la figura 4.3, donde Srepresenta la proyecci´on del orbitador sobre la esfera celeste, Nel nodo de la ´orbita y S0la proyecci´on de Ssobre el plano fundamental Oxy. El tri´angulo esf´erico SNS0contiene toda la informaci´on necesaria para determinar las relaciones buscadas si tenemos en cuenta que la distancia angular SS0representa la latitud planetoc´entrica, mientras que la distancia angular NS0 viene dada por la relaci´on θ+λ−Ω, siendo θ(t) el ´angulo de rotaci´on o tiempo sid´ereo del planeta, como se aprecia en la gr´afica de la derecha de la figura 4.3. r Ω ω+f i N S S0 γ ψ θ+λ−Ω γ N S0 G Ω λ θ Figura 4.3: Relaci´on entre las coordenadas planetoc´entricas y los elementos orbitales. Las coordenadas esf´ericas del orbitador respecto al sistema polar-nodal, (Abad, 2012), {l,n×l,n}, corresponde a la terna ordenada (r, θ+λ−Ω, ψ). Con la aplicaci´on de la matriz de rotaci´on dada por R1(i)R3(ω+f), podemos cambiar al sistema orbital {u,v,n}para obtener las siguientes expresiones:    rcos ψcos(θ+λ−Ω) rcos ψsin(θ+λ−Ω) rsin ψ   =   rcos(ω+f) rsin(ω+f) cos i rsin(ω+f) sin i   ,(4.3) 3Habitualmente se usa la latitud planetogr´afica φen lugar de la planetoc´entrica ψ, sin embargo para determinar las propiedades de los arcos cerrados es preferible usar ´esta ´ultima. 90 Arcos orbitales cerrados o equivalentemente a,              r(t) = a(1 −e2) 1 + ecos f(t), λ(t) = Ω −θ(t)−tan−1(cos[ω+f(t)],sin[ω+f(t)] cos i), ψ(t) = sin−1(sin[ω+f(t)] sin i), (4.4) donde, hemos llamado tan−1(c, s) a la funci´on que calcula el ´angulo αtal que cos α=c, sin α=sy hemos tenido en cuenta que tanto θcomo fson funciones de t. 4.1.2. Algunas propiedades de los arcos keplerianos cerrados Usando las coordenadas esf´ericas planetoc´entricas podemos reescribir la condici´on (4.1) en la forma rT−r0= 0, λT−λ0= 0, ψT−ψ0= 0,(4.5) donde el sub´ındice representa el instante en que se eval´ua la coordenada. Llevando las expresiones (4.4) a las igualdades (4.5) podremos encontrar, de forma anal´ıtica, algunas propiedades de las ´orbitas que poseen los arcos keplerianos cerrados. As´ı, por ejemplo, de la igualdad de distancias en las expresiones (4.5) deducimos que rT−r0= 0,⇐⇒ a(1 −e2) 1 + ecos fT−a(1 −e2) 1 + ecos f0 = 0,⇐⇒ e(cos fT−cos f0) = 0, esto, se cumple para, e= 0,ofT= 2π−f0=−f0,(4.6) de donde se deduce que, o bien la ´orbita es circular, o el punto medio del arco cerrado coincide con el periastro o el apoastro. En relaci´on a la condici´on de la latitud en (4.5) y al sustituir por la tercera expresi´on de la ecuaci´on (4.4), para cada instante de tiempo se tiene que ψT−ψ0= 0,⇐⇒ Arcos keplerianos cerrados 91 sin ψT−sin ψ0= 0,⇐⇒ sin(ω+fT) sin i−sin(ω+f0) sin i= 0,⇐⇒ sin i[sin(ω+fT)−sin(ω+f0)] = 0, condici´on, esta ´ultima, que se cumple cuando i= 0,o sin(ω+fT)−sin(ω−f0) = 0.(4.7) Ya que en la ´ultima expresi´on aparece el argumento de periastro y la anomal´ıa verdadera en los instantes de tiempo inicial y final, y sabiendo que fT=−f0, podremos reemplazar la segunda de las expresiones (4.7) por 2 sin f0cos ω= 0, condici´on que se verifica para los valores f0= 0,oω=nπ 2,3π 2o.(4.8) La condici´on f0= 0, es incompatible con la condici´on fT=−f0, puesto que esto conducir´ıa a que el punto inicial y el punto final del arco son el mismo lo que es absurdo, luego la igualdad de latitudes demuestra la condici´on i= 0,oω=nπ 2,3π 2o,(4.9) lo que indica que la ´orbita con un arco kepleriano cerrado, es ecuatorial o bien el argumento del periastro vale π/2o3π/2. De hecho tendremos uno de los dos valores para el arco cerrado con ´orbita directa y el otro para el de ´orbita retr´ograda. Por otro lado, la tercera de las relaciones (4.4) permite poner f0= sin−1hsin ψ0 sin ii−ω, (4.10) que nos da la relaci´on entre la latitud del v´ertice del arco, la inclinaci´on de la ´orbita y el argumento del periastro (uno de los dos posibles). En este punto puede demostrarse que la inclinaci´on de la ´orbita que contiene un arco kepleriano cerrado depende ´unicamente de la latitud, ψ, del v´ertice y del periodo, T, del arco. Para apreciar este hecho, partiremos de las expresiones (4.2) que representan la posici´on de los dos extremos del arco en el sistema espacial. En esta expresi´on tomaremos, para simplificar los c´alculos y sin p´erdida de generalidad, el valor θ0= 0, con lo que obtendremos x0=   rcos ψcos λ rcos ψsin λ rsin ψ   ,yxT=   rcos ψcos(λ+ωT) rcos ψsin(λ+ωT) rsin ψ   .(4.11) 92 Arcos orbitales cerrados La inclinaci´on, i, de la ´orbita que pasa por esos dos puntos vendr´a dada por la expresi´on cos i=n·e3, donde e3representa la direcci´on del eje Oz del sistema espacial y nviene dado por n=x0×xT kx0×xTk.(4.12) Si sustituimos en la expresi´on (4.12) los vectores dados en (4.11) y nos quedamos con la tercera componente de nobtenemos, despu´es de la correspondiente manipulaci´on trigonom´etrica y simplificaci´on de t´erminos, la expresi´on cos i=cos ψcos(ωT) sin ωT 2p3 + 2 cos(ωT) cos2ψ−cos(2ψ),(4.13) que demuestra que la inclinaci´on ies una funci´on i=i(ψ, T). La figura 4.4 muestra, a la izquierda, las curvas i(T) para diferentes valores de ψ. Si buscamos arcos cerrados, æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à à à à à à à à à ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ô ô ô ô ô ô ô ô ô ô ô ô ô ô ô ô ô ô ô ô ô ô 5 10 15 20 50 100 150 ΤHhL iHéL ô Ψ = 85 é ò Ψ = 55 é ì Ψ = 35 é à Ψ = 15 é æ Ψ = 1é 0 10 20 30 40 50 60 2 4 6 8 10 12 ΨHéL ΤHhL Figura 4.4: Izquierda: curvas i(T) para distintos valores de ψ. Derecha: curva T(ψ) para la inclinaci´on cr´ıtica. con una inclinaci´on dada i0, podemos observar que la relaci´on i0=i(ψ, T) determina una funci´on impl´ıcita T=T(ψ;i0) que da el valor del periodo del arco para un v´ertice de latitud ψy la inclinaci´on deseada. La parte derecha de la figura 4.4 muestra esta funci´on T(ψ;ic) particularizada para la inclinaci´on cr´ıtica4ic. En apartados posteriores buscaremos y estudiaremos las propiedades de los arcos cerrados para distintos valores de uyTo lo que es igual, para distintos valores de (r, λ, ψ, T). Si atendemos a las propiedades vistas hasta aqu´ı y a las relaciones (4.4) podemos restringir m´as el rango de b´usqueda de arcos cerrados. En efecto, observando la segunda de las relaciones (4.4) podemos deducir que la longitud es irrelevante en la b´usqueda de arcos cerrados, puesto que al cambiar ´esta 4La inclinaci´on cr´ıtica ic≈63◦.43 se define a partir de la relaci´on 1 −5 cos2ic= 0. Tiene la importante propiedad de que mantiene constante el periastro para el problema principal del sat´elite. ´ Orbitas cuasiestacionarias 99 O O ρ0ρ zi zi aiai ui ui Figura 4.9: Distancia cenital de un punto del arco visto desde: a) (izquierda) el punto cuyo cenit coincide con el v´ertice del arco; b) (derecha) un punto cualquiera del planeta definido por el vector ρ. indican la separaci´on m´as extrema del arco respecto de su cenit. A este ´angulo que representaremos por la letra αle llamaremos ´angulo del arco kepleriano cerrado, y su valor vendr´a dado por α= max(zi) = max cos−1ai·ρ0 kaikkρ0k:ti∈[0,T].(4.16) Es importante resaltar que este ´angulo es completamente distinto al ´angulo de visi´on desde el sat´elite, (´angulo del cono de visibilidad del planeta visto desde el orbitador), que indica el ´area que cubre sobre la superficie del cuerpo central un sat´elite desde una posici´on arbitraria del mismo. Las propiedades y el c´alculo de este ´angulo, no ser´an tomadas en consideraci´on en esta memoria, pero se pueden encontrar en Miani and Agrawal (2011). El ´angulo αpuede definirse para cualquier arco kepleriano cerrado, independientemente de que ´este pueda o no formar una ´orbita cuasiestacionaria. De hecho, un valor del arco menor que π/2 indicar´a que el arco se encuentra siempre sobre el horizonte mientras que un valor mayor indica que habr´a instantes en los que el orbitador no ser´a visible, luego estableceremos como condici´on para generar una ´orbita cuasiestacionaria que α < π/2. Adem´as cuanto m´as peque˜no sea el ´angulo αm´as “cuasiestacionaria” ser´a la ´orbita, en el sentido de la proximidad de todos sus puntos al v´ertice, o tama˜no de la ventana de visibilidad. Dado un v´ertice u0y un periodo de tiempo Texisten dos arcos keplerianos cerrados para estas condiciones, uno directo y otro retr´ogrado. En lo que sigue y con objeto de simplificar el estudio de las ´orbitas cuasiestacionarias elegiremos, de entre estas dos soluciones, la que tenga menor ´angulo α. Adem´as, si este valor es mayor que π/2 diremos que el arco kepleriano correspondiente no puede dar lugar a una 100 Arcos orbitales cerrados ´orbita cuasiestacionaria. Si queremos establecer las condiciones de visibilidad desde un punto cualquiera de la Tierra de vector ρ(ver parte derecha de la figura 4.9) podremos modificar la expresi´on (4.2) sustituyendo ρpor ρ0de forma que tendremos α(ρ) = max cos−1ai·ρ kaikkρk:ai=u(ti)−ρ, ti∈[0,T].(4.17) 4.2.1. ´ Orbitas cuasiestacionarias en la Luna Las propiedades de las ´orbitas cuasiestacionarias dependen no solo de los elementos orbitales de la ´orbita que nos da un arco cerrado, sino tambi´en de la velocidad de rotaci´on del cuerpo central. Para ilustrar estas ´orbitas, hemos elegido un cuerpo central como la Luna, de rotaci´on lenta, lo que nos permitir´a construirlas con un coste de mantenimiento mucho menor que para el caso de la Tierra. Una representaci´on gr´afica de los arcos keplerianos nos permite tener otra perspectiva de c´omo son ´estos al ser propagados en un periodo de tiempo dado. En la figura 4.10, podemos apreciar cinco diferentes arcos keplerianos cerrados para el caso lunar, los cuales poseen como punto com´un el v´ertice de coordenadas polares (10 rL,0◦,40◦). A partir de esta terna ordenada se calculan las dos ´orbitas, directa y retr´ograda, obtenidas con el m´etodo de Lambert, para el arco cerrado con v´ertice en ese punto y con periodos, en d´ıas, iguales a 4.1,7.5,9.1,18.5 y 20.0. El conjunto de soluciones se muestran en la gr´afica superior izquierda. Estas soluciones presentan formas y geometr´ıas similares, tanto la soluci´on directa como la retr´ograda, aunque en cada caso el ´angulo αes distinto. No siempre la soluci´on directa define la ´orbita cuasiestacionaria pues el valor de α, en algunos casos, es menor para la ´orbita retr´ograda. En la gr´afica superior izquierda de la figura 4.10 se presentan una serie de arcos keplerianos cerrados para el v´ertice y periodos Tanteriormente se˜nalados. No todos estos arcos presentan valores del ´angulo α < π/2, esto es, no todos definen ´orbitas cuasiestacionarias. En la figura superior derecha y en la inferior se muestran los arcos que s´ı forman una ´orbita cuasiestacionaria. En la gr´afica inferior de la figura 4.10, vemos como es la traza respectiva de los ´ Orbitas cuasiestacionarias 101 Figura 4.10: Arcos keplerianos cerrados para la Luna. La gr´afica superior izquierda presenta distintos arcos directos y retr´ogrados para el v´ertice (10 rL,0◦,40◦) y distintos valores de T. La gr´afica superior derecha presenta las soluciones (directas y retr´ogradas) cuyos αson menores a 90◦. La gr´afica inferior central corresponde a la traza o proyecci´on sobre la superficie lunar de las ´orbitas cuasiestacionarias de la parte superior derecha. arcos anteriormente indicados. Esta traza muestra la posibilidad de dise˜nar misiones espaciales de cobertura regional para un cuerpo central cualquiera. Una idea m´as precisa de la forma y la geometr´ıa nos la aportan los elementos orbitales, y para los arcos keplerianos cerrados de la figura 4.10 (derecha superior y central inferior), se muestran en la tabla 4.1 los valores respectivos. La primer columna representa el tipo de ´orbita directa, Od, o retr´ograda, Or, luego est´a el periodo orbital del arco Ten d´ıas, el semieje mayor, a, en radios lunares, la excentricidad, e, la inclinaci´on, i, en grados y el ´angulo α, en grados. Podemos apreciar 102 Arcos orbitales cerrados que los arcos keplerianos cerrados son muy exc´entricos y abarcan distintos ´angulos de visibilidad con distintos periodos orbitales. Tabla 4.1: Elementos orbitales de los arcos keplerianos cerrados para el caso de lunar. Od/OrT(d) a(rL)e i (◦)α(◦) Od4.1 15.3574 0.9701 43.2845 44.7570 Od7.5 22.3501 0.9473 52.2107 76.7305 Od9.1 25.2784 0.9395 59.1742 88.8106 Or18.5 39.8646 0.9660 122.1880 86.8859 Or20.0 41.9179 0.9760 128.4420 75.2432 En lo que sigue se efectuar´a un estudio m´as sistem´atico de las propiedades de las ´orbitas cuasiestacionarias para la Luna. Para ello se ha tenido en cuenta el rango de valores siguiente: La distancia rse tomar´a dentro del intervalo [1,50], tomando el radio lunar, rL, como unidad. Aunque para la Luna el radio de influencia de la misma tiene un valor aproximado de 38 rL, se ha elegido el valor 50 como valor m´aximo de la distancia ya que ´este est´a m´as pr´oximo al valor correspondiente al semieje de las ´orbitas selenos´ıncronas. La longitud λno ser´a tenida en cuenta pues, como se ha dicho en el apartado 4.1.2, variar λequivale a variar el ´angulo del nodo Ω en la misma cantidad. La latitud selenoc´entrica ψse toma en el intervalo [0, π/2]. No se consideran valores negativos de la latitud porque de acuerdo con las relaciones del apartado 4.1.2 la ´orbita de un arco kepleriano cerrado con v´ertice en el punto (r, λ, ψ) y periodo Tse corresponde con el de v´ertice (r, λ, −ψ) cambiando directa por retr´ograda (o viceversa) y el argumento del periastro ωpor ω+π. El periodo Tdel arco cerrado se toma entre 0 y el periodo de rotaci´on de la Luna (unos 27.32 d´ıas). Para adaptarnos a este periodo se usar´a el d´ıa como unidad de tiempo. En primer lugar analizaremos, a partir de las figuras 4.11 y 4.12, el valor del ´angulo del arco kepleriano cerrado para distintos valores de r, ψ yT. Estas figuras muestran en el eje Oy el valor de α, y en el eje Ox la distancia r. Cada curva, representada con un s´ımbolo diferente muestra el ´angulo para un tiempo Tdado. Finalmente, cada una de la cuatro figuras representa un valor de la latitud ψ. ´ Orbitas cuasiestacionarias 103 æ æ æ æ æ æ æ à à à à à à à ì ì ì ì ì ì ì ò ò ò ò ò ò ò ô ô ô ô ô ô ô ç ç ç ç ç ç ç 0 10 20 30 40 50 0 20 40 60 80 100 120 140 rHrLL ΑHéL Ψ = 1é æ Τ = 1 d à Τ = 5 d ì Τ = 10 d ò Τ = 15 d ô Τ = 20 d ç Τ = 25 d æ æ æ æ æ æ æ à à à à à à à ì ì ì ì ì ì ì ò ò ò ò ò ò ò ô ô ô ô ô ô ô ç ç ç ç ç ç ç 0 10 20 30 40 50 0 20 40 60 80 100 120 140 rHrLL ΑHéL Ψ = 25 é Figura 4.11: Variaci´on de αvs rpara distintos Ty las latitudes ψ= 1◦yψ= 25◦. æ æ æ æ æ æ æ à à à à à à à ì ì ì ì ì ì ì ò ò ò ò ò ò ò ô ô ô ô ô ô ô ç ç ç ç ç ç ç 0 10 20 30 40 50 0 20 40 60 80 100 rHrLL ΑHéL Ψ = 45 é æ Τ = 1 d à Τ = 5 d ì Τ = 10 d ò Τ = 15 d ô Τ = 20 d ç Τ = 25 d æ æ æ æ æ æ æ à à à à à à à ì ì ì ì ì ì ì ò ò ò ò ò ò ò ô ô ô ô ô ô ô ç ç ç ç ç ç ç 0 10 20 30 40 50 0 2 4 6 8 10 rHrLL ΑHéL Ψ = 85 é Figura 4.12: Variaci´on de αvs rpara distintos Ty las latitudes ψ= 45◦yψ= 85◦. Observemos que para valores bajos de la latitud ψel ´angulo αdisminuye conforme aumentamos la distancia, y ´este va en aumento al incrementar el periodo T. La izquierda de la figura 4.11, que representa una latitud ψ= 1◦, muestra un valor del ´angulo αque var´ıa entre 70◦y 5◦para un periodo T= 15 d´ıas. Para un periodo mayor el ´angulo αaumenta mientras que para un periodo menor disminuye. Para latitudes m´as altas, el comportamiento es diferente. La figura 4.12, que representa latitudes ψ= 45◦yψ= 85◦la curva, para cada valor de Tpresenta una apariencia horizontal. La realidad es que se presenta una disminuci´on de αpara un Tdado pero ´esta disminuci´on es muy peque˜na y no puede apreciarse en la figura. Al contrario que para latitudes peque˜nas αaumenta al disminuir T, y disminuye cuando la latitud decrece. Por ejemplo, para la latitud ψ= 45◦el valor del ´angulo cuando T= 15 d´ıas es de unos 55◦, mientras que para ψ= 85◦con el mismo Tse obtiene un αde unos 6◦. El cambio en el comportamiento entre latitudes altas y bajas puede observarse mejor en la gr´afica de la derecha de la figura 4.11. En esta gr´afica vemos dos comportamientos distintos de las curvas para un Tdado. Para algunas se observa una forma casi horizontal, como para latitudes m´as altas, pero para otros valores la curva 104 Arcos orbitales cerrados de α(r) va disminuyendo hasta que alcanza un comportamiento casi horizontal. La curva para T= 10 d´ıas muestra claramente este fen´omeno. Para comprender mejor ´esto, calculemos y dibujemos las ´orbitas cuasiestacionarias de v´ertice (r, ψ = 25◦), con r= 1,5,15,25,35,45,50, y para un valor T= 10 d´ıas. La figura 4.13 muestra esas ´orbitas en dos y tres dimensiones. Figura 4.13: ´ Orbitas cuasiestacionarias correspondientes a un arco de v´ertice (r, ψ = 25◦), con r= 1,5,15,25,35,45 y 50, para un valor de T= 10 d´ıas. Los dibujos tridimensionales de las ´orbitas se muestran desde dos puntos de vista diferentes. La imagen de la izquierda muestra como las ´orbitas, conforme nos alejamos de la Luna, son cada vez m´as peque˜nas. La imagen de la derecha muestra las mismas ´orbitas desde una perspectiva m´as alineada con el punto de observaci´on. En esta ´ultima observamos que la visibilidad del arco cerrado se va reduciendo hasta que, a partir de un determinado valor, dicha reducci´on de tama˜no se hace mucho m´as moderada. La figura inferior muestra la traza de las ´orbitas que tiene una forma que muestra el mismo comportamiento. ´ Orbitas cuasiestacionarias 105 Despu´es de analizar el ´angulo del arco kepleriano cerrado pasaremos a estudiar el valor del coste del mantenimiento de la ´orbita en t´erminos de ∆v. La figura 4.14 muestra este coste para el caso de latitud ψ= 45◦y para distintos r∈[1,50]. A la izquierda se da el coste absoluto en km/s de cada maniobra. Como puede observarse el coste es muy grande para ´orbitas de v´ertice pr´oximo a la superficie (valores de r peque˜nos) mientras que disminuye para distancias mayores. Esta medida del coste resulta muy poco representativa pues depende del intervalo de tiempo Tque se tarda en efectuar dos maniobras consecutivas. Por esta raz´on se usa en la parte derecha de la misma figura otro par´ametro que representa el promedio del ∆vdgastado por d´ıa. Para esto usaremos las unidades (km/s)/d. Con este par´ametro ∆vd, que nos da el coste, por d´ıa, de la ´orbita y no de cada maniobra, podemos observar de nuevo la misma tendencia: para que el coste sea peque˜no la distancia rdebe ser grande, como se muestra en la gr´afica de la derecha de la figura 4.14. æ æ æ æ æ æ æ à à à à à à à ì ì ì ì ì ì ì ò ò ò ò ò ò ò ô ô ô ô ô ô ô ç ç ç ç ç ç ç 0 10 20 30 40 50 0 1 2 3 4 5 rHrLL DvHkmsL Ψ = 45 é æ Τ = 1 d à Τ = 5 d ì Τ = 10 d ò Τ = 15 d ô Τ = 20 d ç Τ = 25 d æ æ æ æ æ æ æ à à à à à à ì ì ì ì ì ì ì ò ò ò ò ò ò ò ô ô ô ô ô ô ô ç ç ç ç ç ç ç 0 10 20 30 40 50 0.0 0.1 0.2 0.3 0.4 0.5 0.6 rHrLL DvdHkmsLd Ψ = 45 é Figura 4.14: Valor de ∆v(izquierda) y ∆vd(derecha) para distintos valores de ryT, con una latitud ψ= 45◦. La figura 4.15 muestra el coste por d´ıa de la ´orbita para latitudes ψ= 1◦(izquierda) y ψ= 60◦(derecha). En ambos casos se muestra el coste para valores de rmayores que 25 radios lunares y para distintos valores de T. La figura muestra que el coste es menor cuanto menor sea la latitud, cuanto menor sea el periodo Ty cuanto mayor sea la distancia r. Resulta tambi´en de inter´es en el estudio de las ´orbitas cuasiestacionarias analizar el valor de la excentricidad de la ´orbita de la cual el arco es una secci´on. La figura 4.16 muestra la evoluci´on de esta excentricidad en los casos ψ= 1◦(izquierda) yψ= 60◦(derecha). Como puede verse para una latitud baja, la excentricidad toma cualquier valor entre cero y uno. Tiende a un valor cero para distancias, r, muy altas, mientras que se hace muy grande para distancias peque˜nas. Por otro lado, la excentricidad resulta menor cuando el periodo Taumenta. Observemos 106 Arcos orbitales cerrados æ æ æ æ æ æ à à à à à à ì ì ì ì ì ì ò ò ò ò ò ò ô ô ô ô ô ô 25 30 35 40 45 50 0.000 0.002 0.004 0.006 0.008 0.010 0.012 0.014 rHrLL DvdHkmsLd Ψ = 1é æ Τ = 18 d à Τ = 20 d ì Τ = 22 d ò Τ = 24 d ô Τ = 26 d æ æ æ æ æ æ à à à à à à ì ì ì ì ì ì ò ò ò ò ò ò ô ô ô ô ô ô 25 30 35 40 45 50 0.020 0.025 0.030 0.035 0.040 0.045 rHrLL DvdHkmsLd Ψ = 60 é Figura 4.15: Valor de ∆vdpara distintos valores de ryTy latitudes ψ= 1◦, ψ = 60◦. tambi´en, en la parte derecha de esta figura, que al aumentar la latitud no aparecen excentricidades pr´oximas a cero, manteni´endose, para toda distancia y todo periodo valores grandes de la excentricidad æ æ æ æ æ æ æ à à à à à à à ì ì ì ì ì ì ì ò ò ò ò ò ò ò ô ô ô ô ô ô ô ç ç ç ç ç ç ç 0 10 20 30 40 50 0.0 0.2 0.4 0.6 0.8 1.0 rHrLL e Ψ = 1é æ Τ = 1 d à Τ = 5 d ì Τ = 10 d ò Τ = 15 d ô Τ = 20 d ç Τ = 25 d æ æ æ æ æ æ æ à à à à à à à ì ì ì ì ì ì ì ò ò ò ò ò ò ò ô ô ô ô ô ô ô ç ç ç ç ç ç ç 0 10 20 30 40 50 0.75 0.80 0.85 0.90 0.95 1.00 rHrLL e Ψ = 60 é Figura 4.16: Valor de la excentricidad epara distintos valores de ryTy latitudes ψ= 1◦, ψ = 60◦. Finalmente, la figura 4.17 muestra la evoluci´on de la distancia en el periastro para todo tipo de estas ´orbitas. Valores por debajo de la unidad, en la zona sombreada, nos indican ´orbitas bal´ısticas que si se recorren sali´endose del arco cerrado conducen a una colisi´on con la superficie lunar. æ æ æ æ æ æ æ à à à à à à à ì ì ì ì ì ì ì ò ò ò ò ò ò ò ô ô ô ô ô ô ô ç ç ç ç ç ç ç æ æ æ æ æ æ æ à à à à à à à ì ì ì ì ì ì ì ò ò ò ò ò ò ò ô ô ô ô ô ô ô ç ç ç ç ç ç ç 0 10 20 30 40 50 0 2 4 6 8 10 rHrLL rpHrLL Ψ = 60 é æ Τ = 1 d à Τ = 5 d ì Τ = 10 d ò Τ = 15 d ô Τ = 20 d ç Τ = 25 d æ æ æ æ æ æ à à à à à à ì ì ì ì ì ì ò ò ò ò ò ò ô ô ô ô ô ô 25 30 35 40 45 50 0 1 2 3 4 5 rHrLL rpHrLL Ψ = 60 é æ Τ = 18 d à Τ = 20 d ì Τ = 22 d ò Τ = 24 d ô Τ = 26 d Figura 4.17: Valor de la distancia en el periastro rppara distintos valores de ryT, con una latitud ψ= 60◦. La zona sombreada representa ´orbitas bal´ısticas. Arcos orbitales cerrados: mantenimiento de ´orbitas cuasiestacionarias 107 4.3. Arcos orbitales cerrados: mantenimiento de ´orbitas cuasiestacionarias En el cap´ıtulo anterior se ha distinguido entre arcos keplerianos y arcos orbitales. Los primeros se determinaban por medio del m´etodo de Lambert, mientras que los segundos se obten´ıan aplicando el m´etodo de Lambert generalizado, esto es, mejorando la aproximaci´on obtenida por el m´etodo de Lambert. Esta extensi´on de arcos keplerianos a arcos orbitales se puede realizar tambi´en para los arcos cerrados de manera que, al igual que el m´etodo de Lambert, permite obtener los arcos keplerianos cerrados, el m´etodo de Lambert generalizado permite obtener los arcos orbitales cerrados. El concepto de arco orbital cerrado se hace especialmente ´util cuando se consideran ´orbitas cuasiestacionarias. En efecto, para el mantenimiento del orbitador en estas ´orbitas es preciso la realizaci´on de una maniobra cada periodo de tiempo T para lo cual es necesario que el punto donde se realice la maniobra sea el mismo que el punto inicial, esto es el v´ertice del bucle. Si pretendemos mantener una ´orbita cuasiestacionaria en un contexto orbital con un modelo con perturbaciones un arco kepleriano cerrado se convierte en abierto puesto que las perturbaciones separan el punto inicial del final. Para mantener la ´orbita es preciso calcular el arco orbital cerrado y no el arco kepleriano. En lo que sigue, analizaremos el proceso de mantenimiento de una ´orbita cuasiestacionaria en la Luna. Para ello consideraremos, al igual que en la parte final del cap´ıtulo anterior (apartado 3.2.5), un modelo de ´orbita lunar que tenga en cuenta todo el modelo gravitacional LP165P (165×165) as´ı, como la perturbaci´on producida por la Tierra, cuyo potencial viene dado por V⊕. Hemos establecido cuatro estrategias de mantenimiento de una ´orbita cuasiestacionaria basadas en la correcci´on sistem´atica de la ´orbita cada tiempo T, con un ∆vcalculado de cuatro formas distintas y hemos analizado la evoluci´on de la ´orbita calculando la separaci´on entre el punto inicial y el punto final al cabo de un gran n´umero de vueltas a la ´orbita. Con este an´alisis veremos cual de las cuatro estrategias es mejor y podremos medir la separaci´on final de la ´orbita respecto de su posici´on te´orica. Este estudio sustituye, en este caso, al estudio de la estabilidad de ´orbitas peri´odicas que no puede ser aplicado porque la periodicidad de ´estas ´orbitas se consigue por un proceso de correcci´on de la misma. 108 Arcos orbitales cerrados Las estrategias de correcci´on se pueden resumir en los siguientes puntos: a. La primera, que denotaremos por las siglas AKNA, consiste en calcular una ´unica vez el arco kepleriano cerrado (se usa el m´etodo de Lambert) y aplicar, en cada intervalo de tiempo T, un ∆v=XK 0−XK T, obtenido restando a la velocidad del punto final del arco kepleriano cerrado la del punto inicial. b. La segunda, que denotaremos por las siglas AONA, consiste en calcular una ´unica vez el arco orbital cerrado (se usa el m´etodo de Lambert generalizado) y aplicar, en cada intervalo de tiempo T, un ∆v=XO 0−XO T, obtenido restando a la velocidad del punto final del arco orbital cerrado la del punto inicial. c. La tercera, que denotaremos por las siglas AKA, consiste en calcular el arco kepleriano cerrado cada intervalo de tiempo Ty hacer la correcci´on de ∆v actualizando ´este cada vuelta. Esta estrategia est´a basada en que el punto final siempre est´a ligeramente desplazado respecto de la posici´on inicial por lo que arco kepleriano cerrado ser´a diferente al cambiar el v´ertice, por ello al final de cada vuelta se calcula el nuevo arco y por tanto la nueva correcci´on ∆v. d. La cuarta, que denotaremos por las siglas AOA, consiste en calcular el arco orbital cerrado cada intervalo de tiempo Ty hacer la correcci´on de ∆vactualizando ´este en cada vuelta. Al igual que en la estrategia AKA, ´esta se basa en recalcular el arco orbital debido al desplazamiento inevitable del v´ertice. Comenzaremos el estudio considerando un modelo gravitacional simplificado donde se considera ´unicamente el problema principal del sat´elite, esto es, el efecto del arm´onico zonal (J2). Con este modelo se ha propagado una ´orbita cuasiestacionaria con v´ertice en (r= 25rL, ψ = 35◦) y periodo T= 8 d durante 30 arcos orbitales o periodos (unos 240 d´ıas). En la gr´afica de la figura 4.18 se muestra la evoluci´on del error en la distancia entre el punto inicial y el punto final despu´es de cada vuelta (el eje Ox representa el n´umero de arcos orbitales propagados). En la parte izquierda la figura 4.18 se puede apreciar que para las estrategias AKNA y AKA, el error es de unos 35 m. en el primer arco, cantidad que va aumentando progresivamente llegando a errores de 3 km para la tercera vuelta del modelo AKNA y para la decimocuarta en el modelo AKA. La misma gr´afica muestra que la estrategia AONA presenta un crecimiento del error m´as moderado con un valor de 500 m, tras treinta vueltas, mientras que para la estrategia AOA no puede apreciarse 115 Entre los posibles trabajos futuros, derivados de la presente memoria, podemos mencionar los siguientes: Aplicaci´on de la extensi´on del m´etodo de correcci´on de ´orbitas a la b´usqueda de ´orbitas peri´odicas en el modelo orbital. En particular, podemos mencionar las ´orbitas sim´etricas alrededor de lunas de planetas, ´orbitas congeladas, etc. Extender el m´etodo de Lambert generalizado para m´ultiples periodos orbitales. Extender el an´alisis y la utilidad de las ´orbitas tipo Molniya, y otras ´orbitas que repitan la traza, para ´orbitas en torno a la Luna y a Marte. Adem´as se podr´a relacionar estas ´orbitas con el concepto de constelaciones regionales de sat´elites mencionado en la introduci´on. Bibliograf´ıa Abad, A. (2012). Astrodin´amica. Editor Bubok Publishing S.L. Espa˜na. Abad, A., Barrio, R., Blesa, F., and Rodriguez, M. (2012). Algorithm 924: Tides aTaylor Series Integrator for Differential EquationS.ACM Transactions on Mathematical Software (TOMS), 39(1):1–28. Abad, A., Barrio, R., and Dena, A. (2011). Computing periodic orbits with arbitrary precision. Physical Review E, E84(016701):1–6. Abad, A. and Lacruz, E. (2013). Computing derivatives of a gravity potential by using automatic differentiation. Celestial Mechanics and Dynamical Astronomy, 117(2):187–200. Anderson, P. and Macdonald, M. (2010). Extension of earth objects using low-thrust propulsion. 61th International Astronautical Congress, IAC. Czech Republic. Avanzini, G. (2008). A simple Lambert algorithm. Journal of Guidance, Control, and Dynamics, 31(6):1587–1594. Balmino, G., Barriot, J., Koop, R. Middel, B., Thong, N. C., and Vermeer, M. (1991). Simulation of gravity gradients: a comparison study. Bulletin g´eod´esique, 65(4):218–229. Battin, R. H. (1977). Lambert’s problem revisited.”. AIAA Journal, 15(5):707–713. Battin, R. H. (1999). Introduction to the Mathematics and Methods of Astrophysics. American Institute of Aeronautics and Astronautics, Inc. Reston, Virginia. Battin, R. H., Fill, T. J., and Shepperd, S. W. (1978). A new transformation invariant in the orbital boundary-value problem. Journal of Guidance, Control, and Dynamics, 1(1):50–55. Battin, R. H. and Vaughan, R. M. (1984). An elegant Lambert algorithm. Journal of Guidance, Control, and Dynamics, 7(6):662–670. 118 Bibliograf´ıa Burt, E. G. C. and Elliot, H. (1968). Space science and electrical propulsion. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 308(1493):217–241. Capderou, M. (2005). Satellites - Orbits and missions. SpringerVerlag, France. Casotto, S. and Fantino, E. (2007). Evaluation of methods for spherical harmonic synthesis of the gravitational potential and its gradients. Advances in Space Research, 40(1):69–75. Chandra, R., Dagum, L., Kohr, D., Maydan, D., J., M., and Menon, R. (2001). Parallel Programming in OpenMP. Morgan Kaufmann Publishers. Burlington. Chapman, B., Jost, G., and Van Der Pas, R. (2008). Using OpenMP: Portable Shared Memory Parallel Programming. Scientific and Engineering Computation Series. The MIT Press, Cambridge, Massachusetts. Chia-Chun, G. C. (2005). Applied Orbit Perturbation and Maintenance. The Aerospace Press, AIAA. California, USA. Cunningham, L. E. (1970). On the computation of the spherical harmonic terms needed during the numerical integration of the orbital motion of an artificial satellite. Celestial Mechanics, 2(2):207–216. Curtis, H. D. (2012). Orbital Mechanics for Engineering Students. ButterworthHeinemann. London. Danby, J. M. A. (1988). Fundamentals of Celestial Mechanics. Willmann-Bell, Virginia. Demmel, J. W. (1997). Applied Numerical Linear Algebra. Society for Industrial and Applied Mathematics, Philadelphia. Evans, B. G. (1999). Satellite Communication Systems. The Institution of Electrical Engineers, London. Flohrer, T., Choc, R., and Bastida, B. (2011). Classification of geosynchronous objects. GEN-DB-LOG-00086-OPS-GR, (14). Fukushima, T. (2012). Parallel computation of satellite orbit acceleration. Computers and Geosciences, 49:1–9. Gauss, C. F. (1809). Theoria motus corporum celestiaum in section-ibus conic solem ambientium. (English translation by C.H Davis (1857). Little, Brown and Ca. Boston. Bibliograf´ıa 119 Gooding, R. H. (1990). A procedure for the solution of Lambert’s orbital boundaryvalue problem. Celestial Mechanics and Dynamical Astronomy, 48(2):145–165. Griewank, A. and Walther, A. (2008). Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation. SIAM. Philadelphia USA. Gruber, T., Bode, A., Reigber, C., Schwintzer, R., and Balmino, G. (2000). GRIM5C1:Combination solution of the global gravity field to degree and order 120. Geophysical Research Letters, 27(24):4005–4008. Hairer, E., Nørsett, S. P., and Waner, G. (1993). Solving Ordinary Differential Equations. Non-stiff Problems. Springer Ser. Comput. Math. 8. Springer-Verlag. New York. Heiskanen, V. and Moritz, H. (1967). Theory of Satellite Geodesy. Blaisdell Publ Comp. Waltham, MA. USA. Holmes, S. A. and Featherstone, W. E. (2002). A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions, Journal of Geodesy, 76(5):279–299. Izzo, D. (2005). Lambert’s problem for exponential sinusoids. Journal of Guidance, Control, and Dynamics, 29(5):1242–1245. Jacchia, L. G. (1958). The earth’s gravitational potential as derived from satellites 1957 Beta one and 1958 Beta two. Smithsonian Astronomy Observatory, Special Report, 19(1). Kaula, W. H. (1966). Theory of Satellite Geodesy Applications of Satellites to Geodesy. Blaisdell Publishing Company, New York. King-Hele, D. G. (1958a). Analysis of the orbits of the russian satellites. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 253(1275):529–538. King-Hele, D. G. (1958b). The effect of the earth’s oblateness on the orbit of a near satellite. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 247(1248):49–72. Konopliv, A. S., Asmar, S. W., Carranza, E., Sjogren, W. L., and Yuan, D. N. (2001). Recent gravity models as a result of the Lunar Prospector mission. Icarus, 150(1):1–18. 120 Bibliograf´ıa Konopliv, A. S., Binder, A. B., Hood, L. L., Kucinskas, A. B., Sjogren, W. L., and Williams, J. G. (1998). Improved gravity field of the Moon from Lunar Prospector. Science, 281(5382):1476–1480. Konopliv, A. S. and Yoder, C. F. (1996). Venusian k2 tidal love number from Magellan and PVO tracking data. Geophysical Research Letters, 23(14):1857– 1860. Konopliv, C., Banerdt, W. B., and Sjogren, W. L. (1999). Venus gravity: 180th degree and order model. Icarus, 139(1):3–18. Kozai, Y. (1961). Tesseral harmonics of the potential of the earth as derived from satellite motions. Smithsonian Astronomy Observatory, Special Report, (72). Lacruz, E. and Abad, C. (2008). Automatizaci´on del proceso del c´alculo de las coordenadas geoc´entricas de sat´elites geoestacionarios a trav´es de observaciones astrom´etricas. TL1050 L3, B.I.A.C.I, Universidad de Los Andes, Venezuela. Lancaster, E. R. and Blanchard, R. C. (1969). A unified form of Lambert’s theorem. Technical Note NASA TMX-63355, Goddard Space Flight Center. Maryland. La´ınez, M. D. and Romay, M. M. (2009). Odandt process evolution based on interoperability between different navigation satellite system. Proceeding of the 22nd International Technical Meeting of The Sattellite Division of the Institute of Navigation, (ION GNSS), savannah, GA. La´ınez, M. D., Romay, M. M., and Alcantarilla, I. (2009). Approach for developing an advanced regional navigation satellite system. International Symposium on GPS/GNSS, Korea. Lemoine, F. G., Kenyon, S. C., Factor, J. K., Trimmer, R. G., Pavlis, N. K., Chinn, D. S., and Cox, C. M. (1998). The development of the joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) geopotential model EGM96. NASA Tech. Publi TP-1998-206861, NASA Goddard Space Flight Center, Washington, D. C. Lemoine, F. G., Smith, D. E., Kunz, L., Smith, R., and Pavlis, E. C. (1996). The development of the NASA GSFC and NIMA joint geopotential model.”. Gravity, Geoid and Marine Geodesy, 117:461–469. Lemoine, F. G., Smith, D. E., Rowlands, D. D., Zuber, M. T., Neumann, G. A., Chinn, D. S., and Pavlis, D. E. (2001). An improved solution of the gravity field of Mars (GMM-2B) from Mars Global Surveyor. Journal of Geophysical Research, 106(E10):23359–23376. Bibliograf´ıa 121 Lemoine, F. G., Smith, D. E., Zuber, M. T., Neumann, G. A., and Rowlands, D. D. (1997). A 70th degree lunar gravity model (GLGM-2) from Clementine and other tracking data. Journal of Geophysical Research, 102(E7):16339–16359. Lerch, F. J., Wagner, C. A., Smith, D. E., Sandson, M., Brownd, J. E., and Richardson, J. (1972). NASA Tech. Memo. X-65970. Lerch, F. J., Klosko, S. M., Laubscher, R. E., and Wagner, C. A. (1979). Gravity model improvement using GEOS 3 (GEM 9 and 10). Journal of Geophysical Research, 84(B8):3897–3916. Lerch, F. J., Nerem, R. S., Putney, B. H., Felsentreger, T. L., Sanchez, B. V., Marshall, J. A., Klosko, S., Patel, G., Williamson, R., Chinn, D., Chan, J., Rachlin, K., Chandler, N., Mccarthy, J. J., Luthcke, S., Pavlis, N. K., Pavlis, D. E., Robbins, J., and Kapoor, S. (1994). A geopotential model from satellite tracking, altimeter, and surface gravity data: GEM-T3.Journal of Geophysical Research, 99(B2):2815–2839. Lerch, F. J., Putney, B. H. Wagner, C. A., and Klosko, S. M. (1981). Goddard earth models for oceanographic applications (GEM 10B and IOC). Marine Goedesy, 5(2):145–187. Lorrell, J. (1969). Spherical harmonics applications to geodesy-some frequently used formulas. Jet Propulsion Laboratory, TM 311-112. Lundberg, J. B. and Schutzf, B. E. (1988). Recursion formulas of legendre functions for use with nonsingular geopotential models. Journal of Guidance, Control, and Dynamics, 11(1):31–38. Maral, G. and Bousquet, M. (2009). Satellite Communications Systems: Systems, Techniques and Technology. Jhon Wiley, and Sons, Ltd. UK. McInnes, C. R. (2011). Displaced non-keplerian orbits using impulsive thrust. Celestial Mechanics and Dynamical Astronomy, 110(3):199–215. McInnes, R. C. (1997). The existence and stability of families of displacement twobody orbits. Celestial Mechanics and Dynamical Astronomy, 67(2):167–180. McInnes, R. C. (1999). Solar Sailing, Technology, Dynamics and Mission Applications. Springer Praxis. Germany. McKay, R. J., Bosquillon de Frescheville, F., Vasile, M., McInnes, R. C., and Biggs, J. D. (2009). Non-keplerian orbits using low thrust, high ISP propulsion system. 60th International Astronautical Congress, IAC. Korea. 122 Bibliograf´ıa Metris, G., X, J., and I., W. (1998). Derivatives of the gravity potential with respect to rectangular coordinates. Celestial Mechanics and Dynamical Astronomy, (2):137–151. Miani, A. K. and Agrawal, V. (2011). Satellite Technology. Principles and Applications. Jhon Wiley and Sons, Ltd. Chichester. Milani, A., Nobili, A. M., and Farinella, P. (1987). Non-Gravitational Perturbations and Satellite Geodesy. Adam Hilger. England. Milani, A., Tommei, G., Farnocchia, D., Rossi, A., Schildknecht, T., and Jehn, R. (2011). Correlation and orbit determination of space objects based on sparse optical data. Monthly Notices of the Royal Astronomical Society, 417(3):2094– 2103. Montojo, F. J., L´opez Moratalla, T., and Abad, C. (2011). Astrometric position and orbit determination of geostationary satellites. Advanced in Space Research, 417(3):1043–1053. Montojo, F. J., L´opez Moratalla, T., Abad, C., and Mui˜nos, J. L. Astrometric reduction of geostationary satellites optical observations for orbit determination (PASAGE), journal=Revista Mexicana de Astronom´ıa y Astrof´ısica, year = 2008, volume = 34, pages=45–48. Multon, F. R. (1970). An Introduction to Celestial Mechanics. Dover Publications. New York. Neidinger, R. D. (1992). An efficient method for the numerical evaluation of partial derivatives of arbitrary order. ACM Transactions on Mathematical Software, 18(2):159–173. Neumann, G. A., Lemoine, F. G., Smith, D. E., and Zuber, M. T. (2003). The mars orbiter laser altimeter archive: Final precision experiment data record release and status of radiometry. In Lunar and Planetary Institute Science Conference Abstracts XXXIV, 34:1978. Nock, K. T. (1984). Rendezvous with Saturn’s rings. Anneux des Plan`etes Planetary Rings, IAU. Coloquium, (75). Olvers, F. W. and Smith, J. M. (1983). Associated Legendre functions on the cut. Journal of Computational Physics, (3):502–518. Bibliograf´ıa 123 Pavlis, N. K., Holmes, S. A., Kenyon, S. C., and Factor, J. K. (2012). The development and evaluation of the earth gravitational model 2008 (EGM2008). Journal of Geophysical Research, 117(B4):B04406. Pines, S. (1973). Uniform representation of the gravitational potential and its derivatives. AIAA Journal, (11):1508–1511. Prussing, J. E. (1979). Geometrical interpretation of the angles αand βin Lambert’s problem. Journal of Guidance, Control, and Dynamics, 2(3):442–443. Prussing, J. E. and Conway, B. A. (1993). Orbital Mechanics. Oxford University Press. New York. Racca, G. D. (2003). New challenges to trajectory design by the use of electric propulsion and other new means of wandering in the solar system. Celestial Mechanics and Dynamical Astronomy, 85(1):1–24. Rall, L. B. and Corliss, G. F. (1996). An introduction to automatic differentiation. In: C. Berz Bischof, Griewank (eds). Computational Differentiation: Techniques, Applications, and Tools, pages 1–17. Reigber, C., Balmino, G., Muller, H., Bosch, W., and Moynot, B. (1985). Grim gravity model improvement using LAGEOS (GRIM3-L1). Journal of Geophysical Research, 90(B11):9285–9299. Reigber, C., Schwintzer, P., Barth, W., Massmann, F. H., Riamondo, J. C., and Bode, A. (1993). GRIM4-C1, -C2p: combination solutions of the global earth gravity field. Surveys in Geophysics, 14(4-5):381–393. Smith, D. E., Lerch, F. J., Nerem, R. S., Zuber, M. T., Patel, G. B., Fricke, S. K., and Lemoine, F. G. (1993). An improved gravity model for Mars: Goddard Mars model 1. Journal of Geophysical Research, 98(E11):20871–20889. Soop, E. M. (1994). Handbook of Geostationary Orbits. Kluwer Academic, London and Microcosm, Inc, California. Stoer, J. and Bulirsch, R. (1993). Introduction to Numerical Analysis. SpringerVerlag. Berlin. Tapley, B. D., Watkins, M. M., Ries, J. C., Davis, G. W., Eanes, R. J., Poole, S. R., Rim, H. J., Schutz, B. E., Shum, C. K., Nerem, R. S., Lerch, F. J., Marshall, J. A., Klosko, S. M., Pavlis, N. K., and Williamson, R. G. (1996). The joint gravity model 3, JGM-3.Journal of Geophysical Research, 101(B12):28029–28049. 124 Bibliograf´ıa Torge, W. (2001). Geodesy. De Gruyter. Hannover. Tscherning, C. C. (1976). Computation of the second-order derivatives of the normal potential based on the representation by a legendre series. Manuscr Geod, 1:71–92. Tsukanov, I. and Hall, M. (2003). Data structure and algorithms for fast automatic differentiation. International Journal for Numerical Methods in Engineering, 56(13):1949–1972. Vallado, D. A. (2001). Fundamentals of Astrodynamics and Applications. Microcosm Press. California and Kluwer Academic Publishers. The Netherlands. Wertz, J. R. and Larson, W. J. (2010). Space Mission Analysis and Design. Microcosm Press, Hawthorne, CA and Springer, New York. Wiggins, R. A. and Saito, M. (1971). Evaluation of computational algorithms for the associated legendre polynomials by interval analysis. Bulletin of the Seismological Society of America, (2):375–381.