Mathesis Wiki

Die Koordinaten aller Objekte innerhalb des Sonnensystems werden numerisch für jeden Zeitschritt neu berechnet. Dafür wird für jeden Frame die Krafteinwirkung aller Objekte auf alle Objekte ermittelt und der Position zugerechnet.

Innerhalb dieses Systems befindet sich der Satellit. Um einen energieeffizienten Raumflug zu errechnen wird der Hohmann-Transfer genutzt, welcher die Geschwindigkeitsdifferenzen sowie die Zeitpunkte der Beschleunigungen ausgibt. Nachdem die Beschleunigungen an gegebenen Zeitpunkten durchgeführt wurden, wird der Satellit rein durch die Gravitation der Sonne sein Ziel erreichen.

Der Hohmann-Transfer ist eine energie- und zeiteffiziente Methode, von einer Umlaufbahn um einen Körper auf eine andere zu gelangen. In unserem Fall wollen wir unsere Rakete von einem Planeten, beziehungsweise dessen Umlaufbahn um die Sonne zu einen anderen Planeten lenken. Beim Hohmanntansfer wird die Rakete zweimal beschleunigt, einmal, um von der ursprünglichen Umlaufbahn auf eine Ellipsenbahn zu kommen, und dann, um von der Ellipsenbahn auf die endgültige Umlaufbahn zu kommen. Durch die erste Beschleunigung wird das Aphel der Flugbahn auf die geplante Umlaufbahn angehoben, durch die zweite Beschleunigung, wenn die Rakete das Aphel erreicht hat, wird auch das Perihel angehoben, damit die Rakete die Ellipsenbahn verlässt und auf neuen Umlaufbahn bleibt. Aphel ist der von der Sonne am weitesten entfernteste Punkt und das Perihel der sonnennächste Punkt auf der Ellipse.

Die Geschwindigkeitsänderungen lassen sich mit folgenden
Konstanten berechnen:

$M$ = Masse des Zentralkörpers, in diesem Fall der Sonne
$G$ = Gravitationskonstante
$ra$ = Zielradius
$re$ = Anfangsradius
$a$ = große Halbachse der Ellipse = $(re+ra)/2\\$ $va$ = Anfangsgeschwindigkeit
$vz$ = Bahngeschwindigkeit des Zielorbits

Die Formel für die erste Geschwindigkeitsänderung $dv1$ ist:

$dv1 = sqrt[(M*G*ra)/(re*a)] - va$

Die Formel für die zweite Geschwindigkeitsänderung $dv2$ ist:

$dv2 = vz - sprt[(M*G*re)/(ra*a)]$

Weil wir die Rakete nicht nur von einer Umlaufbahn auf eine andere schicken wollen, sondern sie bei einem Planeten startet und auch bei dem anderen Planeten und nicht nur auf dessen Umlaufbahn ankommen soll, müssen auch die Positionen der Planeten zueinander und deren Geschwindigkeiten berücksichtigt werden. Dafür muss die Dauer des Fluges der Rakete von der einen Umlaufbahn zur anderen berechnet werden und die Rakete muss genau in dem Moment starten, wenn der Zielplanet zu dem Punkt genau so lange braucht wie die Rakete.

def gravitation(obj):
    a = [0, 0, 0]
    for obj_ in objects:  # Einfluss aller anderen Objekte wird berechnet
        if obj != obj_:
            nvector_l = sqrt(
                (obj_.pos[0] - obj.pos[0]) ** 2 + (obj_.pos[1] - obj.pos[1]) ** 2 + (obj_.pos[2] - obj.pos[2]) ** 2)
            nvector = ((obj_.pos[0] - obj.pos[0]) / nvector_l, (obj_.pos[1] - obj.pos[1]) / nvector_l,
                       (obj_.pos[2] - obj.pos[2]) / nvector_l)
            acc = G * (obj_.m / (nvector_l ** 2))
            a[0] += acc * nvector[0]
            a[1] += acc * nvector[1]  # Normalenvektor wird berechnet und mit der Beschleunigung multipliziert
            a[2] += acc * nvector[2]
 
    return a  # nach Durchlauf aller Objekte wird der Beschleunigungsvektor zurückgegeben

In der Gravitationfunktion wird die Kraft aller Objekte auf ein Objekt errechnet. Dazu geht die for-Schleife alle Objekte durch und berechnet einen genormten Vektor (Länge : 1) und die Beschleunigung zu diesem.
Der Wert wird jeden Durchgang der Schleife auf den Beschleunigungsvektor addiert, welcher dann am Ende zurückgegeben wird

def import_objects():
    i = 0
    for id_ in [1, 2, 3, 4, 5, 6, 7, 8, 301, -143205]:  # ids Horizon Objekte
        obj = Horizons(id=id_, location="@sun", epochs=Time(time).jd, id_type='id').vectors()
        objects.append(
            object(sizes[i],
                   names[i], id_, mass[i], obj["x"][0], obj["y"][0], obj["z"][0], obj["vx"][0],
                   obj["vy"][0], obj["vz"][0]))
        i += 1
    objects.append(object(80, "Sun", 0, 1988500E24, 0, 0, 0))
    # objects.append(object(5, "Satelit", 11, 10, objects[2].pos[0] + 0.01, objects[2].pos[1] + 0.01, objects[2].pos[2] + 0.01, objects[2].v[0], objects[2].v[1], objects[2].v[2]))
    # hier lassen sich weitere beliebige Objekte hinzufügen

Um die Planetendaten zu einem gegebenen Zeitpunkt zu erhalten haben wir das Modul Horizons von astroquery genutzt. Dieses greift direkt über telnet auf die HORIZONS-Datenbank der NASA zu.

Jedes Objekt in dieser Datenbank hat eine ID, mit der man sich die gewünschten Daten abgreifen kann. NASA HORIZONS
Es lassen sich auch eigene Objekte, wie die Sonne hinzufügen.

class object():
    def __init__(self, size, name, id, m, x, y, z, vx=0, vy=0, vz=0):
        global objects
        self.size = size
        self.id = id
        self.name = name
        self.m = m
        self.pos = [x * 149597900000, y * 149597900000, z * 149597900000]  # AU in m
        self.v = [vx * 1731460, vy * 1731460, vz * 1731460]  # AU/d in m/s
 
    def update(self):
        dt = 100
        a = gravitation(self)
        self.v[0] += dt * a[0]
        self.v[1] += dt * a[1]
        self.v[2] += dt * a[2]
        self.pos[0] += dt * self.v[0]
        self.pos[1] += dt * self.v[1]
        self.pos[2] += dt * self.v[2]

Die Klasse object() bildet das Kernstück der Simulation. Alles zu Simulierende wird durch das Erstellen einer Instanz dem System hinzugefügt.
Neben den grundlegenden Variablen wie Positions- und Geschwindigkeitsvektor hat die Klasse auch die Funktion update()
Diese ruft bei jeder Aktivierung die Gravitationsfunktion mit sich selbst als Parameter auf.
Nach Durchlauf der Funktion werden die Eigenschaften Geschwindigkeit und Position des Objekts aktualisiert

data = {"x": [],  # Speicherort der Frame
        "y": [],
        "z": [],
        "size": [],
        "Object": [],
        "date": []}
 
 
def get_data():
    i = 1
    j = 365
 
    step = 0  # Tag:86400 Sekunden, pro tag 864 Updates--> Beschleunigung mit Faktor 100
    steps = 864
    while i <= j:
        print("calc frame :", str(i))
        for obj in objects:
            # geht objekte durch und ruft update-funktion auf
            data["x"].append(obj.pos[0] / 149597900000)
            data["y"].append(obj.pos[1] / 149597900000)  # speichert objektvariablen in data in AU
            data["z"].append(obj.pos[2] / 149597900000)
            data["size"].append(obj.size)
            data["Object"].append(obj.name)
            data["date"].append(i)
        while step < steps:
            for obj_ in objects:
                obj_.update()
 
            step += 1
        step = 0
        i += 1

Die Funktion get_data() dient zum Aktualisieren und Speichern der Schritte in dem Dictionary data. Um eine annährend genaue, aber auch relativ schnelle Berechnung zu erhalten haben wir uns für eine Einteilung von 864 Zeitschritten pro Tag à 100 Sekunden entschieden.
Diese Aufteilung ist klein genug um auch schnelle Trabanten mit geringer Orbithöhe zuverlässig zu simulieren.
Das Programm geht die Tage und Zeitschritte der Tage durch und speichert den Zustand aller Objekte.