Wie bestimmt man eine Umlaufbahn eines Satelliten für eine Kollisionserkennung?

Wie kann ich ein Paar Anfangsbedingungen für einen zweiten Satelliten bestimmen, der zu einer Kollision führen würde?

Die Ausbreitung einer Umlaufbahn ist das Lösen eines Satzes von Differentialgleichungen durch numerisches Integrieren über die Zeit. Die Haupteffekte sind Schwerkraft und Luftwiderstand. Solange das Raumfahrzeug einen konstanten Luftwiderstandsbeiwert hat (es nicht taumelt), können diese Gleichungen zeitlich genauso einfach rückwärts ausgeführt werden wie sie vorwärts ausgeführt werden können. Wie unten noch einmal erwähnt, vergessen Sie nicht, das Vorzeichen der Widerstandskraft zu ändern, wenn Sie in der Zeit rückwärts laufen, und wenn Sie andere sphärische Harmonische als die zonalen Harmonischen (J2, J3, J4, ...) verwenden, tun Sie dies nicht Vergessen Sie, die Erde auch rückwärts zu drehen !

Das Verfahren wäre, die Bewegung des ersten Raumfahrzeugs aus Ihrem Anfangszustand zu integrieren$\mathbf{v_0}, \mathbf{x_0}$für einen Tag. Rufen Sie die endgültige Position und Geschwindigkeit auf$\mathbf{x_1}, \mathbf{v_1}$.

Überprüfen Sie Ihre Integration noch einmal, indem Sie das Raumschiff bei starten$\mathbf{x_1}, -\mathbf{v_1}$(die endgültige Position, aber mit der entgegengesetzten Geschwindigkeit) und einen Tag lang propagieren, um sicherzustellen, dass sie ankommt$\mathbf{x_0}, -\mathbf{v_0}$.

Starten Sie nun ein zweites Raumfahrzeug am Kollisionspunkt$\mathbf{x_1}$aber wählen Sie eine andere Geschwindigkeit als $-\mathbf{v_1}$, und propagieren Sie es für einen Tag. Rufen Sie das Ergebnis auf$\mathbf{x_2}, -\mathbf{v_2}$(beachten Sie das Minuszeichen vor der Geschwindigkeit).

Wenn Sie eine andere Geschwindigkeit wählen, vergessen Sie nicht, sicherzustellen, dass sie nicht auch die Erdatmosphäre schneidet. Sobald Sie Ihre anfängliche Umlaufbahn festgelegt haben (siehe unten), besteht eine gute Möglichkeit, dies zu tun, darin, die Radialgeschwindigkeit gleich zu halten und einfach die Tangentialgeschwindigkeit um den Radiusvektor zu drehen. Dies führt zu einer nahezu identisch geformten Umlaufbahn, jedoch mit einer anderen Neigung.

Du bist fertig!

Wenn Sie das erste Raumschiff an starten$\mathbf{v_0}, \mathbf{x_0}$und das zweite Raumschiff bei$\mathbf{x_2}, \mathbf{v_2}$und beide für einen Tag integrieren, sollen beide ankommen$\mathbf{x_1}, \mathbf{v_1}$innerhalb der Genauigkeitsgrenzen Ihres RK-Methodenintegrators und anderer numerischer Artefakte.

Dies sollte für eine zentrale Gravitationskraft mit sphärischen Harmonischen (J2 usw.) funktionieren, aber ohne Sonne oder Mond, plus einem Widerstandsterm, der von der Geschwindigkeit im Quadrat und der Höhe abhängig ist, wobei ein atmosphärisches Modell verwendet wird. (Siehe Atmosphärischer Widerstandseffekt für einige Ideen.) Wenn Sie jedoch Widerstand verwenden, müssen Sie das Vorzeichen Ihrer Widerstandskraft ändern, wenn Sie sich zeitlich rückwärts ausbreiten. Denken Sie daran, dass dies nur für einen konstanten Luftwiderstandsbeiwert und keine pathologischen Probleme wie den atmosphärischen Wiedereintritt funktioniert.

Wenn Sie neben den zonalen Harmonischen (J2, J3, J4, ...) auch andere sphärische Harmonische für die Schwerkraft verwenden, müssen Sie die Erde auch rückwärts drehen .

Sie müssen auch darauf achten, wohin jede Umlaufbahn führt. Ihre Anfangsbedingungen scheinen eine Periapsis von etwa 3660 km zu haben, was einer Höhe von etwa -2710 km entspricht!! die unter der Erdoberfläche liegt.

Unten habe ich mit einem einfachen Python-Skript nur die zentrale Kraft und kein Ziehen verwendet. Sie können sehen, dass das Raumschiff unter die Erdoberfläche eintaucht. Sie können J2 und höhere Harmonische hinzufügen und ziehen, nachdem Sie damit begonnen haben. (Siehe diese Antwort für einige Ideen zum Hinzufügen von J2 und relativistischer Korrektur, aber Sie benötigen ein formelleres Koordinatensystem, um es richtig zu machen.)

Ich habe eine spezifische Energie von -5,42E+06 Joule/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()