7.15. El problema de los dos cuerpos en SPICE

Algunos de los procedimientos descritos en este capítulo han sido implementados en la biblioteca de rutinas de SPICE. Las siguientes son las rutinas disponibles en dicho sistema para realizar tareas relacionadas con el problema de dos cuerpos:

  • oscelt(X,t,mu): Calcula los elementos orbitales de la órbita osculatriz (ver Sección Orbita osculatriz) correspondiente con el vector de estado X. Esta rutina devuelve los siguientes elementos (en ese orden): $q$, $e$, $i$, $\Omega$, $\omega$, $M$, donde $q$ es la distancia al periapsis y $M$ es la anomalía media (elíptica o hiperbólica.) Todos los ángulos son calculados en radianes. La rutina devuelve además el valor de t y de mu, que aunque son los mismos que los provistos en la entrada, permiten tener una salida consistente con la rutina conics descrita más abajo. Esta rutina realiza el mismo trabajo que la rutina estado_a_elementos descrita en la Sección Determinación de la órbita.

  • osceltx(X,t,mu): Esta rutina hace los mismo que oscelt pero devuelve una lista extendida de elementos orbitales, en este orden: $q$, $e$, $i$, $\Omega$, $\omega$, $M$, $t$, $\mu$, $f$ (anomalía verdadera), $a$, $T$ (período orbital, cuando es posible calcularlo.)

  • conics(E,t): Calcula el vector de estado en el tiempo $t$ correspondiente a la cónica con elementos orbitales E: $q$ , $e$, $i$, $\Omega$, $\omega$, $M_0$, $t_0$, $\mu$; donde $t_0$ y $M_0$ son el tiempo y la anomalía media correspondiente a los elementos provistos. La rutina devuelve el vector de estado. Esta rutina realiza, parcialmente, el mismo trabajo que la rutina elementos_a_estado descrita en la Sección Prediccion del vector de estado, con la diferencia que además y a diferencia de nuestra rutina que solo hace una transformación geométrica, conics también realiza la propagación del estado en el tiempo.

  • prop2b(mu,X0,dt): Propaga el estado inicial X0 del vector relativo de un sistema con parámetro gravitacional mu por un tiempo dt=t-t0. Esta rutina utiliza internamente el formalismo de variables universales que hemos esbozado en las últimas secciones. La rutina devuelve el estado relativo en el tiempo t. Esta rutina realiza el mismo trabajo que la rutina propaga_f_g descrita en la Sección Funciones $f$ y $g$.

Nota: Singularidades y errores. El lector más curioso podrá preguntarse si las rutinas de descritas aquí ya habían sido implementadas en distintos apartes en este capítulo, cuál puede ser entonces el interés de introducir las rutinas específicas de SPICE. Hay dos razones fundamentales para hacerlo. La primera es que las rutinas del sistema de NASA tienen una larga historia de desarrollo y pruebas que las hace muy confiables. Podemos usar esas rutinas para comprobar las que desarrollemos por nuestra cuenta.

La segunda razón y en realidad la más importante es que las transformaciones estudiadas en este capítulo (de elementos a vector de estado, de vector de estado a elementos, de vector de estado a vector de estado) tienen algunas limitaciones numéricas. Así por ejemplo, los elementos osculatrices calculados a partir del vector de estado pueden ser muy inciertos cuando la inclinación es cercana a cero o a 180$^\circ$ o cuando la excentricidad es muy cercana a 1. La solución a la ecuación universal de Kepler en $g$ (usada por nuestras rutinas y por las rutinas de SPICE) debe calcularse con precaución para valores grandes del tiempo (en los que la precisión numérica de las series de Stumpff puede estar muy limitada.) Si bien las rutinas de SPICE no necesariamente corrijen todos esos inconvenientes, el sistema viene dotado de mecanismos de control de errores que no hemos implementado en las rutinas desarrolladas en y para el libro.

A continuación mostramos algunos ejemplos del uso de estas rutinas y comparamos sus resultados con las rutinas desarrolladas en este capítulo. Usaremos para ello el sistema de ejemplo introducido en la Sección Un ejemplo numérico.

Para usar las rutinas de SPICE debemos primero construir el vector de estado inicial:

In [1]:
#Sistema
from numpy import array
t0=0
sistema=[
    dict(m=1.0,
         r=array([0.0,0.0,+0.3]),
         v=array([+1.0,0.0,0.5])),
    dict(m=0.5,
         r=array([+1.0,0.0,0.0]),
         v=array([0.0,+1.0,0.0])),
]

#Condiciones iniciales
m1=sistema[0]["m"]
r1_0=sistema[0]["r"]
v1_0=sistema[0]["v"]

m2=sistema[1]["m"]
r2_0=sistema[1]["r"]
v2_0=sistema[1]["v"]

mu=m1+m2 

#Posición y velocidad relativa inicial
rvec0=r1_0-r2_0
vvec0=v1_0-v2_0

#Vector de estado inicial
from numpy import append
X0=append(rvec0,vvec0)
X0 = [-1.   0.   0.3  1.  -1.   0.5]

Calculemos los elementos orbitales correspondientes a este estado:

In [3]:
from spiceypy import oscelt
q,e,i,W,w,M0,t0,mu=oscelt(X0,t0,mu)
q = 0.669944
e = 0.721536
Inclinación = 40.5106 grados
Longitud del nodo ascendente = 159.444 grados
Argumento del periapsis = 107.911 grados
Anomalía media = 347.311 grados

Que coincide con los resultados que obtuvimos en la Sección Un ejemplo numérico. Una buena manera de verificar si los resultados producidos por esta rutina son correctos es invocar su inversa conics:

In [5]:
from spiceypy import conics
X0=conics([q,e,i,W,w,M0,t0,mu],t0)
X0 = [-1.   0.   0.3  1.  -1.   0.5]

Que coincide con el estado inicial. Finalmente podemos propagar este estado hasta $t=10$ como lo hicimos en secciones anterioes:

In [7]:
#Tiempo de propagación
t0=0.0
t=10.0
dt=t-t0

from spiceypy import prop2b
X=prop2b(mu,X0,dt)
Estado final: [-0.0642266  3.2416631 -2.5740625 -0.3095438  0.0535113  0.0500541]

Que de nuevo, coincide con los resultados obtenidos en el Alg. (ejemplo_propaga_fg). La propagación también puede hacerse con la rutina conics:

In [9]:
X=conics([q,e,i,W,w,M0,t0,mu],t)
Estado final: [-0.0642266  3.2416631 -2.5740625 -0.3095438  0.0535113  0.0500541]

Que como era de esperarse coincide con los resultados anteriores.