[PYTHON] Essayez de simuler le mouvement du système solaire

Ceci est le calendrier de l'Avent du Nouvel An envoyé par le "Lehman Sat Project", une organisation qui développe l'espace comme passe-temps. Vous pouvez consulter l'article de synthèse sur ici.

who are you? Je travaille en tant que système de contrôle d'attitude dans l'équipe de développement rsp-01 d'un satellite artificiel ultra-petit cubique de 10 cm visant une mission d'autoportrait avec les bras déployés. Je fais généralement des calculs numériques pour les moteurs de fusée à l'université.

introduction

Des articles traitant de l'orbite de l'ISS et de satellite attitude sont apparus. Dans cet article, je voudrais viser à afficher les mouvements des planètes du système solaire plus éloignées. Vous trouverez ci-dessous un graphique de l'orbite calculé en fonction des éléments orbitaux. Mercure, Vénus, ..., Kaioh sont affichés de l'intérieur. L'unité de la figure est une unité astronomique, qui est la distance moyenne entre le soleil et la terre.

Orbite planétaire

1. Définition des variables nécessaires à la détermination de l'orbite

L'orbite planétaire peut être déterminée comme une forme elliptique et nécessite un demi-grand axe $ a $ et une excentricité $ e $. Les deux paramètres peuvent être déterminés avec Julius Century $ T_ {TDB} $ par rapport à midi le 1er janvier 2000 comme arguments. Pour les coefficients, reportez-vous à "Bases de l'analyse de mission et de conception de trajectoire" et aux éléments de l'annexe. Il est créé comme py.

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

2. Dérivation de l'angle de séparation du point proche par la méthode de Newton Rapson

L'angle moyen de séparation du point proche $ M $, qui indique la position du satellite à chaque instant, est également déterminé par l'équation ci-dessus. L'angle de séparation du point proche $ E $ nécessaire pour déterminer les coordonnées de la planète est calculé en utilisant le taux de divergence $ e $ obtenu jusqu'à présent et l'angle moyen de séparation du point proche $ M $. L'équation de Kepler est utilisée à ce moment.

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

Cette formule utilise une solution numérique car il n'y a pas de solution analytique. Ici, si $ f (E) $ est défini

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

Et

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

En répétant la mise à jour des coordonnées, la solution de l'équation est abordée comme le montre la figure ci-dessous. Lorsque le changement de la valeur de $ E $ devient suffisamment petit, le calcul itératif est arrêté.

3. 3. Conversion en système de coordonnées x, y

Si la distance est définie comme indiqué sur la figure, les coordonnées $ x $ et $ y $ de la planète sont définies par l'équation suivante en utilisant les propriétés de l'ellipse.

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

1. Et 2. En remplaçant $ a $, $ e $ et $ E $ obtenus dans, la position de la planète $ x $, $ y $ dans un certain siècle de Julius $ T_ {TDB} $ a été obtenue.
2. ~ 3. Vous pouvez suivre le mouvement de la planète en exécutant à plusieurs reprises le calcul de.

Cet article a été créé dans le cadre du projet Lehman Sat Calendrier de l'Avent. Le Lehman Sat Project est une organisation privée active sous le slogan "Rassemblons-nous et développons l'espace". Chacun peut être impliqué dans un «projet de développement spatial» qui ne peut être expérimenté nulle part ailleurs. Si vous êtes intéressé, n'hésitez pas à visiter https://www.rymansat.com/join.

Nous recommandons également @ kentani's "Tout le monde peut contrôler à distance un satellite artificiel et une télévision à la maison".

appendice

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()