[PYTHON] 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.

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.

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.

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