2 ecuaciones de movimiento de cuerpo en 3D

$state=[rx,ry,rz,vx,vy,vz]$

$\dot{state}=[vx,vy,vz,ax,ay,az]$

$a=\dfrac{\mu}{r^2}$(Ley de gravitación universal de Newton / aceleración de dos cuerpos en forma escalar)

Los errores de división en cero se evitan utilizando la ecuación vectorial para la aceleración, ya que cualquiera de los componentes de posición puede ser igual a cero, pero la norma será un número positivo (suponiendo que la norma en sí misma no sea 0)

$\vec{a}=\dfrac{-\mu}{|\vec{r}|^3}\vec{r}$(Ley de gravitación universal de Newton / aceleración de dos cuerpos en forma vectorial. Tenga en cuenta que se agregó el signo negativo ya que la aceleración gravitatoria está en la dirección opuesta al vector de posición.

Aquí está la versión de Python de esta ODE:

def two_body_ode( t, state, mu = pd.earth[ 'mu' ] ):
    r = state[ :3 ]
    a = -mu * r / np.linalg.norm( r ) ** 3

    return np.array( [
        state[ 3 ], state[ 4 ], state[ 5 ],
        a[ 0 ], a[ 1 ], a[ 2 ] ] )

Y una trama para divertirse:

Esta respuesta a ¿Cómo resolver el problema de dos cuerpos en el marco ECI a través de la integración numérica? dice:

En lugar de una única ecuación diferencial de segundo orden

$$\ddot{\mathbf{r}} = -\frac{\mu}{r^3}\mathbf{r}$$

Podemos resolver el siguiente par de ecuaciones diferenciales de primer orden en paralelo

$$\dot{\mathbf{v}} = -\frac{\mu}{r^3}\mathbf{r}$$ $$\dot{\mathbf{r}} = \mathbf{v}$$

usando varios métodos simples que incluyen bibliotecas estándar o implementaciones propias de Runge-Kutta, incluido mi RK4/5 simple favorito con tamaño de paso variable .

[...]

A partir de esta respuesta, puede ver una implementación 2D utilizando no solo el monopolo$1/r^2$término de gravedad pero el cuadrupolo adicional$J_2$término para la forma achatada y el campo de la Tierra. Para obtener más información, consulte esta respuesta a Problemas para derivar componentes rectangulares de aceleración del satélite en órbita alrededor de la Tierra con consideración J2.

Aquí$\mu$con el parámetro gravitacional estándar que sería el producto$GM$si conociéramos cada parámetro con precisión. Resulta que para los cuerpos del sistema solar podemos determinar su producto con mucha más precisión que cualquiera de ellos solo. Si le gusta trabajar sin unidades y su cuerpo en órbita tiene una masa despreciable, configure$\mu = 1$y tu periodo sera$2 \pi$.

Entonces la pregunta es:

¿Cómo se escriben las ecuaciones de movimiento de 2 cuerpos en 3D como un sistema de ecuaciones diferenciales en la forma Xdot = f(X)?

Escrito y publicado,$\dot{\mathbf{r}} = \mathbf{v}$y$\dot{\mathbf{v}} = -\frac{\mu}{r^3}\mathbf{r}$convertirse en:

$$\dot{x} = v_x$$ $$\dot{y} = v_y$$ $$\dot{z} = v_z$$ $$\dot{v_x} = -\frac{\mu}{r^3}x$$ $$\dot{v_y} = -\frac{\mu}{r^3}y$$ $$\dot{v_z} = -\frac{\mu}{r^3}z$$ $$r = \left(x^2 + y^2 + z^2 \right)^{1/2}$$

Si:

$$\mathbf{X} = \begin{bmatrix} x \\ y \\ z \\ v_x \\ v_y \\ v_z \\ \end{bmatrix}$$

entonces

$$\mathbf{\dot{X}} = \begin{bmatrix} v_x \\ v_y \\ v_z \\ -\frac{x \ \mu}{(x^2 + y^2 + z^2)^{3/2}} \\ -\frac{y \ \mu}{(x^2 + y^2 + z^2)^{3/2}} \\ -\frac{z \ \mu}{(x^2 + y^2 + z^2)^{3/2}} \\ \end{bmatrix}$$