Try to simulate the movement of the solar system

This is a New Year Advent calendar sent by the "Lehman Sat Project", an organization that develops space as a hobby. You can see the summary article from here.

who are you? I am working as an attitude control system in the rsp-01 development team of a 10 cm cubic ultra-small artificial satellite aiming for a selfie mission with the arms deployed. I usually do numerical calculations for rocket engines at university.

Introduction

Articles dealing with ISS orbit and satellite attitude have popped up. In this article, I would like to aim to display the motion of solar system planets farther away. Below is a plot of the orbit calculated based on the orbital elements. Mercury, Venus, ..., Neptune are displayed from the inside. The unit in the figure is one astronomical unit, which is the average distance between the sun and the earth. orbit_1.png

Planetary orbit

1. 1. Definition of variables required for orbit determination

The planetary orbit can be determined as an elliptical shape and requires a semi-major axis $ a $ and an eccentricity $ e $. The two parameters can be determined with the Julius Century $ T_ {TDB} $ as arguments, relative to noon, January 1, 2000. For the coefficients, refer to "Basics of Mission Analysis and Orbital Design" and the elements of the appendix. It is created like py.

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

2. Derivation of perigee elongation by Newton's method

The average perigee elongation $ M $, which indicates the position of the satellite at each time, is also determined by the above equation. The perigee distance $ E $ required to determine the coordinates of the planet is calculated using the eccentricity $ e $ obtained so far and the average perigee distance $ M $. The Kepler equation is used at this time.

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

This equation uses a numerical solution because there is no analytical solution. Here, if $ f (E) $ is defined

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

And

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

By repeating the update of the coordinates, the solution of the equation is approached as shown in the figure below. The iterative calculation is terminated when the change in the value of $ E $ becomes small enough. image.png

3. 3. Conversion to x, y coordinate system

When the distance is defined as shown in the figure, the $ x $ and $ y $ coordinates of the planet are defined by the following equation using the properties of the ellipse. Elements.png

\begin{align}
x&=a\mathrm{cos}E-ae\\
y&=a\sqrt{1-e^2}\mathrm{sin}E
\end{align}
    1. And 2. By substituting $ a $, $ e $, and $ E $ obtained in, the position of the planet $ x $, $ y $ in a certain Julius century $ T_ {TDB} $ was obtained.
    1. ~ 3. You can follow the movement of the planet by repeatedly executing the calculation of.

This article was created in the Lehman Sat Project Advent Calendar. The Lehman Sat Project is a private organization that works under the slogan "Let's get together and develop space". Everyone can be involved in a "space development project" that cannot be experienced anywhere else. If you are interested, please feel free to visit https://www.rymansat.com/join.

@ kentani's "Artificial satellite with everyone, remote control of home TV by one person" is also recommended.

appendix

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

Try to simulate the movement of the solar system
Try to solve the problems / problems of "Matrix Programmer" (Chapter 1)
Try to estimate the number of likes on Twitter
Try to get the contents of Word with Golang
Try to get the function list of Python> os package
Try to evaluate the performance of machine learning / regression model
Try to evaluate the performance of machine learning / classification model
Try to improve the accuracy of Twitter like number estimation
Try to solve the problems / problems of "Matrix Programmer" (Chapter 0 Functions)
Try to automate the operation of network devices with Python
Try to extract the features of the sensor data with CNN
Try to introduce the theme to Pelican
Cython to try in the shortest
The fastest way to try EfficientNet
Supplement to the explanation of vscode
The easiest way to try PyQtGraph
[Linux] I tried to summarize the command of resource confirmation system
[Note] Let's try to predict the amount of electricity used! (Part 1)
[Cloudian # 9] Try to display the metadata of the object in Python (boto3)
[Python] Try to graph from the image of Ring Fit [OCR]
First python ② Try to write code while examining the features of python
Try to solve the N Queens problem with SA of PyQUBO
Try to model the cumulative return of rollovers in futures trading
Try using Elasticsearch as the foundation of a question answering system
Try to predict the triplet of boat race by ranking learning
The story of trying to reconnect the client
Script to change the description of fasta
10 methods to improve the accuracy of BERT
How to check the version of Django
Try to make a kernel of Jupyter
The story of adding MeCab to ubuntu 16.04
Try to estimate the parameters of the gamma distribution while simply implementing MCMC
Try to image the elevation data of the Geographical Survey Institute with Python
Try to detect fusion movement using AnyMotion
Try to get the road surface condition using big data of road surface management
[Python] Try pydash of the Python version of lodash
Try to prepare each environment of kivy
Python amateurs try to summarize the list ①
Try using n to downgrade the version of Node.js you have installed
The story of pep8 changing to pycodestyle
Try to separate the background and moving object of the video with OpenCV
I made a function to see the movement of a two-dimensional array (Python)
Try to measure the position of the object on the desk (real coordinate system) from the camera image with Python + OpenCV
Using the python class, try to roughly simulate the production behavior of 1 farmer and the consumption behavior of 1 consumer ① ~ Time to explain the conditions ~
Let's simulate the transition of infection rate with respect to population density with python
Until you try to let DNN learn the truth of the image using Colab
Try to display the Fibonacci sequence in various languages in the name of algorithm practice
Organize Python tools to speed up the initial movement of data analysis competitions
I tried to solve the 2020 version of 100 language processing knocks [Chapter 1: Preparatory movement 00-04]
Try to display the railway data of national land numerical information in 3D
[Verification] Try to align the point cloud with the optimization function of pytorch Part 1
The story of switching the Azure App Service web system from Windows to Linux
I tried to solve the 2020 version of 100 language processing knocks [Chapter 1: Preparatory movement 05-09]
How to calculate the volatility of a brand
Try to solve the fizzbuzz problem with Keras
Try installing only the core part of Ubuntu
Terms closely related to the X Window System
Try the python version of emacs-org parser orgparse
How to find the area of the Voronoi diagram
Try adding fisheye lens distortion to the image
Converting the coordinate system to ECEF and geodesy