¿Cómo determinar una órbita de un satélite para una detección de colisión?

¿Cómo puedo determinar un par de condiciones iniciales para un segundo satélite que conduciría a una colisión?

Propagar una órbita es resolver un conjunto de ecuaciones diferenciales integrando numéricamente en el tiempo. Los principales efectos son la gravedad y la resistencia. Siempre que la nave espacial tenga un coeficiente de arrastre constante (no esté dando vueltas), estas ecuaciones se pueden ejecutar hacia atrás en el tiempo con la misma facilidad con la que se pueden ejecutar hacia adelante. Como se indica nuevamente a continuación, cuando retroceda en el tiempo , no olvide cambiar el signo de la fuerza de arrastre y si está utilizando armónicos esféricos distintos de los armónicos zonales (J2, J3, J4, ...), no lo haga . ¡ Olvídate de rotar la Tierra hacia atrás también!

El procedimiento sería integrar el movimiento de la primera nave espacial desde su estado inicial$\mathbf{v_0}, \mathbf{x_0}$por un día. Llame a la posición final y la velocidad$\mathbf{x_1}, \mathbf{v_1}$.

Vuelva a verificar su integración iniciando la nave espacial en$\mathbf{x_1}, -\mathbf{v_1}$(la posición final, pero con la velocidad opuesta) y propagándose durante un día para asegurarse de que llega a$\mathbf{x_0}, -\mathbf{v_0}$.

Ahora inicie una segunda nave espacial en el punto de colisión.$\mathbf{x_1}$pero elija una velocidad diferente a $-\mathbf{v_1}$, y propagarlo por un día. llama al resultado$\mathbf{x_2}, -\mathbf{v_2}$(Observe el signo negativo delante de la velocidad).

Cuando elijas una velocidad diferente, no olvides asegurarte de que no intersecte también la atmósfera terrestre. Una vez que fije su órbita inicial (vea a continuación), una buena manera de hacerlo es mantener la velocidad radial igual y simplemente rotar la velocidad tangencial alrededor del vector de radio. Eso dará como resultado una órbita de forma casi idéntica pero con una inclinación diferente.

¡Estás listo!

Si inicias la primera nave espacial en$\mathbf{v_0}, \mathbf{x_0}$y la segunda nave espacial en$\mathbf{x_2}, \mathbf{v_2}$e integrar ambos por un día, ambos deben llegar a$\mathbf{x_1}, \mathbf{v_1}$dentro de los límites de precisión de su integrador de métodos RK y otros artefactos numéricos.

Esto debería funcionar para una fuerza gravitacional central con armónicos esféricos (J2, etc.) pero sin Sol ni Luna, además de un término de arrastre que depende de la velocidad al cuadrado y de la altitud utilizando algún modelo atmosférico. (Consulte Efecto de arrastre atmosférico para obtener algunas ideas). Sin embargo, si usa el arrastre, debe cambiar el signo de su fuerza de arrastre cuando se propaga hacia atrás en el tiempo. Recordatorio de que esto solo funciona para un coeficiente de arrastre constante y sin problemas patológicos como el reingreso atmosférico.

Además, si está utilizando armónicos esféricos para la gravedad además de los armónicos zonales (J2, J3, J4, ...), también deberá girar la Tierra hacia atrás.

También debe prestar atención a dónde va cada órbita. ¡Sus condiciones iniciales parecen tener un periápside de unos 3660 km, que es una altitud de unos -2710 km! que está debajo de la superficie de la Tierra.

A continuación, he integrado usando solo la fuerza central y sin arrastrar con un simple script de Python. Puedes ver que la nave espacial se sumerge debajo de la superficie de la Tierra. Puede agregar J2 y armónicos superiores y arrastrar después de que esto funcione. (Consulte esta respuesta para obtener algunas ideas sobre cómo agregar J2 y la corrección relativista, pero necesitará un sistema de coordenadas más formal para hacerlo correctamente).

Obtuve una energía específica de -5.42E+06 Joules/kg.

def deriv(X, t):

    x, v = X.reshape(2, -1)
    a    = -GMe * x * ((x**2).sum())**-1.5
    return np.hstack((v, a))

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

# https://space.stackexchange.com/questions/30125/how-to-determine-an-orbit-of-a-satellite-for-a-collision-detection

GMe = 3.986E+14
X0  = np.array([742453.224, -6345418.953, -3394102.502,
                -5142.381, -4487.449, -7264.009])

r0 = np.sqrt((X0[:3]**2).sum())
v0 = np.sqrt((X0[3:]**2).sum())

KE, PE = 0.5 * v0**2,  -GMe/r0
E      = KE + PE
print "E (Joules/kg): ", E

minutes = np.arange(2000.)
time    = 60. * minutes

answer, info = ODEint(deriv, X0, time, full_output=True)

xx, vv = answer.T.reshape(2, 3, -1)
rr     = np.sqrt((xx**2).sum(axis=0))

if True:
    plt.figure()
    plt.subplot(2, 1, 1)
    for thing in xx:
        plt.plot(minutes,  thing/1000.)
    plt.title('x, y, z (km) vs. time (minutes)')
    plt.subplot(2, 1, 2)
    plt.title('radius (km) vs. time (minutes)')
    plt.plot(minutes, rr/1000.)
    plt.plot(minutes, 6378.137 * np.ones_like(rr), '-k')
    plt.show()