[PYTHON] Versuchen Sie, die Bewegung des Sonnensystems zu simulieren

Dies ist der Neujahrs-Adventskalender des "Lehman Sat Project", einer Organisation, die den Weltraum als Hobby entwickelt. Sie können den zusammenfassenden Artikel von [hier] sehen (https://qiita.com/rsp-dgkz/items/497225739b88c817a90c).

who are you? Ich arbeite als Lageregelungssystem im rsp-01-Entwicklungsteam eines 10-cm-kubischen ultrakleinen künstlichen Satelliten, der auf eine Selbstporträt-Mission mit eingesetzten Armen abzielt. Normalerweise mache ich numerische Berechnungen für Raketentriebwerke an der Universität.

Einführung

Artikel über ISS-Umlaufbahn und Satelliteneinstellung sind aufgetaucht. In diesem Artikel möchte ich versuchen, die Bewegungen von Planeten des Sonnensystems weiter entfernt darzustellen. Unten sehen Sie eine grafische Darstellung der Umlaufbahn, die anhand der Umlaufbahnelemente berechnet wurde. Merkur, Venus, ..., Kaioh werden von innen angezeigt. Die Einheit in der Abbildung ist eine astronomische Einheit, dh die durchschnittliche Entfernung zwischen Sonne und Erde.

Planetenbahn

1. 1. Definition der für die Umlaufbahnbestimmung erforderlichen Variablen

Die Planetenbahn kann als elliptische Form bestimmt werden und erfordert eine Semi-Major-Achse $ a $ und eine Exzentrizität $ e $. Die beiden Parameter können mit Julius Century $ T_ {TDB} $ relativ zum Mittag am 1. Januar 2000 als Argumente bestimmt werden. Die Koeffizienten finden Sie unter "Grundlagen der Missionsanalyse und des Flugbahndesigns" und in den Elementen des Anhangs. Es ist wie py erstellt.

\begin{align}
\sum_{k=1}^4 c_{p,k}{T_{TDB}}^{k-1}
\end{align}

2. Ableitung des Nahpunkttrennungswinkels nach der Newton-Rapson-Methode

Der durchschnittliche Nahpunkt-Trennungswinkel $ M $, der die Position des Satelliten zu jedem Zeitpunkt angibt, wird ebenfalls durch die obige Gleichung bestimmt. Der zur Bestimmung der Koordinaten des Planeten erforderliche Nahpunkttrennungswinkel $ E $ wird unter Verwendung der bisher erhaltenen Divergenzrate $ e $ und des durchschnittlichen Nahpunkttrennungswinkels $ M $ berechnet. Zu diesem Zeitpunkt wird die Kepler-Gleichung verwendet.

\begin{align}
M=E-e\mathrm{sin}E
\end{align}

Diese Formel verwendet eine numerische Lösung, da es keine analytische Lösung gibt. Hier, wenn $ f (E) $ definiert ist

\begin{align}
f(E)&=M-E+e\mathrm{sin}E\\
f'(E)&=e\mathrm{cos}E-1
\end{align}

Und

\begin{align}
E'&=E-f(E)/f'(E)\\
\end{align}

Durch Wiederholen der Aktualisierung der Koordinaten wird die Lösung der Gleichung wie in der folgenden Abbildung gezeigt angegangen. Wenn die Änderung des Werts von $ E $ klein genug wird, wird die iterative Berechnung gestoppt.

3. 3. Konvertierung in das x, y-Koordinatensystem

Wenn der Abstand wie in der Abbildung gezeigt definiert ist, werden die $ x $ - und $ y $ -Koordinaten des Planeten durch die folgende Gleichung unter Verwendung der Eigenschaften der Ellipse definiert.

\begin{align}
x&=a\mathrm{cos}E-ae\\
y&=a\sqrt{1-e^2}\mathrm{sin}E
\end{align}

1. Und 2. Durch Ersetzen von $ a $, $ e $ und $ E $, die in erhalten wurden, wurde die Position des Planeten $ x $, $ y $ in einem bestimmten Julius-Jahrhundert $ T_ {TDB} $ erhalten.
1. ~ 3. Sie können die Bewegung des Planeten verfolgen, indem Sie die Berechnung von wiederholt ausführen.

Dieser Artikel wurde im Lehman Sat-Projekt Adventskalender erstellt. Das Lehman Sat Project ist eine private Organisation, die unter dem Motto "Lasst uns zusammenkommen und Raum entwickeln" arbeitet. Jeder kann an einem "Weltraumentwicklungsprojekt" beteiligt sein, das nirgendwo anders erlebt werden kann. Wenn Sie interessiert sind, besuchen Sie bitte https://www.rymansat.com/join.

Wir empfehlen auch @ kentanis "Jeder kann einen künstlichen Satelliten und einen Fernseher zu Hause fernsteuern".

Blinddarm

`elements.py`


import numpy as np

def elems(num, JC):
    coeff = np.empty((6, 4))
    JC_term = np.array([[1.], [JC], [JC ** 2], [JC ** 3]])

    if num == 1: #mercury
        coeff = np.array([[0.38709831, 0, 0, 0], # semimajor axis: a
                          [0.20563175, 0.000020406   , -0.0000000284,  0.000000017], # eccentricity: e
                          [7.00498600, -0.00595160   , 0.00000081000,  0.000000041], # orbital inclination: i (for 3D)
                          [48.3308930, -0.12542290   , -0.0000883300, -0.000000196], # longitude of the ascending node: Omega ( for 3D)
                          [77.4561190, 0.158864300   , -0.0000134300,  0.000000039], # Perihelion longitude: ~omega
                          [252.250906, 149472.6746358, -0.0000053500,  0.000000002]]) # mean latitude: lamda
    elif num == 2: #venus
        coeff = np.array([[0.72332982, 0, 0, 0],
                          [0.00677188, -0.000047766, 0.0000000975,  0.000000044],
                          [3.39466200, -0.000856800, -0.00003244,  0.000000010],
                          [76.6799200, -0.278008000, -0.00014256, -0.000000198],
                          [131.563707, 0.004864600 , -0.00138232, -0.000005332],
                          [181.979801, 58517.815676, 0.000001650, -0.000000002]])
    elif num == 3: #earth
        coeff = np.array([[1.000001018, 0, 0, 0],			
                          [0.01670862,	-0.0000420370,	-0.0000001236,	0.00000000004],
                          [0.00000000,	0.01305460000,	-0.0000093100,	-0.0000000340],
                          [0.00000000,	0.00000000000,	 0.0000000000,	0.00000000000],
                          [102.937348,	0.32255570000,	 0.0001502600,	0.00000047800],
                          [100.466449,	35999.3728519,	-0.0000056800,	0.00000000000]])
    elif num == 4: #mars
        coeff = np.array([[1.523679342, 0, 0, 0],
                          [0.093400620, 0.00009048300, -0.0000000806, -0.00000000035],
                          [1.849726000, -0.0081479000, -0.0000225500, -0.00000002700],
                          [49.55809300, -0.2949846000, -0.0006399300, -0.00000214300],
                          [336.0602340, 0.44388980000, -0.0001732100, 0.000000300000],
                          [355.4332750, 19140.2993313, 0.00000261000, -0.00000000300]])
    elif num == 5: #jupiter
        coeff = np.array([[5.202603191, 0.0000001913,             0,           0],
                          [0.048494850, 0.0001632440, -0.0000004719, -0.000000002],
                          [1.303270000, -0.001987200, 0.0000331800 , 0.0000000920],
                          [100.4644410, 0.1766828000, 0.0009038700 , -0.000007032],
                          [14.33130900, 0.2155525000, 0.0007225200 , -0.000004590],
                          [34.35148400, 3034.9056746, -0.0000850100,  0.000000004]])
    elif num == 6: #saturn
        coeff = np.array([[9.554909596, -0.0000021389, 0, 0],
                          [0.055086200, -0.0003468180, 0.0000006456, 0.0000000034],
                          [2.488878000, 0.00255150000, -0.000049030, 0.0000000180],
                          [113.6655240, -0.2566649000, -0.000183450, 0.0000003570],
                          [93.05678700, 0.56654960000, 0.0005280900, 0.0000048820],
                          [50.07747100, 1222.11379430, -0.000085010, 0.0000000040]])
    elif num == 7: #uranus
        coeff = np.array([[19.218446062, -0.0000000372, 0.00000000098, 0],
                          [0.04629590, -0.000027337, 0.0000000790, 0.00000000025],
                          [0.77319600, -0.001686900, 0.0000034900, 0.00000001600],
                          [74.0059470, 0.0741461000, 0.0004054000, 0.00000010400],
                          [173.005159, 0.0893206000, -0.000094700, 0.00000041430],
                          [314.055005, 428.46699830, -0.000004860, 0.00000000600]])
    elif num == 8: #neptune
        coeff = np.array([[30.110386869, -0.0000001663, 0.00000000069, 0],
                          [0.0089880900, 0.00000640800, -0.0000000008, 0],
                          [1.7699520000, 0.00022570000, 0.00000023000, 0.0000000000],
                          [131.78405700, -0.0061651000, -0.0000021900, -0.000000078],
                          [48.123691000, 0.02915870000, 0.00007051000, 0.0000000000],
                          [304.34866500, 218.486200200, 0.00000059000, -0.000000002]])

    temp =  np.dot(coeff, JC_term)
    elements = np.empty((3, 1))
    elements[0] = temp[0]
    elements[1] = temp[1]
    elements[2] = temp[5] - temp[4] # M = lamda - ~omega
    elements[2] = np.deg2rad(elements[2])
    return elements

`plotResult.py`


import matplotlib.pyplot as plt

def plot_2D(state_x, state_y):
    fig = plt.figure()
    plt.plot(state_x[0, :], state_y[0, :], color = 'skyblue')
    plt.plot(state_x[1, :], state_y[1, :], color = 'yellow')
    plt.plot(state_x[2, :], state_y[2, :], color = 'blue')
    plt.plot(state_x[3, :], state_y[3, :], color = 'red')
    plt.plot(state_x[4, :], state_y[4, :], color = 'orange')
    plt.plot(state_x[5, :], state_y[5, :], color = 'springgreen')
    plt.plot(state_x[6, :], state_y[6, :], color = 'violet')    
    plt.plot(state_x[7, :], state_y[7, :], color = 'purple')
    plt.show()

`main.py`


import numpy as np
import elements as elems
import plotResult as pltResult

def newtonRaphson(e, meanE):
    E = 2 * meanE
    eps = 1e-10
    while True:
        E -= (meanE - E + e * np.sin(E)) / (e * np.cos(E) - 1)
        if abs(meanE - E + e * np.sin(E)) < eps:
            return E

def main():
    JDbase = 2451545 # Julian century := 0 -> Julian day := 2451545
    JY2JD = 36525 # Julian year = 0 -> Julian Day = 36525
    dJD = 1
    totalJD = 65000
    JD = np.arange(JDbase, JDbase + totalJD, dJD)
    imax = len(JD)
    planet = 8
    dimension = 2

    state_x = np.empty((planet, imax))
    state_y = np.empty((planet, imax))

    temp_elem = np.empty(3, ) # a, e, M
    r = np.zeros((dimension, 1))

    for i in range(0, imax):
        JC = (JD[i] - JDbase) / JY2JD
        for planum in range(1, planet + 1):
            temp_elem = elems.elems(planum, JC)
            E = newtonRaphson(temp_elem[1], temp_elem[2])
            state_x[planum - 1][i] = temp_elem[0][0] * (np.cos(E) - temp_elem[1][0]) #x
            state_y[planum - 1][i] = temp_elem[0][0] * np.sqrt(1 - temp_elem[1][0] ** 2) * np.sin(E) #y

    pltResult.plot_2D(state_x, state_y)

if __name__ == '__main__':
    main()