リビジョン 6cba17cc
| learn/python/satellite/orbit/orbit_sorabatake/orbit_odeint.py | ||
|---|---|---|
|
import numpy as np # 数値計算用ライブラリ
|
||
|
from scipy.integrate import odeint # 常微分方程式を解く関数
|
||
|
import matplotlib.pyplot as plt # グラフを描画するライブラリ
|
||
|
|
||
|
|
||
|
def func(x, t):
|
||
|
"""
|
||
|
二体間の万有引力が作用する運動の常微分方程式を数値計算で解く
|
||
|
"""
|
||
|
GM = 398600.4354360959 # [km^3/s^2] 地球の重力定数（地球の質量と万有引力定数の積）
|
||
|
r = np.linalg.norm(x[0:3])
|
||
|
dxdt = [x[3],
|
||
|
x[4],
|
||
|
x[5],
|
||
|
-GM * x[0] / r**3,
|
||
|
-GM * x[1] / r**3,
|
||
|
-GM * x[2] / r**3]
|
||
|
return dxdt
|
||
|
|
||
|
|
||
|
# 微分方程式の初期条件
|
||
|
x0 = [10000, 0, 0, 0, 7, 0] # 位置 (x,y,z) と 速度 (vx, vy, vz)
|
||
|
t = np.linspace(0, 86400, 1000) # 1日分の軌道伝播 1000の区分で計算
|
||
|
|
||
|
# 微分方程式を数値計算で解く
|
||
|
sol = odeint(func, x0, t)
|
||
|
|
||
|
# 軌道を2D(x-y平面）にプロット
|
||
|
plt.plot(sol[:, 0], sol[:, 1], 'b')
|
||
|
plt.grid()
|
||
|
plt.gca().set_aspect('equal')
|
||
|
plt.show()
|
||
| learn/python/satellite/orbit/orbit_sorabatake/requirements.txt | ||
|---|---|---|
|
contourpy==1.3.2
|
||
|
cycler==0.12.1
|
||
|
fonttools==4.58.0
|
||
|
kiwisolver==1.4.8
|
||
|
matplotlib==3.10.3
|
||
|
numpy==2.2.6
|
||
|
packaging==25.0
|
||
|
pillow==11.2.1
|
||
|
pyparsing==3.2.3
|
||
|
python-dateutil==2.9.0.post0
|
||
|
scipy==1.15.3
|
||
|
six==1.17.0
|
||
Add orbit plot by differential equations by following the 'Sorabatake blog'