Ćwiczenie: symulacja mechaniki orbitalnej 2D (python)

Tylko małe zastrzeżenie wcześniej: nigdy nie studiowałem astronomii ani żadnych nauk ścisłych w tym zakresie (nawet IT), więc staram się wypełnić tę lukę samokształceniem. Astronomia jest jednym z obszarów, który zwrócił moją uwagę, a moja idea samokształcenia opiera się na podejściu stosowanym. Od razu do rzeczy - to model symulacji orbitalnej, nad którym od niechcenia pracuję, kiedy mam czas / nastrój. Moim głównym celem jest stworzenie kompletnego układu słonecznego w ruchu i możliwość planowania startów statków kosmicznych na inne planety.

Wszyscy mogą swobodnie wybrać ten projekt w dowolnym momencie i dobrze się bawić eksperymentując!

aktualizacja!!! (10 listopada)

prędkość jest teraz właściwa deltaV i dając dodatkowy ruch oblicza teraz sumaryczny wektor prędkości
możesz umieścić tyle obiektów statycznych, ile chcesz, za każdym razem, gdy obiekt w ruchu sprawdza wektory grawitacyjne ze wszystkich źródeł (i sprawdza kolizję)
znacznie poprawiła wydajność obliczeń
poprawka do konta dla interaktywnego modu w matplotlib. Wygląda na to, że jest to domyślna opcja tylko dla ipython. Zwykły python3 wymaga tego wyrażenia wprost.

Zasadniczo można teraz „wystrzelić” statek kosmiczny z powierzchni Ziemi i zaplanować misję na Księżyc, dokonując korekt wektora deltaV za pomocą giveMotion (). Następnie próbuje zaimplementować globalną zmienną czasową, aby umożliwić jednoczesny ruch, np. Księżyc okrąża Ziemię, podczas gdy statek kosmiczny próbuje manewru wspomagania grawitacyjnego.

Komentarze i sugestie dotyczące ulepszeń są zawsze mile widziane!

Wykonano w Python3 z biblioteką matplotlib

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)

Jak to działa

Wszystko sprowadza się do dwóch rzeczy:

Tworzenie obiektu jak za Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24)pomocą parametrów pozycji na siatce (1 jednostka siatki to domyślnie 1000 km, ale można to również zmienić), promień w jednostkach siatki i masa w kg.
Nadanie obiektowi pewnej deltaV, takiej jak Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000)oczywiście, wymaga Craft = Object(...)utworzenia w pierwszej kolejności, jak wspomniano w poprzednim punkcie. Parametry tutaj są deltaVwm / s (zwróć uwagę, że na razie przyspieszenie jest chwilowe), motionDirectionjest kierunkiem deltaV w stopniach (od aktualnej pozycji wyobraź sobie okrąg o 360 stopni wokół obiektu, więc kierunek jest punktem na tym okręgu), a na koniec parametr timeokreśla liczbę sekund po przejściu deltaV obiekt będzie monitorowany. Kolejne giveMotion()rozpoczęcie od ostatniej pozycji poprzedniej giveMotion().

Pytania:

Czy to prawidłowy algorytm do obliczania orbit?
Jakie oczywiste ulepszenia należy wprowadzić?
Rozważałem zmienną „timeScale”, która zoptymalizuje obliczenia, ponieważ może nie być konieczne ponowne obliczanie wektorów i pozycji na każdą sekundę. Czy są jakieś przemyślenia na temat tego, jak należy go wdrożyć, czy ogólnie jest to dobry pomysł? (utrata dokładności vs poprawiona wydajność)

Zasadniczo moim celem jest rozpoczęcie dyskusji na ten temat i sprawdzenie, dokąd prowadzi. I jeśli to możliwe, naucz się (a jeszcze lepiej - naucz) czegoś nowego i interesującego.

Zapraszam do eksperymentowania!

Spróbuj użyć:

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)

Dzięki dwóm oparzeniom - jednemu postępowi na orbicie Ziemi i jednemu ruchowi wstecz na orbicie Księżyca osiągnąłem stabilną orbitę Księżyca. Czy są one zbliżone do teoretycznie oczekiwanych wartości?

Sugerowane ćwiczenie: Wypróbuj w 3 oparzeniach - stabilna orbita Ziemi od powierzchni Ziemi, spalenie postępowe, aby dotrzeć do Księżyca, spalenie wsteczne, aby ustabilizować orbitę wokół Księżyca. Następnie spróbuj zminimalizować deltaV.

Uwaga: Planuję zaktualizować kod obszernymi komentarzami dla tych, którzy nie znają składni python3.

orbit gravity n-body-simulations

— przestrzeń państwowa
źródło

Bardzo dobry pomysł na samokształcenie! Czy byłoby możliwe podsumowanie twoich formuł dla tych z nas, którzy nie znają składni Pythona?

Pewnie tak. Będę robił obszerniejsze komentarze w kodzie dla tych, którzy chcą go podnieść i podsumować ogólną logikę w samym pytaniu.

— Statespace

Z czubka mojej głowy: rozważ użycie wektora prędkości zamiast traktować inaczej prędkość i kierunek. Kiedy mówisz „F = m * v”, masz na myśli „F = m * a”? Zakładasz, że Ziemia się nie porusza, ponieważ jest o wiele cięższa niż asteroida? Zastanów się, czy nie spojrzeć na github.com/barrycarter/bcapps/blob/master/bc-grav-sim.pl

— barrycarter

Możesz wykonać ruch dowolnego obiektu, w tym Ziemi. Do celów testowych podałem tylko relację obiekt -> Ziemia w głównej pętli. Można łatwo przekonwertować, że każdy obiekt odnosi się do wszystkich innych tworzonych obiektów. I każdy obiekt może mieć własny wektor ruchu. Powód, dla którego tego nie zrobiłem - bardzo wolne obliczenia nawet dla 1 obiektu. Mam nadzieję, że skalowanie jednostek czasu powinno bardzo pomóc, ale nadal nie jestem pewien, jak zrobić to dobrze.

— Statespace

DOBRZE. Myśl: czy wykonasz symulację dla dwóch rzeczywistych obiektów (np. Ziemia / Księżyc lub Ziemia / Słońce) i porównujesz swoje wyniki z ssd.jpl.nasa.gov/?horizons w celu uzyskania dokładności? Nie będzie to idealne z powodu zaburzeń pochodzących z innych źródeł, ale da ci pojęcie o dokładności?

— barrycarter

$m_1, m_2$

F = m a

$F = ma$

a

$a$

F_{21} = \frac{G m_{1} m_{2}}{| r_{21} |^{3}} r_{21}

$F_{21} = \frac{G m_1 m_2}{|r_{21}|^3}r_{21}$

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

r_{21} = (\begin{matrix} x_{1} - x_{2} \\ y_{1} - y_{2} \end{matrix}) .

$r_{21} = \begin{pmatrix} x_1-x_2 \\ y_1-y_2 \end{pmatrix}.$

| r | = \sqrt{(x_{1} - x_{2})^{2} + (y_{1} - y_{2})^{2}} .

$|r| = \sqrt{(x_1-x_2)^2 + (y_1-y_2)^2}.$

a = F / m

$a=F/m$

\begin{aligned} x_{1}^{″} (t) & = \frac{G m_{2} (x_{2} - x_{1})}{| r |^{3}} \\ y_{1}^{″} (t) & = \frac{G m_{2} (y_{2} - y_{1})}{| r |^{3}} \\ x_{2}^{″} (t) & = \frac{G m_{1} (x_{1} - x_{2})}{| r |^{3}} \\ y_{2}^{″} (t) & = \frac{G m_{1} (y_{1} - y_{2})}{| r |^{3}} . \end{aligned}

$\begin{align} x_1''(t) & = \frac{G m_2 (x_2-x_1)}{|r|^3} \\ y_1''(t) & = \frac{G m_2 (y_2-y_1)}{|r|^3} \\ x_2''(t) & = \frac{G m_1 (x_1-x_2)}{|r|^3} \\ y_2''(t) & = \frac{G m_1 (y_1-y_2)}{|r|^3}. \end{align}$

Układ równań różniczkowych zwyczajnych (ODE) wraz z położeniami początkowymi i prędkościami stanowi problem wartości początkowej. Zwykle podejście polega na napisaniu tego jako układu 8 równań pierwszego rzędu i zastosowaniu metody Runge-Kutty lub metody wieloetapowej w celu ich rozwiązania.

Jeśli zastosujesz coś prostego, na przykład Euler do przodu lub Euler do tyłu, zobaczysz, że Ziemia spiralnie przechodzi odpowiednio w nieskończoność lub w stronę słońca, ale jest to efekt błędów numerycznych. Jeśli użyjesz dokładniejszej metody, takiej jak klasyczna metoda Runge-Kutta czwartego rzędu, okaże się, że pozostaje ona przez pewien czas blisko prawdziwej orbity, ale ostatecznie ostatecznie przechodzi w nieskończoność. Właściwe podejście polega na zastosowaniu metody symplektycznej, która utrzyma Ziemię na właściwej orbicie - chociaż jej faza nadal będzie wyłączona z powodu błędów numerycznych.

W przypadku problemu dwóch ciał można uzyskać prostszy układ, opierając układ współrzędnych wokół środka masy. Ale myślę, że powyższe sformułowanie jest jaśniejsze, jeśli jest dla ciebie nowe.

— David Ketcheson
źródło

Potrwa to trochę czasu.

— Statespace,

Wciąż trawię. Zbyt wiele nieznanych mi słów, ale jakoś czuję, że w pewnym momencie się tam dostanę. Na razie mój własny algorytm jest wystarczająco dobry, aby wszystko po prostu działało. Ale kiedy podłączę jednoczesny ruch - będę zmuszony dostać się do literatury i poczytać o odpowiednich algorytmach. Biorąc pod uwagę, że ograniczenia współczesnego sprzętu są znacznie luźniejsze, mogę sobie pozwolić na wygłupianie się za pomocą prostych równań. Nie boję się jednak długo.

— Statespace,

Rzeczywiście metody symplektyczne są zdecydowanie najdokładniejsze, ale myślę, że trudno jest je wdrożyć komuś bez doświadczenia w nauce. Zamiast tego możesz użyć bardzo prostej metody Eulera wraz z korektą Feynmana. Nie sądzę, że potrzebujesz czegoś bardziej złożonego niż to do celów samokształcenia.

— chrispap