Determinación de la posición orbital en un punto futuro en el tiempo

Planteamiento del problema

El problema que quieres resolver se llama problema de Kepler . En su formulación del problema, está comenzando con los vectores de estado orbital cartesianos (también llamados elementos cartesianos ): es decir, la posición inicial y la velocidad.

Como ha descubierto, la única forma de propagar los elementos cartesianos en el tiempo es mediante integración numérica. Esto funciona bien, pero puede ser lento si desea una alta precisión, y hay algunos problemas numéricos (los errores causados ​​por el redondeo [se acumulan lentamente y muchos integradores provocan una deriva de energía ). Puede solucionar algunos de estos problemas mediante el uso de un integrador de orden superior ( Runge-Kutta es uno popular) que le permite dar pasos más grandes para el mismo nivel de precisión u obtener una mejor precisión para el mismo tamaño de paso. Sin embargo, esto es un poco excesivo para una simulación simple.

Si su simulación se puede tratar como un problema de dos cuerpos , entonces las cosas se simplifican dramáticamente. El problema de los dos cuerpos es una buena simplificación si los objetos de simulación se ven afectados principalmente por un solo objeto grande. Por ejemplo, la Tierra que viaja alrededor del Sol o una nave espacial que viaja en una órbita terrestre baja están bien modelados como un problema de dos cuerpos; sin embargo, una nave espacial que viaja de la Tierra a la Luna no lo es (más sobre esto más adelante).

Dado que está tratando de modelar las posiciones de los planetas con una precisión media, la reducción al problema de dos cuerpos debería funcionar para usted.

Definición de términos

La solución tradicional al problema de los dos cuerpos implica una forma diferente de representar la posición del cuerpo en órbita, llamados elementos orbitales keplerianos (también llamados simplemente elementos orbitales ). En lugar de especificar la posición y la velocidad, especifican seis parámetros diferentes de la órbita (si solo desea acceder al código, puede omitir esta parte):

1. Semieje mayor,$a$: La mitad del diámetro máximo de la órbita elíptica, ( = radio del ciclo si la órbita es circular). La energía y el período de la órbita dependen únicamente de$a$. El recto semi-latus $\ell$, el "ancho" de la órbita, puede ser una mejor opción para órbitas cercanas a las parabólicas (como los asteroides) o que cambian de elípticas a hiperbólicas (como las naves espaciales interplanetarias). Los dos están relacionados por$\ell=a(1-e^2)$.

2. Excentricidad,$e$: La "puntalidad* de la órbita. Va desde$e=0$para una órbita perfectamente circular, a$e=1$para una órbita parabólica, a$e>1$para órbitas hiperbólicas. Mercurio es el planeta más excéntrico con$e \approx 0.2$. Las naves espaciales en órbita terrestre suelen tener$e<0.01$.

Aparte de$e$y$a$podemos determinar los puntos más lejanos y más cercanos de la órbita, el apoapsis y el periapsis ( absides juntos ):

$$r_a = a(1+e) \\ r_p = a(1-e)$$ _{La denominación de estos puntos es un poco divertida: apoapsis y periapsis son términos genéricos, pero las órbitas alrededor de cuerpos particulares tienen términos específicos: una nave espacial alrededor de la Tierra tiene apogeo y perigeo , mientras que la Tierra (en órbita alrededor del Sol) tiene un afelio y perihelio .}

Los dos parámetros$a$y$e$son suficientes para determinar la forma de la órbita. Los siguientes tres parámetros definen la orientación de la órbita relativa a un sistema de coordenadas que consiste en un plano de referencia y una dirección de referencia (paralela al plano).

Para casi todas las órbitas del sistema solar, el sistema de coordenadas utilizado es el sistema de coordenadas de la eclíptica . El plano de referencia es el plano de la eclíptica , el plano de la órbita de la Tierra alrededor del Sol. La dirección de referencia es el punto del equinoccio de primavera , la dirección de la Tierra al Sol en el momento del equinoccio de primavera. Dado que ambas referencias se desplazan lentamente con el tiempo, debemos especificar un momento particular en el que se definen estas referencias, llamado época . El más común es J2000 , mediodía del 1 de enero de 2000 (UTC).

Las órbitas centradas en la Tierra suelen utilizar el sistema de coordenadas ecuatoriales , cuyo plano de referencia es el ecuador de la Tierra. La situación con la época es un poco complicada, así que no entraré en eso aquí.

Los siguientes parámetros ubican la órbita con respecto a la órbita terrestre:

3. Inclinación,$i$: el ángulo entre el plano de la órbita y el plano de referencia. La inclinación entre 90 y 180 grados se refiere a una órbita retrógrada , que orbita "hacia atrás" desde la dirección habitual.

4. Longitud del nodo ascendente,$\Omega$: el nodo ascendente es donde la órbita cruza desde debajo del plano de referencia hacia arriba. (Está en la intersección entre el plano orbital y el plano de referencia)$\Omega$es el ángulo entre este punto y la dirección de referencia, medido en sentido antihorario.

5. Argumento del periapsis,$\omega$: el ángulo entre el nodo ascendente y el periapsis (el punto más bajo de la órbita). Para órbitas con muy poca inclinación donde la ubicación del nodo ascendente es difícil de determinar (ya que es la intersección entre dos planos casi paralelos), en su lugar usamos la longitud del periapsis $\varpi = \Omega + \omega$.

El sexto parámetro define la posición del objeto en su órbita. Hay un par de opciones diferentes, pero la más común es:

• Anomalía media,$M$: un ángulo "imaginario" que es cero en el periapsis y aumenta a un ritmo constante de 360 ​​grados por órbita.

La tasa a la que$M$cambios se llama el movimiento medio ,$n$, igual a$2\pi/T$. Normalmente tienes una medida de$M$en una epoca particular$t_0$, llamado (como era de esperar) la anomalía media en epoch ,$M_0$.

Al igual que el argumento del periapsis, para órbitas de baja inclinación usamos un valor relacionado, la longitud media ,$L=\varpi + M$.

El ángulo real entre el cuerpo en órbita y el periapsis se denomina anomalía verdadera .$\nu$. Este es el ángulo que necesitamos para calcular la posición del cuerpo. Desafortunadamente no hay manera de calcular directamente$\nu$de$M$. En su lugar, primero resolvemos la anomalía excéntrica $E$:

$$M = E - e \sin E$$

Esto se llama ecuación de Kepler y no se puede resolver analíticamente. Una vez que tengamos$E$sin embargo, hay una expresión relativamente simple para$\nu$.

Cálculo de la posición a partir de elementos orbitales

Realizaremos este cálculo en tres pasos: primero, resolveremos la ecuación de Kepler. Segundo, calcularemos la posición 2d del cuerpo en el plano orbital. Por último, rotaremos nuestra posición 2D en coordenadas 3D. Daré algo de "pseudocódigo" en Javascript para la mayoría de estas tareas.

Asumiré que está utilizando un conjunto de elementos como estos del sitio web de JPL . Estos utilizan$L$y$\varpi$en vez de$M$y$\omega$. La tabla da dos valores para cada uno de los elementos; la segunda es la derivada del tiempo. Si usa los valores de esta tabla, también debe usar los derivados.

Calcular el tiempo$t$en siglos desde J2000:

// month is zero-indexed, so 0 is January
var tMillisFromJ2000 = Date.now() - Date.UTC(2000, 0, 1, 12, 0, 0);
var tCenturiesFromJ2000 = tMillisFromJ2000 / (1000*60*60*24*365.25*100);

Ahora calculamos los valores actuales de cada uno de los parámetros orbitales. Por ejemplo, el semieje mayor de la Tierra, utilizando los valores de la Tabla 1 (válidos desde 1800 hasta 2500):

// a0 = 1.00000261; adot = 0.00000562
var a = a0 + adot * tCenturiesFromJ2000;

(Tenga en cuenta que los valores en realidad se dan para "EM Barycenter ", el centro de masa del sistema Tierra-Luna. La Tierra está a unos 4600 kilómetros del baricentro en dirección opuesta a la Luna. Si desea corregir esto imprecisión, también necesitará simular el movimiento de la Luna, pero eso probablemente sea excesivo).

La Tabla 2a proporciona elementos que son precisos desde el 3000 a. C. hasta el 3000 d. C.; sin embargo, si usa los elementos de la tabla 2a, debe complementarlos con correcciones a$L$de la Tabla 2b! Por ejemplo, aquí está calculando la longitud de Saturno:

// L0 = 34.33479152; Ldot = 3034.90371757
// b = -0.00012452          
// c =  0.06064060
// s = -0.35635438
// f = 38.35125000
var L = L0 + Ldot * tCenturiesFromJ2000
           + b * Math.pow(tCenturiesFromJ2000, 2)
           + c * Math.cos(f * tCenturiesFromJ2000)
           + s * Math.sin(f * tCenturiesFromJ2000);

No necesitamos calcular explícitamente el movimiento medio y sumarlo a$L$, ya que ambas tablas lo incluyen en$\dot L$.

Ahora estamos listos para calcular$M$y$\omega$( w):

var M = L - p \\ p is the longitude of periapsis
var w = p - W \\ W is the longitude of the ascending node

En el paso 2: necesitamos resolver la ecuación de Kepler:

$$M = E - e \sin E$$

Podemos resolver esto numéricamente usando el método de Newton . Resolver la ecuación de Kepler es equivalente a encontrar las raíces de$f(E) = E - e \sin E - M$. Dado$E_i$, una estimación de$E$, podemos usar el método de Newton para encontrar una mejor estimación:

$$E_{i+1} = E_i - f(E_i) / f'(E_i) \\ f'(E) = 1 - e \cos E$$

Dado que la parte no lineal$e \sin E$es muy pequeño, podemos empezar con la estimación$E=M$. Nuestro código se parece a esto:

E = M;
while(true) {
  var dE = (E - e * Math.sin(E) - M)/(1 - e * Math.cos(E));
  E -= dE;
  if( Math.abs(dE) < 1e-6 ) break;
}

Ahora hay dos formas de calcular la posición a partir de la anomalía excéntrica. Primero podemos calcular la anomalía y el radio verdaderos (la posición del objeto en coordenadas polares) y luego convertirlos a coordenadas rectangulares; sin embargo, si aplicamos un poco de geometría, podemos calcular las coordenadas directamente desde$E$:

var P = a * (Math.cos(E) - e);
var Q = a * Math.sin(E) * Math.sqrt(1 - Math.pow(e, 2));

( Py Qformar un sistema de coordenadas 2d en el plano de la órbita, +Papuntando hacia el periapsis).

Finalmente, podemos rotar estas coordenadas en el sistema de coordenadas 3D completo:

// rotate by argument of periapsis
var x = Math.cos(w) * P - Math.sin(w) * Q;
var y = Math.sin(w) * P + Math.cos(w) * Q;
// rotate by inclination
var z = Math.sin(i) * y;
    y = Math.cos(i) * y;
// rotate by longitude of ascending node
var xtemp = x;
x = Math.cos(W) * xtemp - Math.sin(W) * y;
y = Math.sin(W) * xtemp + Math.cos(W) * y;

( x, y, y zestarán en unidades de AU.)

¡Y tu estas listo!

---

Algunos consejos:

• Si también desea calcular la velocidad, puede hacerlo al mismo tiempo que calcula$P$y$Q$, luego gírelo de la misma manera.

  $$\dot M = n = \dot L \\ \dot M = \dot E - e (\cos E) \dot E \\ \dot E = \dot M / (1 - e \cos E) \\ \dot P = -a (\sin E) \dot E \qquad \dot Q = a (\cos E) \dot E \sqrt{1 - e^2}$$ Tenga en cuenta que no incluyo ninguno de los derivados (excepto$\dot L$) en este cálculo, ya que no afectan mucho el resultado. Podrías codificar esto como:

   var vP = - a * Math.sin(E) * Ldot / (1 - e * Math.cos(E));
   var vQ = a * Math.cos(E) * Math.sqrt(1 - e*e) * Ldot / (1 - e * Math.cos(E));
  

  Tenga en cuenta que las velocidades estarán en AU por siglo.

• Si está actualizando las posiciones con mucha frecuencia, puede usar el valor anterior de$E$para sembrar el método de Newton, y hacer un número fijo de iteraciones (probablemente solo una sería suficiente). Sin embargo, tenga en cuenta que debe mantener ese valor de$E$local para cada objeto!

• También puede usar un número fijo de iteraciones para la solución inicial. Incluso para$e=0.2$, después de tres iteraciones el error en$E$se trata solo de$10^{-13}$, y después de cuatro iteraciones el error es más pequeño que el error de redondeo de un IEEE se duplica hasta$e=0.42$.

---

Si desea obtener más información, puede buscar en línea, pero si está realmente interesado, debe leer un texto introductorio sobre mecánica orbital. Yo personalmente recomiendo Fundamentals of Astrodynamics de Bate, Mueller y White (pdf) . Mi papá usó este libro cuando estaba en la universidad y me pareció más fácil de leer que mi libro de texto de la universidad. Le interesaría el Capítulo 4, Posición y velocidad como función del tiempo.

Dado que es solo un juego, ¿estarías contento con las órbitas circulares y las órbitas de los planetas solo afectadas por el cuerpo central? En ese caso, la propagación es bastante simple. En el plano de la órbita con el cuerpo central en (0,0), la posición en función del tiempo es:

$$x(t)=a\cos\left({2\pi \left(t-t_0\right)\over T}\right)$$

$$y(t)=a\sin\left({2\pi \left(t-t_0\right)\over T}\right)$$

dónde$a$es el eje semi-mayor, o realmente solo el radio de la órbita en este caso,$T$es el período de la órbita, y$t_0$determina la fase de la órbita, donde en$t=t_0$, el planeta está en el eje x en el lado positivo.

Para que las órbitas de los diferentes planetas sean consistentes entre sí, solo necesita definir el$GM$del cuerpo central, que llamaremos$\mu$. Entonces para cualquier radio de órbita$a$, el período de la órbita está relacionado con$a$por:

$$T=2\pi\sqrt{a^3\over\mu}$$

Si bien ya ha habido una respuesta aceptada de alta calidad durante años, aquí hay algunos antecedentes adicionales, algunos recursos particularmente útiles y consejos adicionales para la propagación en órbita por primera vez.

---

Si no está haciendo física de N-cuerpos, por lo que los planetas no interactúan, entonces puede usar soluciones analíticas para el problema de Kepler. Eventualmente, te darás cuenta de que también necesitas resolver órbitas hiperbólicas en algún momento. Eso lo llevará a formulaciones de variables universales para resolver el problema de Kepler.

Las mejores soluciones a eso probablemente serán el método de Goodyear:

W. Goodyear, "Solución de forma cerrada completamente general para coordenadas y derivadas parciales del problema de dos cuerpos", The Astronomical Journal, vol. 70, No. 3, 1965, pp. 189–192 (o el documento NTRS TD de la NASA sobre el mismo material )

Método de Shepherd:

Shepperd, SW Celestial Mechanics (1985) 35: 129. https://doi.org/10.1007/BF01227666

O Danby-Stumpff:

Danby, JMA Mecánica Celestial (1987) 40: 303. https://doi.org/10.1007/BF01235847

Aquí hay un código de MATLAB que podría ser útil (y mucho más accesible), aunque los fragmentos de código aleatorios en matlabcentral están lejos de garantizarse que estén libres de errores y parece que este código puede carecer de una normalización útil de sus entradas (generalmente vas querer normalizar a la escala de su problema para que haga matemáticas en unidades donde r0-bar = 1.0 y mu-bar = 1.0 y donde v-bar = 1 es la velocidad en una órbita circular en r0 o algo así) .

Si va a hacer una integración de N-cuerpos del movimiento planetario, entonces creo que tendrá que usar la integración numérica. Runge-Kutta violará la conservación de la energía, por lo que es probable que desee utilizar la integración simpléctica . El integrador simpléctico de cuarto orden en ese artículo no es tan difícil de codificar, aunque eso lo deja con la dificultad de adivinar el paso de tiempo correcto (nuevamente, la normalización ayuda porque una órbita planetaria circular y un LEO circular son el mismo problema solo que con diferentes escalas de distancia ) y con interpolación de los puntos interiores (y debe tener cuidado con el fenómeno de Runge , pero no he luchado con eso, así que no sé qué enfoque tomar allí).

Si va a utilizar Runge-Kutta, entonces Dormand-Prince con lado de paso dinámico y su interpolador de tercer orden será muy conveniente, y es lo que utiliza Matlab en su solucionador ode45.

Probablemente recomendaría comenzar con la implementación de runge-kutta más simple basada en la facilidad de codificación, pero si está haciendo runge-kutta en cada marca de física para avanzar un paso, eso es bastante brutal y los errores eventualmente se sumarán, pero podrías hacer un prototipo de esa manera. En algún momento, querrá ir a un sistema en el que resuelva el problema durante muchos pasos de tiempo en el futuro, y luego use una función de interpolación para elegir la solución en pasos de tiempo intermedios (que es el punto de mi mención Dormand- Prince y su función interpoladora).