Resolviendo computacionalmente el problema de dos cuerpos para la aproximación de cónicas parcheadas

Probablemente ya encontraste una respuesta cuando publicaste esta pregunta, pero todavía doy una respuesta en caso de que te ayude.

Usar los elementos orbitales es la "manera correcta" de hacerlo. Hay 6 de ellos, descritos en este artículo de Wikipedia :

• excentricidad$e$
• semieje mayor$a$
• inclinación$i$
• longitud del nodo ascendente$\Omega$
• argumento del periapsis$\omega$
• verdadera anomalía$\theta$(que depende del tiempo)

Los primeros cinco describen completamente la geometría de una órbita en 3D. La verdadera anomalía es un ángulo que se refiere a la posición del cuerpo/nave espacial en esa órbita. Puede calcular todos estos parámetros a partir de vectores de estado inicial (posición y velocidad), conociendo el parámetro gravitatorio estándar del cuerpo atractor. René Schwarz describe las matemáticas en este PDF sobre cómo convertir vectores de estado orbitales cartesianos en elementos orbitales.

Una vez que tenga los elementos de una órbita, puede usarlos para calcular los vectores de estado (vector de posición/radio y vector de velocidad) en cualquier momento siguiendo este otro PDF sobre cómo convertir elementos orbitales en vectores de estado cartesianos. Sin embargo, algunas fórmulas descritas en el documento solo funcionan para órbitas elípticas. Para órbitas hiperbólicas (cuando$e > 1$) tendrás que usar otras fórmulas para la anomalía excéntrica y verdadera (ver este artículo de Wikipedia y la segunda respuesta en esta publicación ), así como adaptar la fórmula para el vector de velocidad del documento (ver esta publicación ).

Las matrices de rotación de aspecto "complicado" que establecen la posición y la velocidad de la nave espacial/cuerpo en el espacio 3D se pueden desarrollar en estos pasos:

• Elija los vectores de referencia para "derecha" y "arriba", por ejemplo, la dirección x (1, 0, 0) y la dirección y (0, 1, 0), los llamo$\vec{u_x}$y$\vec{u_y}$.
• Calcule el vector de nodo ascendente$\vec{u_{asc}}$: es$\vec{u_x}$girado alrededor$\vec{u_y}$por un ángulo$\Omega$;
• Establecer la posición$\vec{r} = \vec{o}$y la velocidad$\vec{v} = \dot{\vec{o}}$como se describe en el documento (con la fórmula modificada en caso$e > 1$);
• Girar$\vec{r}$y$\vec{v}$alrededor$\vec{u_y}$de un ángulo$\Omega + \omega$;
• Girar de nuevo$\vec{r}$y$\vec{v}$alrededor$\vec{u_{asc}}$de un ángulo$i$.

(Tenga en cuenta que en los documentos de René Schwarz, considera el eje z como el eje "arriba", en lugar del eje y)

El caso de una órbita parabólica ($e = 1$) puede ser, como dijo notovny en un comentario, cortésmente ignorado (aplicar el algoritmo de avestruz ).

Sin embargo, tenga en cuenta que es posible que deba tener en cuenta los casos extremos. es decir, cuando se trata de órbitas que son muy cercanas a las circulares o tienen una inclinación nula. En estos casos, las fórmulas para ángulos, como la anomalía verdadera en el primer documento, no funcionarán (darán valores de NaN debido a la precisión del punto flotante y los redondeos). Por lo tanto, para estos casos particulares, deberá establecer vectores predeterminados, por ejemplo, el vector excéntrico, y calcular los ángulos a partir de él.

En cuanto a las unidades, las fórmulas de hecho esperan unidades SI. (Sin embargo, puede usar cualquier múltiplo de estas unidades siempre que lo tenga en cuenta en las ecuaciones).