[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. orbit_1.png

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é. image.png

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. Elements.png

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

Recommended Posts

Essayez de simuler le mouvement du système solaire
Essayez de résoudre les problèmes / problèmes du "programmeur matriciel" (Chapitre 1)
Essayez d'estimer le nombre de likes sur Twitter
Essayez d'obtenir le contenu de Word avec Golang
Essayez d'obtenir la liste des fonctions du paquet Python> os
Essayez d'évaluer les performances du modèle d'apprentissage automatique / de régression
Essayez d'évaluer les performances du modèle d'apprentissage automatique / de classification
Essayez d'améliorer la précision de l'estimation du nombre de Twitter
Essayez de résoudre les problèmes / problèmes du "programmeur matriciel" (fonction du chapitre 0)
Essayez d'automatiser le fonctionnement des périphériques réseau avec Python
Essayez d'extraire les caractéristiques des données de capteur avec CNN
Essayez d'introduire le thème sur Pelican
Essayez Cython dans les plus brefs délais
Le moyen le plus rapide d'essayer EfficientNet
Supplément à l'explication de vscode
La façon la plus simple d'essayer PyQtGraph
[Linux] J'ai essayé de résumer les commandes de confirmation des ressources
[Note] Essayons de prédire la quantité d'électricité utilisée! (Partie 1)
Premier python ② Essayez d'écrire du code tout en examinant les fonctionnalités de python
Essayez de résoudre le problème N Queen avec SA de PyQUBO
Essayez de modéliser le rendement cumulatif du roulement dans le trading à terme
Essayez d'utiliser Elasticsearch comme base de votre système de questions et réponses
Essayez de prédire le triplet de la course de bateaux en classant l'apprentissage
L'histoire d'essayer de reconnecter le client
Script pour changer la description de fasta
10 méthodes pour améliorer la précision de BERT
Comment vérifier la version de Django
Faisons un noyau jupyter
L'histoire de la mise en place de MeCab dans Ubuntu 16.04
Essayez d'estimer les paramètres de la distribution gamma tout en implémentant simplement MCMC
Essayez d'imaginer les données d'élévation du National Land Research Institute avec Python
Essayez de détecter les mouvements de fusion en utilisant AnyMotion
Essayez d'obtenir l'état de la surface de la route en utilisant de grandes données de gestion de la surface de la route
[Python] Essayez pydash de la version Python de lodash
Essayez de préparer chaque environnement de kivy
Python amateur tente de résumer la liste ①
Essayez d'utiliser n pour rétrograder la version de Node.js que vous avez installée
L'histoire du changement de pep8 en pycodestyle
Essayez de séparer l'arrière-plan et l'objet en mouvement de la vidéo avec OpenCV
J'ai créé une fonction pour voir le mouvement d'un tableau à deux dimensions (Python)
Essayez de mesurer la position d'un objet sur le bureau (système de coordonnées réel) à partir de l'image de la caméra avec Python + OpenCV
En utilisant la classe python, essayez de simuler grossièrement le comportement de production d'un agriculteur et le comportement de consommation d'un consommateur ① ~ Il est temps d'expliquer les conditions ~
Simulons la transition du taux d'infection par rapport à la densité de population avec python
Jusqu'à ce que vous essayiez de laisser DNN apprendre la vérité de l'image en utilisant Colab
Essayez d'afficher la séquence de Fibonacci dans différentes langues comme pratique d'algorithme
Organiser les outils Python pour accélérer le mouvement initial des compétitions d'analyse de données
J'ai essayé de résoudre la version 2020 de 100 coups de traitement de langue [Chapitre 1: Mouvement préparatoire 00-04]
Essayez d'afficher les données ferroviaires des informations numériques des terres nationales en 3D
[Vérification] Essayez d'aligner le groupe de points avec la fonction d'optimisation de pytorch Partie 1
L'histoire du passage du système Web Azure App Service de Windows à Linux
J'ai essayé de résoudre la version 2020 de 100 traitements linguistiques [Chapitre 1: Mouvement préparatoire 05-09]
Comment calculer la volatilité d'une marque
Essayez de résoudre le problème du fizzbuzz avec Keras
Essayez d'installer uniquement la partie principale d'Ubuntu
Termes étroitement liés au système X Window
Essai du parseur d'emacs-org orgparse pour python
Comment trouver la zone du diagramme de Boronoi
Essayez d'ajouter la distorsion de l'objectif fisheye à l'image
Conversion du système de coordonnées en ECEF et géodésique