Elementos orbitales y conversión estado ↔ elementos

1. Motivación

Un vector de estado $(\mathbf r, \mathbf v)$ fija la órbita, pero sus seis números cambian a cada instante y no dicen nada a simple vista. Los elementos orbitales son otra parametrización de las mismas seis constantes: cinco describen la geometría de la cónica y su orientación en el espacio, y son constantes en el problema de dos cuerpos; el sexto dice dónde está el cuerpo. Las efemérides del JPL, los catálogos de asteroides y los TLE de los satélites se publican en elementos; convertir en ambos sentidos es la operación más básica de la astrodinámica.

2. Los seis elementos keplerianos

Fijamos un sistema inercial $(\hat{\mathbf x}, \hat{\mathbf y}, \hat{\mathbf z})$: el plano $xy$ es el plano de referencia (la eclíptica J2000 para órbitas heliocéntricas, el ecuador ICRF para satélites terrestres) y $\hat{\mathbf x}$ apunta al equinoccio vernal $\gamma$ de J2000. Con $\mathbf h = \mathbf r\times\mathbf v$ y $\mathbf A = \mathbf v\times\mathbf h - \mu\hat{\mathbf r}$ (lección 0):

El vector $\mathbf n = \hat{\mathbf z}\times\mathbf h = (-h_y, h_x, 0)$ está en la intersección de los dos planos y apunta hacia el nodo ascendente, el punto en que el cuerpo cruza el plano de referencia hacia $z > 0$ (Problema 2). Una órbita es prograda si $i < 90°$ y retrógrada si $i > 90°$. Los ángulos $\Omega$ y $\omega$ se miden en el sentido del movimiento, en $[0, 2\pi)$. Como alternativa a $\omega$ y $\nu$ se usan la longitud del periapsis $\varpi = \Omega + \omega$ y la longitud media $\lambda = \varpi + M$, que siguen definidas cuando $i\to0$.

3. El sistema perifocal y las matrices de rotación

En el sistema perifocal $(\hat{\mathbf p}, \hat{\mathbf q}, \hat{\mathbf w})$, con $\hat{\mathbf p} = \mathbf A/A$ hacia el periapsis, $\hat{\mathbf w} = \mathbf h/h$ y $\hat{\mathbf q} = \hat{\mathbf w}\times\hat{\mathbf p}$, la órbita es plana y su descripción es la de la lección 0:

(la segunda es la hodógrafa del Problema 3 de la lección 0). Para pasar al sistema inercial hay que deshacer tres giros. Definimos las rotaciones pasivas (giran los ejes, no el vector) alrededor de $x$ y de $z$:

Del sistema inercial al perifocal se llega con la secuencia de Euler 3-1-3: girar $\Omega$ alrededor de $\hat{\mathbf z}$ (ahora $\hat{\mathbf x}'$ apunta al nodo), girar $i$ alrededor de la línea de nodos (ahora $\hat{\mathbf z}'' = \hat{\mathbf w}$) y girar $\omega$ alrededor de $\hat{\mathbf w}$ (ahora $\hat{\mathbf x}''' = \hat{\mathbf p}$). En componentes, $\mathbf r_{pf} = R_3(\omega)R_1(i)R_3(\Omega)\,\mathbf r$. Como cada $R_k$ es ortogonal, la inversa es la transpuesta, $R_k(\theta)^{-1} = R_k(\theta)^T = R_k(-\theta)$, y

Multiplicando (Problema 1) se obtiene la matriz explícita, con $c_\Omega = \cos\Omega$, $s_\Omega = \sin\Omega$, etc.:

Sus columnas son $\hat{\mathbf p}$, $\hat{\mathbf q}$, $\hat{\mathbf w}$ expresados en el sistema inercial; en particular la tercera columna es $\mathbf h/h$, lo que da una comprobación inmediata de cualquier implementación (el test de elementos.ts lo usa). La misma $Q$ transforma velocidades, porque la rotación es constante en el tiempo.

4. Algoritmo estado → elementos

Dados $\mathbf r$, $\mathbf v$, $\mu$ (Curtis, algoritmo 4.2; Vallado, algoritmo 9):

1. $\mathbf h = \mathbf r\times\mathbf v$, $h = \lVert\mathbf h\rVert$; $\ \mathbf n = \hat{\mathbf z}\times\mathbf h$, $n = \lVert\mathbf n\rVert$.
2. $\mathbf e = (\mathbf v\times\mathbf h)/\mu - \mathbf r/r$, $e = \lVert\mathbf e\rVert$; $\ \varepsilon = v^2/2 - \mu/r$, $a = -\mu/(2\varepsilon)$ (si $e \ne 1$), $p = h^2/\mu$.
3. $i = \arccos(h_z/h)$. Como $h > 0$ y $i\in[0,\pi]$, no hay ambigüedad.
4. $\Omega = \operatorname{atan2}(n_y, n_x)$, llevado a $[0, 2\pi)$. (Equivale a $\arccos(n_x/n)$ con la corrección $\Omega \to 2\pi - \Omega$ si $n_y < 0$.)
5. $\omega = \arccos\big((\mathbf n\cdot\mathbf e)/(n e)\big)$, y $\omega \to 2\pi - \omega$ si $e_z < 0$ (el periapsis está bajo el plano de referencia).
6. $\nu = \arccos\big((\mathbf e\cdot\mathbf r)/(e r)\big)$, y $\nu \to 2\pi - \nu$ si $\mathbf r\cdot\mathbf v < 0$ (el cuerpo se acerca al periapsis).
7. Si se quiere $M$: $E = 2\operatorname{atan2}\big(\sqrt{1-e}\,\sin\tfrac\nu2,\ \sqrt{1+e}\,\cos\tfrac\nu2\big)$ y $M = E - e\sin E$; el paso por el periapsis es $\tau = t - M/n$.

Casos degenerados. Si $i = 0$ o $\pi$, $\mathbf n = \mathbf 0$ y $\Omega$ no existe: se toma $\Omega = 0$ y $\omega$ se mide desde $\hat{\mathbf x}$ (es la longitud del periapsis $\varpi$). Si $e = 0$, $\omega$ no existe: se toma $\omega = 0$ y $\nu$ se mide desde el nodo (argumento de la latitud $u = \omega + \nu$). Si ambas cosas ocurren, $\nu$ se mide desde $\hat{\mathbf x}$ (longitud verdadera). En la implementación hay que decidir con una tolerancia ($e < 10^{-10}$, $n < 10^{-10}h$) y documentarla: la sección 7 explica por qué.

5. Algoritmo elementos → estado

Dados $a, e, i, \Omega, \omega$ y $\nu$ (o $M$, resolviendo la ecuación de Kepler para obtener $E$ y luego $\nu$):

1. $p = a(1-e^2)$, $r = p/(1+e\cos\nu)$.
2. $\mathbf r_{pf}$ y $\mathbf v_{pf}$ como en la sección 3.
3. $\mathbf r = Q\,\mathbf r_{pf}$, $\mathbf v = Q\,\mathbf v_{pf}$ con $Q = R_3(-\Omega)R_1(-i)R_3(-\omega)$.

Para propagar la órbita un intervalo $\Delta t$ basta cambiar $M \to M + n\,\Delta t$ y rehacer el paso 3: ningún otro elemento varía en el problema de dos cuerpos. Esta es la forma más barata y más precisa de propagar una órbita kepleriana (no hay error de integración, solo el de resolver la ecuación de Kepler, que es de redondeo).

6. Verificación: Marte en J2000

Partimos del vector de Horizons (lección 0, sección 6) y aplicamos el algoritmo de la sección 4. Comparamos con la tabla ELEMENTS de Horizons para la misma época y el mismo marco (eclíptica J2000, centro Sol, $\mu = 1.3271248287\times10^{11}$ km³/s²):

La concordancia es la del redondeo de los datos de entrada (16 cifras). El paso por el perihelio, $\tau = 2451545.0 - M/n$ con $n = 6.0652706\times10^{-6}$ °/s $= 1.0585894\times10^{-7}$ rad/s, cae el 25 de noviembre de 1999. La conversión inversa devuelve el vector original con error menor que 1 cm en posición y $10^{-9}$ km/s en velocidad (test «ida y vuelta»).

Osculador vs. medio. Estos son elementos osculadores: los de la cónica tangente a la trayectoria real en ese instante. Un mes después Horizons da valores ligeramente distintos porque Júpiter y los demás planetas perturban a Marte. Para la Tierra el efecto es aún más visible: su vector en J2000 da $a = 1.49665\times10^{8}$ km (no 1 UA) e $i = 4.2\times10^{-4}$ grados (no 0), porque el centro de la Tierra orbita el baricentro Tierra–Luna (lección 0, Problema 6).

7. Análisis numérico: condicionamiento

Cuadrantes. Las fórmulas con $\arccos$ devuelven ángulos en $[0,\pi]$; sin las correcciones de signo de los pasos 4–6 se obtiene el ángulo equivocado la mitad de las veces. Con $\operatorname{atan2}$ se evita el problema de raíz: por ejemplo $\nu = \operatorname{atan2}(e\sin\nu,\ e\cos\nu)$ con $e\cos\nu = p/r - 1$ y $e\sin\nu = (\mathbf r\cdot\mathbf v)\,h/(\mu r)$, que además es robusto cuando $e\to0$.

Singularidades. $\Omega$ se obtiene de $\mathbf n$, cuyo módulo es $h\sin i$; para $i\to0$ el cálculo divide por un número pequeño y el error en $\Omega$ crece como $1/\sin i$: para la Tierra ($i = 4.18\times10^{-4}$°) el factor de amplificación es $1/\sin i \approx 1.4\times10^{5}$. Del mismo modo, $\omega$ y $\nu$ se obtienen de $\mathbf e$ y su error crece como $1/e$. Las sumas $\varpi = \Omega+\omega$ y $\lambda = \varpi + M$ están, en cambio, bien definidas y bien condicionadas: los errores opuestos de $\Omega$ y $\omega$ se cancelan. Esto motiva los elementos equinocciales (Broucke & Cefola 1972):

junto con $a$: son continuos y no singulares para $e = 0$ e $i = 0$ (solo fallan en $i = 180°$), y son los que usan los propagadores semianalíticos modernos.

Tolerancias. La decisión «$e = 0$» con un umbral de $10^{-10}$ introduce una discontinuidad: para $e = 1.1\times10^{-10}$ se calcula $\omega$ (con error relativo $\epsilon_{\text{máq}}/e \sim 10^{-6}$) y para $e = 0.9\times10^{-10}$ se fija $\omega = 0$. Las magnitudes físicas ($\mathbf r$, $\mathbf v$, $\varpi+\nu$) son continuas a través del umbral; solo la partición en $\omega$ y $\nu$ salta. Un buen test unitario comprueba precisamente eso.

8. Simulador

El simulador trabaja en el plano perifocal (el eje $x$ apunta al periapsis). Ajusta $a = 1.5237$ UA y $e = 0.0933$ para ver la órbita de Marte en su propio plano; lo que falta para colocarla en el espacio son exactamente $i$, $\Omega$ y $\omega$: la matriz $Q$. Comprueba que con $e = 0$ la posición del periapsis (y por tanto $\omega$) deja de tener sentido.

Problemas

1. (Demostración) Multiplica $R_3(-\Omega)R_1(-i)R_3(-\omega)$ y obtén la matriz $Q$ de la sección 3. Comprueba que $Q^TQ = I$, que $\det Q = 1$ y que su tercera columna es $(\sin\Omega\sin i,\ -\cos\Omega\sin i,\ \cos i)$. Interpreta esta última.

2. (Demostración) Demuestra que $\mathbf n = \hat{\mathbf z}\times\mathbf h$ apunta al nodo ascendente y no al descendente.

3. (Cuadrantes) Justifica las dos correcciones de signo: $\omega\to2\pi-\omega$ si $e_z < 0$ y $\nu\to2\pi-\nu$ si $\mathbf r\cdot\mathbf v < 0$.

4. (Propagación) Con los elementos de Marte en J2000, calcula su anomalía verdadera y su distancia al Sol 100 días después, en el modelo de dos cuerpos. (Horizons: $n = 6.0652706\times10^{-6}$ °/s, $M_0 = 19.356483°$.)

5. (Equinocciales) Expresa $e$, $i$, $\Omega$ y $\varpi$ en función de $(h_e, k_e, p_e, q_e)$ y demuestra que la aplicación $(e,\varpi)\mapsto(h_e,k_e)$ es un difeomorfismo del disco $e<1$ que incluye el origen, mientras que $(e,\varpi)$ como coordenadas polares no lo son en $e = 0$.

6. (Computacional) Implementa ambos algoritmos (o usa src/lib/math/elementos.ts) y: (a) reproduce la tabla de la sección 6; (b) comprueba la ida y vuelta para $10^4$ estados aleatorios con $e\in[0,0.99]$, $i\in[0,\pi]$ y mide el error máximo; (c) repite con $e = 10^{-k}$, $k = 2,\dots,12$, y representa el error en $\omega$ y en $\varpi + \nu$ frente a $e$. ¿Qué pendiente logarítmica esperas en cada caso?

import numpy as np
def estado_a_elementos(r, v, mu):
    h = np.cross(r, v); hn = np.linalg.norm(h); n = np.array([-h[1], h[0], 0.0])
    ev = np.cross(v, h)/mu - r/np.linalg.norm(r); e = np.linalg.norm(ev)
    a = -mu/(2*(np.dot(v,v)/2 - mu/np.linalg.norm(r)))
    i = np.arccos(h[2]/hn); Om = np.arctan2(n[1], n[0]) % (2*np.pi)
    w = np.arccos(np.dot(n, ev)/(np.linalg.norm(n)*e)); w = 2*np.pi - w if ev[2] < 0 else w
    nu = np.arccos(np.dot(ev, r)/(e*np.linalg.norm(r))); nu = 2*np.pi - nu if np.dot(r, v) < 0 else nu
    return a, e, i, Om, w, nu

7. (Satélite terrestre) Un satélite tiene, respecto al ecuador ICRF, $\mathbf r = (-6045, -3490, 2500)$ km y $\mathbf v = (-3.457, 6.618, 2.533)$ km/s con $\mu = 398\,600$ km³/s² (Curtis, ejemplo 4.3). Calcula sus elementos y di si la órbita es prograda o retrógrada y si está sobre o bajo el ecuador en ese instante.

Autoevaluación

Lecturas

• Curtis (2020), cap. 4: la presentación más directa de los algoritmos 4.2 (estado → elementos) y 4.5 (elementos → estado), con ejemplos numéricos 4.3 y 4.7 que reproducen nuestros tests.
• Vallado (2013), §2.4: algoritmos 9 y 10 con todos los casos degenerados, y §2.4.3 sobre elementos equinocciales y otros conjuntos no singulares.
• Murray & Dermott (1999), §2.8: la conversión en el contexto de la dinámica del Sistema Solar, con la notación $(\varpi, \lambda)$.
• Broucke & Cefola (1972), DOI 10.1007/BF01228432: el artículo original de los elementos equinocciales.
• JPL Horizons, documentación de las tablas ELEMENTS (convenciones de $\Omega$, $\omega$, $T_p$ y «Keplerian GM»).