Integrador EVORB

(version 15, Julio de 2009)
Integracion orbital numerica de sistema de M particulas sin masa en el campo gravitacional de N planetas masivos y una estrella central.
VIRTUDES
  • Es un integrador orbital muy simple de manejar y no por ello menos preciso: se requiere un unico archivo de entrada: "datos.dat"
  • Esta pensado para quien desea iniciarse en la experimentacion de evolucion orbital en sistemas planetarios.
  • El programa maneja muy bien los encuentros de las particulas (sin masa, por ejemplo asteroides y cometas) con los planetas (masivos).
  • Considera efectos relativistas (opcional)
  • Integra hacia el futuro o hacia el pasado
  • Sistema Operativo: Este programa viene en 2 sabores, Windows o Linux,
  • Instalación: No necesita instalacion, simplemente descomprima evorb15.zip y copie los archivos a una carpeta.
    LIMITACIONES
  • Maximo de objetos que integra: 500 (modificable).
  • No integra un sistema estelar binario
  • No integra orbitas hiperbolicas, al detectarse una orbita (heliocentrica) hiperbolica el objeto es eliminado.
  • No puede haber encuentros proximos entre objetos masivos (planetas) pues el programa pierde precision.
  • No considera achatamiento de estrella central (queda para futuras ediciones).
    Este programa determina la evolucion de los elementos orbitales heliocentricos (a, e, i, longitud del nodo ascendente, argumento del perihelio, anomalia media) de un sistema de planetas (con cierta masa en masas solares y cierto radio en km) y particulas (sin radio ni masa como asteroides y cometas) bajo el efecto gravitacional de un cuerpo central (que puede ser el Sol u otra estrella) y de los demas cuerpos masivos. Se necesita al menos un planeta y una particula. Este programa se distribuye para fines puramente experimentales para el curso de Mecanica Celeste de Facultad de Ciencias. Si su nave espacial se aparta de la trayectoria indicada por este programa no nos hacemos responsables. Si necesita apoyo celeste consulte (gallardo@fisica.edu.uy) o lea este Manual del Mecánico Celeste (de dudosa ayuda). Utiliza el integrador original de A. Brunini (UNLP) y la adaptacion es de T. Gallardo (FC-UdelaR). Es un leapfrog tipo avace-impulso-avance. Referencias: En publicaciones cientificas este integrador puede ser referido como Fernández, J. A., Gallardo, T., and Brunini, A.: 2002, "Are There Many Inactive Jupiter Family Comets Among the Near-Earth Asteroid Population?", Icarus 159, 358-368. Este integrador con mayores o menores variantes ha sido usado en varios articulos cientificos sobre dinamica de asteroides, cometas y objetos transneptunianos. 
    ESQUEMA DE USO:

   1                 2
datos.dat ------> evorb15.exe -------> orbitas.sal  (elementos orbitales para cada intervalo de escritura)
                                          xyzv.sal  (posiciones y velocidades para cada intervalo de escritura)
                                       encuent.sal  (registro de encuentros de particulas con planetas)
                                       finales.sal  (elementos finales de los cuerpos eliminados)
                                       colisio.sal  (elementos de cuerpos que colisionaron)
                                       energia.sal  (control de precision)
                                        condin.xxx  (para continuar la integracion)

                                      3
encuent.sal + orbitas.sal ------> sacando.exe -----> evoxxxx.dat (evolucion del cuerpo xxxx)
                                                     encxxxx.dat (encuentros del cuerpo xxxx)

                                      3
                 xyzv.sal ------> sacandox.exe -----> xyzxxxx.dat (evolucion del cuerpo xxxx
                                                                   en coordenadas cartesianas)

    1) PRIMER PASO: preparar datos.dat

    En datos.dat (editarlo con bloc de notas o TextPad, atencion: no alterar el orden de las lineas ni introducir nuevas) se deben indicar los detalles de la integracion numerica:
    1. indice para considerar o no efectos relativistas (0, 1, 2)
    2. indice para guardar o no archivo con posiciones y velocidades (0, 1)
    3. distancia minima y maxima a la estrella permitidas (por si no interesa lo que ocurra con cuerpos muy proximos a la estrella o muy lejanos)
    4. masa central en masas solares y radio en kms (el radio para simular colisiones con la estrella)
    5. tiempo total de la integracion en años julianos (de 365.25 dias). Si desea integrar hacia el pasado basta poner un tiempo negativo.
    6. intervalo de escritura de datos (si es muy chico se llenara el disco duro...). DEBE SER MULTIPLO DEL PASO.
    7. paso: si es grande va rapido pero acumula errores, se recomienda hacer algunas pruebas con paso del orden de la cincuenta-ava parte del menor periodo orbital de los objetos intervinientes
    8. indice: poniendo 0 empieza una nueva integracion y poniendo 1 continua con una anterior
    9. numero de planetas que intervienen, pueden ser reales o ficticios (minimo uno)
    10. epoca para los elementos de esos planetas
    11. los elementos de esos planetas, una linea por cada planeta que interviene
    12. numero de particulas que intervienen (minimo una)
    13. epoca para los elementos de esas particulas (este sera el instante inicial de la integracion)
    14. los elementos de esas particulas, una linea por cada particula que interviene (maximo 500 particulas+planetas)
    Ejemplo, si desea la evolucion de un asteroide o cometa obtenga sus elementos orbitales de alguno de los lugares citados abajo o del Astronomical Almanac (seccion G), pongalo como particula con un numero identificatorio arbitrario (menor a 10000) y corra el programa usando todos los (8) planetas. A veces es posible despreciar Mercurio y agregar su masa al Sol.

    Importante: no elimine los comentarios (*) de datos.dat pues el programa asume que estan alli. Luego de la lectura de la ultima linea de datos de las particulas el programa no leera nada mas de datos.dat asi que puede escribir y guardar lo que quiera alli.

    Un practico conversor de fechas a JD se encuentra en http://aa.usno.navy.mil/AA/data/docs/JulianDate.html

    Para saber el JD exacto de un evento en la integracion que ocurre en el tiempo t basta calcular: JD(t) = JD(EpocaParticulas) + t*365.25



    Ejemplo de datos.dat: 
    	*INDICE CORRECCION RELATIVISTA (0 NO, 1 modelo Venturini-Gallardo, 2 formula exacta)
	0
	*INDICE PARA GUARDAR (1) O NO (0) ARCHIVO CON (X,Y,Z,Vx,Vy,Vz)
	1
	*LIMITES INFERIOR Y SUPERIOR EN UAs PARA ELIMINACION DE CUERPOS
	0.005  1000.0
	*MASA CENTRAL EN MASAS SOLARES Y RADIO EN kms
	1.0  696000.0
	*TIEMPO MAXIMO DE LA SIMULACION EN AÑOS:
	1000.0
	*INTERVALO DE ESCRITURA EN AÑOS (multiplo del PASO):
	100.0
	*PASO DE INTEGRACION en años:
	0.005
	*INDICE PARA RE-ARRANQUE (0 ES UNA CORRIDA NUEVA, 1 ES CONTINUACION)
	0
	*NUMERO DE PLANETAS (minimo uno)
	8
	*JD epoca de planetas
	2451400.5
	*elementos J2000 de los planetas: N a e i nodo argper anomedia masa radio(km)
	1  0.38709897  0.205622188  7.0050053  48.3307302  29.1239378 303.4539666 1.66013679D-07  2440.0
	2  0.72333793  0.006771606  3.3945712  76.6821802  54.8201213 178.9704486 2.44783833D-06  6052.0
	3  0.99998836  0.016736826  0.0001213  27.6464539  75.3942093 215.0034486 3.04043264D-06  6371.0
	4  1.52367503  0.093290001  1.8498920  49.5616397 286.5090192 303.6562023 3.22715144D-07  3390.0
	5  5.20399661  0.048719684  1.3046775 100.4829904 275.0927233   6.8013147 9.54791938D-04 69911.0
	6  9.58195033  0.055272778  2.4852156 113.6436321 335.5040835 315.9800041 2.85885980D-04 58232.0
	7 19.24744286  0.043713914  0.7727272  74.0128993  97.4115217 140.3579592 4.36624404D-05 25362.0
	8 30.14125054  0.011166945  1.7676331 131.7957546 259.2549938 273.2914807 5.15138902D-05 24624.0
	*NUMERO DE PARTICULAS (minimo una)
	4
	*JD epoca de particulas
	2452880.5
	*elementos J2000 de ESAS particulas: N a e i nodo argper anomedia
	10 6.927713 0.8156838 29.0741  218.0082  171.1487  177.1351
	11 6.001769 0.7876825 18.70094  22.21415 257.2521    0.0513757
	12 5.671776 0.8240879 54.69256 270.5487  206.7027  244.5646
	13 5.132881 0.7743644  4.87674  14.41821  22.32553 196.9951
    Con estos datos estariamos integrando sin considerar efectos relativistas por 1000 años con salida cada 100 años y con paso 0.005 años a los 8 planetas mas 4 objetos que en este caso son los cometas P/McNaught-Russell(1994X1) , 66P/duToit , 8P/Tuttle y 85P/Boethin. Se guardaran tambien los datos de posiciones y velocidades y todo objeto que se aproxime 0.005 UA del Sol o se aleje mas de 1000.0 UA del Sol sera eliminado.

    El numero entero identificatorio de cada objeto no debe repetirse y debe estar comprendido entre 1 y 9999. Es sano adjudicar numeros bajos (1, 2,...) a los planetas y altos a las particulas (10, 11, ...).

    2) SEGUNDO PASO: ejecutar evorb15.exe

    Ejecute evorb15.exe (haciendo doble clic en Windows o ./evorb15 en Linux) y espere.... Una vez terminada la integracion aparecerán los siguientes archivos que pueden abrirse con bloc de notas o TexPad:
    • En "orbitas.sal" esta el detalle de evolucion de todos los cuerpos. En general este es el archivo que estamos buscando.
    • En "xyzv.sal" esta el detalle de evolucion de todos los cuerpos en coordenadas rectangulares (yo esto lo uso poco).
    • En "encuent.sal" habra un listado y detalles de los encuentros a menos de 2 radios de Hill: instante, distancia minima en radios planetarios (si es menor que 1 hubo colision), velocidad de encuentro en pericentro, energia planetocentrica (si es negativa estaba temporariamente capturado), planeta, particula.
    • En "finales.sal" estan las condiciones finales de los objetos que fueron eliminados por eyeccion del sistema o por colision. Se guardan los elementos orbitales anteriores al impulso que eliminó a la particula.
    • En "colisio.sal" estan los datos de los objetos que colisionaron con un planeta (eventos raros) o con el Sol (indice 0).
    • En "energia.sal" tenemos para control la evolucion de la energia del sistema de cuerpos masivos, que debe mantenerse constante con el tiempo. Se incluye tambien la energia cinetica y la potencial a efectos de observar su oscilacion.
    • No intente abrir el archivo "condin.xxx" pues esta en binario y contiene los datos de los objetos de la ultima integracion por si ud desea continuarla luego (por ejemplo a causa de un corte de energia) colocando previamente el indice 1 en datos.dat.
    • Si hay un solo planeta, automaticamente se generara un archivo "cjacobi.sal" conteniendo la evolución de la cte de Jacobi para la primer particula. Se utiliza para probar la precision del integrador colocando un planeta con e=0, i=0.

    3) TERCER PASO: ejecutar sacando.exe y/o sacandox.exe

    Para obtener un archivo con la evolucion para un objeto particular corra el programa "sacando.exe" haciendo doble clic. Ud introducira el numero del objeto y se generan 2 archivos para el objeto conteniendo su evolucion orbital y su historia de encuentros. Estos archivos son facilmente graficables con algun programa (como gnuplot). Por ejemplo, en esta figura graficamos la evolucion de las excentricidades de las orbitas de Venus, Tierra (verde y azul, vea que están acopladas) y Marte (rojo, vea que no siempre será tan excentrico como en la actualidad) en los próximos 5 millones de años. Si le interesa la evolucion de la posicion y velocidad de algun objeto debera ejecutar "sacandox.exe". Esta ultima opcion suele emplearse para ver el detalle del encuentro entre 2 cuerpos o para graficar la trayectoria espacial de un objeto (se requerirá un intervalo de escritura bastante pequeño).     		?

    Continuacion De Integracion Anterior

    Si desea continuar una integracion terminada o cortada por un corte de energia tiene 2 opciones:
    1. en datos.dat ponga el nuevo tiempo final (que puede ser el original si se corto la integracion) y sustituya el 0 por un 1 en el indice de re-arranque. El archivo "condin.xxx" debe estar disponible. Ejecute el programa. Los nuevos resultados generados de la integracion se agregan a los archivos anteriores (si ya existen) o se crean nuevos archivos (si no existen).
    2. otra opcion (menos precisa) es tomar los elementos de orbitas.sal y ponerlos en datos.dat agregando masas y radios de planetas y poniendo la misma epoca (cualquiera, es irrelevante, pero la misma para planetas y particulas) para planetas y particulas. Previamente guarde en lugar seguro los archivos generados por la corrida anterior.

    Advertencias

    • Si el paso es el adecuado y no hay encuentros la integracion es muy precisa
    • Si hay encuentros los resultados para las particulas son precisos para los primeros encuentros luego los resultados seran cualitativamente validos
    • Orbitas con periodo menor que paso*50 y que ademas tienen excentricidad mayor que 0.98 comienzan a perder precision en la integracion. La solucion es disminuir el paso.
    • Los encuentros con los planetas cuyo periodo sea < paso*50 estaran mal calculados. Ejemplo: si se integra a Mercurio con particulas que lo encuentran sera necesario usar paso = 0.005 años de lo contrario apareceran eliminaciones espureas por mal manejo de encuentros.
    • La correccion relativista (colocando indice 1) se introduce con una perturbacion que emula perfectamente el efecto secular en el argumento del perihelio (Venturini-Gallardo 2010). Esta perturbacion genera una pequeña peturbacion espurea en la anomalia media pero esto no tiene ningun efecto dinamico importante. Es posible aplicar la correccion relativista exacta estandar debida al cuerpo central (ver Benitez & Gallardo 2008, CMDA) colocando indice 2, pero si el paso no es suficientemente pequeño se introduciran errores (a su propio riesgo).
    • Originalmente EVORB fue simpletico pero actualmente en esta version es un simple leapfrog con una muy buena precision. Si duda chequee con Mercury por ejemplo, lo cual es muy saludable siempre.
    • La version linux es nueva, no ha sido muy chequeada aun...

    Bugs

    En la version 15 se introducen ajustes para calcular con mayor precision orbitas de muy alta excentricidad (0.99999). En la version 14 se usa una correccion relativista (Venturini-Gallardo 2010) que no genera errores en caso de orbitas de perihelio muy pequeño como ocurria en la version 13. En la version 13 se agregan 2 lineas al inicio de "datos.dat" para hacer opcional la correccion relativista. Se corrige conteo del tiempo en pantalla (en integraciones largas se acumulaba un error de redondeo). A partir de la version 12.R se consideran efectos relativistas provocados por el cuerpo central segun Quinn et al. (1991, AJ, 101, 2287). Cuando la distancia heliocentrica es muy pequeña la correccion relativista introduce un error notorio. Por esta razon en la version 12.RC limitamos la aplicacion de la correccion relativista a distancias mayores que un cierto limite (r**1.5 > 20*paso). La version 12.0 sustituye a las anteriores y supuestamente hemos corregido los errores de las versiones anteriores. En la version 11.5 corregimos error en el calculo de la energia total del sistema (para control) que aparecia en caso de sistemas donde la estrella central tuviera masa diferente de 1. Ahora ademas se guardan la energia cinetica y la potencial para los curiosos. En la version 11.4, gracias a Esmeralda Mallada, corregimos un error en el calculo de la distancia de un encuentro que se manifiesta en caso de encuentros muy alejados (mas de un radio de Hill), por causa de este error aparecian encuentros muy proximos en "encuent.sal" cuando en realidad eran muy alejados. En la version 11.2 corregimos un pequeño error en el calculo del instante de los encuentros proximos. En alguna compilacion anterior el programa abortaba por faltar una libreria que el usuario desprevenido seguramente no tendria. Ahora se soluciono en esta version compilada con g77.

    Temas Relacionados

    Elementos Orbitales

    NOTA: Ud puede estar obsesionado con utilizar elementos orbitales superprecisos. No le servira de nada si no usa en el integrador el mismo modelo dinamico utilizado para la determinacion de esos elementos orbitales. Por ejemplo, si usa elementos determinados con un modelo relativista debera usar la opcion de integracion relativista. PROGRAMAS PARA SELECCION Y VISUALIZACION DE ORBITAS. Para seleccionar orbitas de MPCORB: BROMPCTM. Si le interesa la trayectoria de un objeto respecto a un sistema rotante Sol-Planeta: NEAPLOT.

    PROGRAMA PARA DETERMINACION DE ORBITAS. A partir de datos observacionales es posible determinar orbitas en forma sencilla utilizando por ejemplo FindOrb. Tambien es posible hacerlo de forma un poquito mas complicada pero precisa con OrbFit.

    Otros integradores

    Fortran

    Compilador y Editor FORCE 2.0 (libre) http://force.lepsch.com/

    Accesorios para manejo de datos

    Un práctico y poderoso graficador de datos lo encuentra en GNUPLOT. Un editor muy recomendable con sorting y manejo de columnas y bloques de datos es TextPad.