リビジョン 4a887337
| learn/python/satellite/orbit/orbit_sorabatake/orbit_odeint_3d.py | ||
|---|---|---|
|
import numpy as np # 数値計算用ライブラリ
|
||
|
from scipy.integrate import odeint # 常微分方程式を解く関数
|
||
|
import matplotlib.pyplot as plt # グラフを描画するライブラリ
|
||
|
|
||
|
|
||
|
GM = 398600.4354360959 # [km^3/s^2] 地球の重力定数（地球の質量と万有引力定数の積）
|
||
|
|
||
|
|
||
|
def func(pv_state, t):
|
||
|
"""
|
||
|
二体間の万有引力が作用する運動の常微分方程式を数値計算で解く
|
||
|
"""
|
||
|
x, y, z, vx, vy, vz = pv_state
|
||
|
r = np.linalg.norm([x, y, z])
|
||
|
|
||
|
dpdt = [vx,
|
||
|
vy,
|
||
|
vz,
|
||
|
-GM * x / r**3,
|
||
|
-GM * y / r**3,
|
||
|
-GM * z / r**3]
|
||
|
return dpdt
|
||
|
|
||
|
|
||
|
# 微分方程式の初期条件
|
||
|
pv0 = [10000, 0, 0, 0, 6, 3.6] # 位置 (x,y,z) と 速度 (vx, vy, vz)
|
||
|
t = np.linspace(0, 86400, 1000) # 1日分の軌道伝播 1000の区分で計算
|
||
|
|
||
|
# 微分方程式を数値計算で解く
|
||
|
sol = odeint(func, pv0, t)
|
||
|
|
||
|
# 軌道を3Dにプロット
|
||
|
fig = plt.figure()
|
||
|
ax = fig.add_subplot(projection='3d')
|
||
|
ax.plot(sol[:, 0], sol[:, 1], sol[:, 2], 'b')
|
||
|
ax.set_aspect('equal')
|
||
|
plt.show()
|
||
Add orbit 3D plot by differential equations