Übung: 2D-Simulation der Orbitalmechanik (Python)

Nur ein kleiner Hinweis: Ich habe noch nie Astronomie oder irgendwelche exakten Wissenschaften studiert (nicht einmal IT), also versuche ich, diese Lücke durch Autodidaktik zu füllen. Die Astronomie ist einer der Bereiche, auf die ich aufmerksam geworden bin, und meine Vorstellung von Selbstbildung beruht auf einem angewandten Ansatz. Also gleich zur Sache - das ist ein Orbital-Simulationsmodell, an dem ich nebenbei arbeite, wenn ich Zeit / Stimmung habe. Mein Hauptziel ist die Schaffung eines vollständigen Sonnensystems in Bewegung und die Fähigkeit, Raumschiffstarts auf andere Planeten zu planen.

Es steht Ihnen allen frei, dieses Projekt jederzeit in die Hand zu nehmen und Spaß beim Experimentieren zu haben!

aktualisieren!!! (10. Nov.)

• Geschwindigkeit ist jetzt richtig deltaV und zusätzliche Bewegung berechnet jetzt den Summenvektor der Geschwindigkeit
• Sie können beliebig viele statische Objekte auf jedes in Bewegung befindliche Einheitsobjekt setzen und nach Schwerkraftvektoren aus allen Quellen suchen (und nach Kollisionen suchen).
• erheblich verbessert die Leistung von Berechnungen
• Ein Fix für den interaktiven Mod in Matplotlib. Sieht so aus, als wäre dies die Standardoption nur für ipython. Für reguläres python3 ist diese Anweisung explizit erforderlich.

Grundsätzlich ist es jetzt möglich, ein Raumschiff von der Erdoberfläche aus zu "starten" und eine Mission zum Mond zu zeichnen, indem DeltaV-Vektorkorrekturen über giveMotion () vorgenommen werden. Als nächstes wird versucht, eine globale Zeitvariable zu implementieren, um eine gleichzeitige Bewegung zu ermöglichen, z. B. umkreist der Mond die Erde, während das Raumschiff ein Manöver zur Unterstützung der Schwerkraft ausprobiert.

Kommentare und Verbesserungsvorschläge sind immer willkommen!

In Python3 mit der matplotlib-Bibliothek erledigt

import matplotlib.pyplot as plt
import math
plt.ion()

G = 6.673e-11  # gravity constant
gridArea = [0, 200, 0, 200]  # margins of the coordinate grid
gridScale = 1000000  # 1 unit of grid equals 1000000m or 1000km

plt.clf()  # clear plot area
plt.axis(gridArea)  # create new coordinate grid
plt.grid(b="on")  # place grid

class Object:
    _instances = []
    def __init__(self, name, position, radius, mass):
        self.name = name
        self.position = position
        self.radius = radius  # in grid values
        self.mass = mass
        self.placeObject()
        self.velocity = 0
        Object._instances.append(self)

    def placeObject(self):
        drawObject = plt.Circle(self.position, radius=self.radius, fill=False, color="black")
        plt.gca().add_patch(drawObject)
        plt.show()

    def giveMotion(self, deltaV, motionDirection, time):
        if self.velocity != 0:
            x_comp = math.sin(math.radians(self.motionDirection))*self.velocity
            y_comp = math.cos(math.radians(self.motionDirection))*self.velocity
            x_comp += math.sin(math.radians(motionDirection))*deltaV
            y_comp += math.cos(math.radians(motionDirection))*deltaV
            self.velocity = math.sqrt((x_comp**2)+(y_comp**2))

            if x_comp > 0 and y_comp > 0:  # calculate degrees depending on the coordinate quadrant
                self.motionDirection = math.degrees(math.asin(abs(x_comp)/self.velocity))  # update motion direction
            elif x_comp > 0 and y_comp < 0:
                self.motionDirection = math.degrees(math.asin(abs(y_comp)/self.velocity)) + 90
            elif x_comp < 0 and y_comp < 0:
                self.motionDirection = math.degrees(math.asin(abs(x_comp)/self.velocity)) + 180
            else:
                self.motionDirection = math.degrees(math.asin(abs(y_comp)/self.velocity)) + 270

        else:
            self.velocity = self.velocity + deltaV  # in m/s
            self.motionDirection = motionDirection  # degrees
        self.time = time  # in seconds
        self.vectorUpdate()

    def vectorUpdate(self):
        self.placeObject()
        data = []

        for t in range(self.time):
            motionForce = self.mass * self.velocity  # F = m * v
            x_net = 0
            y_net = 0
            for x in [y for y in Object._instances if y is not self]:
                distance = math.sqrt(((self.position[0]-x.position[0])**2) +
                             (self.position[1]-x.position[1])**2)
                gravityForce = G*(self.mass * x.mass)/((distance*gridScale)**2)

                x_pos = self.position[0] - x.position[0]
                y_pos = self.position[1] - x.position[1]

                if x_pos <= 0 and y_pos > 0:  # calculate degrees depending on the coordinate quadrant
                    gravityDirection = math.degrees(math.asin(abs(y_pos)/distance))+90

                elif x_pos > 0 and y_pos >= 0:
                    gravityDirection = math.degrees(math.asin(abs(x_pos)/distance))+180

                elif x_pos >= 0 and y_pos < 0:
                    gravityDirection = math.degrees(math.asin(abs(y_pos)/distance))+270

                else:
                    gravityDirection = math.degrees(math.asin(abs(x_pos)/distance))

                x_gF = gravityForce * math.sin(math.radians(gravityDirection))  # x component of vector
                y_gF = gravityForce * math.cos(math.radians(gravityDirection))  # y component of vector

                x_net += x_gF
                y_net += y_gF

            x_mF = motionForce * math.sin(math.radians(self.motionDirection))
            y_mF = motionForce * math.cos(math.radians(self.motionDirection))
            x_net += x_mF
            y_net += y_mF
            netForce = math.sqrt((x_net**2)+(y_net**2))

            if x_net > 0 and y_net > 0:  # calculate degrees depending on the coordinate quadrant
                self.motionDirection = math.degrees(math.asin(abs(x_net)/netForce))  # update motion direction
            elif x_net > 0 and y_net < 0:
                self.motionDirection = math.degrees(math.asin(abs(y_net)/netForce)) + 90
            elif x_net < 0 and y_net < 0:
                self.motionDirection = math.degrees(math.asin(abs(x_net)/netForce)) + 180
            else:
                self.motionDirection = math.degrees(math.asin(abs(y_net)/netForce)) + 270

            self.velocity = netForce/self.mass  # update velocity
            traveled = self.velocity/gridScale  # grid distance traveled per 1 sec
            self.position = (self.position[0] + math.sin(math.radians(self.motionDirection))*traveled,
                             self.position[1] + math.cos(math.radians(self.motionDirection))*traveled)  # update pos
            data.append([self.position[0], self.position[1]])

            collision = 0
            for x in [y for y in Object._instances if y is not self]:
                if (self.position[0] - x.position[0])**2 + (self.position[1] - x.position[1])**2 <= x.radius**2:
                    collision = 1
                    break
            if collision != 0:
                print("Collision!")
                break

        plt.plot([x[0] for x in data], [x[1] for x in data])

Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24)
Moon = Object(name="Moon", position=(100.0, 100.0), radius=1.737, mass = 7.347e22)  # position not to real scale
Craft = Object(name="SpaceCraft", position=(49.0, 40.0), radius=1, mass=1.0e4)

Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000)
Craft.giveMotion(deltaV=2000.0, motionDirection=90, time=60000)
plt.show(block=True)

Wie es funktioniert

Alles läuft auf zwei Dinge hinaus:

1. Erstellen von Objekten wie Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24)mit Positionsparametern auf dem Gitter (1 Einheit des Gitters ist standardmäßig 1000 km, aber dies kann auch geändert werden), Radius in Gittereinheiten und Masse in kg.
2. Geben Sie dem Objekt ein DeltaV, wie es Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000)offensichtlich Craft = Object(...)an erster Stelle erstellt werden muss, wie im vorherigen Punkt erwähnt. Die Parameter hier sind deltaVin m / s (beachten Sie, dass die Beschleunigung momentan erfolgt), motionDirectiondie Richtung von deltaV in Grad (von der aktuellen Position ausgehend stellen Sie sich einen 360-Grad-Kreis um das Objekt vor, die Richtung ist also ein Punkt auf diesem Kreis) und der Parameter schließlich, timewie viele Sekunden Nach dem DeltaV wird die Push-Trajektorie des Objekts überwacht. Anschließender giveMotion()Start von der letzten Position der vorherigen giveMotion().

Fragen:

1. Ist dies ein gültiger Algorithmus zur Berechnung der Umlaufbahnen?
2. Was sind die offensichtlichen Verbesserungen?
3. Ich habe die Variable "timeScale" in Betracht gezogen, um die Berechnungen zu optimieren, da es möglicherweise nicht erforderlich ist, Vektoren und Positionen für jede Sekunde neu zu berechnen. Überlegen Sie, wie es umgesetzt werden soll, oder ist es im Allgemeinen eine gute Idee? (Genauigkeitsverlust gegenüber verbesserter Leistung)

Grundsätzlich ist es mein Ziel, eine Diskussion über das Thema zu beginnen und zu sehen, wohin es führt. Und wenn möglich etwas Neues und Interessantes lernen (oder noch besser lehren).

Fühlen Sie sich frei zu experimentieren!

Versuchen Sie es mit:

Earth = Object(name="Earth", position=(50.0, 100.0), radius=6.371, mass=5.972e24)
Moon = Object(name="Moon", position=(434.0, 100.0), radius=1.737, mass = 7.347e22)
Craft = Object(name="SpaceCraft", position=(43.0, 100.0), radius=1, mass=1.0e4)

Craft.giveMotion(deltaV=10575.0, motionDirection=180, time=322000)
Craft.giveMotion(deltaV=400.0, motionDirection=180, time=50000)

Mit zwei Verbrennungen - einer auf der Erdumlaufbahn und einer auf der Mondumlaufbahn - erreichte ich eine stabile Mondumlaufbahn. Liegen diese nahe an den theoretisch erwarteten Werten?

Empfohlene Übung: Versuchen Sie es in 3 Verbrennungen - stabile Erdumlaufbahn von der Erdoberfläche, progressive Verbrennung, um den Mond zu erreichen, retrograde Verbrennung, um die Umlaufbahn um den Mond zu stabilisieren. Versuchen Sie dann, DeltaV zu minimieren.

Hinweis: Ich plane, den Code mit umfangreichen Kommentaren für diejenigen zu aktualisieren, die nicht mit der Python3-Syntax vertraut sind.

$m_1, m_2$

$r_{21}$$F_{12}=-F_{21}$$r_{12}=-r_{21}$$(x_1,y_1)$$(x_2,y_2)$

und

Dieses System der gewöhnlichen Differentialgleichungen (ODEs) bildet zusammen mit den Ausgangspositionen und -geschwindigkeiten ein Anfangswertproblem. Der übliche Ansatz ist, dies als ein System erster Ordnung aus 8 Gleichungen zu schreiben und eine Runge-Kutta- oder Mehrschrittmethode anzuwenden, um sie zu lösen.

Wenn Sie etwas Einfaches wie Vorwärts-Euler oder Rückwärts-Euler anwenden, sehen Sie, wie sich die Erde zur Sonne hin ins Unendliche oder hinein windet, aber das ist ein Effekt der numerischen Fehler. Wenn Sie eine genauere Methode verwenden, wie die klassische Runge-Kutta-Methode 4. Ordnung, werden Sie feststellen, dass sie für eine Weile in der Nähe einer echten Umlaufbahn bleibt, sich aber schließlich bis ins Unendliche fortsetzt. Der richtige Ansatz ist die Verwendung einer symplektischen Methode, die die Erde in der richtigen Umlaufbahn hält - obwohl ihre Phase aufgrund von numerischen Fehlern immer noch nicht stimmt.

Für das 2-Körper-Problem ist es möglich, ein einfacheres System abzuleiten, indem Sie Ihr Koordinatensystem um den Schwerpunkt legen. Aber ich denke, die obige Formulierung ist klarer, wenn dies für Sie neu ist.