Uso de la posición, la velocidad y la masa en 2D para determinar las ecuaciones de posición paramétricas de un cuerpo en órbita

Esta respuesta es un complemento de la respuesta de Chris White.

En primer lugar, no hay ecuaciones explícitas para la posición de un objeto siguiendo una órbita de Kepler en función del tiempo. Sin embargo, cuando se conocen las condiciones iniciales, se puede encontrar el camino que seguirá el objeto, así como la velocidad, aceleración, etc. en cualquier posición dada. Este camino puede ser descrito por la siguiente ecuación:

$$r=\frac{a(1-e^2)}{1+e\cos{\theta}}$$ dónde $r$es la distancia entre el objeto (pequeño) y el foco de la órbita (que es igual a la posición del objeto fijo grande o el centro de masa de los dos si el segundo objeto no está fijo), $a$es el semieje mayor, $e$es la excentricidad y $\theta$es la anomalía verdadera (el ángulo entre la posición del objeto, el foco y la posición del objeto en el acercamiento más cercano/periapsis).

El semieje mayor y la excentricidad se pueden deducir de las condiciones iniciales:

$$a=\frac{\mu r}{2\mu-rv^2}$$ $$e=\sqrt{1+\frac{rv_{\theta}^2}{\mu}\left(\frac{rv^2}{\mu}-2\right)}$$ dónde $\mu$es el parámetro gravitacional del gran objeto fijo ( $\mu=GM$), $v$es la velocidad del objeto y $v_r$es la componente tangencial de la velocidad.

Esta componente tangencial de la velocidad se puede encontrar usando algo de álgebra lineal (el equivalente resultante del producto vectorial en 2D):

$$v_{\theta}=\frac{xv_y-yv_x}{\sqrt{x^2+y^2}}$$

La posición del periapsis puede tener cualquier ángulo relativo al foco y al eje en el que representas las posiciones. Por lo tanto, deberá compensar/rotar la ruta que ha encontrado en un ángulo específico para que coincida con las condiciones iniciales.

Este ángulo de compensación de$\theta$, llamémoslo$\Delta\theta$, se puede calcular de la siguiente manera:

$$\Delta\theta={\rm sign}\left(v_{\theta}v_r\right)\cos^{-1}\left(\frac{a(1-e^2)-\sqrt{x^2+y^2}}{e\sqrt{x^2+y^2}}\right)-{\rm atan2}\left(y,x\right)$$ donde la funcion ${\rm sign}(x)$se refiere a si $x$es positivo o negativo $\left(\frac{x}{|x|}\right)$y $v_r$es la velocidad radial, que se puede encontrar mediante: $$v_r=\frac{xv_x+yv_y}{\sqrt{x^2+y^2}}$$ Este ángulo de compensación se elige de manera que la ruta resultante esté definida por: $$x_{path}=r(\theta)\cos\left(\theta+\Delta\theta\right)$$ $$y_{path}=-r(\theta)\sin\left(\theta+\Delta\theta\right)$$

Tengo una nota final sobre la elección de la gama de$\theta$, desde cuando$e\ge1$entonces el camino ya no será una elipse, sino una parábola o una hipérbola, que se extiende hasta el infinito. Para evitar intentar dibujar esto, puede limitar el rango de$\theta$. Por ejemplo, eligiendo un radio máximo,$r_{max}$. El rango de$\theta$será entonces:

$$\theta_{max}={\rm real}\left(\cos^{-1}\left(\frac{a(1-e^2)-r_{max}}{er_{max}}\right)\right)$$ $$\theta\in[-\theta_{max},\theta_{max}]$$

Parece que ya has hecho la parte difícil, que es evolucionar la posición del objeto en función del tiempo. Y además, la simulación parece estable en varias órbitas. (Pero eventualmente las cosas comienzan a salir mal; es posible que desee ver una respuesta que escribí a ¿ Cuál es la forma correcta de integrarse en simulaciones de astronomía? )

Entonces, entiendo que todo lo que realmente necesita es la órbita completa trazada. En ese caso, no hay necesidad de convertirlo en una función del tiempo real (eso es difícil). En su lugar, podemos recurrir a la Primera Ley de Kepler , que dice que la separación$r$entre los cuerpos obedece

$$r = \frac{r_\mathrm{max}(1-e)}{1+e\cos(\theta)}.$$ $r_\mathrm{max}$es la separación máxima, que creo que es la separación inicial en su simulación particular. la excentricidad $e$(fórmula aquí ) está dada por $$e^2 = 1 + \frac{2Er^4\dot{\theta}^2}{G^2M^2m},$$ dónde $E$es la energía orbital, $G$es la constante gravitacional , $M$es la masa que es tan pesada que esencialmente no se mueve, y $m$es la masa más ligera. Esto se puede reescribir más convenientemente $$e^2 = 1 + 2 \left(\frac{rv}{GM}\right)^2 \left(\frac{v^2}{2} - \frac{GM}{r}\right),$$ dónde $v$es la velocidad de la partícula en movimiento.

Tenga en cuenta que la convención es para$\theta$para medir el ángulo desde el acercamiento más cercano (es decir,$\theta = 0$es el negativo$y$-eje en su simulación). Si quieres$\theta$aumente con el tiempo, entonces la transformación (no estándar) a coordenadas cartesianas es

$$\begin{align} x & = -r \sin(\theta) \\ y & = -r \cos(\theta), \end{align}$$ asumiendo $r = 0$corresponde al objeto masivo. En cualquier caso, para trazar la órbita todo lo que necesita hacer es calcular la constante $e$en un punto de la órbita, muestra $\theta$con tantos puntos como quieras, calcula las separaciones $r(\theta)$, y convertir el $(r, \theta)$pares a $(x, y)$.

Las simulaciones orbitales se pueden manejar usando las siguientes relaciones:

$$\begin{eqnarray} \mathbf F&=&m\mathbf a=m\frac{d^2\mathbf x}{dt^2}\tag{a} \\ \mathbf v&=&\frac{d\mathbf x}{dt}\tag{b}\\ \mathbf a&=&\frac{d\mathbf v}{dt}\tag{c} \end{eqnarray}$$ La fuerza que actúa sobre dos cuerpos, masa $M$y $m$viene dada por la ley gravitatoria de Newton $$\mathbf F=G\frac{mM}{r^2}\hat{\mathbf r}=G\frac{mM}{r^3}\left(x\hat{\mathbf x}+y\hat{\mathbf y}\right)\tag{d}$$ con $r=\sqrt{x^2+y^2}$. Entonces el algoritmo queda así:

1. Calcular$\mathbf F$de$\mathbf x$
2. Determinar la nueva velocidad$\mathbf v$a través de$\mathbf v^{n+1}=\mathbf v^n+\mathbf a\,dt$
3. Determinar la nueva posición$\mathbf x$a través de$\mathbf x^{n+1}=\mathbf x^n+\mathbf v\,dt$
4. Vuelva al Paso 1.

Esto se puede hacer usando las ecuaciones exactas dadas arriba, pero no será numéricamente estable por mucho tiempo. Los métodos de orden superior como Runge-Kutta 3/4/5 mantendrán su estabilidad durante períodos de tiempo más prolongados, pero no serán perfectos.

Si esta no es la información que necesitas, házmelo saber e intentaré ayudarte en lo que pueda, pero me parece que solo necesitabas saber la Ecuación (d) en relación con las demás.