¿Cómo calcular programáticamente elementos orbitales usando vectores de posición/velocidad?

Dados vectores de posición y velocidad inerciales (ECI) centrados en la Tierra$\vec{r}$y$\vec{v}$, puedes resolver directamente los elementos orbitales clásicos$(a,e,i,\Omega,\omega,\nu)$de la siguiente manera (algoritmos primero, seguidos de pseudocódigo en la parte inferior):

Primero, resuelve el momento angular

$$\vec{h}=\vec{r} \times \vec{v}$$

entonces el vector de nodo

$$\hat{n}=\hat{K} \times \vec{h}$$ que se utilizará más adelante.

El vector de excentricidad es entonces

$$\vec{e} = \frac{(v^2-\mu /r)\vec{r}-(\vec{r} \cdot \vec{v})\vec{v}}{\mu}$$

y$e=|\vec{e}|$.

La energía mecánica específica es

$$E=\frac{v^2}{2}-\frac{\mu}{r}$$

Si$e\neq 1$, después

$$a = -\frac{\mu}{2E}$$ $$p=a(1-e^2)$$ De lo contrario, $$p=\frac{h^2}{\mu}$$ $$a=\infty$$

Ahora,

$$i=\cos^{-1}{\frac{h_K}{h}}$$ $$\Omega=\cos^{-1}{\frac{n_I}{n}}$$ $$\omega=\cos^{-1}{\frac{\vec{n}\cdot\vec{e}}{ne}}$$ $$\nu=\cos^{-1}{\frac{\vec{e}\cdot\vec{r}}{er}}$$

Y deberá realizar las siguientes comprobaciones: Si$n_J<0$, después$\Omega=360^{\circ}-\Omega$,

Si$e_K<0$, después$\omega=360^{\circ}-\omega$, y

Si$\vec{r}\cdot\vec{v}<0$, después$\nu=360^{\circ}-\nu$.

Tenga en cuenta que se encontrará con problemas (singularidades) para ciertos casos: órbitas circulares ($e\approx 0$), y órbitas ecuatoriales ($i\approx 0$), particularmente. En estos casos, normalmente introduce una variable nueva y menos problemática, como la longitud media o la longitud real del perigeo.

h=cross(r,v)
nhat=cross([0 0 1],h)

evec = ((mag(v)^2-mu/mag(r))*r-dot(r,v)*v)/mu
e = mag(evec)

energy = mag(v)^2/2-mu/mag(r)

if abs(e-1.0)>eps
   a = -mu/(2*energy)
   p = a*(1-e^2)
else
   p = mag(h)^2/mu
   a = inf

i = acos(h(3)/mag(h))

Omega = acos(n(1)/mag(n))

if n(2)<0
   Omega = 360-Omega

argp = acos(dot(n,evec)/(mag(n)*e))

if e(3)<0
   argp = 360-argp

nu = acos(dot(evec,r)/(e*mag(r))

if dot(r,v)<0
   nu = 360 - nu

Nota : esto se deriva del método establecido en Fundamentos de astrodinámica y aplicaciones , de Vallado, 2007.

La primera clave para resolver esto es obtener el sistema de coordenadas correcto. Hay dos sistemas de coordenadas de uso común para tales cosas. Son los marcos Earth Centered Earth Fixed (ECEF) y Earth Centered Interial (ECI) . A medianoche, estos dos se alinean exactamente, pero divergen en otros momentos, según la rotación de la Tierra. ECEF funciona mejor para las cosas en la Tierra (si no se está moviendo, debe tener una velocidad de 0. ECEF tiene esto en cuenta, la velocidad de ECI hará que se mueva con la rotación de la Tierra), ECI funciona mejor para las cosas en órbita (en órbita). a los objetos no les importa la rotación de la Tierra, al menos, a la física no le importa). ¡Asegúrate de que los sistemas de coordenadas sean correctos!

Bien, entonces tienes una posición y velocidad en coordenadas ECI, ¿qué haces? Hay un documento excelente que describe todo el proceso, del cual copiaré las fórmulas finales aquí. También hay algunas buenas fuentes aquí , aquí y aquí . Recomiendo encarecidamente leerlos detenidamente. La incertidumbre es mucho más difícil, así que supongamos que tiene un conocimiento perfecto de la velocidad y la posición. Específicamente, los 6 Elementos Keplarianos clásicos son excentricidad (e), inclinación (i), ascensión recta del nodo ascendente ($\Omega$), argumento del perigeo ($\omega$), semieje mayor (a), y tiempo de paso del perigeo ($T_O$).

Debo mencionar que estoy siguiendo principalmente el método de determinación orbital de Laplace , existe una metodología competidora conocida como el método de Gauss. Pero finalmente esto se redujo a descifrar el código de Matlab .

Semieje mayor

$W_s = \frac{1}{2}*v^2s - \text{mus}./r;$

$a = -mus/2./W_s$; %semieje mayor

Excentricidad

L = [rs(2,:).*vs(3,:) - rs(3,:).*vs(2,:);...
     rs(3,:).*vs(1,:) - rs(1,:).*vs(3,:);...
     rs(1,:).*vs(2,:) - rs(2,:).*vs(1,:)]; %angular momentum

$p = \sum{L^2}./mus;$% recto semilato

$e = \sqrt{1 - p/a}; %eccentricity$

Inclinación

$I = atan(\frac{\sqrt{L(1,:)^2 + L(2,:)^2}}{L(3,:)});$

Argumentos del pericentro

$\omega = atan2(\frac{(vs(1,:).*L(2,:) - vs(2,:).*L(1,:))./mus - rs(3,:)./r)./(e.*sin(I))}{((\sqrt{L2s}.*vs(3,:))./mus - (L(1,:).*rs(2,:) - L(2,:).*rs(1,:))./(\sqrt{L2s}.*r))./(e.*sin(I)))}$

Longitud del nodo ascendente

$\Omega = atan2(-L(2,:),L(1,:));$

Tiempo de paso del perigeo:

$T_0 = -(E - e.*sin(E))./\sqrt{mus.*a.^-3}$

OrbitalPy tiene una elements_from_state_vectorfunción útil que hace exactamente eso:

https://github.com/RazerM/orbital/blob/0.7.0/orbital/utilities.py#L252

Puede verificar que las matemáticas sean las mismas que en la respuesta del usuario 29.