Trabajo de Grado Final

FACULTAD DE INGENIERÍAS ´ EN ASTROINGENIERÍA ALFA ORION ´ GRUPO DE INVESTIGACION

TRABAJO DE GRADO:

´ DE LOS PARAMETROS ´ DETERMINACION ORBITALES DE ASTEROIDES Y COMETAS USANDO LOS INVARIANTES DEL MOVIMIENTO

Un trabajo de grado presentado como requisito parcial para optar por el t´ıtulo de Ingeniero F´ısico de la Universidad Tecnológica de Pereira Presentado por: Santiago Jim´ enez Villarraga

Asesorado por: Msc. Edwin Andrés Quintero

“When I consider the short duration of my life, swallowed up in the eternity before and after, the little space which I fill, and even can see, engulfed in the infinite immensity of spaces which I am ignorant, and which know me not , I am frightened and am astonished at being here rather than there; for there is no reason why here rather than there, why now rather than then... The eternal silence of these infinite spaces frightens me.” Blas Pascal

I

II

Agradecimientos Quiero, en primera instancia, agradecer a m´ıs padres que siempre me brindaron su apoyo durante m´ıs estudios. También deseo agradecer al Ingeniero Edwin Quintero, por su valiosa asesor´ıa y paciencia en la realización de este trabajo. Finalmente quiero agradecer a m´ıs compa˜ neros del grupo de investigación Alfa Orión por su ayuda y acompa˜ namiento durante todo este tiempo.

III

IV

AGRADECIMIENTOS

Resumen Determinar los parámetros orbitales de asteroides y cometas es importante para establecer probabilidades de impactos futuros con la Tierra y para determinar poblaciones de estos cuerpos y sus distribuciones, de tal forma que se puedan validar o generar teor´ıas sobre la formación y evolución inicial del Sistema Solar. En este trabajo se estudia el problema consistente en determinar un conjunto preliminar de parámetros orbitales para cuerpos menores del Sistema Solar a partir de observaciones de sus posiciones relativas en la esfera celeste. Se describe y deducen los elementos principales de un método basado en los invariantes del movimiento para el problema de los dos cuerpos en una variante geocéntrica y en una topocéntrica. También se muestra detalladamente la forma en que se implementó computacionalmente. El algoritmo se prueba inicialmente con cuerpos pertenecientes a las familias de asteroides MBA (main-belt asteroids) y NEO (near earth objects). Se analizaron los resultados y errores obtenidos para ambas versiones del algoritmo. Para tal efecto, se usó dos parámetros que dan estimaciones del error en la forma y en la orientación de la o´rbita. Los resultados iniciales mostraron que el método implementado es adecuado para obtener o´rbitas preliminares. También se obtuvo que, en general, la versión más precisa es la topocéntrica, por lo cual se optó por usar esta versión para el resto del trabajo. Finalmente, el método se aplica a observaciones realizadas desde el Observatorio Astronómico de la Universidad Tecnológica de Pereira (OAUTP). Se describen los procedimientos experimentales que fueron necesarios para obtener estas mediciones y se muestran los resultados de la aplicación del método con estos datos. Los resultados obtenidos muestran que el algoritmo implementado es eficiente y práctico para determinar órbitas iniciales con datos tomados desde el OAUTP, puesto que los errores obtenidos están dentro del margen de errores cuando se usa el problema de los dos cuerpos como modelo dinámico.

V

VI

RESUMEN

Contenido Resumen

V

Índice de figuras

XI

Índice de tablas 1 Introducci´ on 1.1 Introducción . . . 1.2 Antecedentes . . 1.3 Objetivos . . . . 1.3.1 General . 1.3.2 Espec´ıficos

XIII

. . . . .

1 1 3 7 7 7

. . . . . .

9 9 9 11 12 15 16

. . . . .

19 19 20 21 25 26

4 Teor´ıa de eliminaci´ on algebraica 4.1 Introducción a la eliminación algebraica . . . . . . . . . . . . . . . 4.2 Los Teoremas de eliminación y extensión . . . . . . . . . . . . . . 4.3 Complejidad computacional y el Resultante . . . . . . . . . . . .

31 31 36 38

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

2 Teor´ıa b´ asica de ´ orbitas 2.1 El problema de los dos cuerpos . . 2.1.1 La ecuación de movimiento 2.1.2 Solución geométrica . . . . . 2.1.3 La o´rbita en el tiempo . . . 2.2 Sistema de coordenadas . . . . . . 2.3 Parámetros orbitales . . . . . . . . 3 M´ etodos de inversi´ on de ´ orbitas 3.1 El problema de inversión de órbitas 3.2 Método de Laplace . . . . . . . . . 3.2.1 Versión geocéntrica . . . . . 3.2.2 Versión topocéntrica . . . . 3.3 El método de Gauss . . . . . . . .

VII

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

. . . . . .

. . . . .

VIII 5 Fundamento te´ orico del m´ etodo 5.1 Vectores atribu´ıbles . . . . . . . . 5.2 Las integrales del movimiento . . 5.3 Las ecuaciones del método . . . . 5.4 Solución a las ecuaciones . . . . . 5.4.1 Selección de las soluciones 5.5 Implementación numérica . . . . 5.6 Resultados iniciales . . . . . . . . 5.7 Análisis de errores . . . . . . . . 6 Resultados experimentales 6.1 Planeamiento de las observaciones 6.2 Observaciones . . . . . . . . . . . 6.3 Reducción de datos . . . . . . . . 6.4 Resultados observacionales . . . . 6.5 Resultados del algoritmo . . . . .

CONTENIDO

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

41 41 42 43 45 46 47 49 51

. . . . .

55 55 57 62 66 67

7 Conclusiones y trabajo futuro

73

A C´ odigo Fuente de funciones b´ asicas A.1 Convertir UT a JD . . . . . . . . . . . . . . . . . A.2 Convertir JD a UT . . . . . . . . . . . . . . . . . A.3 Transforma la ascensión recta a radianes . . . . . A.4 Transformar la declinación a radianes . . . . . . A.5 Resuelve la ecuación de Kepler . . . . . . . . . . . A.6 Calcula el vector de estado . . . . . . . . . . . . . A.7 Calcula el vector de estado de la Tierra . . . . . . A.8 Parámetros orbitales a partir del vector de estado A.9 Obtiene el vector ecl´ıptico ecuatorial de la Tierra. A.10 Obtiene el vector atribu´ıble . . . . . . . . . . . . A.11 Determina errores orbitales . . . . . . . . . . . . . A.12 Grafica la trayectoria . . . . . . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

77 77 77 78 78 79 81 82 83 85 86 87 87

´ B C´ odigo Fuente para Algebra de polinomios B.1 Suma de polinomios . . . . . . . . . . . . . . B.2 Producto de polinomios . . . . . . . . . . . B.3 Simplificación de polinomios . . . . . . . . . B.4 Evaluación de polinomios en una variable . . B.5 Evaluación de polinomios de dos variables . B.6 Algoritmo FFT . . . . . . . . . . . . . . . . B.7 Algoritmo IDFT . . . . . . . . . . . . . . . B.8 Cálculo del resultante . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

91 91 91 92 92 92 93 93 94

C C´ odigo Fuente m´ etodo principal

. . . . . . . .

. . . . . . . .

. . . . . . . .

97

CONTENIDO

IX

Bibliograf´ıa

101

X

CONTENIDO

Índice de figuras Figura 2.1 Parámetros orbitales de una o´rbita en el problema de los dos cuerpos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 Figura 2.2 Relación entre los parámetros orbitales y el vector de estado. 17 Figura 3.1 Movimiento de un asteroide alrededor del Sol y observado desde la Tierra. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figura 5.1 Diagrama de flujo del algoritmo principal . . . . . . . . . . Figura 5.2 Trayectorias heliocéntricas para los cuerpos de estudio de la familia MBA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figura 5.3 Trayectorias heliocéntricas obtenidas para los cuerpos NEO. Figura 5.4 Errores obtenidos para cada cuerpo. En la parte izquierda se muestra el error en la forma de la órbita, y en la derecha, el error en la orientación de la o´rbita. . . . . . . . . . . . . . . . . . . . . Figura 5.5 Combinación de ambos errores en el plano complejo. En azul se muestra los errores obtenidos con la versión geocéntrica y en rojo los de la versión topocéntrica. El error es menor entre más cerca se esté del origen. . . . . . . . . . . . . . . . . . . . . . . . . Figura 6.1 Ventana de b´ usqueda del programa TOOLKIT FOR CCD ASTROMETRY. . . . . . . . . . . . . . . . . . . . . . . . . . . . Figura 6.2 Lista de asteroides que satisfacen las condiciones de la ventana de b´ usqueda. . . . . . . . . . . . . . . . . . . . . . . . . . . . Figura 6.3 Camino en el programa para encontrar campos de referencia. Figura 6.4 Campo de referencia obtenido del survey ESO-DSS I/II. . Figura 6.5 Imágen obtenida del asteroide 675. . . . . . . . . . . . . . Figura 6.6 Imágenes tomadas del asteroide 675 . . . . . . . . . . . . . Figura 6.7 Imágenes tomadas del asteroide 654 . . . . . . . . . . . . . Figura 6.8 Configuración de la sección CCD. . . . . . . . . . . . . . . Figura 6.9 Configuración de la sección Program. . . . . . . . . . . . . Figura 6.10 Ajuste manual de los parámetros instrumentales. . . . . . Figura 6.11 Trayectoria calculada para el asteroide 675 . . . . . . . . . Figura 6.12 Trayectoria calculada para el asteroide 654 . . . . . . . . . XI

19 48 51 52

53

54 56 57 59 59 60 61 61 63 64 66 69 70

XII

ÍNDICE DE FIGURAS Figura A.1 Diagrama de flujo de la solución a la ecuación de Kepler . Figura A.2 Diagrama de flujo del algoritmo que obtiene el vector de estado a partir de los parámetros orbitales. . . . . . . . . . . . . . Figura A.3 Diagrama de flujo del algoritmo que obtiene los parámetros orbitales a partir del vector de estado. . . . . . . . . . . . . . . .

79 81 84

Índice de tablas Tabla 5.1 Tiempos y lugares de observación de los asteroides utilizados como objeto de estudio. . . . . . . . . . . . . . . . . . . . . . . . Tabla 5.2 Comparación entre los parámetros orbitales. . . . . . . . . .

49 50

Tabla Tabla Tabla Tabla Tabla Tabla Tabla

67 67 68 68 69 70 71

6.1 6.2 6.3 6.4 6.5 6.6 6.7

Datos astrométricos obtenidos para 675. . . . . Datos astrométricos obtenidos para 654. . . . . Datos astrométricos adicionales para 675. . . . Parámetros orbitales de 675. . . . . . . . . . . Datos astrométricos adicionales para 654. . . . Parámetros orbitales de 654. . . . . . . . . . . Errores para los parámetros orbitales obtenidos

XIII

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

XIV

ÍNDICE DE TABLAS

Cap´ıtulo 1 Definici´ on del Problema y objetivos del trabajo 1.1

Introducci´ on

Conocer de forma precisa la o´rbita de asteroides y cometas es una tarea fundamental para el ser humano, en tanto que es muy importante la anticipación de colisiones futuras con el planeta. Existen evidencias que sugieren que colisiones en el pasado han provocado incluso, la extinción de especies. El evento registrado más reciente de una colisión fue en el 2013, cuando el Meteorito Chelyabinsk impactó sobre Rusia, [1]. El proyecto WISE (Wide-field infrared Survey Explorer) de la NASA estima que existen aproximadamente 19500 cuerpos cercanos a la Tierra NEO (Near-Earth objects)1 , de los cuales la mitad a´ un no han sido descubiertos por lo comunidad cient´ıfica. Piazi en 1801 descubrió Ceres [2], aunque este solo fue capaz de seguir el cuerpo durante pocas noches, puesto que no exist´ıa forma de predecir su movimiento. El problema de inversión de o´rbitas consiste en obtener los parámetros orbitales que describen el movimiento de un cuerpo a partir de información astrométrica. A los resultados obtenidos se le asocia una incertidumbre y una confiabilidad. Los primeros en proponer una solución efectiva a este problema fueron Gauss [3], [4] y Laplace [5, 6], quienes desarrollaron métodos que permiten conocer de forma aproximada la órbita del cuerpo usando observaciones astrométricas distintas de este. Su solución, si bien aproximada, fue el primer paso para conocer de forma precisa el movimiento de Ceres cuando más observaciones estuvieran disponibles. En general, el proceso de obtener una órbita confiable se divide en tres etapas distintas [7, 8]: la o´rbita inicial (o preliminar), órbita de m´ınimos cuadrados y 1

Entre estos est´ an los PHA, asteroides potencialmente peligrosos, los cuales tienen órbitas que se acercan bastante a la Tierra. Hoy en d´ıa se han detectado 1588, pero se estima que a´ un faltan muchos PHA’s por descubrir. Fuente: http:\\neo.jpl.nasa.gov/faq/#neo .

1

´ CAPÍTULO 1. INTRODUCCION

2

control de calidad. Desde un punto de vista matemático, lo descrito antes es un problema de m´ınimos cuadrados no lineal, en el cual la o´rbita preliminar es una primera estimación. Por la no-linealidad del problema se hace necesario recurrir a alg´ un método numérico para resolver el problema de m´ınimos cuadrados. El más com´ un de estos métodos en la literatura [6, 5] es una variación del método de Newton-Raphson conocido como método de corrección diferencial. Como es bien sabido, para que estos métodos converjan la ra´ız inicial 2 debe estar dentro del intervalo de conergencia del método, razón por la cual no deben ahorrarse esfuerzos para que la órbita preliminar sea tan buena como sea posible. Los métodos clásicos de Gauss y Laplace son considerados aceptables en el pro´ posito de encontrar una o´rbita preliminar del cuerpo estudiado [9], sin embargo, estos procedimientos fueron concebidos hace casi dos siglos ya, en tiempos en los que adquirir información astrométrica era bastante complicado y en la que poco o nada se conoc´ıa sobre la estructura del Sistema Solar [7]. El avance de la tecnolog´ıa ha provocado cambios fundamentales en el problema de inversión de o´rbitas. La mayor capacidad de los telescopios y la aparición de las camáras CCD [10] que permiten tomar una cantidad considerable de imágenes del cielo en cuestión de minutos ha generado una revolución en esta área en las u ´ltimas decadas [11]. Se pueden mencionar dos razones fundamentales por las que se han generado cambios y creado problemas modernos en esta a´rea: Los métodos clásicos presentan una convergencia muy débil [2, 8] si la información utilizada representa un arco con una curvatura muy peque˜ na de la o´rbita. Y, en general, esto es lo que ocurre normalmente: Se identifica que el cuerpo se está moviendo en la esfera celeste mediante la toma de una serie de imágenes (entre 3 y 6) en el rango de una hora de observación y se espera utilizar todas estas imágenes con el objetivo de obtener la órbita. Este tren de observaciones cercanas se conoce en la literatura como un VSA (a very short arc, [12]) si de este tren se puede conseguir una o´rbita preliminar y se conoce como un TSA (A too short arc) si no es posible obtenerla. Luego ha surgido la cuestión de cómo mejorar los métodos clásicos para aumentar su intervalo de convergencia o qué métodos nuevos pueden ser utlizados en lugar de los clásicos y que permitan utilizar toda la información disponible. Trabajos en cada forma de enfrentar el problema han aparecido en los u ´ltimos a˜ nos, [8, 13]. La gran cantidad de informaci´ on en intervalos tan cortos de tiempo ha hecho que la confiabilidad sea un reto. Esto se da por dos razones: en la época de Gauss y de Laplace no se conoc´ıa mucho sobre los distintos tipos de cuerpos menores del sistema solar, y una suposición fundamental para que 2

En este caso la ra´ız inical es la ´ orbita preliminar.

1.2. ANTECEDENTES

3

el método de Laplace sea aceptable es considerar que las observaciones son geocéntricas y que esto no afecta demasiado el resultado. No obstante, esto no es cierto para todas las familias de asteroides que se han descubierto en la actualidad. En general, para los asteroides del cinturón de asteroides (MBA, por Main Belt asteroids) s´ı se puede obtener resultados aceptables con esta suposición. Sin embargo, para objetos cercanos a la Tierra (NEO’s) y objetos tras-neptunianos (TNO’s) estas consideraciones fallan. En definitiva puede decirse que los métodos clásicos presentan problemas en la actualidad por razones que estaban fuera del alcance de los autores clásicos, puesto que en aquella época ni se conoc´ıa las caracter´ısticas de los cuerpos observables ni se contaba con la poderosa herramienta de los computadores. Ahora mismo se cuenta con la capacidad computacional suficiente para que el enfoque no sea la simplicidad sino la exactitud y además se conocen las limitaciones de los métodos clásicos, lo cual permite que se intenten resolver estas limitaciones ya sea haciendo mejores de estos o proponiendo métodos completamente alternativos.

1.2

Antecedentes

En la antig¨ uedad, obtener datos astrométricos de cuerpos celestes era una tarea dif´ıcil e incluso tediosa, debido a que requer´ıa seguir el cuerpo durante varias noches diferentes y las medidas se hac´ıan de forma poco automatizada, comparando la posición del cuerpo con estrellas vecinas. Hoy en d´ıa, por el contrario, la adquisición de datos astrométricos se ha convertido en un procedimiento relativamente sencillo, en el cuál no se requiere de equipos muy sofisticados para efectuarse. La aparición de las cámaras CCD y la mejora en las caracter´ısticas ópticas de los telescopios ha hecho posible obtener fotos de un mismo cuerpo varias veces en una misma noche, conllevando esto a un incremento exponencial de la cantidad de datos disponibles. Un cuerpo en movimiento, como un asteroide o un cometa, se identifica al comparar fotograf´ıas del lugar de la esfera celeste que se está estudiando y notar cambios de este con respecto a las estrellas fijas. Esto se puede realizar varias veces en una misma noche, y las observaciones entregan lo que se conoce como un arco demasiado corto (A too short arc) [14]. Esto, sin embargo, ha creado problemas adicionales en el proceso de determinación de o´rbitas, uno de los cuales se conoce como el problema de linkage ([8], [2], [15]). Este problema consiste en determinar con exactitud si un conjunto de datos, usualmente conocidos como atribuibles, pertenecen al mismo cuerpo. Este problema aparece porque, como se dijo antes, la información obtenida por mediciones es bastante elevada, y es frecuente que en un mismo grupo de datos se encuentre información sobre más de un cuerpo. Otro problema propio del alto n´ umero de datos es la obtención de métodos más efectivos, rápidos y precisos para obtener de manera preliminar

4


la órbita del cuerpo en estudio que permita el uso de más información para su consecución, puesto que los métodos clásicos no son efectivos, generalmente, con la utilización de este tipo de datos. Para la solución de estos problemas se han propuesto varios métodos, que en adelante se describirán. En [16] se estudia el problema de la identificación de asteroides, es decir, usar los datos insuficientes obtenidos del cuerpo para calcular una región de futura aparición. Para ello se establecen algor´ıtmos semi-lineales para calcular la región de confianza que de mejores predicciones que los algoritmos clásicos. El método se probó para varios cuerpos como los asteroides 4161 PL, 92B4 y 9076P2 para los cuales se encontró que las observaciones de recuperación se encuentran dentro de la región admisible calculada con σ = 3 (desviación estándar) en la mayor´ıa de los casos, y en algunos hasta σ = 6. El método es todav´ıa muy teórico y necesita ser probado para casos más generales. En [17] y [12] se contin´ uan estudiando los resultados de [16] y además se fijan conceptos nuevos e importantes que surgen en el problema moderno de la determinación de órbitas. El problema analizado en [17] es el de atribución, es decir, determinar numéricamente si un arco detectado corresponde a un objeto previamente descubierto. El algoritmo consiste en realizar 3 filtros distintos a los datos en los cuales se tienen en cuenta las incertidumbres en los datos calculados. Usando este algoritmo y en un a˜ no de prueba se consiguieron realizar 2577 atribuciones de las cuales 1502 se realizaron utilizando solo el método propuesto en el art´ıculo. Son factibles de mejoras tanto el procedimiento para seleccionar posibles atribuciones como los métodos para resolverlas numéricamente. La importancia de estos trabajos sobre el que aqu´ı se propone es que en estos se introducen nuevas definiciones teóricas, procedimientos numéricos y técnicas experimentales que se deben tener en cuenta para la determinación de o´rbitas de asteroides y cometas. Otro problema que ha sido objeto de alto interés es el de obtener mejores métodos para calcular la órbita preliminar directamente o en su defecto, realizar una extensión o mejora a los métodos clásicos existentes, de tal forma que se puedan aplicar sobre los conjuntos de datos modernos. En el problema de determinar los parámetros orbitales de un cuerpo, la o´rbita preliminar es el primer paso, y su importancia radica en la posibilidad de efectuar un seguimiento al cuerpo estudiado, lo cual permite obtener más datos del mismo y de esta forma, mejorar la precisión con la que se conocen los parámetros orbitales. Las aproximaciones modernas a este tema usan el problema de los dos cuerpos en su gran mayor´ıa, pero a la vez hacen uso de los invariantes de movimiento, la energ´ıa del sistema y el momento angular.

1.2. ANTECEDENTES

5

Mirtorabi, [13], describe un método para determinación de órbitas que ofrece más exactitud que los clásicos. En este art´ıculo es extiende el método de Gauss de 3 a N observaciones y se usan las invariantes del movimiento: la energ´ıa total y el momentum angular para mejorar en un proceso iterativo los datos calculados. Se prueba el método con datos virtuales obteniéndose resultados con errores inferiores a 1%. El procedimiento se debe probar sobre datos reales y se deben estudiar casos con errores que siguen una distribución distinta a la gaussiana, además se debe analizar la convergencia del método. Una forma alternativa de proceder, es la que se usa en [18]. Aqu´ı se enfrenta el problema de determinar una primera o´rbita para conjuntos de datos peque˜ nos, con una longitud de arco corta, con los cuales muchos de los métodos clásicos fallar´ıan. En estos trabajos se parte también de la conservación del momentum angular y la energ´ıa del sistema, se propone escribir sistemas de ecuaciones que se resuelven utilizando algoritmos numéricos tradicionales, como el método multidimensional de Newton-Raphson y modificaciones al método de Lagrange y de Laplace. Estas técnicas se aplican en estos trabajos en la descripción del movimiento de satélites artificiales en su órbita alrededor de la Tierra, con una o varias observaciones, en la que tienen que los datos calculados poseen errores en la mayor´ıa de los casos, inferiores a 1%. Una conclusión importante también sobre este procedimiento, es la rapidez con la que se efect´ ua, permitiendo una excelente estimación de datos con muy pocas observaciones, lo que reduce, en principio, problemas de pérdidas de los cuerpos estudiados. Vale la pena aclarar que en la forma en que se encuentra construido este procedimiento, los resultados siempre serán locales y por lo tanto pueden encontrarse ocasiones en que soluciones admisibles sean descartadas, lo cual podr´ıa afectar la calidad de los datos obtenidos. En [19] y [20] también se utiliza la conservación de estas magnitudes f´ısicas para escribir un sistema de ecuaciones para la distancia topocéntrica y la velocidad radial del cuerpo para dos épocas distintas. Sin embargo, en lugar de resolver estas por métodos numéricos como los anteriores, se utiliza a´lgebra computacional para reducir el sistema de ecuaciones a una sola ecuación que se resuelve después utilizando el algoritmo DFT o un procedimiento de eliminación de variables. Al ejecutar esto se utiliza también el hecho de que se tienen más datos de los que son estrictamente necesarios, lo cual permite dise˜ nar un algoritmo para resolver el problema de linkage. Además también se puede determinar el error del procedimiento mediante la construcción de una matriz de covarianzas. Ambos métodos de solución se probaron con datos del asteroide 1999 NR23 para las épocas t1=5400 y t2=54109 (tiempos en MJD). Se encontró que en ambos casos, el error entre los datos calculados y los reales fue del orden de 3 exp(-5) AU

6


para las distancias topocéntricas y que el error en los parámetros orbitales es del mismo orden que el esperado para el modelo de los dos cuerpos, que se usa como modelo de fuerza en estos problemas, es decir, el método de solución propuesto y su implementación arrojan unas excelentes primeras aproximaciones. Si bien en estos trabajos se usó el algoritmo DFT y una eliminación de variables para reducir el sistema de ecuaciones, también pueden buscarse métodos alternativos que reduzcan a´ un más la complejidad computacional de estos algoritmos. En [21] se aplican estos métodos teóricos a la basura espacial. El objetivo principal de este art´ıculo es investigar la posibilidad de estructurar una red global en la cual mediante el uso de sensores y observaciones se pueda formar un catálogo completo de la basura espacial, de tal forma que se puedan prevenir colisiones y fragmentaciones. En este documento se utiliza el método de las integrales de Kepler para obtener de forma inicial, los parámetros o´rbitales que describen las órbitas de LEO’s (lowEarth objects) principalmente. Los datos usados son simulaciones de estos cuerpos, y el algoritmo se prueba para datos distribuidos durante una semana y para datos distribuidos durante tres semanas. Se encuentra que el procedimiento para una semana no es concluyente, puesto que los cuerpos dentro de una envoltura de exactitud menor al 1 % no superan el 80 %, en el mejor de los casos. Sin embargo, se encontró que para los datos de tres semanas los resultados si pueden ser concluyentes puesto que en este caso la proporción de cuerpos dentro del 1% de exactitud siempre es mayor al 94.3 %. En este punto se evidencia la importancia que reviste la realización del proyecto que aqu´ı se plantea, puesto que la mayor´ıa de las aplicaciones que se encuentran hasta hace muy poco tiempo están en una fase de desarrollo y simulación, haciendo importante la prueba experimental de los mismos desde un observatorio astronómico más peque˜ no. A nivel nacional se han realizado estudios en el a´rea de astrometr´ıa de asteroides y cometas. Estos estudios se han realizado principalmente desde el observatorio astronómico de la Universidad de Nari˜ no. Por ejemplo, en [22], se realiza astrometr´ıa al asteroide 2003Q0104. Se realizan 72 observaciones distintas del cuerpo y mediante el método de Gauss se calculan los parámetros órbitales del objeto. Sin embargo, los métodos utilizados en este documento son los tradicionales, por lo que se concluye nivel nacional no se conocen estudios realizados en la dirección que aqu´ı se propone, lo cual hace que cobre particular relevancia el inicio de cualquier investigación en esta a´rea.

1.3. OBJETIVOS

1.3 1.3.1

7

Objetivos General

Implementar un algoritmo computacional basado en el método de las Integrales de Kepler para el cálculo de parámetros órbitales de asteroides y cometas y verificarlo experimentalmente con datos tomados desde el Observatorio Astronómico de la Universidad Tecnológica de Pereira.

1.3.2

Espec´ıficos

1. Desarrollar un algoritmo numérico para la determinación de parámetros orbitales usando los métodos teóricos estudiados en el proyecto. 2. Implementar el algoritmo numérico en una herramienta computacional y determinar los errores que aparecen en el mismo. 3. Aplicar técnicas instrumentales para la toma de las fotograf´ıas del cuerpo de estudio y para la extracción de los datos necesarios para el algoritmo implementado. 4. Evaluar los resultados de aplicar el algoritmo implementado a los datos obtenidos desde el observatorio astronómico mediante la comparación de los mismos con respecto a los datos entregados por el MPC y la JPL3 .

3

La JPL (Jet Propulsion Laboratory) a través de su programa JPL Small-Body Database Browser y la IAU (International Astronomical Union), a través del MPC (Minor Planet Center) proveen bases de datos de todos los asteroides y cometas conocidos, que al ser construida por datos tomados por muchos observatorios alrededor del mundo, se trata de valores muy exactos.

8


Cap´ıtulo 2 Teor´ıa b´ asica de ´ orbitas Se hará una breve descripción de la formulación técnica del problema que se va a estudiar. La información aqu´ı plasmada proviene de libros como [6, 5, 23]. Lo primero que se hará es discutir el modelo dinámico que se empleará en las secciones subsiguientes. Como el propósito del trabajo es determinar o´rbitas iniciales, un modelo apropiado por simplicidad y eficiencia es el de los dos cuerpos [7, 13] o problema de fuerzas centrales 1 .

2.1 2.1.1

El problema de los dos cuerpos La ecuaci´ on de movimiento

El problema que se quiere resolver es el movimiento de una masa puntual en el espacio2 , cuando está sometida a una fuerza que depende solo de la distancia (en este caso inversamente proporcional al cuadrado de la distancia) y cuya l´ınea de acción pasa en todo instante por un punto fijo, que se considera el origen del sistema de coordenadas. La simetr´ıa existente en la definición del problema hace suponer que un buen sistema de coordenadas es el esférico. As´ı, si la part´ıcula estudiada tiene masa m, las energ´ıas cinética y potencial están dadas por: Z 1 2 2 2 2 2 ˙2 V (r) = − F (r) dr (2.1) T = m r˙ + r sin θϕ˙ + r θ , 2 Con: F (r) = − 1

µ r2

(2.2)

En este caso se puede entender como un problema de fuerza central (como lo entiende [24]) debido a que se considera el movimiento de cuerpos menores (con masas muy peque˜ nas) alrededor del Sol, que para todos los efectos está fijo en el origen de un sistema de coordenadas adecuado [23, 25] 2 Esta es otra aproximaci´ on inicial. Se considera que los cuerpos son esféricos y de densidad uniforme, por lo cual se pueden aproximar como masas puntuales.

9

´ ´ CAPÍTULO 2. TEORÍA BASICA DE ORBITAS

10

Donde µ es una constante gravitacional 3 . Con estas expresiones se escribe el lagrangiano del sistema y utilizando las ecuaciones de Euler-Lagrange [24, 26] se obtienen las ecuaciones de movimiento: d (mr) ˙ − mr sin2 θϕ˙ 2 − mrθ˙2 = F (r) (2.3) dt d 2 ˙ (2.4) mr θ − mr2 ϕ˙ 2 sin θ cos θ = 0 dt d mr2 sin2 θϕ˙ = 0 (2.5) dt Como se nota en (2.1), la coordenada ϕ es una coordenada ignorable [24], lo cual implica que: ∂T = mr2 sin2 θϕ˙ = K = cons pϕ = ∂ ϕ˙ Esta ecuación podr´ıa utilizarse para eliminar ϕ˙ del resto de ecuaciones. No obstante, si para describir el movimiento solo son fundamentales r y θ esto quiere decir que el movimiento es bidimensional, el cual entonces se presenta en ϕ = const, lo que implica que: ϕ˙ = 0 ⇒ K = 0 En definitiva, la energ´ıa cinética se puede escribir como: 1 T = m r˙ 2 + r2 θ˙2 2

(2.6)

Y las ecuaciones de movimiento son: d (mr) ˙ − mrθ˙2 = F (r) dt d 2 ˙ mr θ = 0 dt Esta u ´ltima genera otra conservación de momento, puesto que: pθ =

(2.7) (2.8)

∂T = mr2 θ˙ = k = cons ∂ θ˙

De la que se sigue que: dθ h = 2 (2.9) dt r Donde h es el momento angular por unidad de masa. Además del momento, como la fuerza es potencial, la energ´ıa también se conserva: µ 1 2 m r˙ + r2 θ˙2 − = E 0 2 r 3

M´ as adelante se definir´ a y dar´ a el valor utilizado en los cálculos numéricos.

2.1. EL PROBLEMA DE LOS DOS CUERPOS

11

Donde E 0 es la energ´ıa del movimiento. Esta expresión se puede reescribir, utilizando (2.9) y simplificando, como: h2 µ r˙ + 2 − 2 = 2E r r 2

(2.10)

Donde E es la energ´ıa total por unidad de masa. Esta es la ecuación cuya solución permite conocer la posición del cuerpo con respecto al origen, en cualquier instante dado de tiempo.

2.1.2

Soluci´ on geom´ etrica

Resolver (2.10) no es tan sencillo puesto que es una ecuación no lineal. Por tal motivo, primero se resuelve para r en función de θ, i.e., r (θ). Para ello, notése que: dr dr dθ = dt dθ dt Utilizando (2.9): h dr d h dr = 2 =− (2.11) dt r dθ dθ r Reemplazando en (2.10) y ordenando un poco, se obtiene: s r 2 2 2µ h µ2 h 2µ µ2 d h = 2E + − 2 = − − + 2 − dθ r r r h2 r2 r h Haciendo: α2 = 2E +

µ2 h2

(2.12)

La anterior ecuación queda como: d − dθ

s 2 h h µ = α2 − − r r h

Definiendo una nueva variable, ψ como: ψ=

h µ − r h

La anterior ecuación se convierte en: −

dψ p 2 = α − ψ2 dθ

La cual es una ecuación diferencial con solución: ψ = α cos (θ + λ)

(2.13)


12

Donde λ es la cosntante de integración. Usando (2.12) y (2.13) se obtiene: r h µ µ2 − = + 2E cos (θ + λ) r h h2 Despejando de esta ecuación a r: r= µ h2

1 q 2 + h1 2E + µh2 cos (θ + λ)

(2.14)

La cual se conoce como la ecuación de trayectorias porque indica el valor de r en función de θ. En astrodinámica, θ recibe el nombre de anomal´ıa verdadera [6, 5]. Lo interesante de esta ecuación es que corresponde a la ecuación generalizada de las cónicas en coordenadas polares [23]. Esto implica que no solo el movimiento se restringe a un plano, sino que además la trayectoria es alguna de las secciones cónicas, con la forma precisa determinada por las constantes h, E, λ. Es importante ver que este resultado es una generalización de las leyes de Kepler para el movimiento planetario. Como los cuerpos de estudio en este trabajo son asteroides y cometas, cuyas o´rbitas (al menos las encontradas hasta ahora [8]) son el´ıpticas (a pesar de que los cometas siguen órbitas muy excéntricas), el desarrollo que sigue se hará para este caso. No obstante, como puede consultarse en [23, 6] los otros casos son similares. Comparando (2.14) con la ecuación de una elipse [23] se llega a: s 2 2Eh2 h = a 1 − e2 , e= + 1, λ=0 µ µ2

(2.15)

De la segunda de estas ecuaciones se deduce que: 1 − e2 = −

2Eh2 µ2

Y usando esta expresión en la primera se obtiene: E=−

µ 2a

(2.16)

Lo que implica que la energ´ıa para una o´rbita el´ıptica es negativa.

2.1.3

La o ´rbita en el tiempo

Se acaba de obtener una expresión que relaciona θ (la anomal´ıa verdadera) con r. El propósito ahora es obtener la relación entre θ y el tiempo t, de tal manera que

2.1. EL PROBLEMA DE LOS DOS CUERPOS

13

impl´ıcitamente se tenga la relación r (t). Utilizando (2.15) en (2.14) se obtiene: r=

a (1 − e2 ) 1 + e cos θ

(2.17)

Usando la ecuación (2.9): dθ h = 2 dt r Y la primera ecuación de (2.15) y (2.17), esta ecuación se transforma en: p aµ (1 − e2 ) (1 + e cos θ)2 dθ = dt a2 (1 − e2 )2 La cual se puede reescribir como: √ µdt dθ 3 2 = 3 (1 + e cos θ) a 2 (1 − e2 ) 2

(2.18)

Como la integral de la izquierda es una integral el´ıptica (siempre que e 6= 1), se debe recurrir a encontrar primero la expresión para r (t) para de esta forma hallar la solución para θ. (2.10) se puede escribir, utilizando (2.16) como: r 2µ µ h2 dr = − − 2 dt r a r Multiplicando esta ecuación por r se obtiene: r µ dr r = 2µr − r2 − h2 dt a Definiendo el movimiento medio n de la trayectoria como: µ = n 2 a3

(2.19)

Con esta ecuación y con la primera ecuación de (2.15) se puede reemplazar para obtener: p dr r = na 2ar − r2 − a2 (1 − e2 ) dt Separando variables y factorizando se llega entonces a: rdr q = nadt a2 e2 − (r − a)2

(2.20)

14


Para continuar, es conveniente introducir una nueva variable ξ, que se llama anomal´ıa excéntrica [6]. Está definida mediante la expresión: r = a (1 − e cos ξ)

(2.21)

De donde: dr = ae sin ξdξ Reemplazando esta ecuación y (2.21) en (2.20) se obtiene: a2 e (1 − e cos ξ) sin ξdξ p = nadt a2 e2 − a2 e2 cos2 ξ La cual, luego de simplificar se reduce a: (1 − e cos ξ) dξ = ndt Si se define la condición inicial de que en el instante t0 en el que la part´ıcula pasa por el pericentro de la o´rbita se cumpla ξ = 0, la ecuación anterior, luego de ser integrada se convierte en: ξ − e sin ξ = n (t − t0 )

(2.22)

Que se conoce como ecuación de Kepler. (2.21) es la ecuación que otorga r (t), puesto que aunque la ecuación de Kepler es no-lineal, esta siempre se puede resolver numéricamente en función del tiempo t. Para encontrar θ (t) basta con igualar las expresiones que definen a r, es decir, (2.17) y (2.21): a (1 − e2 ) = a (1 − e cos ξ) 1 + e cos θ De donde se deduce cosξ − e cos θ = (2.23) 1 − e cos ξ Después de un proceso algebraico sencillo (ver [6, 23]) se obtiene la expresión: r θ 1+e ξ tan = tan (2.24) 2 1−e 2 La idea con la anterior disgresión en este tema básico es mostrar que en un problema de dos cuerpos, el conocimiento de unos pocos parámetros (a, e, t0 ) permite determinar la posición de la part´ıcula en el plano de su órbita para cualquier instante dado de tiempo. Ahora bien, todos los cuerpos orbitando el Sol se mueven en un plano, pero no necesariamente en el mismo, razón por la cual se hace necesario adoptar un sistema de cooordenadas fijo y estándar, en el que se pueda comparar posiciones entre asteroides y planetas, por ejemplo.

2.2. SISTEMA DE COORDENADAS

2.2

15

Sistema de coordenadas

El sistema de coordenadas más utilizado en esta área se conoce como sistema ecl´ıptico (o ecuatorial) heliocéntrico y recibe este nombre por su definición. El origen del sistema es el Sol. El plano x − y es el plano de la ecl´ıptica4 (o el plano del ecuador celeste, para el caso ecuatorial). El eje x va en dirección al punto vernal [23] y es un sistema dextrógiro. Es un sistema cartesiano, donde la unidad de longitud usual es la AU , que se define como: 1AU = 149597870700m En lo que sigue, este será el sistema de coordenadas utilizado en la descripción del movimiento de los asteroides. Sin embargo, no debe confundirse este sistema con los sistemas de coordenadas que se usan sobre la bóveda celeste para indicar la posición aparente de los cuerpos con respecto al observador en la Tierra. Estos sistemas de coordenadas utilizados sobre la bóveda celeste son sistemas angulares. Como su propósito es describir la posición aparente de los cuerpos (más no su distancia) son necesarios solo dos coordenadas. Existen muchos de estos sistemas, como el altazimutal, ecuatorial horario, ecuatorial absoluto, galáctico, etc. [25]. No obstante, en los catálogos astronómicos y en especial en los catálogos con información astrométrica se siguen las siguientes convenciones que se adoptarán en este trabajo: La unidad de tiempo es el d´ıa Juliano (JD) [25]. Esto implica que se hizo necesario la implementación de código que haga la conversión entre (JD) y tiempo universal (UT), que es la medida de tiempo a la que el ser humano está acostumbrado, y viceversa. Los códigos se muestran en los apéndices A.1 y A.2. El sistema de coordenadas utilizado es el ecuatorial absoluto (α, δ) con respecto a una fecha fija (J2000) donde α es la ascensión recta y δ la declinación. El formato utilizado para estos datos es el del MPC, que es el siguiente:

Para las ascensión recta, el formato es: HH MM SS.ddd (horas, minutos y segundos ) y para la declinación es sDD MM SS.dd (s es el signo, DD grados, MM minutos, SS.dd segundos). Las funciones que convierten estos datos a radianes, de tal forma que se pueda trabajar con ellos se muestra en los apéndices A.3 y A.4 4

La ecl´ıptica es el plano de la ´ orbita de la Tierra en su movimiento alrededor del Sol, [5].

16

2.3


Par´ ametros orbitales

Como se mencionó antes, el conocimiento de (a, e, t0 ) permite conocer la posición del cuerpo en su o´rbita para cualquier instante dado de tiempo. Si se quiere conocer la posición en el sistema de coordenadas ecl´ıptico heliocéntrico, deben incluirse 3 parámetros adicionales que junto a los anteriores se denominan parámetros orbitales. Son muy importantes porque en el problema de los dos cuerpos son constantes en el tiempo y porque permiten conocer con precisión la posición y velocidad (es decir, el vector de estado) del cuerpo estudiado. Por claridad, estos son los parámetros orbitales: la excentricidad e, el semi-eje

Figura 2.1: Parámetros orbitales de una órbita en el problema de los dos cuerpos. mayor a (o la distancia al pericentro q), la inclinación de la o´rbita con respecto al plano de referencia i, la longitud del nodo ascendente Ω, el argumento de latitud del pericentro ω, y el instante de paso por el pericentro t0 . En la Fig. 2.1 se muestra el significado geométrico de estos parámetros. Ya se definió que son a, e, t0 . La inclinación es el ángulo que forma el plano de la o´rbita con respecto al plano fundamental. El nodo ascendente es el punto sobre la o´rbita donde el cuerpo pasa de tener coordenada z negativa a positiva. La longitud del nodo ascendente Ω es el ańgulo medido sobre el plano fundamental, entre del eje x y el nodo ascendente, medido en sentido directo. El argumento de latitud del pericentro es el ángulo medido sobre el plano de la o´rbita, entre el nodo ascendente y la l´ınea de apsides (es decir, entre el nodo ascendente y el pericentro).

´ 2.3. PARAMETROS ORBITALES

17

Lo anterior quiere decir que el conocimiento de los parámetros orbitales es suficiente para determinar el movimiento de asteroides y cometas que es, en u ´ltimas, el objetivo que se persigue. En [25, 5, 23] se puede consultar las ecuaciones que hacen esto posible 5 . También es muy fácil determinar que el conocimiento del vector de estado (r, r˙ ) para un instante dado de tiempo permite hallar los parámetros orbitales y por consiguiente, el vector de estado en cualquier otro instante dado de tiempo. No se entrará en detalle en este aspecto, pero esto quiere decir que los parámetros orbitales y el vector de estado son dos representaciones diferentes de un mismo espacio: el espacio de o´rbitas (orbital space, [8])6 . Lo anterior se puede resumir en el diagrama mostrado en la Fig.2.2.

Figura 2.2: Representación esquemática de la relación entre los parámetros orbitales y el vector de estado. Ya antes se hab´ıa indicado que la ecuación de Kepler era el puente que permite pasar de los parámetros orbitales al vector de estado (que es la parte inferior del diagrama) y se nota ahora que la energ´ıa total y el momento angular son el puente que permite hacer el proceso inverso. La conclusión es bastante interesante: el objetivo es conocer los parámetros orbitales (puesto que con estos se conoce la trayectoria) pero para hacerlo se puede hacer el paso intermedio de conocer primero el vector estado para un instante en particular. En el cap´ıtulo siguiente se mostrará que los métodos de determinación de órbitas hacen esto para obtener los parámetros orbitales. Entendido lo anterior, se hizo entonces necesario escribir códigos que hagan el cambio de representación: de parámetros a vector de estado y viceversa. Para esto se implementó un algoritmo para resolver numéricamente la ecuación de Kepler. Se usó un método descrito en [6] que posee convergencia cuártica . En el apéndice A.5 se muestra el diagrama de flujo de este método por claridad y también se muestra la implementación computacional. Las funciones que transforman entre ambas representaciones se detallan en los apéndices A.6 y A.8, con el caso especial de la tierra implementado en el apéndice A.7. 5

No se incluir´ a aqu´ı esta información, no obstante en el apéndice se muestra la implementaci´ on numérica. 6 Por esta raz´ on, los par´ ametros orbitales son 6, igual al n´ umero de datos en los dos vectores tridimensionales.

18


Cap´ıtulo 3 M´ etodos de inversi´ on de o ´rbitas El propósito de este cap´ıtulo es describir cuál es el problema de inversión de o´rbitas, explicar brevemente los métodos clásicos de Laplace y Gauss que solucionan tal problema y por u ´ltimo, introducir algunos de los términos y nociones que se usarán en los cap´ıtulos posteriores.

3.1

El problema de inversi´ on de o ´rbitas

Figura 3.1: Movimiento de un asteroide alrededor del Sol y observado desde la Tierra. Tomado de [13]. Un cuerpo menor (en este caso espec´ıfico, un asteroide) P se mueve alrededor del Sol con radio vector r. Este a su vez es observado desde la Tierra E que posee radio vector R, ver Fig. 3.1. El problema de inversión de órbitas consiste en el cómo a partir de información astrométrica,1 - que provee solo coordenadas 1

En cap´ıtulos siguientes se explica en detalle qué es Astrometr´ıa y se describe el proceso para

19

20

´ ´ DE ORBITAS ´ CAPÍTULO 3. METODOS DE INVERSION

angulares de la posición relativa (usualmente ascensión recta y declinación)- se puede obtener los paramétros orbitales de P. Como se estableció en el cap´ıtulo anterior, para esto es suficiente con determinar el vector de estado para un instante dado. Para determinar el vector de estado, primero notése en la Fig. 3.1 que el vector posición está dado por: r = R + ρˆ ρ (3.1) Además, R es el vector posición de la Tierra que se conoce en cualquier instante dado de tiempo. Luego en la expresión para r, la incógnita es ρ, que es la distancia topocéntrica entre el cuerpo y la Tierra. El hecho de que a partir de astrometr´ıa no pueda obtenerse información concerniente con esta distancia, hace que el problema no tenga solución cerrada y que se deban utilizar esquemas iterativos para intentar aproximar esta cantidad. Con este panorama puede entonces decirse que los métodos de inversión de o´rbitas son procedimientos basados en algunas consideraciones f´ısicas y matemáticas que intentan obtener numéricamente el valor del vector de estado. A continuación entonces se describen dos de los principales métodos clásicos que resuelven este problema.

3.2

M´ etodo de Laplace

Def´ınase el vector posición topocéntrica del cuerpo como ρ = ρˆ ρ. Este indica la posición con respecto al observador del cuerpo estudiado. La derivada con respecto al tiempo de este vector se obtiene como: d dρ = ρˆ ˙ ρ + ρ (ˆ ρ) dt dt

(3.2)

Utilizando la longitud de arco de la curva sobre la esfere celeste s se obtiene [2]: d dˆ ρ ds (ˆ ρ) = = η νˆ (3.3) dt ds dt En la cual se ha definido el movimiento medio η y un vector unitario νˆ como: η=

ds , dt

νˆ =

dˆ ρ ds

(3.4)

ˆ , puesto que este es de magnitud Ahora bien, dtd (ˆ ρ) es un vector perpendicular a ρ unitaria. As´ı, de (3.3) se deduce que: ˆ=0 νˆ · ρ

(3.5)

obtener tales mediciones desde el Observatorio Astronómico de la Universidad Tecnológica de Pereira.

´ 3.2. METODO DE LAPLACE

21

Lo cual permite definir una base tridimensional si se agrega el vector: ˆ =ρ ˆ × νˆ n

(3.6)

Por las relaciones anteriores, es fácil verificar que: d dˆ ρ dˆ ν ˆ= ·ρ (ˆ ρ · νˆ ) − νˆ · = −1 ds ds ds dˆ ν 1 d · νˆ = kˆ ν k2 = 0 ds 2 ds ˆ Con estas propiedades de ddsν , se puede escribir este vector como [8]:

dˆ ν ˆ = −ˆ ρ + κn ds

(3.7)

Donde κ se conoce como curvatura geodésica y es una función escalar desconocida. Reemplazando (3.3) en (3.2) se obtiene que: dρ = ρˆ ˙ ρ + ρη νˆ dt Derivando esta expresión una vez más con respecto al tiempo, utilizando (3.7) y factorizando se llega a: d2 ρ 2 2 ˆ ˆ ˆ n = ρ ¨ − ρη ρ + (2 ρη ˙ + ρ η) ˙ ν + ρη κ dt2

3.2.1

(3.8)

Versi´ on geoc´ entrica

En el cap´ıtulo anterior se discutió el problema de los dos cuerpos en el cual se hace la suposición fuerte2 de que el cuerpo de estudio solo se ve afectado gravitacionalmente por el Sol y que no hay ninguna otra fuerza relevante para la trayectoria. En este caso se desarrollará el método con otra aproximación importante: se considerará que las observaciones se realizan desde el centro de la Tierra, [5, 12]. Esto se conoce como aproximación geocéntrica. Usando la Ley de gravitación de Newton, se tiene las expresiones para las aceleraciones heliocéntricas de la Tierra y del cuerpo estudiado respectivamente como: ¨ = − µ R, R R3

¨r = −

µ r r3

(3.9)

De (3.1): ¨ = −µ ¨ = ¨r − R ρ 2

R+ρ R − 3 3 r R

Pero v´ alida para deteminar ´ orbitas iniciales con buena precisión.


22

Multiplicando escalarmente esta expresión con νˆ : R+ρ R ¨ · νˆ −µ − 3 · νˆ = ρ r3 R Usando (3.8) en esta ecuación: 1 1 − µ 3 − 3 (R · νˆ ) = 2ρη ˙ + ρη˙ r R

(3.10)

En esta expresión se concluye que si se conocen los valores de ρ y r, se pueden encontrar los valores de ρ˙ 3 . ˆ Haciendo el producto punto con n: ¨·n ˆ = −µ ρ

R R+ρ − 3 3 r R

ˆ ·n

ˆ esta expresión se reduce a: De (3.8) y la definición de n 1 1 ˆ = ρη 2 κ − µ 3 − 3 (R · n) r R

(3.11)

La cual se puede escribir como: 1 1 − 3 = −ρ 3 r R

η2κ ˆ µ (R · n)

Finalmente, multiplicando a ambos lados por −R3 se obtiene: 1−

ρ R3 = C r3 R

(3.12)

η 2 κR3 ˆ ·n ˆ µ R

(3.13)

Donde C=

La ecuación (3.12) se conoce en la literatura [6] como ecuación dinámica. En la Fig. 3.1 se aprecia la relación geométrica entre los vectores R, r y ρ. En particular, notése que por la ley de cosenos se puede escribir: r2 = R2 + ρ2 + 2Rρ cos 3

(3.14)

ˆ . En la Siempre y cuando sea posible determinar, de alguna manera, los valores de η, η˙ y ν siguiente subsecci´ on se mostrar´ a c´ omo, a partir de las observaciones, pueden hacerse estimaciones de estos par´ ametros.


23

ˆ ·ρ ˆ . De (3.12) se obtiene que: Donde cos = R R R3 R2 2R3 R6 2 ρ= 1− 3 , ρ = 2 1− 3 + 6 C r C r r Al reemplazar estos valores en la ecuación geométrica (3.14) se obtiene: R2 2R3 R6 2R2 cos 3 2 2 3 r − R r =R + 2 1− 3 + 6 + C r r Cr3 Multiplicando por C 2 r6 y reordenando se obtiene: C 2 r8 − R2 r6 C 2 + 2C cos + 1 + 2R5 r3 (1 + C cos ) − R8 = 0

(3.15)

La cual es una ecuación de orden 8 para r. C´ alculo de los coeficientes El propósito ahora es describir cómo se pueden obtener todos los coeficientes que aparecen en las ecuaciones anteriores de tal manera que se pueda encontrar las variables de intéres: r, r, ˙ ρ, ρ. ˙ Como se requiere determinar seis constantes del movimiento y como de cada observación se obtiene dos datos angulares diferentes, esto quiere decir que se requieren al menos 3 observaciones distintas del asteroide para determinar sus parámetros orbitales [5, 6]. No obstante, usualmente se tiene más de estas 3 necesarias: ˆ i = (cos αi cos δi , sin αi cos δi , sin δi ) , ρ i≥3 Con estos datos de ascensión recta αi y declinación δi puede hacerse una inter˙ δ¨ polación cuadrática de tal forma que se puedan obtener los valores de α, ˙ α ¨ y δ, 4 para un tiempo medio t¯ . Con estos datos se puede obtener el valor del movimiento medio η, puesto que como se deduce fácilmente de la definición: p (3.16) η = α˙ 2 cos2 δ + δ˙ 2 La curvatura κ se puede obtener a través de la ecuación [2, 8]: 1 ¨ κ= 3 δ α˙ − α ¨ δ˙ cos δ + α˙ η 2 + δ˙ 2 sin δ η

(3.17)

Y, finalmente: η˙ = 4

1 α ¨ α˙ cos2 δ + δ¨δ˙ − α˙ 2 δ˙ cos δ sin δ η

(3.18)

En el cap´ıtulo siguiente se ver´ a que en la terminolog´ıa moderna, esto se conoce como el vector atribu´ıble del arco de observaciones [12].

24


ˆ En efecto, Con estos datos se puede cálcular los vectores unitarios νˆ , n. νˆ =

ˆ dδ ∂ ρ ˆ ˆ ˆ dα ∂ ρ dt ∂ ρ dt ∂ ρ dˆ ρ = + = α˙ + δ˙ ds ds ∂α ds ∂δ ds ∂α ds ∂δ

Lo que implica que: 1 νˆ = η

ˆ ˆ ∂ρ ∂ρ ˙ α˙ +δ ∂α ∂δ

(3.19)

ˆ se tiene que: Para n ˆ =ρ ˆ × νˆ = ρ ˆ× n

ˆ dδ ∂ ρ ˆ dα ∂ ρ + ds ∂α ds ∂δ

De donde se obtiene que, [2] ˆ= n

1 η

δ˙ ˆ α˙ cos δˆ ρδ − ρ cos δ α

! (3.20)

Con estas ecuaciones se pueden obtener todos los coeficientes y vectores que son relevantes para el método de Laplace. Obteniendo los par´ ametros orbitales Con lo anterior, se puede resumir los pasos en los que el método de Laplace resuelve el problema de inversión de o´rbitas. ˙ δ¨ para 1. A partir de las observaciones se obtiene los valores de α, α, ˙ α ¨ y δ, δ, un tiempo medio t¯. 2. Se obtiene el valores de las constantes (inluyendo C y cos ) a través de las ecuaciones (3.19), (3.20), (3.18),(3.17),(3.16),(3.13). 3. Se resuelve (3.15) para r. 4. Con cada solución r se utiliza (3.14) para obtener los valores admisibles de ρ para la r correspondiente. 5. Con (3.10) se obtiene los respectivos valores para la velocidad radial ρ. ˙ 6. Con estos valores es suficiente para determinar el vector de estado para el tiempo medio t¯. 7. Finalmente, se utiliza la función A.8 para obtener los parámetros orbitales.


3.2.2

25

Versi´ on topoc´ entrica

La forma en la que se describe el método en el apartado anterior es una versión moderna de las ideas originales de Laplace. En la época de este cient´ıfico, la suposición de que las observaciones se hacen desde el centro de la Tierra funcionaba bien puesto que los cuerpos observables para aquella época se mueven de tal forma que esta aproximación no afecta considerablemente el resultado [7]. No obstante, hoy en d´ıa se conocen cuerpos que se mueven muy cerca a la Tierra (y por consiguiente muy rápido) y para estos casos esta aproximación puede producir errores imporantes [19]. Se explicará brevemente cómo se puede corregir esta suposición en el método anterior, pero solo se mostrará el cambio en las ecuaciones principales, ya que el método completo es análogo. Como antes (ecuación (3.1)) se tiene que: ρ=r−R La cuestión ahora es considerar que R es el vector posición del observador sobre la superficie terrestre, que se relaciona con la posición del centro de la Tierra como: R = R⊕ + P

(3.21)

En la cual R⊕ es la posición heliocéntrica del centro de la Tierra y P es la posición geocéntrica del observador. De esta expresión y del modelo dinámico se obtiene la ecuación: ¨=− ρ

µr µR⊕ ¨ + 3 −P r3 R⊕

(3.22)

ˆ y utilizando (3.8): Haciendo el producto escalar con n " # ˆ⊕ ·n ˆ⊕ ·n ˆ ·n ˆ ˆ ˆ d2 ρ R R P ¨ ·n ˆ = ρη 2 κ = µ R⊕ ˆ ·n − R⊕ −P 3 −P 3 dt2 R⊕ r3 r Ahora bien, de esta expresión puede eliminarse, por tener una magnitud muy peque˜ na, el término5 ˆ ·n ˆ P P 3 r As´ı, se obtiene la ecuación dinámica para la versión topocéntrica: C 5

R3 ρ = (1 − Λn ) − 3⊕ R⊕ r

Esto es admisible puesto que, en general, r > R⊕ .

(3.23)


26 En la que:

C=

Λn =

3.3

3 ηκR⊕ ˆ µRˆ⊕ · n

3 ¨ ˆ R⊕ P·n = ˆ⊕ ·n ˆ µR

¨ ·n ˆ P µ ˆ Rˆ⊕ · n 2

(3.24)

(3.25)

R⊕

El m´ etodo de Gauss

Fue propuesto por el gran matemático alemán Carl Friedrich Gauss en 1801, [4, 3]. Este es, sin duda alguna, el método clásico de determinación de o´rbitas más importante debido a la propiedad que posee de facilitar el uso de observaciones topocéntricas [8]. Hoy en d´ıa a´ un es usual que muchos observatorios alrededor del mundo utilicen este método en versiones modernas para estimar o´rbitas preliminares [13, 9]. La primer gran diferencia de este método comparado con el de Laplace es que en este caso se usan solo tres observaciones distintas del cuerpo de estudio para determinar la o´rbita inicial. Si más observaciones están disponibles, lo que se hace es mejorar los parámetros inicialmente cálculados usando el método de min´ımos cuadrados [6, 25], procedimiento que se conoce usualmente como método de corrección diferencial. Dadas tres observaciones distintas (αi , δi ) realizadas en instantes ti se tiene, como antes: ri = Ri + ρi (3.26) El método de Laplace estima el vector de estado para el tiempo medio t¯ que resulta de interpolar las observaciones. En este caso, el método se propone estimar el vector de estado para la observación de la mitad. Es por esta razón que el método se conoce en muchas referencias como un método central6 [6]. Por notación, sean los tiempos de observación t1 < t2 < t3 . Como se dijo antes (en la sección 2.1.1), el movimiento del asteroide alrededor del Sol se da en un plano fijo, por lo tanto se puede escribir: λ1 r1 + λ3 r3 = r2

(3.27)

En la cual λ1 , λ3 ∈ R. Haciendo el producto cruz con r1 (y con r3 ) y, posteriorˆ (que es el vector unitario que apunta en mente multiplicando escalarmente por c 6

Notése que el método de Laplace también es central si se ejecuta con solo 3 observaciones distintas.

´ 3.3. EL METODO DE GAUSS

27

la misma dirección que el momento angular por unidad de masa) se obtiene: λ1 =

ˆ r2 × r3 · c , ˆ r1 × r3 · c

λ3 =

ˆ r1 × r2 · c ˆ r1 × r3 · c

(3.28)

Reemplazando (3.26) en (3.27): λ1 (R1 + ρ1 ) + λ3 (R3 + ρ3 ) = R2 + ρ2 ˆ1 × ρ ˆ 3 y utilizando la ortogonalidad de estos Haciendo el producto punto con ρ vectores, la anterior ecuación se puede reescribir como: ˆ1 × ρ ˆ 3 · [λ1 R1 − R2 + λ3 R3 ] = ρ2 (ˆ ˆ3 · ρ ˆ 2) ρ ρ1 × ρ

(3.29)

Se puede aproximar los valores de r1 y de r3 a partir de los valores de r2 y de r˙2 a partir de [27, 28, 6]: ri = fi r2 + gi r˙ 2 (3.30) Donde fi y gi son las llamadas series f y g: µt2i2 3 + O ∆t 2r23 µt2i2 gi = ti2 1 − 3 + O ∆t4 6r2 fi = 1 −

(3.31)

(3.32)

Con tij = ti − tj . De estas ecuaciones es fácil deducir que: ri × r2 = −gi c,

r1 × r3 = (f1 g3 − f3 g1 ) c

(3.33)

Reemplazando en (3.28) se obtiene: λ1 =

g3 , f1 g3 − f3 g1

λ3 =

−g1 f1 g3 − f3 g1

(3.34)

El denominador de estas expresiones f1 g3 − f3 g1 , puede expandirse usando las definiciones de las series f y g (3.31) (3.32) para obtener, luego de no considerar términos de orden superior a 3: µt231 (3.35) f1 g3 − f3 g1 = t31 1 − 3 + O ∆t4 6r2 Usando esta ecuación y (3.31) (3.32) en (3.34) se encuentra que: t32 µ 2 2 λ1 = 1 + 3 t31 − t32 + O ∆t3 t31 6r2 t21 µ 2 2 λ3 = 1 + 3 t31 − t21 + O ∆t3 t31 6r2

(3.36)

(3.37)


28

Por simplificar la notación, sea: V = ρ1 × ρ2 · ρ3

(3.38)

Usando las expresiones: t231 − t232 = t21 (t31 + t32 ) ,

t231 − t221 = t32 (t31 + t21 )

(3.39)

Y reemplazando (3.37) y (3.36) en (3.29) se llega a: ˆ1 × ρ ˆ 3 · (t32 R1 − t31 R2 + t21 R3 ) − V ρ2 t31 = ρ µ ˆ1 × ρ ˆ 3 · [t32 t21 (t31 + t32 ) R1 + t32 t21 (t31 + t21 ) R3 ] + O ∆t4 + 3ρ 6r2

(3.40)

Si no se consideran los términos O (∆t4 ) se observa que en la anterior ecuación el factor de r13 es: 2

B (R1 , R3 ) =

µ ˆ1 × ρ ˆ 3 · [(t31 + t32 ) R1 + (t31 + t21 ) R3 ] t32 t21 ρ 6

As´ı, si en (3.26) se multiplica por −

R23 B(R1 ,R3 )

(3.41)

se llega a:

V ρ2 t31 R3 A (R1 , R2 , R3 ) R23 = 32 + B (R1 , R3 ) r2 B (R1 , R3 )

(3.42)

En la cual: ˆ1 × ρ ˆ 3 · [t32 R1 − t31 R2 + t21 R3 ] A (R1 , R2 , R3 ) = R23 ρ Haciendo C2 =

V t31 R24 , B (R1 , R3 )

γ2 =

A (R1 , R2 , R3 ) B (R1 , R3 )

(3.43)

(3.44)

Utilizando esta notación, (3.42) se convierte en: C2

ρ2 R3 = γ2 − 32 R2 r2

(3.45)

Que se conoce como la ecuación dinámica del método de Gauus [8]. La ecuación geométrica es la misma que la del método de Laplace. Luego con ambas ecuaciones pueden obtenerse todos los valores posibles de r2 . Esta es una descripción básica de la componente fundamental del método de Gauss. No obstante, para obtener el vector r˙ 2 (de tal forma que se puedan obtener los parámetros orbitales) debe buscarse por un método adicional. Existen muchos de estos [28, 6], no obstante solo se describirá brevemente la fórmula de Gibbs.

´ 3.3. EL METODO DE GAUSS

29

El primer paso para utilizar la fórmula de Gibbs consiste en cálcular valores para r1 y r3 . Usando (3.36) y (3.37) se obtiene los valores de λ1 y λ3 . Haciendo el ˆ1 × ρ ˆ 2 se obtiene una ecuación lineal de la cual se producto escalar de (3.27) con ρ ˆ2 × ρ ˆ 3 se encuentra puede hallar ρ3 y del producto escalar de esta ecuación con ρ ρ1 . Con estos datos se tiene entonces valores para r1 ,r2 y r3 . La función de Gibss es entonces [28]: r˙ 2 = d2 r2 − d1 r1 + d3 r3 (3.46) Donde: di = Gi + Hi ri−3 , G1 =

i = 1, 2, 3.

t232 t21 t32 t31

(3.47) (3.48)

t21 , G2 = G1 − G3 (3.49) t32 t31 Finalmente, después de todo el procedimiento descrito antes, se ha obtenido el vector de estado para la observación central. Con esto es suficiente para dar una estimación de los parámetros orbitales. Sin embargo, algunas personas han propuesto un procedimiento iterativo para obtener mejores resultados [8]. El método consiste en utilizar un propagador de dos cuerpos para obtener los vectores r1 y r3 . Con estos valores y (3.27) se puede obtener un mejor valor para ρ2 . El proceso se contin´ ua hasta que se obtenga convergencia. Otras personas han estudiado casos en los que este procedimiento puede derivar en divergencias importantes [29]. G2 =

30


Cap´ıtulo 4 Teor´ıa de eliminaci´ on algebraica El propósito de esta parte de la tésis es describir en detalle el tema central del trabajo, para lo cual este cap´ıtulo se centra en los aspectos estrictamente matemáticos, de tal forma que en los cap´ıtulos subsiguientes se pueda hacer uso de los resultados principales de esta a´rea.

4.1

Introducci´ on a la eliminaci´ on algebraica

La teor´ıa de eliminación es una rama de la geometr´ıa algebraica [30, 31, 32] que estudia la solución a sistemas de ecuaciones polinómicas no-lineales. Como muchos aspectos de la Ciencia e Ingenier´ıa pueden ser modelados a través de ecuaciones polinómicas [33], esta a´rea ha recibido mucha atención durante las u ´ltimas décadas. Puede decirse que, en un sentido muy profundo, la eliminación algebraica busca y establece generalizaciones al método de Gauss-Jordan para sistemas de ecuaciones lineales. Se hará una breve descripción de los elementos principales de esta teor´ıa que se utilizarán más adelante en este trabajo. Un hecho notable de la Geometr´ıa algebraica consiste en la idea de transformar un problema geométrico (como lo es la solución a un sistema de ecuaciones) en un problema algebraico análogo. Uno de los elementos más importantes del a´lgebra moderna es el estudio de estructuras y las relaciones entre estas. En este caso en particular, la estructura básica es el anillo de polinomios en n variables k [x1 , . . . , xn ], que es el conjunto que agrupa todos los polinomios en n variables xi , con coeficientes en el campo k 1 . k [x1 , . . . , xn ] satisface las propiedades para ser un anillo conmutativo con unidad. A continuación se presentan las definiciones y teoremas más importantes, cuyas demostraciones pueden ser encontradas en [30, 31, 32]. 1

Esto es bastante general. No obstante, la eliminación algebraica requiere que este campo sea algebraicamente cerrado. En lo subsiguiente, k se entenderá como C.

31

´ ALGEBRAICA CAPÍTULO 4. TEORÍA DE ELIMINACION

32

Definici´ on 4.1.1 Sea I ⊆ k [x1 , . . . , xn ]. I se denomina un Ideal de k [x1 , . . . , xn ] si satisface: 1. 0 ∈ I. 2. Si f, g ∈ I ⇒ f + g ∈ I. 3. Si f ∈ I y g ∈ k [x1 , . . . , xn ] ⇒ f g ∈ I. Los ideales juegan un papel central en el a´lgebra conmutativa, y son especialmente importantes para el anillo de polinomios por el siguiente teorema que provó Hilbert [30]: Teorema 4.1.1 (Teorema de la base de Hilbert). Todo ideal I ⊂ k [x1 , . . . , xn ] tiene un conjunto generador finito. Es decir, ∃ g1 , . . . , gt ∈ k [x1 , . . . , xn ] tal que: I = Donde I =< g1 , . . . , gt > significa que dado f ∈ I existen h1 , . . . , ht ∈ k [x1 , . . . , xn ] tal que: X f = h1 g1 + . . . ht gt = gj hj j

Notése que esta noción es análoga a la de base para espacios vectoriales. Ya se introdujo el objeto algebraico principal. Ahora se definirá el objeto geométrico de intéres. Definici´ on 4.1.2 Se k un campo y sean f1 , . . . fs ∈ k [x1 , . . . , xn ]. Se define: V (f1 , . . . , fs ) = {(a1 , . . . an ) : fi (a1 , . . . an ) = 0 ∀ 1 ≤ i ≤ s} Se llama a V (f1 , . . . , fs ) la variedad af´ın definida por f1 , . . . , fs . Es importante notar que la variedad es el conjunto de puntos que satisfacen el sistema de ecuaciones: f1 (x1 , . . . , xn ) = 0 f2 (x1 , . . . , xn ) = 0 .. .. . . fs (x1 , . . . , xn ) = 0 Esta observación es bastante interesante, puesto que surge la idea de que si se quiere resolver un sistema de ecuaciones polinómicas: f1 (x1 , . . . , xn ) = c1 f2 (x1 , . . . , xn ) = c2 .. .. . . fs (x1 , . . . , xn ) = cs

´ A LA ELIMINACION ´ ALGEBRAICA 4.1. INTRODUCCION

33

Primero se puede formar un objeto algebraico: I =

(4.1)

Y determinar la variedad V (I), en donde esta variedad se define como: V (I) = {(a1 , . . . an ) : f (a1 , . . . an ) = 0 ∀ f ∈ I} Ahora bien, esto tiene sentido si y sólo si: V (I) = V () es igual a V (f1 , . . . , fs ). Este hecho es bastante fácil de establecer. Teorema 4.1.2 V (I) es una variedad af´ın. Además, si I =< f1 , . . . fs > entonces V (I) = V (f1 , . . . , fs ) La conclusión aqu´ı es bastante importante. Si se quiere resolver un sistema de ecuaciones polinómicas puede intentarse primero formar el ideal I de la ecuación (4.1) y determinar su variedad V (I), que gracias al Teorema anterior, coincidirá con las soluciones al sistema. La idea anterior solo es valiosa si el cálculo de V (I) resulta más fácil que el de V (f1 , . . . , fs ), porque de lo contrario no habr´ıa beneficio. Pues bien, en general s´ı es posible obtener de una manera mucho más sencilla a V (I). Primero se notará ´ que, como en Algebra lineal, se pueden hacer cambios de base a un ideal sin alterar la variedad que define. Lema 4.1.1 Sea I ⊂ k [x1 , . . . , xn ] un ideal y sean f1 , . . . , fs ∈ k [x1 , . . . , xn ]. Los siguientes enunciados son equivalentes. 1. f1 , . . . , fs ∈ I. 2. ⊂ I . El anterior lema permite establecer cuando un ideal es un subcojunto de otro. Con esto puede entonces definirse la igualdad de ideales. Definici´ on 4.1.3 Dos ideales I =< f1 , . . . , fs > y J =< g1 , . . . , gt > son iguales ⇔ I ⊂ J y J ⊂ I. Como se dijo anteriormente, f1 , . . . , fs y g1 , . . . , gt son dos bases diferentes para el mismo ideal. Ahora se mostrará que el cambio de base no altera la variedad que define un ideal. Teorema 4.1.3 Si f1 , . . . , fs y g1 , . . . , gt son bases del mismo ideal en k [x1 , . . . , xn ], tal que: < f1 , . . . , fs >=< g1 , . . . , gt >, entonces se tiene que: V (f1 , . . . , fs ) = V (g1 , . . . , gt )

34


Conociendo que cambiar la base de un ideal no altera la variedad que define, se puede preguntar si entre las infinitas bases diferentes para un ideal, existe alguna en especial que facilite el cálculo de V (I). Pues bien, Buchberger en 1960 [34] e Hironoka [35] de manera independiente en 1959 descubrieron una base especial para ideales que cumplen este propósito. Buchberger las llamó bases de Groebner [33, 30, 31, 32] en honor a su director de tésis. Antes de introducir su definición es importante mencionar un elemento fundamental en esta a´rea: los ordenamientos de monomios. Para tal efecto se introduce la siguiente convención. Definici´ on 4.1.4 Un polinomio f ∈ k [x1 , . . . , xn ] se puede escribir de manera abstracta como: X f= xα α α

xα1 1 ...xαnn

se denomina un monomio del polinomio y α se En la cual cada x = conoce como el vector de exponentes, [30]. Se puede reconstruir un monomio xα = xα1 1 ...xαnn de la n-tupla de exponentes α = (α1 , ..., αn ) ∈ Zn≥0 . Con estas definiciones ahora se puede definir la manera en la que se ordenarán monomios, lo cual permitirá definir el término “más grande” del polinomio. Definici´ on 4.1.5 Un orden de monomios en k [x1 , ..., xn ] es cualquier orden total n > en Z≥0 que satisface: α ≥ 0 ∀α ∈ Zn≥0 . Si α > β y γ ∈ Zn≥0 , entonces α + γ > β + γ

Cualquier orden que cumpla estas condiciones puede utilizarse. No obstante, los tradicionales y más empleados son los o´rdenes lexicográfico (lex ), lexicográfico graduado (grlex ), el lexicográfico inverso graduado (grevlex ) o combinaciones de estos [32]. Definici´ on 4.1.6 El término lider LT (f ) de un polinomio f se define como el monomio que es mayor con respecto al ordenamiento de monomios utilizado. El coeficiente lider es el coeficiente en k del término lider. Por ejemplo, usando orden lex si f (x, y, z) = 2x2 yz 4 +12xy 8 −4 entonces LT (f ) = 2x2 yz 4 . LM (f ) = x2 yz 4 . LC (f ) = 2 Con el orden lexicográfico graduado: f (x, y, z) = 2x2 yz 4 + 12xy 8 − 4 entonces LT (f ) = 12xy 8 . Con estos conceptos ahora s´ı puede definirse la noción de base de Groebner. Es importante aclarar que existen muchas maneras equivalentes a la siguiente de hacer esta definición [32].

´ A LA ELIMINACION ´ ALGEBRAICA 4.1. INTRODUCCION

35

Definici´ on 4.1.7 (Base de Groebner). F´ıjese un ordenamiento de monomios. G = {g1 , ..., gs } ⊆ I se denomina una base de Groebner (o base estándar) para I si: < LT (g1 ) , ..., LT (gs ) >=< LT (I) > Un hecho importante que aqu´ı no se describirá para no alejarse del tema central, es que todo ideal de polinomios tiene una base de Groebner [30]. Si a la definición se le agrega una condición adicional, se adquiere una condición de unicidad. Definici´ on 4.1.8 Una base de Groebner reducida para el ideal I es una base de Groebner en la que: LC (p) = 1 ∀p ∈ G. ∀p ∈ G, ning´ un monomio de p pertenece a < LT (G − {p}) >

La unicidad consiste en que: Teorema 4.1.4 F´ıjese un ordenamiento de monomios en un ideal I. Entonces I tiene una base de Groebner reducida u ńica. Es necesario saber cómo discernir si una base dada es una base de Groebner. Este criterio, conocido hoy en d´ıa como criterio de Buchberger, representa uno de los ´ resultados más importante del Algebra computacional. Teorema 4.1.5 Sea I un ideal. Una base G = {g1 , ..., gs } para I es una base de Groebner si y sólo si: G

S (gi , gj ) = 0 para todos los pares i 6= j, en la que: S (f, g) =

xγ xγ f− g LT (f ) LT (g)

Donde : xγ = LCM (LM (f ) , LM (g)) y donde S¯G denota el residuo de la división de S por G. Antes de discutir las propiedades de estas bases, debe verse cómo, dado una base cualquiera para un ideal, se puede obtener la base de Groebner con cualquier ordenamiento de monomios utilizado. Fue el propio Buchberger el que propuso un algoritmo para lograr este objetivo [33]. El algoritmo básico se basa en el teorema anterior:


36

Entrada: Un conjunto {f1 , . . . , fs } Salida: Una base de Groebner G = {g1 , . . . , gt } 0 G := F , G = φ ; 0 while G 6= G do 0 G = G; 0 for {p, q} ∈ G , p 6= q do G

0

S = S (p, q) ; if S 6= 0 then G = G ∪ {S} end end end Que si bien es matemáticamente correcto, era poco práctico por los grandes tiempos de ejecución que requer´ıa. Desde aquella época el mismo Buchberger y otros autores han propuesto métodos para optimizar el algoritmo. Las bases de Groebner tienen muchas implicaciones teóricas y aplicaciones prácticas. No obstante, como se ha mostrado desde el principio, el intéres en este caso es en las propiedades de resolver ecuaciones polinómicas. En la sección que sigue se mostrará cómo las bases de Groebner hacen que el cálculo de V (I) sea mucho más sencillo.

4.2

Los Teoremas de eliminaci´ on y extensi´ on

Antes de presentar los resultados abstractos, se mostrará un breve ejemplo para clarificar la teor´ıa. Supóngase que se quiere resolver el sistema: x2 + 2y 2 + xy = 2,

4x2 − y 2 = 1

Siguiendo las ideas expuestas anteriormente, una buena idea ser´ıa formar el ideal I = Y determinar V (I). También se dijo que las bases de Groebner eran especiales porque simplificaban esta tarea. Para ver el porqué, utilizando cualquier libreria ´ de Algebra simbólica como las funciones dentro del Kernel de Mathematica® , se obtiene la base de Groebner para este ideal utilizando orden lex. El resultado obtenido es: I =<28x − 67y + 77y 3 , 77y 4 − 130y 2 + 49> Notése que el polinomio de la derecha solo tiene la variable y, lo que implica que se pueden encontrar su ra´ıces con alguno de los métodos para polinomios en una variable y para cada solución de y se reemplaza en el primer polinomio para

´ Y EXTENSION ´ 4.2. LOS TEOREMAS DE ELIMINACION

37

obtener la respectiva solución en x. Como la variedad no cambia por el cambio de base que se hizo, las soluciones as´ı encontradas son soluciones al sistema original. Este ejemplo, si bien sencillo, ilustra el poder de la teor´ıa descrita antes en la solución de sistemas de ecuaciones polinómicas. Es importante recalcar lo similar de este procedimiento con el método de Gauss-Jordan para sistemas lineales. De hecho, es muy fácil ver que si el sistema de ecuaciones es lineal, el resultado de obtener la base de Groebner para el orden lex es la forma escalonada reducida de este método [33]. En este ejemplo puede diferenciarse dos pasos claros: el primero que consiste en obtener la base de Groebner para el ideal, en el cual se elimina variables y el segundo cuando se reemplazan soluciones de una variable para encontrar las soluciones para la otra variable. El primero se llama paso de eliminación y el segundo paso de extensión y las condiciones en las que pueden hacerse están soportadas por los teoremas que se mostrarán en esta sección, pero antes de hacerlo se deben construir ciertas definiciones. Definici´ on 4.2.1 Dado un ideal I =< f1 , . . . , fs >, el l-ésimo ideal de eliminación Il es el ideal de k [xl+1 , . . . , xn ] definido por: Il = I ∩ k [xl+1 , . . . , xn ] Es sencillo probar que Il es de hecho un ideal. También es muy claro de esta definición que el l-ésimo ideal de eliminación está formado por todos aquellos polinomios que pertenecen a I que no tienen ninguna de las variables x1 , . . . , xl . En estos términos, se puede decir que eliminar las variables x1 , . . . , xn de un sistema de ecuaciones significa encontrar polinomios que pertenezcan al l-ésimo ideal de eliminación. El ejemplo anterior mostró que obtener la base de Groebner es una manera sistemática de encontrar estos polinomios. Sin embargo, para que esto se cumpla deben satisfacerse ciertas condiciones que se enmarcan en el famoso teorema de eliminaci´ on [30]. Teorema 4.2.1 (Teorema de Eliminación). Sea I ⊂ k [x1 , . . . , xn ] un ideal y sea G una base de Groebner para I con respecto al orden lex, con x1 > . . . > xn . Entonces, para cada 0 ≤ l ≤ n el conjunto: Gl = G ∩ k [xl+1 , . . . , xn ] Es una base de Groebner para el l-ésimo ideal de eliminación Il . Resolviendo el paso de eliminación a través del cálculo de la base de Groebner correspondiente se elimina una, dos y más variables, y de esta manera se encuentran los polinomios en los ideales de eliminación que pueden utilizarse para encontrar las ra´ıces de estas variables. Estas soluciones se conocen como soluciones


38

parciales. La idea es extender estas soluciones parciales a soluciones completas mediante el reemplazo en los polinomios restantes. El paso de extensi´ on es más s´ util que el de eliminación, en especial, en este caso se requiere que el campo sobre el que se trabaja sea algebraicamente cerrado2 . En este caso en particular, el campo de intéres es C. Además, se mostrará solo el caso en el que se quiere eliminar una variable x1 , ya que esté será el caso en el que se aplique más adelante. Teorema 4.2.2 (Teorema de extensi´ on) Sea I =< f1 , . . . , fs > un ideal de C [x1 , . . . , xn ] y sea I1 el primer ideal de eliminación de I. Para cada 1 ≤ i ≤ s, las fi se pueden escribir de la forma: i fi = gi [x2 , . . . , xn ] xN 1 + hi (x1 , . . . , xn )

Donde Ni ≥ 0, gi 6= 0 ∈ C [x2 , . . . , xn ] y las hi (x1 , . . . , xn ) son funciones en las que x1 tiene grado menor a Ni . Supóngase que se tiene una solución parcial (a2 , . . . , an ) ∈ V (I1 ). Si (a2 , . . . , an ) ∈ / V (g1 , . . . , gs ) ⇒ ∃ a1 ∈ C tal que (a1 , a2 , . . . , an ) ∈ V (I). Estos son los dos teoremas principales de esta área. En principio, permiten encontrar soluciones a cualquier sistema de ecuaciones polinómicas, de manera cerrada o aproximada. No obstante, un problema importante en el paso de eliminación ha hecho que muchas personas piensen en maneras alternativas de proceder.

4.3

Complejidad computacional y el Resultante

El problema consiste en que el algortimo de Buchberger mostrado antes tiene una complejidad computacional alta. Este algoritmo exige grandes espacios en memoria y tiempos de ejecución muy altos. Y la complejidad aumenta considerablemente si el algoritmo se implementa utilizando el orden lex 3 [36, 37]. Esto ha llevado a que muchas personas se cuestionen sobre la utilidad práctica de esta teor´ıa, ya que incluso con ejemplos relativamente sencillos puede requerirse muchos recursos computacionales. Han aparecido dos formas de intentar solventar estas dificultades. Una consiste en hacer mejoras tanto teóricas como prácticas a los algoritmos [38, 39, 40, 41] y la otra, que podr´ıa llamarse la persepectiva semi-simbólica, consiste en buscar una combinación entre las herramientas teóricas mostradas antes y métodos numéricos, [42, 43]. Como a´ un se realiza demasiada investigación sobre las algoritmos y su complejidad y además su implementación es compleja, se opta por 2

Porque de otra manera algunas soluciones podr´ıan perderse [30]. Notése que es con este orden que se enuncia en teorema de eliminación. No cualquier orden es admisible para este teorema. Los o´rdenes que funcionan se denominan, com´ unmente, órdenes de eliminaci´ on. El orden lex es el ejemplo más com´ un de estos. 3

4.3. COMPLEJIDAD COMPUTACIONAL Y EL RESULTANTE

39

seguir las persepectiva semi-simbólica. La idea es buscar una manera alternativa y menos costosa de encontrar polinomios en el ideal de eliminación. El cálculo de la base de Groebner se hace precisamente porque es una manera directa de encontrar estos polinomios. Existe una manera diferente de acceder al ideal de eliminación. Este consiste en formar un polinomio, que se denomina el Resultante, el cual es especial porque está dentro del ideal de eliminación [30]. Se mostrará brevemente cuál es la definición, sin embargo se hará para el caso de solo dos ecuaciones polinómicas, puesto que este será el caso que se necesitará más adelante. P P Definici´ on 4.3.1 Dados dos polinomios f = α xα y g = β xβ , el resultante de f y g, con respecto a x2 , . . . , xn se define como: a0 b0 a1 a0 b1 b0 ... a a ... b b 2 1 2 1 . .. ... ... .. a . b 0 0 .. .. Resf,g (x2 , . . . , xn ) = (4.2) . a1 . b1 a1 bm . . .. .. al bm ... ... al bm En donde los ai y los bi son los coeficientes polinómicos de f y g, si estos polinomios se ven como polinomios en la variable x1 : f = a0 xl1 + . . . + al ,

a0 6= 0

g = b0 x m 1 + . . . + bm ,

b0 6= 0

Como se puede ver claramente, el resultante es un polinomio que no contiene la variable x1 . El siguiente teorema garantiza que este polinomio está en el primer ideal de eliminación I1 . P P Teorema 4.3.1 Sean f = α xα y g = β xβ dos polinomios en k [x1 , . . . , xn ]. El resultante de f y g, con respecto a la variable x1 , Resf,g (x2 , . . . , xn ) ∈ I1 .

El c´ alculo del Resultante El resultante podr´ıa obtenerse haciendo el determinante de manera simbólica. Sin embargo, este procedimiento puede resultar, como antes, en tiempos de ejecución muy altos. Por tal razón, la forma más efectiva de calcularlo es a través de una evaluación e interpolación numérica [19, 20].

40


Lo que se hace es evaluar el resultante en unos puntos seleccionados, utilizando la propiedad de que: det (Resf,g (x2 , . . . , xn )) |x=xk = det (Resf,g (x2 , . . . , xn ) |x=xk )

(4.3)

Lo cual implica que se evaluan cada uno de los ai y bi de la ecuación (4.2) y posteriormente se cálcula el determinante numéricamente. Con los puntos evaluados del resultante se puede realizar una interpolación, lo cual permite conocer con una precisión aceptable el valor de los coeficientes numéricos de este polinomio. Con esto ya se pueden obtener sus ra´ıces. Ahora bien, existen muchas formas de hacer la evaluación e interpolación de polinomios. No obstante, la evaluación e interpolación se hará utilizando la transormada discreta de Fourier y su inversa, ya que es un método mucho más eficiente que los tradicionales [44]. En el cap´ıtulo siguiente se mostrará esto de manera más expl´ıcita. Para ver la diferencia en eficiencia, por ejemplo, una evaluación com´ un 2 tiene una complejidad de O (n ), mientras que si se utiliza DFT la complejidad se reduce a O (n log n). En este punto ya se puede hacer un breve resumen de la aproximación que se empleará en los capitulos posteriores para resolver ecuaciones polinómicas no lineales. Se utilizará eliminación algebraica, sin embargo, el paso de eliminación se hará por razones de eficiencia mediante el resultante, el cual se calculará utilizando DFT e IDFT. Para hacer uso de esta teor´ıa, se hizo necesario la implementación de rutinas que ejecuten procedimientos estándares entre polinomios como la suma B.1, producto B.2, simplificación B.3, etc. Para tal efecto, cada polinomio se representará en Matlab a través de una matriz, en la que cada fila representa los monomios. Por ejemplo, el polinomio f = 4x2 y 3 − 8y + 2x4 y + 1 se representa como:   4 2 3  −8 0 1   f =  2 4 1  1 0 0 Se implementaron rutinas para evaluar polinomios en una B.4 y dos variables B.5, también se implementó el algoritmo FFT B.6 y el IDFT B.7 y una rutina para obtener el resultante de manera numérica B.8.

Cap´ıtulo 5 Fundamento te´ orico del m´ etodo 5.1

Vectores atribu´ıbles

En este cap´ıtulo se introducirá y describirá en detalle el método basado en los invariantes del movimiento [8, 19, 20]. Como se mencionó antes, la idea es poder usar toda la información contenida en un arco de observaciones. Este se define como un conjunto de observaciones astrométricas: (ti , αi , δi )

i = 1, ..., m

Donde, como antes, αi representa mediciones de ascensión recta, δi las declinaciones y m ≥ 3 es el n´ umero total de observaciones disponibles. El propósito es al igual que en los métodos clásicos, obtener una estimación del vector de estado para un tiempo dado. Recordando del cap´ıtulo 3 que estos vectores están dados por: ˙ + ρˆ ˆ α α˙ + ρ ˆ δ δ˙ r = R + ρˆ ρ, r˙ = R ˙ρ + ρ ρ (5.1) ˙ es el vector de estado de la Tierra (que son datos conocidos), α, Donde R, R ˙ δ˙ son las derivadas con respecto al tiempo de las coordenadas angulares. Notése que ˙ de acuerdo a lo visto en el cap´ıtulo 3, si se pueden determinar los valores de α, ˙ δ, las u ńicas incógnitas en (5.1) son la distancia topocéntrica ρ y la velocidad radial ρ. ˙ Aqu´ı entonces se puede establecer que los métodos basados en los invariantes del movimiento tienen como objetivo encontrar ecuaciones polinómicas en las incógnitas ρ, ρ˙ de tal forma que se pueda resolver numéricamente para encontrar el vector de estado[18, 19, 7]. Por lo anterior, los datos se entienden de una manera radicalmente diferente. En lugar de utilizar cada dato de manera individual, lo que se hace es definir un vector, que se conoce en la literatura como vector atribu´ıble [12], que representa 41

´ ´ CAPÍTULO 5. FUNDAMENTO TEORICO DEL METODO

42

un arco completo. Este vector se define como: π π Λ = α, δ, α, ˙ δ˙ ∈ [−π, π) × − , × R2 2 2

(5.2)

Y se obtiene numéricamente a partir de un arco de observaciones (ti , αi , δi ) , i = 1, ..., m mediante una interpolación lineal (si m = 2) o cuadrática (si m > 2) para el tiempo medio t¯ [8]. Al proceso de interpolación para encontrar el atribu´ıble Λ se le puede asociar una matriz de covarianzas ΓΛ que es importante para estimar la incertidumbre de la órbita calculada [19]. En el apéndice A.9 se cálcula el vector posición de la tierra y en el apéndice A.10 se muestra el código de la funci´ on ˙ para implementada que obtiene el vector atribuible y el vector de estado R, R el tiempo medio t¯. Interpolaci´ on de Poincar´ e ˙ es el vector de estado de la Tierra. Sin embargo, aqu´ı hay Como se dijo, R, R que mencionar dos formas diferentes de obtener estos vectores. En la aproximación geocéntrica, estos vectores representan el centro de la Tierra, puesto que se considera que las observaciones se realizan desde este punto. La perspectiva topocéntrica, por el contrario, asume que el vector posición debe ser calculado mediante la suma de la posición heliocéntrica del centro de la Tierra más la posición geocéntrica del observador: R = R + Robs Y de manera análoga para el vector velocidad, para lo cual es necesario conocer las constantes de paralaje del lugar de observación1 . Además de esto, Poincaré sugirió en 1909 [8, 19], que la posición y velocidad geocéntrica deben ser obtenidos mediante la misma interpolación que se utiliza para obtener los vectores atribuibles y no utilizando fórmulas exactas. Estas aproximaciones entregan valores ligera˙ mente diferentes del vector R, R , y como estos datos se usan en las ecuaciones posteriores, los parámetros orbitales también serán diferentes para cada caso. En las secciones posteriores se muestran los resultados obtenidos con ambas aproximaciones.

5.2

Las integrales del movimiento

Ahora se usa la conservación del momento angular y la energ´ıa del sistema para escribir ecuaciones que se puedan resolver para ρ y ρ, ˙ puesto que seg´ un lo anterior esto es suficiente para obtener el vector estado del cuerpo en estudio. 1

El MPC lista el valor de estas constantes para cada observatorio numerado, ver: http://www.minorplanetcenter.net/iau/lists/ObsCodesF.html

´ 5.3. LAS ECUACIONES DEL METODO

43

El momento angular por unidad de masa es el producto r × r˙ , el cual, utilizando las ecuaciones (5.1) se puede escribir como: c (ρ, ρ) ˙ = r × r˙ = Dρ˙ + Eρ2 + Fρ + G

(5.3)

En la cual: ˆ D=R×ρ

(5.4)

˙ρ × ρ ˆ α + δˆ ˆ δ = ηˆ E = αˆ ˙ρ × ρ n ˙ ×ρ ˙ ˆ α + δR ˆδ + ρ ˆ×R F = αR ˙ ×ρ

(5.5)

˙ G=R×R

(5.7)

(5.6)

Notése que en esta expresión, las u ńicas incógnitas son ρ y ρ, ˙ puesto que los demás datos se pueden extraer del atribu´ıble. También se puede escribir la ecuación polinómica para la energ´ıa de tal forma que dependa solo de ρ y ρ˙ [19] : 2k 2 p 2E (ρ, ρ) ˙ = ρ˙ + c1 ρ˙ + c2 ρ + c3 ρ + c4 − ρ 2 + c5 ρ + c0 2

2

(5.8)

En la cual k es la constante de Gauss y: ˙ ·ρ ˆ c2 = η 2 c0 = |R|2 c1 = 2R 2 ˙ ˙ ρδ ˙ · αˆ ˆ c4 = R c3 = 2R ˙ ρα + δˆ c5 = 2R · ρ

(5.9) (5.10)

La ecuación de energ´ıa (5.8) se puede escribir como: 2k 2 √ 2E (ρ, ρ) ˙ =F− G Donde: F = ρ˙ 2 + c1 ρ˙ + c2 ρ2 + c3 ρ + c4

(5.11)

G = ρ 2 + c5 ρ + c0

(5.12)

Ahora bien, un atribuible y estas ecuaciones no es suficiente para resolver el problema de encontrar ρ y ρ, ˙ puesto que no se conocen los cantidades conservadas E ni c. Por esta razón, si se quiere resolver el problema se deben usar dos atribuibles distintos de tal forma que se puedan eliminar estas incógnitas.

5.3

Las ecuaciones del m´ etodo

Para escribir un sistema deecuaciones en incógnitas ρ y ρ˙ deben utilizarse dos ˙ atribuibles distintos Λ1 = α1 , δ1 , α˙ 1 , δ1 y Λ2 = α2 , δ2 , α˙ 2 , δ˙2 obtenidos en tiempos t¯1 y t¯2 . Por consisitencia y confiabilidad de los resultados, cada atribuible


44

deber´ıa ser obtenido de un mismo observatorio [45, 17]. Como en principio se supone que los atribuibles son del mismo cuerpo, el momento angular en ambos instantes debe ser el mismo. Denotando con un sub´ındice 1 las cantidades referentes al tiempo t¯1 y con 2 al tiempo t2 , esto quiere decir que: c1 = c2 O, utilizando (5.3): D1 ρ˙ 1 + E1 ρ21 + F1 ρ1 + G1 = D2 ρ˙ 2 + E2 ρ22 + F2 ρ2 + G2 De esta ecuación se deduce que: D1 ρ˙ 1 − D2 ρ˙ 2 = E2 ρ22 + F2 ρ2 + G2 − E1 ρ21 − F1 ρ1 − G1

(5.13)

D = D1 × D2

(5.14)

Definiendo: Y haciendo el producto de (5.13) con D: D · E2 ρ22 + F2 ρ2 + G2 − E1 ρ21 − F1 ρ1 − G1 = 0 Esta ecuación de orden dos, con incógnitas ρ1 y ρ2 define la primera ecuación del método: q (ρ1 , ρ2 ) ≡ (D · E2 ) ρ22 + (D · F2 ) ρ2 − (D · E1 ) ρ21 − (D · F1 ) ρ1 + D · (G2 − G1 )

(5.15)

También se puede ver que dado el valor de ρ1 y ρ2 , es sencillo obtener los valores de ρ˙ 1 y ρ˙ 2 a partir de (5.13). En efecto, si en (5.13) se hace el producto cruz con D1 y con D2 y luego se multiplca escalarmente por D se obtiene: ρ˙ 1 =

(E2 ρ22 + F2 ρ2 + G2 − E1 ρ21 − F1 ρ1 − G1 ) × D2 · D D2

(5.16)

(E2 ρ22 + F2 ρ2 + G2 − E1 ρ21 − F1 ρ1 − G1 ) × D1 · D (5.17) D2 Para encontrar otra ecuación en variables ρ1 y ρ2 , se utiliza la conservación de la energ´ıa. Si se utiliza (5.16) y (5.17) en (5.11), la ecuación para la energ´ıa queda con incógnitas ρ1 y ρ2 . As´ı, igualando la energ´ıa para los dos atribu´ıbles: ρ˙ 2 =

2k 2 2k 2 = F2 (ρ1 , ρ2 ) − p F1 (ρ1 , ρ2 ) − p G1 (ρ1 ) G2 (ρ2 ) √ Multiplicando a ambos lados por G1 G2 y elevando al cuadrado: p (F1 − F2 )2 G1 G2 − 4k 4 (G1 + G2 ) = −8k 4 G1 G2

´ A LAS ECUACIONES 5.4. SOLUCION

45

Elevando una vez más al cuadrado esta expresión se llega a la ecuación polinómica: 2 p (ρ1 , ρ2 ) ≡ (F1 − F2 )2 G1 G2 − 4k 4 (G1 + G2 ) − 64k 8 G1 G2 (5.18) Cuyo grado total es 24, [2] De esta manera, se han encontrado dos ecuaciones: (5.15) y (5.18) en las incógnitas ρ1 y ρ2 . Si este sistema de ecuaciones puede ser resuelto, se conocerán los valores posibles de todos los pares (ρ1 , ρ2 ). Con estos datos y (5.16) y (5.17) puede encontrarse los valores posibles de (ρ˙ 1 , ρ˙ 2 ). Ello implica que para cada tiempo t¯1 y t¯2 se tiene los siguientes datos: α1 , δ1 , α˙ 1 , δ˙1 , ρ1 , ρ˙ 1 , α2 , δ2 , α˙ 2 , δ˙2 , ρ2 , ρ˙ 2 (5.19) Usando estos datos puede obtenerse, a través de (5.1), el vector de estado para cualquiera de los dos tiempos y, por consiguiente, los parámetros orbitales. Los datos de la forma (5.19) se conocen en la literatura como elementos orbitales atribuibles,[8]. Como se mencionó anteriormente, estos datos son una representación diferente del espacio de órbitas. La dificultad recae ahora en la solución del sistema de ecuaciones (5.15) y (5.18), puesto que es un sistema 2 × 2 no lineal de orden total 48. La solución de este problema matemático requiere la utilización de las herramientas sofisticadas mostradas en el cap´ıtulo anterior.

5.4

Soluci´ on a las ecuaciones

Utilizando la teor´ıa del cap´ıtulo anterior, se tiene que resulta de las ecuaciones (5.15) y (5.18) es: a20 0 b2 0 a19 a20 b1 b2 .. .. . . b0 b1 Resp,q (ρ2 ) = .. .. . 0 b0 . a0 a1 ... ... 0 a0 0 0

que para este caso el resultante · · · · · · 0 0 · · · 0 . b2 · · · .. . b1 · · · .. .. . b0 b1 0 0 b0

(5.20)

Donde los ai = ai (ρ2 ) y bi = bi (ρ2 ) son los coeficientes de (5.18) y (5.15) respectivamente, si se ven estos polinomios como polinomios en la variable ρ1 . Para utilizar el método de la transformada DFT e IDFT lo que se hace primero, para el caso DFT, es evaluar los polinomios ai y bi (que aparecen como entradas k en el resultante) en las 64 ra´ıces de la unidad ωk = e2πi 64 , con k = 1, . . . 64 2 . Al 2

Se evaluan 64 puntos puesto que se espera utilizar el algoritmo FFT y el grado máximo del resultante es 48, [19]


46

hacer estas evaluaciones, se puede obtener fácilmente el valor del resultante para cada punto, puesto que como se mencionó en el cap´ıtulo anterior: det (Resp,q (ρ2 )) |ωk = det (Resp,q (ρ2 ) |ωk ) Finalmente, se aplica el algoritmo IDFT a los puntos evaluados y la respuesta a este proceso serán valores aproximados a los coeficientes de (4.2). Con esto ya se puede aplicar alguno de los métodos para resolver ecuaciones en una variable.

5.4.1

Selecci´ on de las soluciones

Solucionar (5.15) y (5.18) produce un conjunto de pares (ρ1 , ρ2 ). Cada par de soluciones permite obtener (ρ˙ 1 , ρ˙ 2 ) a través de (5.16) y (5.17). En principio, con estos datos es suficiente para obtener los parámetros orbitales. No obstante, al solucionar estas ecuaciones se encontrarán pares (ρ1 , ρ2 ) que no tendrán significado f´ısico y que aparecen como consecuencia del proceso de manipulación algebraica. Se deben seleccionar como soluciones admisibles solo aquellos pares que satisfagan las ecuaciones iniciales (5.8) y (5.3). En realidad, como las soluciones se obtienen de manera aproximada, ninguna de las soluciones satisfacerá de manera exacta estas ecuaciones. Se aceptan como posibles soluciones solo aquellos pares que satisfagan (5.3) (5.8) con un nivel muy alto de exactitud. Sin embargo, existen muchos casos en que dos o más de las soluciones pasan por este filtro, por lo tanto, es necesario establecer un criterio adicional para seleccionar la solución correcta.3 Sea (ρ1 , ρ2 , ρ˙ 1 , ρ˙ 2 ) una solución. Con estos datos se puede determinar los parámetros orbitales para cualquiera de las dos épocas t¯1 o t¯2 . De hecho, como se usó la conservación de la energ´ıa y momento angular para obtener la solución, los parámetros a, e, i, Ω serán los mismos para ambas épocas. Sin embargo habrán divergencias en los parámetros ω y l: ∆1,2 (Λ1 , Λ2 ) = (ω1 − ω2 , ∆l ) Con ∆l = l1 − (l2 + n (t¯1 − t¯2 )) ,

3

n = ka− 2

Donde ∆1,2 se conoce como condiciones de compatibilidad. Notése que se hizo expl´ıcita la dependencia de esta expresión en ambos atribuibles. Definiendo: A = (Λ1 , Λ2 ) Por la ley de propagación de covarianzas, se tiene la matriz de covarianzas marginal para las condiciones de compatibilidad: T ∂∆1,2 ∂∆1,2 Γ∆1,2 = ΓA (5.21) ∂A ∂A 3

Se asume que siempre habr´ a solo una solución para los casos estudiados. Sin embargo, esta suposici´ on no siempre es cierta, [46]

´ NUMERICA ´ 5.5. IMPLEMENTACION

47

Haciendo C = Γ−1 on [8] como: ∆1,2 , se define la norma de identificaci´ k∆1,2 k2∗ = ∆1,2 C∆T1,2

(5.22)

El criterio adicional consiste en seleccionar la solución (ρ1 , ρ2 , ρ˙ 1 , ρ˙ 2 ) que arroje el valor más peque˜ no para la norma de identificación.

5.5

Implementaci´ on num´ erica

Descrita la teor´ıa del método basado en las integrales de movimiento, y con lo estudiado en cap´ıtulos anteriores, se implementó un conjunto de subrutinas en Matlab ® que buscan obtener los parámetros orbitales a partir de astrometr´ıa. Este incluye rutinas para el acondicionamiento de los datos, transformación entre sistemas de coordenadas cartesianos y esféricos, rutinas para la transformación entre el espacio de fases a los parámetros orbitales, rutinas para el manejo algebraico de ecuaciones polinómicas y, finalmente, una rutina principal que implementa el procedimiento descrito en las secciones precedentes, cuya implementación se muestra en el apéndice C. Si la convergencia se consigue4 , este procedimiento está dise˜ nado para obtener un conjunto de parámetros orbitales para la información suministrada como entrada. Sin embargo, aparecieron algunos problemas con las primeras pruebas de este código. En general, no se estaban obtenido respuestas aceptables para los parámetros orbitales, sin importar la calidad ni cantidad de los datos que se usaran como entrada. Los problemas están asociados a que los coeficientes numéricos de los polinomios p y q (5.15) y (5.18) son valores extremadamente peque˜ nos. Además, no es adecuado hacer un escalamiento de las cantidades (mediante un cambio de unidades, por ejemplo), puesto que como menciona Gronchi (comunicación personal, Junio 5, 2015), el problema es que los coeficientes no son uniformemente peque˜ nos. Ante este panorama, se hizo necesario usar las funciones de ar´ıtmetica de alta precisión del toolbox de matemática simbólica 5 para obtener resultados confiables. En promedio, los tiempos de ejecución están alrededor de TE ≈ 6s en un computador personal. Para mayor claridad, a continuación en la Fig. 5.1 se presenta el diagrama de flujo del método principal. 4

Con convergencia, se quiere decir que el algoritmo llegue hasta el u ´ltimo proceso de ejecución. Fuente: http://www.mathworks.com/products/symbolic/features.html#variable-precisionarithmetic. 5

48


Constantes de paralaje del observatorio (ρ cos ϕ, ρ sin ϕ)

(αj,i , δj,i , tj,i ), con j = 1, . . . , n.

Interpolación de Poincar´ e ˙ para hallar R, R

Cálcular Λi para i=1,2. y sus respectivas matrices de covarianza ΓΛ

Obtener los coeficientes de p y de q

Determinar los polinomios ai y bi y evaluarlos en las ra´ıces de la unidad k ωk = e2πi 64 , con k = 1, . . . 64. Se cálculan las normas de identificación. Se selecciona los parámetros orbitales con menor valor para este parámetro

Evaluar el resultante para cada k = 1, . . . , 64 utilizando el hecho de que: det (Resp,q (ρ2 )) |ωk = det (Resp,q (ρ2 ) |ωk )

Se obtienen los parámetros orbitales

Aplicar IDFT para obtener los coeficientes del resultante

Obtener el vector de estado (r, r˙ )

Encontrar las soluciones reales ρ2,m a este polinomio, donde m es cuando mucho 48.

Se eliminan los pares (ρ1 , ρ2 ) que no satisfacen las ecuaciones de energ´ıa y momento.

Para cada m, encontrar ρ1,m mediante la ecuación q (ρ1 , ρ2,m ) = 0. Se obtiene un par ρ1,m y ρ∗1,m

Se selecciona la solución que de∗ el menor valor de |p (ρ1,m , ρ2,m )| y p ρ1,m , ρ2,m

Es ρ1,m > 0 o ρ∗1,m > 0?

s´ı

Figura 5.1: Diagrama de flujo del algoritmo principal

5.6. RESULTADOS INICIALES

5.6

49

Resultados iniciales

Para probar la correctitud del algoritmo antes de usar datos medidos desde el observatorio astronómico de la Universidad Tecnológica de Pereira y además estudiar la influencia de las aproximaciones geocéntricas y topocéntricas, se probó el método con cuerpos pertenecientes a las familias MBA y NEO, cuyos datos fueron medidos por otros observatorios. En la Tabla 5.1 se listan todos los objetos utilizados como prueba. Además se muestra la procedencia de los datos usados para cálcular el par de atribuibles para cada cuerpo. La información astrométrica se obtuvo del motor de b´ usqueda de observaciones del MPC6 en el cual se muestran todos los datos tomados por observatorios que poseen código MPC. Con t1 se indica el tiempo inicial del arco de observaciones medido por el observatorio correspondiente, que se identifica con el código entregado por el MPC. Asimismo, t2 representa el tiempo inicial del segundo atribuible.

Tabla 5.1: Tiempos y lugares de observación de los asteroides utilizados como objeto de estudio. Asteroide MBA (101878) MBA (675) MBA (914) MBA (1238) MBA (1816) NEO (99942) NEO (11500) NEO (15745) NEO (25330) NOE (65679)

t1 (JD)

Obs.

2454108.63919 G96 2456878.31659 L33 2451713.47464 117 2457022.01582 G45 2457127.78775 703 2453175.67015 695 2451109.46667 561 2451716.43194 610 2451907.58476 704 2447872.54673 675

t2 (JD) 2453999.03472 2456901.72601 2454223.858414 2457105.8128 2457187.76007 2453357.92318 2451151.94163 2451747.68264 2452007.84594 2447828.76856

Obs. 568 703 196 G45 703 E12 428 739 649 010

Procedencia de los datos utilizados. Se muestra el tiempo inicial de cada arco de observación y se se˜ nala el código MPC de los observatorios que realizaron las mediciones. El algoritmo se ejecutó entonces para estos datos y se obtuvo convergencia en cada caso. En la Tabla 2 se muestran los parámetros orbitales nominales (los entregados por el MPC) y el resultado de aplicar el algoritmo en la versión geocéntrica (Geoc.) y en la versión topocéntrica (Topo.). 6

Ver: http://www.minorplanetcenter.netdb search

50

´ ´ CAPÍTULO 5. FUNDAMENTO TEORICO DEL METODO Tabla 5.2: Comparación entre los parámetros orbitales.

Asteroide

101878

675

914

1238

1816

99942

11500

15745

25330

65679

Datos

a(AU )

Par´ ametros orbitales e i() Ω ()

ω ()

MPC Geoc. Topo.

2.2380734 2.2906675 2.2613360

0.1834232 0.2137229 0.2006867

0.60223 0.61409 0.61380

156.26944 156.87946 157.10655

145.11560 145.78970 144.24619

MPC Geoc. Topo.

2.7704278 2.7933342 2.7751505

0.2007596 0.2103692 0.2089904

9.78383 9.77293 9.78112

263.26851 263.21830 263.28802

152.10953 151.59498 151.50585

MPC Geoc. Topo.

2.4577207 2.4527812 2.4548271

0.2146690 0.2105730 0.2111170

25.20561 25.21534 25.22273

255.80617 255.81562 255.83729

49.16424 48.52473 48.46589

MPC Geoc. Topo.

2.6658816 2.6966733 2.6498104

0.1416187 0.1615664 0.1448617

12.15504 12.07645 12.05276

51.95447 51.18613 52.72733

91.86598 139.73490 96.97644

MPC Geoc. Topo.

2.3386797 2.2967588 2.3178117

0.2179828 0.2030637 0.2102513

26.13863 26.11064 26.19727

153.36043 153.23408 153.67579

340.78454 334.29417 339.60845

MPC Geoc. Topo.

0.9221174 0.8919441 0.9218615

0.1912151 0.1987364 0.1913915

3.33062 2.94893 3.33999

204.20352 209.75485 204.52951

126.45688 115.03299 126.27401

MPC Geoc. Topo.

1.0798741 0.9985168 1.0689727

0.3558841 0.1716986 0.3327292

10.30928 4.68421 9.69861

234.45954 236.11125 234.83437

289.41972 290.97121 289.21493

MPC Geoc. Topo.

1.7197574 1.7460786 1.7074082

0.2550697 0.2624876 0.2519446

14.42786 14.75843 14.24257

132.64179 132.70280 132.85465

140.54859 141.08076 139.85713

MPC Geoc. Topo.

1.5405768 1.6319274 1.5678689

0.3707986 0.3991019 0.3795623

14.32802 14.74747 14.40693

50.61908 49.95111 50.54001

85.97887 83.44736 87.20927

MPC Geoc. Topo.

0.9150868 0.8934026 0.9135664

0.2647216 0.3618796 0.2759978

1.29952 1.83796 1.34754

178.23836 181.51654 178.92088

14.89192 14.40954 14.79512

Resultados obtenidos al ejecutar el método. Se comparan los valores nominales de los parámetros orbitales (MPC) con respecto a los obtenidos en el caso geocéntrico (Geoc.) y topocéntrico.

´ 5.7. ANALISIS DE ERRORES

51

Con los parámetros orbitales (los topocéntricos) se puede realizar una propagación de la órbita para encontrar la trayectoria heliocéntrica de cada cuerpo. En la Fig. 5.2 se muestran las trayectorias obtenidas para los cuerpos MBA y en la Fig. 5.3 se muestra las trayectorias para los objetos NEO.

Figura 5.2: Trayectorias heliocéntricas para los cuerpos de estudio de la familia MBA.

5.7

An´ alisis de errores

Analizar estos resultados no es una tarea sencilla. Esto se debe en parte a que realizar comparaciones entre tantos parámetros puede ser una tarea poco clara [47]. Una manera eficiente de analizar estos resultados es reducir el n´ umero de parámetros de error de 6 a 2. Estos son el parámetro de error en la forma de la o´rbita d y en la orientación de la o´rbita Φ [48]. Si los parámetros nominales son (a, e, i, Ω, ω, θ) y los parámetros aproximados son (a∗ , e∗ , i∗ , Ω∗ , ω ∗ , θ∗ ), los errores se estimarán mediante: q d = (a − a∗ )2 + (b − b∗ )2 (5.23) √ Donde b = a 1 − e2 es el semieje menor. El parámetro de error en la orientación de la o´rbita Φ se define como: 1 cos Φ = tr CC ∗T − 1 (5.24) 2


52

Figura 5.3: Trayectorias heliocéntricas obtenidas para los cuerpos NEO. En la cual: h C = R (ω + θ) P (i) R (Ω) = ˆr

ˆ × ˆr h

iT

ˆ h

(5.25)

En la que P , R son la matrices de rotación sobre los ejes x e z, respectivamente y los vectores en la matriz al lado derecho se entienden como vectores columna. Notése que el ańgulo Φ es el ańgulo principal de la matriz correctiva entre las ˆ matrices C y C. Este análisis se implementa en el apéndice A.11 para lo cual fue necesario implementar en A.6 el cálculo del vector de estado para cualquier cuerpo y en A.7 para la Tierra. Al aplicar estas ecuaciones a los resultados obtenidos, se obtuvieron los errores mostrados en la Fig. 5.4. Primero notése en la gráfica del error en la forma que, a excepción de los cuerpos 11500 y 25330, los errores son inferiores a 65 × 10−3 AU , que son valores esperados al utilizar el problema de los dos cuerpos como modelo dinámico [19]. Otra observación importante se obtiene de analizar esta gráfica: en general, tanto para el error en la forma como para el error en la orientación, el uso de la versión topocéntrica del algoritmo produjo errores inferiores con respecto a la geocéntrica, y esto se nota para los ejemplos de ambas familias de asteroides. Sin embargo, también es importante notar que para los asteriores de la familia NEO, las diferencias entre ambas aproximaciones es mayor. Por ejemplo, si se calcula el promedio de las diferencias entre los errores de cada aproximación dG − dT se obtiene un valor de 0.02711AU para los cuerpos MBA y de 0.04950AU para los NEO, lo cual

´ 5.7. ANALISIS DE ERRORES

53

Figura 5.4: Errores obtenidos para cada cuerpo. En la parte izquierda se muestra el error en la forma de la órbita, y en la derecha, el error en la orientación de la órbita. es casi el doble que para el caso MBA. Para observar con mayor claridad los errores obtenidos, se puede aprovechar que el error en la forma es una distancia y que el error en la orientación es un ángulo, para hacer una gráfica en el plano complejo, en la que el error de cada ejemplo se representa exclusivamente por el punto: = deiΦ

(5.26)

En esta representación, el error será inferior para puntos más cercanos al origen y con menor ańgulo de inclinación con respecto al eje x positivo. Los resultados se muestran en la Fig. 5.5. En esta se puede apreciar claramente que los mejores resultados se obtienen con al aproximación topocéntrica. Lo anterior sugiere dos cosas importantes: la primera es que como muchas personas han discutido en el pasado [7, 8], para cuerpos que pertenecen a la familia MBA, la versión geocéntrica produce resultados aceptables que permiten hacer seguimiento futuro al objeto y por lo tanto, un refinamiento de la o´rbita. Esto también es cierto para los métodos clásicos de Laplace y Gauss. Segundo: para objetos pertenecientes a la familia NEO, los errores generados al usar la versión geocéntrica pueden ser muy superiores a los producidos por la versión topocéntrica, y esto puede llegar a ser determinante en la imposibilidad de realizar seguimiento al cuerpo y/o en la imposibilidad de efectuar un refinamiento de los parámetros orbitales a través de la corrección diferencial estándar [6]. Esto se debe a que muchos objetos pertenecientes a esta familia son muy tenues y se mueven muy rápido (i.e, poseen un movimiento diario muy alto), lo cual demanda mayor presición instrumental. Por lo tanto, la versión del algoritmo utilizada para

54


Figura 5.5: Combinación de ambos errores en el plano complejo. En azul se muestra los errores obtenidos con la versión geocéntrica y en rojo los de la versión topocéntrica. El error es menor entre más cerca se esté del origen. estos cuerpos deber´ıa ser siempre la topocéntrica. De aqu´ı en adelante se opta por usar la versión topocéntrica del algoritmo para las pruebas posteriores. El propósito del cap´ıtulo siguiente es mostrar resultados experimentales, con datos tomados desde el observatorio astronómico de la Universidad Tecnológica de Pereira (OAUTP).

Cap´ıtulo 6 Resultados experimentales En este cap´ıtulo se describen los procedimientos aplicados para obtener las mediciones de cada asteroide usado para probar el algoritmo. Además de esto se muestran los resultados de aplicar el método del cap´ıtulo anerior a tales datos. Obtener datos astrométricos es un proceso que puede dividirse en tres pasos que deben ejecutarse de forma adecuada para que la información obtenida tenga la calidad necesaria para obtener o´rbitas preliminares confiables. Se describirá cómo se llevaron a cabo estos pasos durante la realización de este trabajo.

6.1

Planeamiento de las observaciones

El primer paso para realizar observaciones astronómicas es el planeamiento. Este consiste en construir una lista de cuerpos que se quieren observar y sus ubicaciones en el cielo. Para hacer esta lista se debe tener en cuenta, además del intéres propio que represente la observación de estos cuerpos, la posición del cielo en el que la visibilidad es posible. En general, las regiones más adecuadas por la ubicación del OAUTP son desde el zenit hasta el sur y sur-este. En principio, y mientras se adquiere la experiencia necesaria, se elegirán cuerpos con magnitud aparante m ≤ 16 y se seleccionarán cuerpos con identificadores del MPC superiores a 4001 . Una vez analizada la condición del cielo, se utiliza el software TOOLKIT FOR CCD ASTROMETRY que forma parte del proyecto CLEA, el cual es financiado por la NSF, NASA y Gettysburg College2 . Este software está disponible para uso libre en situaciones acádemicas. 1

As´ı lo sugiere el propio MPC a observatorios que empiezan a funcionar y desean obtener el c´ odigo del observatorio. 2 Fuente: http://www3.gettysburg.edu/ marschal/clea/cleahome.html

55

CAPÍTULO 6. RESULTADOS EXPERIMENTALES

56

Después de instalar el programa, descargado e instalado las bases de datos del MPC y del observatorio Lowell, se está en capacidad de hacer b´ usquedas de cuerpos. Para esto se da clic en el botón Orbit search y luego se selecciona la base de datos. Al hacer esto aparecerá una ventana con opciones de b´ usqueda, como se muestra en la Fig. 6.1. Esta ventana, además de la fecha de intéres, ofrece varias

Figura 6.1: Ventana de b´ usqueda del programa TOOLKIT FOR CCD ASTROMETRY.

opciones para limitar o restringir los resultados. Se elegirá limitar la b´ usqueda por regiones del cielo (mediante el establecimiento de l´ımites a los valores de ascensión recta y declinación) y también por rango de magnitudes. Al presionar el botón OK, el programa automáticamente consulta en la base de datos predefinida3 y arroja una lista como resultado (ver Fig. 6.2.). Esta lista es importante porque además de mostrar los cuerpos que cumplen con las condiciones de b´ usqueda, también entrega las coordenadas de sus posiciones, las magnitudes aparentes y además ofrece información adicional como los cuerpos que poseen o´rbitas mal condicionadas, etc. De esta lista se seleccionan unos cuantos asteroides (el n´ umero por jornada de observación es de más o menos 8) que se convierten en los candidatos a ser observados en el cielo. En este punto se pasa a la siguiente etapa del proceso. 3

Siempre se trabajar´ a con la base de datos del MPC.

6.2. OBSERVACIONES

57

Figura 6.2: Lista de asteroides que satisfacen las condiciones de la ventana de b´ usqueda.

6.2

Observaciones

Los datos observacionales se obtienen con los equipos con los que cuenta el OAUTP. El telescopio es un MEADE4 f /10 LX 200 GPS de 1600 con tecnolog´ıa UHTC (Ultra-High Transmission Coatings) y viene integrado con el software AutoStar Suite® , el cual permite controlar el telescopio de manera precisa a través de sus sensores GPS y de movimiento. Este software es además bastante u ´til porque posee bases de datos de muchos catálogos estelares y además de esto, también es posible integrarle el catálogo de cuerpos menores del sistema solar. Esto implica que para centrar los asteroides seleccionados en la etapa de planeación, solo hay que buscarlos por el nombre en el software y llevar el telescopio hasta las coordenadas apropiadas. Sin embargo, para que lo anterior funcione de manera correcta el telescopio debe estar calibrado. Esto quiere decir que los motores y sensores deben estar en perfecto funcionamiento (para cual es necesario realizar varios procedimientos periódicos como el PEC, etc. 5 ). Además de esto el telescopio debe estar apropiadamente alineado (lo cual garantiza que la posición en la que apunta el telescopio coincida con la del software) y antes de cada jornada de observación debe realizarse un proceso de sincronización con algunas estrellas de referencia. La realización de todas estas actividades de mantenimiento es fundamental para preservar las cualidades del equipo. En el OAUTP este trabajo es realizado por los estudiantes y profesores investigadores, pero principalmente por el profesional que tiene como función principal esta tarea. 4

Para informaci´ on detallada del telescopio, consultar: http://www.meade.com/lx200-acf-16with-super-giant-field-tripod.html 5 Consultar el manual del equipo para ver en detalle la naturaleza de estos procedimientos.


58

El OAUTP cuenta con varias cámaras CCD6 que permiten tomar fotograf´ıas de alta calidad. Para la realización de este trabajo se utilizó la cámara SBIG ST2000 XM7 , la cual posee un sensor CCD con p´ıxeles cuadrados de 7.4µm. Antes de cada jornada de observación y después de hacer la sincronización del telescopio, se acopla la cámara y esta se enfoca mediante el movimiento del espejo del telescopio y el uso del microenfocador. Una vez realizado todos estos procedimientos se puede enviar el telescopio mediante el software al sitio de intéres. La cámara se puede configurar para que tome fotos con un tiempo de exposición variable, desde unos cuantos microsegundos hasta minutos. El tiempo de exposición dependerá de dos factores. El primero es la magnitud de los objetos en el campo. Para hacer astrometr´ıa, en el campo de la imágen capturada deben verse, además del asteroide objetivo, varias estrellas que se usan como cuerpos de referencia. Estas estrellas deben estar bien definidas (con un PSF (point-spread function) y SNR (signal to noise ratio) alto, [49, 10]), por lo cual el tiempo de exposición debe ser el suficiente para que as´ı suceda. Segundo, los asteriodes se mueven relativamente rápido en la imágen. El tiempo de exposición no puede ser tan alto como para que el cuerpo aparezca difuminado, porque esto también afecta la calidad de los datos. Esto quiere decir que el tiempo de exposición se deje elegir cuidadosamente al balancear estos factores. En este punto hay un reto instrumental adicional: como es obvio, el asteroide que se quiere medir debe quedar dentro de la imágen obtenida. Por lo tanto es necesario realizar una comparación del campo de la cámara con uno de referencia para garantizar esto, antes de empezar a tomar fotograf´ıas. Existen muchas formas de hacer esta comparación. En este trabajo y después de mucho probar, se decidió usar el software libre SAOImage DS9 8 . Con el programa instalado y al ejecutarse se sigue camino Analysis → Image Servers → ESO-DSS I/II, como se muestra en la Fig. 6.3. Con esto se le está indicando al programa que se desea buscar un campo de referencia en el survey ESO-DSS I/II. Al dar clic se mostrará una ventana en la cual se deben dar las coordenadas del centro de la imagen deseada y del tama˜ no rectangular de la imágen. Como la SBIG-ST 2000 XM produce campos de 10 × 7.5 arcominutos, el tama˜ no para comparar se hace coincidir.

6

Para ver fotograf´ıas y una lista con los equipos del OAUTP, consultar en : http://observatorioastronomico.utp.edu.co/instrumentacion.html 7 Para las especificaciones detalladas, ver: http://archive.sbig.com/sbwhtmls/st2000xm new.htm. 8 Ver: http://ds9.si.edu/site/Home.html.

6.2. OBSERVACIONES

Figura 6.3: Camino en el programa para encontrar campos de referencia.

Figura 6.4: Campo de referencia obtenido del survey ESO-DSS I/II.

59


60

Figura 6.5: Imágen obtenida del asteroide 675. En la Fig. 6.4 se muestra el campo que arroja el programa para una de las observaciones llevadas a cabo. El propósito es confirmar que se esté en el campo correcto para observar el asteroide 675. En la Fig. 6.5 se muestra la imágen tomada en el OAUTP. Notése que después de un análisis cuidadoso es fácil identificar varios asterismos 9 com´ unes en ambas imágenes, lo cual hace concluir que ambas imágenes corresponden al mismo campo. El asteroide se muestra en color rosa en la imágen obtenida desde el OAUTP. Nótese que en la imágen del campo de referencia no se puede apreciar este cuerpo, lo cual es una indicación de que esta fue tomada mucho antes de que la trayectoria del asteroide pasara por ese lugar. Una vez se tiene certeza de que el campo observado es el correcto, se procede a tomar una secuencia de imágenes. De cada una de estas se podrá extraer, si la calidad as´ı lo permite, un solo dato de la forma (α, δ, t). En la terminolog´ıa del cap´ıtulo anterior, la serie de imágenes permite extrar un arco muy peque˜ no de la órbita (αi , δi , ti ). Para probar el algoritmo se hicieron observaciones de los asteroides 675 y 654, los cuales son asteroides que pertenecen al cinturón principal (MBA) con dimensiones de aproximadamente 67.6km y 127.8km, respectivamente. En la Fig. 6.6 se muestran dos de las imágenes tomadas del asteroide 675, en las que se alcanza a apreciar el desplazamiento del asteroide en el campo. 9

Un asterismo es un patr´ on formado por varias estrellas.

6.2. OBSERVACIONES

61

De hecho, para el proceso posterior de reducción de datos es importante ver este movimiento porque de esta manera se puede identificar el o los asteroides dentro del campo con seguridad. En la Fig. 6.7 se muestran dos imágenes capturadas del asteroide 654. En esta imágen también es posible apreciar el movimiento del asteroide con respecto a las estrellas de referencia.

Figura 6.6: Imágenes tomadas del asteroide 675. En rojo se muestra el asteroide. Notése que se alcanza a apreciar el movimiento del asteroide con respecto a las estrellas de fondo.

Figura 6.7: Imágenes tomadas del asteroide 654. En rojo se muestra el asteroide. Notése que se alcanza a apreciar el movimiento del asteroide con respecto a las estrellas de fondo. Para obtener los datos astrométricos a partir de estas imágenes es importante notar que las posiciones que se pueden medir directamente sobre la imágen son la ubicación en p´ıxeles, que son coordenadas rectangulares que resultan de la proyección de la esfera celeste sobre la imágen. Es por esta razón que es importante tener un buen n´ umero de estrellas de referencia bien definidas. Es posible


62

medir las coordenadas en p´ıxeles de estas estrellas y también buscar sus coordenadas angulares en catálogos estelares con información de posición precisa, tales como UCAC-4. Con estos dos conjuntos de datos se puede ajustar un polinomio mediante un procedimiento de m´ınimos cuadrados, de tal manera que se obtenga un función que relacione posición en p´ıxeles con coordenadas angulares. Con esta función se estiman valores de (α, δ) para cualquier punto sobre la imágen, en especial, para el asteroide encontrado. Este procedimiento lo realiza Astrometrica® , del cual se tiene licencia al interior del grupo de investigación. A continuación se muestra en detalles el uso de estre programa, el cual representa la tercera etapa necesaria para obtener datos astrométricos.

6.3

Reducci´ on de datos

La idea es entonces obtener las coordenadas angulares del asteroide desde su ubicación en la imágen. Este procedimiento se ejecutó en Astrometrica, el cuál es un programa más complejo y por lo tanto a continuación se explica en detalle el procedimiento necesario para que funcione corrrectamente. El proceso de configuración inicial del programa de acuerdo con las caracter´ısticas instrumentales con las que se va a usar es de fundamental importancia, puesto que una mala configuración puede impedir el funcionamiento del mismo.

Determinaci´ on de los par´ ametros instrumentales Lo primero que debe hacerse es ingresar los datos instrumentales. Como se verá más adelante, distintas circunstancias hacen que los valores nominales difieran de los valores reales de estos parámetros. Se da clic en F ile > settings. Al hacer esto se abre un recuadro con un conjunto de botones como Observing Site, CCD, etc., que deben llenarse. Observing Site En el recuadro Observing Site se ingresan los datos geográficos y de contacto del observatorio, aclarando que si el observatorio no tiene código MPC la opción MPC CODE debe llenarse con XXX. CCD En este recuadro se ingresan los valores nominales de la cámara con la que se toman las fotos que se quieren procesar. En la Fig 6.8. Se muestra cómo llenar es-

´ DE DATOS 6.3. REDUCCION

63

Figura 6.8: Configuración de la sección CCD. tos parámetros para el caso de la cámara CCD SBIG ST-2000 XM10 . Las opciones Flip Horizontal y Flip vertical permiten organizar la orientación de la imágen que se carga al programa si es necesario. Ninguna de las dos opciones fue requerida para orientar la imágen para las imágenes obtenidas. Los parámetros restantes de esta sección conciernen al tiempo de las observaciones. Debe revisarse en los manuales de la cámara utilizada, por ejemplo, si el tiempo se guarda con el tiempo inicial, en la mitad o al final de la exposición. Program Esta es, sin duda, la sección más importante de la configuración. En el recuadro Object detection se configuran los parámetros que le indican al programa cómo eliminar datos at´ıpicos y como identificar las estrellas y cuerpos menores en la imágen. El parámetro Aperture Radius se refiere al tama˜ no de las estrellas: entre 10

N´ otese que estos par´ ametros difieren de los nominales. Sin embargo, la primera vez que se ingresen estos valores, deben ser los nominales. Más adelante se indicará porqué y cómo varian estos par´ ametros.

64


Figura 6.9: Configuración de la sección Program.

más alto sea el valor de este parámetro el programa buscará por regiones más grandes en la imágen que puedan ser estrellas. Las demás opciones en este recuadro se refieren a la calidad de la imágen. Debe consultarse en cualquier libro o documento sobre imágenes astronómicas para entender el significado de estos términos[10]. Inicialmente, se va a suponer que la imágen no es de muy buena calidad, para lo cual se ingresa los siguientes valores: Detection limit = 6. Minimun HWHM =0.7 y PSF Fit-RMS = 0.2. Esto se refina después de analizar la calidad final de los resultados. En la Fig. 6.9 se muestra la configuración final de esta sección. En la sección de catálogo estelar, se selecciona UCAC-4. Los l´ımites superior e inferior indican qué tan tenues y qué tan brillantes son las estrellas que el programa intenta identificar. En este punto se pone un l´ımite inferior alto de tal forma que se evite el usar estrellas muy tenues en el proceso.

´ DE DATOS 6.3. REDUCCION

65

En el recuadro de Reference star matching se pone 0 la primera vez que se usa el programa. Esto le indica al programa que en la primera ocasión se quiere hacer el proceso de manera manual. En alignment area se ingresa el n´ umero de estrellas con las que se busca realizar la alineación de imágenes. Un valor usual de este parámetro es 30.

Refinando los par´ ametros Para refinar los parámetros se siguieron los siguientes pasos. Primero, se da clic en F ile > Load images y se elige alguna de las imágenes disponibles en formato .FITS . Para hacer astrometr´ıa siempre se recomienda usar más de una imágen del mismo cuerpo, por lo cual se repite el proceso y se cargan todas11 . Posteriormente, se da clic en Astrometry > data reduction. Al hacer esto aparece un recuadro en el que se debe indicar las coordenadas del centro de la imágen o la opción de ingresar cuál es el cuerpo observado. En este caso se está usando información de 2326, 675 y 654, por lo cual se da clic en el botón de b´ usqueda y luego se da ok al cuerpo encontrado. Después de esto aparecerá sobre la imágen desplegada un conjunto de c´ırculos rojos que indican las estrellas de referencia que el programa encontró para ese campo. También aparece, como se observa en la Fig. 6.10, un recuadro con unas direcciones y las opciones magnitude, focal length y Position angle que deben manipularse hasta que los c´ırculos rojos concuerden con las estrellas de la imágen tan preciso como sea posible. En este punto es donde se hace eviente que la longitud focal y el ańgulo no corresponden a los nominales, puesto que deben ser modificados, generalmente, para que se logre un buen acuerdo entre la imágen y los c´ırculos. Cuando se considere que se obtuvo un buen acuerdo entre estrellas y c´ırculos se la clic al botón ok. el programa procederá a realizar la astrometr´ıa y mostrará un cuadro adicional con los datos reducidos, si el proceso está correcto, es decir, si los parámetros ingresados previamente son apropiados para la calidad de la imágen y si la alineación se hizo de forma adecuada. Para finalizar la configuración inicial, se observa los valores que toman estos parámetros (Local length, position angle) y se devuelve a la sección anterior de settings para ingresar los nuevos valores de estos parámetros. Sin embargo, si el procedimiento de fijar los parámetros del programa no se realiza de manera cuidadosa, o si la calidad de la imágen es baja, o si no se hizo una buena alineación, después de presionar el botón OK el programa hará una de dos 11

Obsérvese que al cargar cada imágen aparece un recuadro que indica las coordenadas y tiempo de la im´ agen. Estos datos provienen del archivo .FITS guardado por el software de la c´ amara y por tal raz´ on se considera como el correcto.

66


Figura 6.10: Ajuste manual de los parámetros instrumentales.

cosas: no avanzará de este punto puesto que se espera que el usuario haga una mejor alineación o aparecerá un error diciendo que hubo alguna divergencia en el proceso de fijar las constantes. En ambos casos se sugiere analizar detalladamente los parámetros de calidad del programa, estudiar si la imágen tiene la calidad suficiente para estos, y luego modificar la configuración de tal forma que se pueda usar los datos disponibles. Después de que este procedimiento se realiza correctamente, no es necesario repetirlo para cada imágen siempre y cuando la configuración instrumental permanezca igual: misma cámara, mismo telescopio, etc.

6.4

Resultados observacionales

Al ejecutar la reducción de datos para cada imágen se obtuvieron los resultados que se muestran a continuación. En la Tabla 6.1 se muestran los resultados para 675 y en la Tabla 6.2 para 654. Los datos se muestran en el formato del Minor Planet Center (MCP).

6.5. RESULTADOS DEL ALGORITMO

67

Tabla 6.1: Datos astrométricos obtenidos para 675. Cuerpo 675 675 675 675 675 675 675 675

t(UT) 2014 2014 2014 2014 2014 2014 2014 2014

09 09 09 09 09 09 09 09

16.16787 16.16840 16.16895 16.16948 16.17057 16.20616 16.20778 16.20831

α 22 22 22 22 22 22 22 22

41 41 41 41 41 41 41 41

02.44 02.42 02.38 02.36 02.30 00.50 00.40 00.38

δ +09 +09 +09 +09 +09 +09 +09 +09

16 16 16 16 16 16 16 16

41.4 41.2 40.9 40.7 40.4 29.0 28.2 28.2

Información obtenida para 675.

Tabla 6.2: Datos astrométricos obtenidos para 654. Cuerpo 654 654 654 654 654 654 654 654 654 654

t(UT) 2014 2014 2014 2014 2014 2014 2014 2014 2014 2014

09 09 09 09 09 09 09 09 09 09

16.22959 16.23013 16.23067 16.23120 16.23175 16.23228 16.23282 16.23337 16.23444 16.25144

α 21 21 21 21 21 21 21 21 21 21

23 23 23 23 23 23 23 23 23 23

27.20 27.26 27.21 27.28 27.18 27.28 27.15 27.12 27.02 26.36

δ +09 +09 +09 +09 +09 +09 +09 +09 +09 +09

08 08 08 08 08 08 08 08 08 08

31.2 31.4 30.9 30.1 32.5 31.7 30.1 30.6 29.5 23.0

Información obtenida para 654.

6.5

Resultados del algoritmo

Asteroide 675 Es importante notar que de estos datos puede extraerse un atribuible para cada cuerpo. Para obtener un segundo atribuible se utiliza el motor de b´ usqueda del MPC. En la Tabla 6.3 se muestran los datos usados del MPC para obtener un segundo atribuible para 675.


68

Tabla 6.3: Datos astrométricos adicionales para 675. Cuerpo 675 675 675 675

t(UT) 2014 2014 2014 2014

10 10 10 10

13.20413 13.21267 13.22344 13.23197

α 22 22 22 22

25 25 25 25

58.44 58.31 58.14 58.01

δ +06 +06 +06 +06

26 26 26 26

46.8 43.6 39.7 36.7

Información adicional para 675. Estos fueron datos fueron obtenidos por el observatorio MPC 703. El código se ejecutó para estos datos y en la Tabla 6.4 se muestran los resultados obtenidos. En la Fig. 6.11 se muestra la trayectoria encontrada para este cuerpo.

Tabla 6.4: Parámetros orbitales de 675. Par´ ametro

Obtenido

MPC

a(U A) e i(◦ ) Ω(◦ ) ω(◦ )

2.7976033 0.1994317 10.10851 263.74402 148.68496

2.7704278 0.2007596 9.78383 263.26851 152.10953

Comparación entre los datos obtenidos y los entregados por el MPC para 675.


69

Figura 6.11: Se muestra la trayectoria del asteroide 675 en el sistema ecl´ıptico heliocéntrico, en el cual el Sol está en el origen de coordenadas.

Asteroide 654 En la tabla 6.5 se muestran las datos adicionales que se utilizaron para obtener un segundo atribuible para 654. En la tabla 6.6 se muestran los parámetros orbitales obtenidos y la comparación con los entregados por el MPC y en la Fig 6.12 se observa la órbita obtenida. Tabla 6.5: Datos astrométricos adicionales para 654. Cuerpo 654 654 654 654 654 654 654 654 654

t(UT) 2014 2014 2014 2014 2014 2014 2014 2014 2014

08 08 08 08 08 08 08 08 08

08.81354 08.83865 08.85953 09.83154 09.90981 09.95659 10.81447 10.83525 10.85918

α 22 22 22 21 21 21 21 21 21

00 00 00 59 59 59 58 58 58

36.78 35.28 34.03 36.74 32.01 29.16 38.16 36.89 35.43

δ +10 +10 +10 +10 +10 +10 +10 +10 +10

48 48 48 49 49 49 50 50 50

15.8 17.8 19.6 31.9 37.3 40.5 31.9 33.1 34.4

Información adicional para 654. Estos fueron datos fueron obtenidos por el observatorio MPC l33.


70

Tabla 6.6: Parámetros orbitales de 654. Par´ ametro

Obtenido

MPC

a(U A) e i(◦ ) Ω(◦ ) ω(◦ )

2.3695859 0.2453443 14.82003 270.62349 215.84125

2.2967431 0.2313217 18.12709 278.47430 214.02028

Comparación entre los datos obtenidos y los entregados por el MPC para 654. Como en el cap´ıtulo anterior, considerando los parámetros nominales tomados del

Figura 6.12: Se muestra la trayectoria del asteroide 654 en el sistema ecl´ıptico heliocéntrico, en el cual el Sol está en el origen de coordenadas. MPC como valores verdaderos, se puede determinar los errores de los parámetros orbitales calculados. En la Tabla 6.7 se muestran los errores obtenidos. En este caso, los errores son similares a los obtenidos en el cap´ıtulo anterior, en especial, es importante notar que los errores en la forma son del orden de 10−3 AU , que son valores aceptables como órbitas iniciales. Esto muestra que los datos medidos desde el OAUTP cuentan con la calidad necesaria para obtener o´rbitas iniciales. También se nota que el método implementado en este trabajo es adecuado para la determinación de órbitas iniciales.


71

Tabla 6.7: Errores para los parámetros orbitales obtenidos Asteroide 675 654

Errores d(AU ) 0.03857 0.05192

Φ 0.096119 0.121481

Errores obtenidos para los dos cuerpos de estudio utilizados.

72


Cap´ıtulo 7 Conclusiones y trabajo futuro En este trabajo se estudia el problema de inversión de o´rbitas. Se describe detalladamente un método basado en las integrales de movimiento para el problema de los dos cuerpos, y se discute la forma en que se implementó computacionalmente, en una versión geocéntrica y en una topocéntrica. El código se escribe en una serie de subrutinas auxiliares y una rutina principal, que recibe dos conjuntos de datos observacionales y las coordenadas de los lugares de observación. Se realizó una prueba inicial del código con datos tomados del motor de b´ usqueda del MPC. El propósito de esta prueba fue analizar la correctitud del código y además determinar la versión más adecuada para los datos tomados por el Observatorio Astronómica de la Universidad Tecnológica de Pereira (OAUTP). Se encontró que el algoritmo implementado es adecuado para obtener o´rbitas iniciales y se encontró que la versión más precisa es la topocéntrica, lo cual es cierto para asteroides de la familia MBA y NEO. Por lo tanto, la versión final del método usada fue la topocéntrica. Posteriormente se describe la metodolog´ıa que es necesario implementar para obtener información astrométrica precisa con los instrumentos con los que cuenta el OAUTP. Se realizaron mediciones de los asteroides 675 y 654 y se aplica la reducción de datos para probar el algoritmo implementado. Los parámetros orbitales obtenidos son cercanos a los valores nominales entregados por el MPC. Esto se comprueba mediante el análisis de errores realizado, que muestra que los errores están dentro de los rangos esperables para o´rbitas iniciales. Es importante realizar pruebas de este algoritmo con una muestra mayor de cuerpos, para lo cual la estrategia más adecuada podr´ıa ser llevar a cabo una simulación de observaciones (mediante el método Monte Carlo, por ejemplo), en la que se pueda estudiar más detalladamente el comportamiento del método cuando se utiliza un conjunto grande de datos. También es importante ahondar en el 73

74

CAPÍTULO 7. CONCLUSIONES Y TRABAJO FUTURO

comportamiento del algoritmo cuando se var´ıa el tiempo transcurrido entre ambos arcos de observaciones. Por u ´ltimo, es importante realizar comparaciones de este método con respecto a otros para establecer ventajas y desventajas, de tal manera que el proceso de inversión de o´rbitas sea más efectivo para los datos tomados desde el OAUTP. Finalmente, en este trabajo se ha avanzado en el a´rea de Astrometr´ıa y Astronom´ıa de posición, puesto que se han podido realizar mediciones de asteroides con la calidad suficiente para obtener en un futuro cercano el código de observatorio, que son otorgados por el MPC. Asimismo, la teor´ıa aqu´ı estudiada permitirá al grupo de investigación en Astroingenier´ıa Alfa Orión abordar otros problemas modernos relacionados con la inversión de o´rbitas, en los que el grupo se encuentra en la capacidad de hacer aportes propios y de realizar colaboraciones con otros grupos y observatorios.

Anexos

75

Ap´ endice A C´ odigo Fuente de funciones b´ asicas A.1

Convertir UT a JD

function res= UTtoJD(X) %X debe tener el formato del %minor planet center en UT: Yyyy mm dd.sss A=X(:,1); M=X(:,2); d=X(:,3); m=length(X(:,1)); for i=1:m if M(i)<=2 A(i)=A(i)-1; M(i)=M(i)+12; end res(i)=floor(365.25*(A(i)+4716))+floor(30.6001*(M(i)+1))-floor(A(i)/100); res(i)=res(i) +floor(floor(A(i)/100)/4)+d(i)-1522.5; end end

A.2

Convertir JD a UT

function res= JDtoUT(t) %T=(t-2415019.5)/(365.25); T=(t-2451545.5)/(365.25); Yr=2000+floor(T); LeapYear=floor((Yr-2000-1)*0.25); days=(t-2451544.5)-((Yr-2000)*365+LeapYear); if days<1 Yr=Yr-1; end

77

78

´ ´ ´ APENDICE A. CODIGO FUENTE DE FUNCIONES BASICAS

LeapYear=floor((Yr-2000-1)*0.25); days=(t-2451544.5)-((Yr-2000)*365+LeapYear); d=floor(days); aux=[31 28 31 30 31 30 31 31 30 31 30 31]; if mod(Yr,4)==0 aux(2)=29; end a=0; i=1; while(a
A.3

Transforma la ascensi´ on recta a radianes

function res=RA(A) %Convierte la ascensi n rescta del format MPC %a radianes format long res=zeros(length(A(:,1)),1); for i=1:length(A(:,1)) res(i,1)=15*(A(i,1)+A(i,2)/60+A(i,3)/3600)*pi/180; if res(i,1)>=pi res(i,:)= -1*(2*pi-res(i,:)); end end end

A.4

Transformar la declinaci´ on a radianes

function res=DEC(A)%Convierte la declinaci n del format MPC a radianes format longg res=zeros(length(A(:,1)),1); for i=1:length(A(:,1)) if A(i,1)~=0 res(i,1)=sign(A(i,1))*(abs(A(i,1))+A(i,2)/60+A(i,3)/3600)*pi/180; else res(i,1)=(abs(A(i,1))+A(i,2)/60+A(i,3)/3600)*pi/180; end end end

´ DE KEPLER A.5. RESUELVE LA ECUACION

A.5

79

Resuelve la ecuaci´ on de Kepler

En la Fig. A.1 se muestra el diagrama de flujo del método numérico utilizado para resolver la ecuación de Kepler y posteriormente se muestra el código implementado.

Se introducen f, f 0 , f 00 , f 00 , el valor del error , δ, la excéntricidad e y M .

se introduce la estimación inicial E = E0 = M + 0.85e ∗ sgn (sin M ) n=1 δ1

=

− ff0(E) (E)

δ2 = − f 0 (E)+f (E) 1 δ1 f 00 (E) 2

∆x =

f (E) f 0 (E)+ 12 δ2 f 00 (E)+ 61 δ22 f 000 (E)

E0

= E + ∆x

k∆xk ≤ δ o kf (E)k ≤

s´ı

Ouput=E 0

E = E0 n = n+1

no

n ≥ 30

s´ı

No converge

Figura A.1: Diagrama de flujo de la solución a la ecuación de Kepler.

80


function res= E(e,M) %Resuelve la ecuaci n de Kepler format longg function res= f(x) res=x-e*sin(x)-M; end function res= f1(x) res=1-e*cos(x); end function res= f2(x) res=e*sin(x); end function res=f3(x) res=e*cos(x); end a=0; M=M-floor(M/2*pi)*2*pi; sigma=sign(sin(M)); x=M+0.85*e*sigma; error=1e-15; k=1; while (a==0) d1=(-1*f(x))/(f1(x)); d2=(-1*f(x))/(f1(x)+0.5*d1*f2(x)); d3=(-1*f(x))/(f1(x)+0.5*d2*f2(x)+0.1666666666666667*d2*d2*f3(x)); aux=x+d3; if (abs(x-aux)=30 error('myApp:argChk', 'eL M todo no converge') end end end

A.6. CALCULA EL VECTOR DE ESTADO

A.6

81

Calcula el vector de estado

En la Fig. A.2 se muestra el diagrama del flujo del código implementado en esta sección. Entrada: parámetros orbitales (a, e, i, Ω, ω, t0 ), la constante gravitacional µ y el tiempo t.

t = t0

s´ı

M = M (t0 )

no ∆t = 86400 (t − t0 ) p M = M0 + ∆t aµ3

Se resuelve la ecuación de Kepler: Ec = E (e, M ) Se obtiene la anomal´ıa verdadera: √ 1+e sin Ec θ = 2 arctan √1−e cos E2c y 2

r = a (1 − e cos Ec )

r = r (cos θ, sin θ, 0) √ µa − sin Ec , 1 − e2 cos Ec , 0 r

√

r˙ =

{r, r˙ } = Rz (−Ω) Rx (−i) Rz (−ω) {r, r˙ } Figura A.2: Diagrama de flujo del algoritmo que obtiene el vector de estado a partir de los parámetros orbitales.

82


function [pos vel]=phasevector(par,t)%Obtiene el vector de %estado a partir de los par metros orbitales para el tiempo t. %Las unidades de entrada deben ser [AU] y [ ] y los %de salida [AU] y [AU/dia] a=par(1)*(149597870700); e=par(2); ic=cspice convrt(par(3),'degrees','radians'); Omega=cspice convrt(par(4),'degrees','radians'); omega=cspice convrt(par(5),'degrees','radians'); t0=par(6); dt=86400*(t-t0); mu=1.32712440041e20; M=dt*sqrt((mu)/(aˆ3)); M=wrapTo2Pi(M); Et=E(e,M); Et=wrapTo2Pi(Et); theta=2*atan2(sqrt(1+e)*sin((Et/2)),sqrt(1-e)*cos((Et/2))); r=a*(1-e*cos(Et)); pos=[r*cos(theta) r*sin(theta) 0]; vel=[-1*sin(Et) sqrt(1-eˆ2)*cos(Et) 0]; vel=((sqrt(mu*a))/(r))*vel; Omega=-1*Omega; omega=-1*omega; ic=-1*ic; pos=[cos(Omega) sin(Omega) 0; -1*sin(Omega) cos(Omega) 0;0 0 1]*..., [1 0 0; 0 cos(ic) sin(ic);0 -1*sin(ic) cos(ic)]*..., [cos(omega) sin(omega) 0; -1*sin(omega) cos(omega) 0;0 0 1]*pos'; pos=(1/(1.49597870691e11))*pos'; vel=[cos(Omega) sin(Omega) 0; -1*sin(Omega) cos(Omega) 0;0 0 1]*..., [1 0 0; 0 cos(ic) sin(ic);0 -1*sin(ic) cos(ic)]*..., [cos(omega) sin(omega) 0; -1*sin(omega) cos(omega) 0;0 0 1]*vel'; vel=(86400/(1.49597870691e11))*vel'; end

A.7

Calcula el vector de estado de la Tierra

function [Tp Tv]=earthphase(t) a=1.00000261; mu=1.32712440041e20; adot=0.00000562; e=0.01671123; edot=-0.00004392; I=-0.00001531; Idot=-0.01294668; L=100.46457166; Ldot=35999.37244981; Lperi=102.93768193; Lperidot=0.32327364;

´ A.8. PARAMETROS ORBITALES A PARTIR DEL VECTOR DE ESTADO83 Omega=0; T=(t-2451545)/(36525); a=a+adot*T; a=a*(149597870700); e=e+edot*T; I=I+Idot*T; L=L+Ldot*T; Lperi=Lperi+Lperidot*T; omega=Lperi; M=L-Lperi; I=cspice convrt(I,'degrees','radians'); omega=cspice convrt(omega,'degrees','radians'); M=cspice convrt(M,'degrees','radians'); M=wrapTo2Pi(M); Et=E(e,M); Et=wrapTo2Pi(Et); theta=2*atan2(sqrt(1+e)*sin((Et/2)),..., sqrt(1-e)*cos((Et/2))); r=a*(1-e*cos(Et)); Tp=[r*cos(theta) r*sin(theta) 0]; Tv=[-1*sin(Et) sqrt(1-eˆ2)*cos(Et) 0]; Tv=((sqrt(mu*a))/(r))*Tv; Omega=-1*Omega; omega=-1*omega; ic=-1*I; Tp=[cos(Omega) sin(Omega) 0; ..., -1*sin(Omega) cos(Omega) 0;0 0 1]*..., [1 0 0; 0 cos(ic) sin(ic);..., 0 -1*sin(ic) cos(ic)]*..., [cos(omega) sin(omega) 0; ..., -1*sin(omega) cos(omega) 0;0 0 1]*Tp'; Tp=(1/(1.49597870691e11))*Tp'; Tv=[cos(Omega) sin(Omega) 0;..., -1*sin(Omega) cos(Omega) 0;0 0 1]*..., [1 0 0; 0 cos(ic) sin(ic);0 -1*sin(ic) cos(ic)]*..., [cos(omega) sin(omega) 0;..., -1*sin(omega) cos(omega) 0;0 0 1]*Tv'; Tv=(86400/(1.49597870691e11))*Tv'; end

A.8

C´ alcula los par´ ametros orbitales a partir del vector de estado

En la Fig. A.3 se muestra el diagrama de flujo del algoritmo que obtiene los parámetros orbitales si como se entrada se tiene el vector de estado y la constante gravitacional.

84

´ ´ ´ APENDICE A. CODIGO FUENTE DE FUNCIONES BASICAS Entrada: parámetros orbitales (a, e,Entrada: i, Ω, ω, t0 ), (r,lar˙ )constante y µ. gravitacional µ y el tiempo t.

h = r × r˙ , ˙ r×h µ

e=

−

r krk

hz i = arccos khk

θ = 2π − e·r arccos kekkrk

no

e

r · r˙ ≥ 0

=

s´ı

e·r θ = arccos kekkrk

kek, tan

θ

E = 2 arctan q 1+e2 , 1−e

n = (0, 0, 1)T × h

Ω = 2π − nx arccos knk

no

ω = 2π − n·e arccos knkkek

no

ny ≥ 0

ez ≥ 0

s´ı

nx Ω = arccos knk

s´ı

n·e ω = arccos knkkek

M = E − e sin E, a

=

1 ˙ 2 krk 2 − µ krk

Figura A.3: Diagrama de flujo del algoritmo que obtiene los parámetros orbitales a partir del vector de estado.

A.9. OBTIENE EL VECTOR ECLÍPTICO ECUATORIAL DE LA TIERRA.85 function res=stateTopar(r,rdot,a,t) k=0.01720209895; ext=dms2degrees([23 27 08.26]); ext=cspice convrt(ext,'degrees','radians'); Tr=[1 0 0; 0 cos(ext) sin(ext);0 -1*sin(ext) cos(ext)]; r=Tr*r'; r=r'; rdot=Tr*rdot'; rdot=rdot'; h=cross(r,rdot); nh=cross([0 0 1],h); Omega =acos(nh(1)/norm(nh)); if nh(2)<0 Omega = 2*pi-Omega; end I=atan2((sqrt(h(1)ˆ2+h(2)ˆ2)),h(3)); e=(1/kˆ2)*cross(rdot,h)-r/(norm(r)); aux1=e(1)*cos(Omega)+e(2)*sin(Omega); aux2=(e(3))/(sin(I)); w=atan2(aux2,aux1); if (1-(norm(r)/a))/(norm(e)) <1 E=wrapTo2Pi(acos((1-(norm(r)/a))/(norm(e)))); T=t-(E-norm(e)*sin(E))*sqrt((aˆ3)/(kˆ2)); medio=k*aˆ(-1*3/2); M=medio*(t-T); tmjd=t-2400000.5; res=[a norm(e) I*(180/pi) Omega*(180/pi)]; res=[res, wrapTo360(w*(180/pi)) wrapTo360(M*(180/pi)) tmjd]; else res=zeros(1,7); end end

A.9

Obtiene el vector ecl´ıptico ecuatorial de la Tierra.

function res= tierra(t)%Entrega las coordenadas cartesianas de la tierra en % el sistema ecl\'iptico ecuatorial d=t - 2451543.5; w=282.9404+(4.70935e-5)*d; w=w*(pi/180); w=wrapTo2Pi(w); a=1.00000261; e=0.016709-(1.151e-9)*d; M=356.0470+(0.9856002585)*d; ext=23.4393- (3.563e-7)*d; ext=ext*(pi/180);

86


M=wrapTo2Pi((M*(pi/180))); anoex=wrapTo2Pi(E(e,M)); x=a*(cos(anoex)-e); y=a*sqrt(1-e*e)*sin(anoex); r = sqrt(x*x + y*y); v=wrapTo2Pi(atan2(y,x)); l=wrapTo2Pi(v+w); x=r*cos(l); y=r*sin(l); z=0; Tr1=[1 0 0;0 cos(ext) -1*sin(ext); 0 sin(ext) cos(ext)]; res=Tr1*[x y z]'; res=-1*res'; end

A.10

Obtiene el vector atribu´ıble

Utiliza A.9 para obtener el vector tierra en el instante t. function [atri Tpos Tvel t]=atribuible(X,a) format longg ext=dms2degrees([23 27 08.26]); ext=cspice convrt(ext,'degrees','radians'); Tr1=[1 0 0; 0 cos(ext) -1*sin(ext);0 sin(ext) cos(ext)]; t=X(:,1); tm=mean(t); for i=1:length(t) Tierra(i,:)=vpa(sym(tierra(t(i))), 128); geoc(i,:)=vpa(sym(topovector(t(i),a)), 128); end ar=vpa(sym(X(:,2)), 128); de=vpa(sym(X(:,3)), 128); m=length(X(:,1)); if m<=2 error('myApp:argChk', 'el n mero de observaciones m debe ser mayor o igual a 3. S\'i m=2, int ntese un modelo lineal') end %C lculo de la matriz de dise o (design matrix) B=zeros(m,3); %El modelo es una funci n cuadr tica for i=1:m B(i,:)=[-1 tm-t(i) -0.5*(t(i)-tm)*(t(i)-tm)]; end C=B'*B; Ar=-1*C\B'*ar; De=-1*C\B'*de; Tx=-1*C\B'*Tierra(:,1); Ty=-1*C\B'*Tierra(:,2); Tz=-1*C\B'*Tierra(:,3); gx=-1*C\B'*geoc(:,1); gy=-1*C\B'*geoc(:,2);

A.11. DETERMINA ERRORES ORBITALES gz=-1*C\B'*geoc(:,3); atri=[Ar(1) De(1) Ar(2) De(2)]; a=tierra(tm); Tpos=[a(1)+gx(1) a(2)+gy(1) a(3)+gz(1)]; Tvel=[Tx(2)+gx(2) Ty(2)+gy(2) Tz(2)+gz(2)]; t=tm; end

A.11

Determina errores orbitales

function [d phi]=orbitalerrors(X) %X es una matriz 2x7 %donde la primera fila son los par metros %verdaderos y la segunda los estimados k=0.01720209895; ext=dms2degrees([23 27 08.26]); ext=cspice convrt(ext,'degrees','radians'); par=X(1,1:6); t=X(1,7); [rv rdotv]=phasevector(par,t); par=X(2,1:6); t=X(2,7); [ra rdota]=phasevector(par,t); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% b1=X(1,1)*sqrt(1-(X(1,2))ˆ2); b2=X(2,1)*sqrt(1-(X(2,2))ˆ2); d=sqrt((X(1,1)-X(2,1))ˆ2+(b1-b2)ˆ2); h=cross(rv,rdotv); h=(1/norm(h))*h; runi=(1/norm(rv))*(rv); auxs=cross(h,runi); C=[runi; auxs; h]; %%%%%%%%%%%%% h=cross(ra,rdota); h=(1/norm(h))*h; runi=(1/norm(ra))*(ra); auxs=cross(h,runi); Cs=[runi; auxs; h]; phi=0.5*(trace(C*transpose(Cs))-1); phi=acos(phi); end

A.12

Grafica la trayectoria

function res= grafica(par,t) res=0;

87

88


k=0.01720209895; e=par(2); a=cspice convrt(par(1),'AU','KM'); inc=cspice convrt(par(3),'degrees','radians'); lnode=cspice convrt(par(4),'degrees','radians'); argp=cspice convrt(par(5),'degrees','radians'); T=par(6); par=[a e inc lnode argp T]; ts=t; for i=1:2000 t=ts+i; aux=grupoao par2vecsta(par,t); r(i,1:3)=cspice convrt(aux(1,1:3),'KM','AU'); Tierra(i,:)=vector tierra(t); end set(figure,'Color',[1 1 1]); subplot(2,2,1) h=plot3(r(:,1),r(:,2),r(:,3),Tierra(:,1),Tierra(:,2),Tierra(:,3)); grid on jk=title('Trayectoria del asteroide 654'); set(jk, 'FontSize', 12) set(jk,'FontWeight','bold') legend('654','Tierra',1); set(h,'linewidth',1.5) xh=xlabel('x(AU)') ; set(xh, 'FontSize', 12) set(xh,'FontWeight','bold') yh=ylabel('y(AU)') ; set(yh, 'FontSize', 12) set(yh,'FontWeight','bold') zh=zlabel('z(AU)') ; set(zh, 'FontSize', 12) set(zh,'FontWeight','bold') subplot(2,2,2) h=plot(r(:,1),r(:,2),Tierra(:,1),Tierra(:,2)); grid on jk=title('Proyecci n sobre el plano x-y'); set(jk, 'FontSize', 12) set(jk,'FontWeight','bold') legend('654','Tierra',1); set(h,'linewidth',1.5) xh=xlabel('x(AU)') ; set(xh, 'FontSize', 12) set(xh,'FontWeight','bold') yh=ylabel('y(AU)') ; set(yh, 'FontSize', 12) set(yh,'FontWeight','bold') subplot(2,2,3) h=plot(r(:,2),r(:,3),Tierra(:,2),Tierra(:,3)); grid on

A.12. GRAFICA LA TRAYECTORIA jk=title('Proyecci n sobre el plano y-z'); set(jk, 'FontSize', 12) set(jk,'FontWeight','bold') legend('654','Tierra',1); set(h,'linewidth',1.5) xh=xlabel('y(AU)') ; set(xh, 'FontSize', 12) set(xh,'FontWeight','bold') yh=ylabel('z(AU)') ; set(yh, 'FontSize', 12) set(yh,'FontWeight','bold') subplot(2,2,4) h=plot(r(:,1),r(:,3),Tierra(:,1),Tierra(:,3)); grid on jk=title('Proyecci n sobre el plano x-z'); set(jk, 'FontSize', 12) set(jk,'FontWeight','bold') legend('654','Tierra',1); set(h,'linewidth',1.5) xh=xlabel('x(AU)') ; set(xh, 'FontSize', 12) set(xh,'FontWeight','bold') yh=ylabel('z(AU)') ; set(yh, 'FontSize', 12) set(yh,'FontWeight','bold') end

89

90


Ap´ endice B ´ C´ odigo Fuente para Algebra de polinomios B.1

Suma de polinomios

Ejecuta la operación X + (c) Y , donde X, Y son polinomios y c ∈ R. function res= suma(X,Y,c) %Hace la operaci n X+(c)Y [N1 N2]=size(Y); aux=zeros(N1,N2); for i=1:N1 aux(i,:)=[c*Y(i,1) Y(i,2:end)]; end res=red([X;aux]); end

B.2

Producto de polinomios

function res= producto(X,Y) N1=length(X(:,1)); N2=length(Y(:,1)); k=0; for i=1:N1 for j=1:N2 aux(j+k,:)= [X(i,1)*Y(j,1) X(i,2:end)+Y(j,2:end)]; end k=i*N2; end res=red(aux); end

91

´ ´ ´ 92 APENDICE B. CODIGO FUENTE PARA ALGEBRA DE POLINOMIOS

B.3

Simplificaci´ on de polinomios

function res= red(a)%simplifica polinomios axu1=expo(a); aux=a; N=length(a(:,1)); res=zeros(1,length(a(1,:))); j=1; r=1; row=true; while (N>=1) coe=aux(1,1); for i=2:N if (aux(1,2:end)-aux(i,2:end)==0) coe=coe+aux(i,1); row=false; r=[r i]; end end if row && (coe~=0) res(j,:)=aux(1,:); j=j+1; elseif (coe~=0) res(j,:)=[coe, aux(1,2:end)]; j=j+1; end aux=removerows(aux,'ind',r); N=length(aux(:,1)); r=1; end end

B.4

Evaluaci´ on de polinomios en una variable

function res= polyevaluas(A,x) %Evalua x en el polinomio de una variable A res=0; for i=1:length(A(:,1)) aux=A(i,2); res=res+A(i,1)*xâux; end end

B.5

Evaluaci´ on de polinomios de dos variables

function res= polyevalua2d(A,x,y) %Evalua x,y en el polinomio de dos variable A

B.6. ALGORITMO FFT res=0; for i=1:length(A(:,1)) aux1=A(i,2); aux2=A(i,3); res=res+(A(i,1))*(xâux1)*(yâux2); end end

B.6

Algoritmo FFT

function res= fftpoly(X) if mod(length(X),2)~=0 X=[X 0] ; end N=length(X); m=1; n=1; for i=1:N if mod(i,2)==0 A(m)=X(i); m=m+1; else B(n)=X(i); n=n+1; end end Af=dfts(A); Bf=dfts(B); Wn=exp((-1*2*pi*1i)/N); w=1; for k=0:N/2-1 res(k+1)=Bf(k+1)+w*Af(k+1); d=k+1+N/2; res(d)=Bf(k+1)-w*Af(k+1); w=w*Wn; end end

B.7

Algoritmo IDFT

function res= invfftpoly(X) if mod(length(X),2)~=0 X=[X 0] ; end N=length(X); m=1; n=1;

93

´ ´ ´ 94 APENDICE B. CODIGO FUENTE PARA ALGEBRA DE POLINOMIOS for i=1:N if mod(i,2)==0 A(m)=X(i); m=m+1; else B(n)=X(i); n=n+1; end end Af=invdfts(A); Bf=invdfts(B); Wn=exp((2*pi*1i)/N); w=1; for k=0:N/2-1 res(k+1)=Bf(k+1)+w*Af(k+1); d=k+1+N/2; res(d)=Bf(k+1)-w*Af(k+1); w=w*Wn; end res=(1/N)*res; for i=1 :length(res) if imag(res(i))<=1e-15 res(i)=real(res(i)); end end end

B.8

C´ alculo del resultante

function res= detsylvester(A,B,x) format longg function res= polyevalua(A,x) res=0; for i=1:length(A(:,1)) aux=A(i,2); res=res+A(i,1)*xâux; end end for j=1:21 p(j)=polyevalua(A{j},x); end for j=1:3 q(j)=polyevalua(B{j},x); end aux1=[p(21) p(20) p(19) p(18) p(17) p(16) p(15) p(14) p(13) p(12)]; aux1=[aux1,[p(11) p(10) p(9) p(8) p(7) p(6) p(5) p(4) p(3) p(2) p(1) 0]; aux2=[0 p(21) p(20) p(19) p(18) p(17) p(16) p(15) p(14) p(13) p(12)]; aux2=[aux2,[p(11) p(10) p(9) p(8) p(7) p(6) p(5) p(4) p(3) p(2) p(1)]]; aux3=[q(3) q(2) q(1) zeros(1,19)];

´ B.8. CALCULO DEL RESULTANTE

95

aux4=[0 q(3) q(2) q(1) zeros(1,18)]; aux5=[0 0 q(3) q(2) q(1) zeros(1,17)]; aux6=[0 0 0 q(3) q(2) q(1) zeros(1,16)]; aux7=[0 0 0 0 q(3) q(2) q(1) zeros(1,15)]; aux8=[0 0 0 0 0 q(3) q(2) q(1) zeros(1,14)]; aux9=[0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,13)]; aux10=[0 0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,12)]; aux11=[0 0 0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,11)]; aux12=[0 0 0 0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,10)]; aux13=[0 0 0 0 0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,9)]; aux14=[0 0 0 0 0 0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,8)]; aux15=[0 0 0 0 0 0 0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,7)]; aux16=[0 0 0 0 0 0 0 0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,6)]; aux17=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,5)]; aux18=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,4)]; aux19=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,3)]; aux20=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,2)]; aux21=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 q(3) q(2) q(1) zeros(1,1)]; aux22=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 q(3) q(2) q(1)]; A=[aux1' aux2' aux3' aux4' aux5' aux6' aux7' aux8' aux9' aux10']; A=[A,aux11' aux12' aux13' aux14' aux15' aux16' aux17' aux18']; A=[A,[aux19' aux20' aux21' aux22']]; res=det(A); end

´ ´ ´ 96 APENDICE B. CODIGO FUENTE PARA ALGEBRA DE POLINOMIOS

Ap´ endice C C´ odigo Fuente m´ etodo principal

function res=metodo(X,Y,x,y)%X es el atribuible 1 y Y el 2. x e y son las %constantes de paralaje de X e Y, respectivamente. format longg digits(128) [at2 Tp2 Tv2 t2]=atribuible(X,x); [at1 Tp1 Tv1 t1]=atribuible(Y,y); at1=vpa(sym(at1), 128); at2=vpa(sym(at2), 128); Tp1=vpa(sym(Tp1), 128); Tp2=vpa(sym(Tp2), 128); Tv1=vpa(sym(Tv1), 128); Tv2=vpa(sym(Tv2), 128); k=0.01720209895; c0=Tp1(1)*Tp1(1)+Tp1(2)*Tp1(2)+Tp1(3)*Tp1(3); d0=Tp2(1)*Tp2(1)+Tp2(2)*Tp2(2)+Tp2(3)*Tp2(3); rho1=[cos(at1(1))*cos(at(2)) sin(at1(1))*cos(at1(2)) sin(at1(2)) ]; rho2=[cos(at2(1))*cos(at2(2)) sin(at2(1))*cos(at(2)) sin(at2(2))]; pa1=[-1*sin(at1(1))*cos(at1(2)) cos(at1(1))*cos(at1(2)) 0]; pa2=[-1*sin(at2(1))*cos(at2(2)) cos(at2(1))*cos(at2(2)) 0]; pd1=[-1*cos(at1(1))*sin(at1(2)) -1*sin(at1(1))*sin(at1(2)) cos(at1(2))]; pd2=[-1*cos(at(1))*sin(atr(2)) -1*sin(at2(1))*sin(at2(2)) cos(at2(2))]; c1=2*dot(Tv1,rho1); d1=2*dot(Tv2,rho2); c2=(cos(at1(2)))*(cos(at1(2)))*(at1(3))ˆ2+(at1(4))ˆ2; d2=(cos(at2(2)))*(cos(at2(2)))*(at(3))ˆ2+(at2(4))ˆ2; c3=2*at1(3)*dot(Tv1,pa1)+2*at1(4)*dot(Tv1,pd1); d3=2*at2(3)*dot(Tv2,pa2)+2*at2(4)*dot(Tv2,pd2); c4=Tv1(1)*Tv1(1)+Tv1(2)*Tv1(2)+Tv1(3)*Tv1(3); d4=Tv2(1)*Tv2(1)+Tv2(2)*Tv2(2)+Tv2(3)*Tv2(3); c5=2*dot(Tp1,rho1); d5=2*dot(Tp2,rho2); D1=cross(Tp1,rho1); D2=cross(Tp2,rho2); E1=at(3)*cross(rho1,pa1)+at1(4)*cross(rho1,pd1);

97

98

´ ´ ´ APENDICE C. CODIGO FUENTE METODO PRINCIPAL

E2=at2(3)*cross(rho2,pa2)+at2(4)*cross(rho2,pd2); F1=at1(3)*cross(Tp1,pa1)+at1(4)*cross(Tp1,pd1)+cross(rho1,Tv1); F2=at2(3)*cross(Tp2,pa2)+at2(4)*cross(Tp2,pd2)+cross(rho2,Tv2); G1=cross(Tp1,Tv1); G2=cross(Tp2,Tv2); POL1=suma([dot(G2-G1,cross(D1,D2)) 0 0],suma([dot(F2,cross(D1,D2)) 0 1],... suma([-1*dot(F1,cross(D1,D2)) 1 0],suma([-1*dot(E1,cross(D1,D2)) 2 0],... [dot(E2,cross(D1,D2)) 0 2],1),1),1),1); POL1=red(POL1); beta=cross(D1,D2); nbeta=beta(1)*beta(1)+beta(2)*beta(2)+beta(3)*beta(3); l1=(dot(beta,cross(E2,D2)))/nbeta; l2=(-1*dot(beta,cross(E1,D2)))/nbeta; l3=(dot(beta,cross(F2,D2)))/nbeta; l4=(-1*dot(beta,cross(F1,D2)))/nbeta; l5=(dot(beta,cross(G2-G1,D2)))/nbeta; h1=(dot(beta,cross(E2,D1)))/nbeta; h2=(-1*dot(beta,cross(E1,D1)))/nbeta; h3=(dot(beta,cross(F2,D1)))/nbeta; h4=(-1*dot(beta,cross(F1,D1)))/nbeta; h5=(dot(beta,cross(G2-G1,D1)))/nbeta; xdot=suma([l5 0 0],suma([l4 1 0],suma([l3 0 1],... suma([l1 0 2],[l2 2 0],1),1),1),1); ydot=suma([h5 0 0],suma([h4 1 0],suma([h3 0 1],... suma([h1 0 2],[h2 2 0],1),1),1),1); xdot2=producto(xdot,xdot); ydot2=producto(ydot,ydot); FH1=suma(suma(suma(suma(xdot2,xdot,c1),[1 2 0],c2),... [1 1 0],c3),[1 0 0],c4); FH2=suma(suma(suma(suma(ydot2,ydot,d1),[1 0 2],d2),... [1 0 1],d3),[1 0 0],d4); GH1=suma(suma([1 2 0],[1 1 0],c5),[1 0 0],c0); GH2=suma(suma([1 0 2],[1 0 1],d5),[1 0 0],d0); EQ=producto(suma(FH1,FH2,-1),suma(FH1,FH2,-1)); EQ=producto(producto(EQ,GH1),GH2); EQ=suma(EQ,suma(GH1,GH2,1),-4*kˆ4); EQ=producto(EQ,EQ); EQ=suma(EQ,producto(GH1,GH2),-64*kˆ8); EQ=red(EQ); p=0; aux1=zeros(1,2); a=cell(21,1); for i=1:21 for j=1:length(EQ(:,1)) if EQ(j,2)==(i-1) aux1=[aux1 ; [EQ(j,1) EQ(j,3)]]; p=1; end end a{i}=red(aux1); aux1=zeros(1,2); end

99 aux2=zeros(1,2); b=cell(3,1); for i=1:3 for j=1:length(POL1(:,1)) if POL1(j,2)==(i-1) aux2=[aux2; [POL1(j,1) POL1(j,3)]]; p=1; end end b{i}=red(aux2); aux2=zeros(1,2); end N=64; aux3=zeros(N,1); for j=0:N-1 auxs=exp((2*pi*1i*j)/N); aux3(j+1)= detsylvester(a,b,auxs); end T=invfftpoly(aux3); T=(1e94)*T(1:49); for i=1:49 r(i)=T(end-(i-1)); end Rai=roots(r'); m=1; sol=0; for i=1:length(Rai(:,1)) aux=imag(Rai(i)); aux2=real(Rai(i)); if (aux==0)&&(aux2>0) sol=[sol; real(Rai(i))]; end end sol=sol(2:end); solpar=zeros(1,2); for i=1:length(sol(:,1)) aux=roots([polyevaluas(b{3},sol(i)) polyevaluas(b{2},sol(i))... polyevaluas(b{1},sol(i))]); auxa=polyevalua2d(EQ,aux(1),sol(i)); auxb=polyevalua2d(EQ,aux(2),sol(i)); if (abs(auxa)>abs(auxb))&&(aux(2)>0)&&(imag(aux(2))==0) solpar=[solpar; [aux(2) sol(i)]]; elseif (aux(1)>0)&&(imag(aux(2))==0) solpar=[solpar; [aux(1) sol(i)]]; end end solpar=solpar(2:end,:); l=1; for i=1:length(solpar(:,1)) spu(i,1)=polyevalua2d(FH1,solpar(i,1),solpar(i,2))-... polyevalua2d(FH2,solpar(i,1),solpar(i,2)); spu(i,1)=spu(i,1)-(2*k*k)/(sqrt(polyevaluas(GH1,solpar(i,1))))...

100


+2*k*k/(sqrt(polyevaluas(GH2,solpar(i,2)))); POLA=suma(producto(producto(suma(FH1,FH2,-1),suma(FH1,FH2,-1)),... producto(GH1,GH2)),suma(GH1,GH2,1),-4*kˆ4); POLB=producto(GH1,GH2); spu(i,2)=polyevalua2d(POLA,solpar(i,1),solpar(i,2))+... (8*kˆ4)*sqrt(polyevalua2d(POLB,solpar(i,1),solpar(i,2))); spuf(i)=abs(spu(i,1))+abs(spu(i,2)); end for i=1:length(spuf) if (abs(spuf(i))==min(abs(spuf)))||(abs(spuf(i))<=1e-6) solf(l,:)=solpar(i,:); if (abs(spuf(i))==min(abs(spuf))) rsl=l; end l=l+1; end end solf=solpar; a=zeros(1,length(solf(:,1))); ext=dms2degrees([23 27 08.26]); ext=cspice convrt(ext,'degrees','radians'); s=1; for i=1:length(solf(:,1)) aux1=solf(i,1); aux2=polyevalua2d(xdot,solf(i,1),solf(i,2)); aux3=solf(i,2); aux4=polyevalua2d(ydot,solf(i,1),solf(i,2)); orb1(i,:)=[double(at1(1)), double(at1(2)),... double(at1(3)),double(at1(4)), aux1, aux2]; orb2(i,:)=[double(at2(1)), double(at2(2)),... double(at2(3)),double(at2(4)), aux3, aux4]; rs1=double(Tp1+rho1*orb1(i,5)); rdot1=double(Tv1+rho1*orb1(i,6)+orb1(i,5)*... orb1(i,3)*pa1+orb1(i,5)*orb1(i,4)*pd1); rs2=double(Tp2+rho2*orb2(i,5)); rdot2=double(Tv2+rho2*orb2(i,6)+orb2(i,5)*... orb2(i,3)*pa2+orb2(i,5)*orb2(i,4)*pd2); p=orb1(i,5); r=orb1(i,6); energia1= double((rˆ2+c1*r+c2*pˆ2+c3*p+c4-(2*kˆ2)/(sqrt(pˆ2+c5*p+c0)))/2); a(i)=-(kˆ2)/(2*energia1); sp=(299792458*60*60*24)/(149597870700); t1=t1-(orb1(i,5)/sp); t2=t2-(orb2(i,5)/sp); if a(i)>0 parametros1(i,:)=stateTopar(rs1,rdot1,a(i),t1); parametros2(i,:)=stateTopar(rs2,rdot2,a(i),t2); penalty(i,1)=parametros1(i,5)-parametros2(i,5); penalty(i,2)=parametros2(i,6)+k*a(i)ˆ(-3/2)*(t1-t2)-parametros2(i,6); identinorm(i)=sqrt(penalty(i,1)ˆ2+penalty(i,1)ˆ2); if identinorm(i)<1e-4 orbitafinal(s,:)=parametros1(i,:);

101 s=s+1; end end end res=orbitafinal(rsl,:); end

102


Bibliograf´ıa [1] Jorge I Zuluaga and Ignacio Ferrin. A preliminary reconstruction of the orbit of the chelyabinsk meteoroid. arXiv preprint arXiv:1302.5377, 2013. [2] Linda Dimare. Problems of Celestial Mechanics. PhD thesis, Sapienza, Universitá Di Roma, Roma, 2010. [3] Donald Teets and Karen Whitehead. The discovery of ceres: How gauss became famous. Mathematics magazine, pages 83–93, 1999. [4] C.F. Gauss. Theoria motus corporum coelestium in sectionibus conicis solem ambientium. 1809. [5] H. Curtis. Orbital Mechanics for Engineering Students. Aerospace Engineering. Elsevier Science, 2013. [6] J.M.A. Danby. Fundamentals of Celestial Mechanics. Willmann-Bell, 1988. [7] Andrea Milani, Giovanni F Gronchi, Davide Farnocchia, Zoran Kneˇzević, Robert Jedicke, Larry Denneau, and Francesco Pierfederici. Topocentric orbit determination: algorithms for the next generation surveys. Icarus, 195(1):474–492, 2008. [8] A. Milani and G. Gronchi. Theory of Orbit Determination. Theory of Orbit Determination. Cambridge University Press, 2010. [9] BG Marsden. Initial orbit determination-the pragmatist’s point of view. The Astronomical Journal, 90:1541–1547, 1985. [10] J. Kovalevsky and P.K. Seidelmann. Fundamentals of Astrometry. Cambridge University Press, 2004. [11] Andrea Milani and Zoran Kneˇzević. From astrometry to celestial mechanics: orbit determination with very short arcs. Celestial Mechanics and Dynamical Astronomy, 92(1-3):1–18, 2005. [12] Andrea Milani, Giovanni F Gronchi, and Zoran Kneˇzević. New definition of discovery for solar system objects. Earth, Moon, and Planets, 100(1-2):83– 116, 2007. 103

104

BIBLIOGRAFÍA

[13] Taghi Mirtorabi. A simple procedure to extend the gauss method of determining orbital parameters from three to n points. Astrophysics and Space Science, 349(1):137–141, 2014. [14] K Muinonen. Orbital covariance eigenproblem for asteroids and comets. Monthly Notices of the Royal Astronomical Society, 280(4):1235–1238, 1996. [15] E Bowell. A new protocol for the operation of the minor planet center. In AAS/Division for Planetary Sciences Meeting Abstracts# 31, volume 31, 1999. [16] Andrea Milani. The asteroid identification problem: I. recovery of lost asteroids. Icarus, 137(2):269–292, 1999. [17] Andrea Milani, Maria Eugenia Sansaturio, and Steven R Chesley. The asteroid identification problem iv: Attributions. Icarus, 151(2):150–159, 2001. [18] LG Taff and DL Hall. The use of angles and angular rates. Celestial mechanics, 16(4):481–488, 1977. [19] Giovanni Federico Gronchi, Linda Dimare, and Andrea Milani. Orbit determination with the two-body integrals. Celestial Mechanics and Dynamical Astronomy, 107(3):299–318, 2010. [20] Giovanni F Gronchi, Davide Farnocchia, and Linda Dimare. Orbit determination with the two-body integrals. ii. Celestial Mechanics and Dynamical Astronomy, 110(3):257–270, 2011. [21] Andrea Milani, Davide Farnocchia, Linda Dimare, Alessandro Rossi, and Fabrizio Bernardi. Innovative observing strategy and orbit determination for low earth orbit space debris. Planetary and Space Science, 62(1):10–22, 2012. [22] Alberto Quijano, M Rojas, and L Chaves. Cálculo de los parámetros orbitales del asteroide 2003qo104. Revista Colombiana de F´ısica, 42(2):4, 2010. [23] J.G.P. Barbosa. Elementos de astronomia de posicion. Coleccion textos. Universidad Nacional de Colombia. Facultad de Ciencias, 2009. [24] L.A. Pars. A treatise on analytical dynamics. London, 1968. [25] D.A. Vallado and W.D. McClain. Fundamentals of Astrodynamics and Applications. Managing Forest Ecosystems. Springer Netherlands, 2001. [26] E.T. Whittaker. A Treatise on the Analytical Dynamics of Particles and Rigid Bodies. Cambridge Mathematical Library. Cambridge University Press, 1988. [27] Edgar Everhart and Edward T Pitkin. Universal variables in the two-body problem. American Journal of Physics, 51(8):712–717, 1983.

BIBLIOGRAFÍA

105

[28] S. Herrick. Astrodynamics: Orbit determination, space navigation, celestial mechanics. Astrodynamics. Van Nostrand Reinhold Co., 1971. [29] Alessandra Celletti and Gabriella Pinzari. Dependence on the observational time intervals and domain of convergence of orbital determination methods. In Periodic, Quasi-Periodic and Chaotic Motions in Celestial Mechanics: Theory and Applications, pages 327–344. Springer, 2006. [30] D.A. Cox, J. Little, and D. O’Shea. Using Algebraic Geometry. Graduate Texts in Mathematics. Springer New York, 2013. [31] D.A. Cox, J. Little, and D. O’Shea. Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra. Undergraduate Texts in Mathematics. Springer International Publishing, 2015. [32] T. Becker, H. Kredel, and V. Weispfenning. Gröbner Bases: A Computational Approach to Commutative Algebra. Graduate Texts in Mathematics. Springer New York, 2012. [33] Bruno Buchberger. Gröbner bases: An algorithmic method in polynomial ideal theory. In Multidimensional Systems Theory–Progress, Directions and Open Problems in Multidimensional Systems, pages 184–232. 1985. [34] B Buchberger. An Algorithm for Finding a Basis for the Residue Class Ring of a Zero-dimensional Polynomial Ideal (German.) Ph. D. PhD thesis, Thesis, Math. Inst., Univ. of Insbruck, Austria, 1965. [35] Heisuke Hironaka. Resolution of singularities of an algebraic variety over a field of characteristic zero: Ii. Annals of Mathematics, pages 205–326, 1964. [36] Ernst W Mayr and Albert R Meyer. The complexity of the word problems for commutative semigroups and polynomial ideals. Advances in mathematics, 46(3):305–329, 1982. [37] Alessandro Giovini, Teo Mora, Gianfranco Niesi, Lorenzo Robbiano, and Carlo Traverso. one sugar cube, please or selection strategies in the buchberger algorithm. In Proceedings of the 1991 international symposium on Symbolic and algebraic computation, pages 49–54. ACM, 1991. [38] Jean-Charles Faugere, Patrizia Gianni, Daniel Lazard, and Teo Mora. Efficient computation of zero-dimensional gröbner bases by change of ordering. Journal of Symbolic Computation, 16(4):329–344, 1993. [39] Jean-Charles Faugere. A new efficient algorithm for computing gröbner bases (f 4). Journal of pure and applied algebra, 139(1):61–88, 1999.

106

BIBLIOGRAFÍA

[40] Elizabeth A Arnold. Modular algorithms for computing gröbner bases. Journal of Symbolic Computation, 35(4):403–419, 2003. [41] Shuhong Gao, Yinhua Guan, and Frank Volny IV. A new incremental algorithm for computing gröbner bases. In Proceedings of the 2010 International Symposium on Symbolic and Algebraic Computation, pages 13–19. ACM, 2010. [42] Ronald N Goldman, Thomas W Sederberg, and David C Anderson. Vector elimination: A technique for the implicitization, inversion, and intersection of planar parametric rational polynomial curves. Computer Aided Geometric Design, 1(4):327–356, 1984. [43] Dinesh Manocha. Solving systems of polynomial equations. Computer Graphics and Applications, IEEE, 14(2):46–55, 1994. [44] D. Bini and V. Pan. Polynomial and Matrix Computations: Fundamental Algorithms. Progress in Theoretical Computer Science. Birkhäuser Boston, 2012. [45] Andrea Milani, Giovanni F Gronchi, Mattia deMichieli Vitturi, and Zoran Kneˇzević. Orbit determination with very short arcs. i admissible regions. Celestial Mechanics and Dynamical Astronomy, 90(1-2):57–85, 2004. [46] Giovanni Federico Gronchi. Multiple solutions in preliminary orbit determination from three observations. Celestial Mechanics and Dynamical Astronomy, 103(4):301–326, 2009. [47] Daniele Mortari, Sante R Scuro, and Christian Bruccoleri. Attitude and orbit error in n-dimensional spaces. The Journal of the Astronautical Sciences, 54(3-4):467–484, 2006. [48] Andrew Vernon Schaeperkoetter. A comprehensive comparison between angles-only initial orbit determination techniques. PhD thesis, Texas A&M University, 2011. [49] Herbert Raab. Detecting and measuring faint point sources with a ccd. In Proceedings of the Meeting on Asteroid and Comets in Europe, pages 1–4, 2003. [50] R. Rosenberg. Analytical Dynamics of Discrete Systems. Mathematical Concepts and Methods in Science and Engineering. Springer US, 2012. [51] S.J. Orfanidis. Optimum signal processing: An introduction. McGraw-Hill international editions. Electrical engineering series. Mcgraw-Hill Book Comp., 1990.

BIBLIOGRAFÍA

107

[52] Mart´ın Avendano, Verónica Mart´ın-Molina, and Jorge Ortigas-Galindo. Solving keplers equation via smales\ alpha-theory. Celestial Mechanics and Dynamical Astronomy, 119(1):27–44, 2014. [53] Reza Raymond Karimi and Daniele Mortari. Initial orbit determination using multiple observations. Celestial Mechanics and Dynamical Astronomy, 109(2):167–180, 2011.

Trabajo de Grado Final

Recommend Documents