衛星軌道を運動方程式を数値的に解く¶
はじめに¶
次の記事で解説されている、ニュートンの運動方程式を数値的に解くことで衛星軌道を計算することを習ってみます。
Pythonを使って人工衛星の軌道を表現する～軌道６要素、TLE～（宙畑）
この記事では、運動方程式に万有引力をさらっと記載して、すぐにPythonコード（常微分方程式を数値計算）に落としています。
しかし、なぜそうなるのかすぐに理解できなかったので、展開を考えながら追ってみました。
運動方程式をPythonで解く¶
運動方程式を衛星に適用¶
人工衛星はニュートンの運動方程式に従って飛行（運動）します。ニュートンの運動方程式は、次の式で表現されます。
$ m\boldsymbol{a} = \boldsymbol{F} $
- mは人工衛星の質量、 a は加速度、 F は人工衛星に作用する力で、太字はベクトル量を示します。
運動方程式の加速度を位置ベクトルで表現¶
加速度は、次の式で位置ベクトルで置き換えることができます。
$ \boldsymbol{a} = \frac{ d^2 \boldsymbol{r}}{dt^2} $
位置ベクトルで運動方程式を表現すると、次の式となります。
$ m \frac{ d^2 \boldsymbol{r}}{dt^2} = \boldsymbol{F} $
万有引力の大きさ¶
人工衛星に作用する力は、地球と人工衛星との2体間に働く万有引力です。地球の重心を原点として、地球重心と人工衛星重心との距離を r、地球の質量を M 、人工衛星の質量を m 、重力定数（万有引力定数）を G とすると、万有引力の法則から
$ F = G\frac{Mm}{r^2} $
と表せます。
- 大きさ（絶対値）を表わすので、F はベクトル量（太字）ではなくスカラー量で記載
- $ G = 6.67430(15) \times 10^{-11} [Nm^2/kg^2] $
CODATA 2022より（国際学術会議 ISC のデータに関する委員会で4年ごとに物理定数について最新の改定値を公開）。(15)は標準不確かさで $ 0.00015 \times 10^{-11} $ となる。 - $ M = 5.972 \times 10^{24} [kg] $
理科年表より
人工衛星に作用する力（ベクトル量）¶
万有引力による人工衛星に作用する力のベクトル量 F は、原点（地球）向きに働くので、衛星の位置ベクトル量 r を用いると、
万有引力の式に $ - \frac{\boldsymbol{r}}{|r|} $ をかけた次の式になると思います。
$ \boldsymbol{F} = -G\frac{Mm}{r^2} \frac{\boldsymbol{r}}{|r|} $ $ = -G\frac{Mm}{r^3}\boldsymbol{r} $
運動方程式を常微分方程式で表現¶
記事では、Pythonのライブラリ SciPy のodeintライブラリ を使って運動方程式を数値的に解いています。
上述の運動方程式を直接コードに記述してはおらず、2元連立1階微分方程式を立てて解いています。
odeintは、2階微分方程式は解けないので、1階微分方程式を2つ連立して解きます。
$ \frac{d\boldsymbol{r}}{dt} = \boldsymbol{v} $
$ \frac{d\boldsymbol{v}}{dt} = -G\frac{M}{r^3} \boldsymbol{r} $
- 前者は、位置 r の1階微分が速度 v となる式
- 後者は、万有引力で表現した運動方程式 $ m\boldsymbol{a} = -G\frac{Mm}{r^3}\boldsymbol{r} $ の両辺に登場する m を割った残りの式で、加速度 a を速度の微分に置き換えた式
Pythonで数値的に解く¶
SciPyのodeintを使って常微分方程式を数値的に解きます。
利用するライブラリ¶
ベクトルや数値生成などを行う numpy ライブラリ、微分方程式を解く scipy ライブラリ（の中のodeintライブラリ）、グラフをプロットして表示する matplotlib ライブラリ（の中のpyplotライブラリ）をimportします。
import numpy as np # 数値計算ライブラリ
from scipy.integrate import odeint # 常微分方程式を解くライブラリ
import matplotlib.pyplot as plt # グラフ描画ライブラリ
- odeintは、Ordinary Differential Equation の略
微分方程式を漸化式で表現¶
微分方程式を、前の時刻の解から次の時刻の解を求める漸化式で表現（離散化）します。
odeintでは、前の時刻の解と時刻を引数として入力し、次の時刻の解を戻り値として返却する関数として定義します。
# 二体問題の運動方程式
def func(x, t):
GM = 398600.4354360959 # 地心重力定数
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
入力のxは、前の時刻における衛星の位置・速度で[x, y, z, vx, vy, yz] の6つのベクトル（配列）です。
戻り値は、前の時刻の値から算出した今の時刻の衛星の位置・速度です。
戻り値の前半の3つは、$ \frac{d\boldsymbol{r}}{dt} = \boldsymbol{v} $ の右辺値の計算結果に対応し、ここでは速度 v が結果となるので、入力xのvx, vy, vzをそのまま入れます。
戻り値の後半の3つは、$ \frac{d\boldsymbol{v}}{dt} = -G\frac{M}{r^3} \boldsymbol{r} $ の右辺値の計算結果に対応し、入力xのx,y,zから計算した結果を入れます。
初期値を定義し微分方程式を解く¶
# 微分方程式の初期条件
x0 = [10000, 0, 0, 0, 7, 0] # 位置(x,y,z) + 速度(vx, vy, yz)
t = np.linspace(0, 86400, 1000) # 1日分の軌道伝播
# 微分方程式の数値計算
sol = odeint(func, x0, t)
初期位置 (10000, 0, 0)、初期速度 (7, 0, 0) としています。
時刻は、開始=0、終了=86400で、1000個の時刻データを生成しています。
odeint関数に、微分方程式の漸化式を定義した関数と、初期状態、時刻を渡します。
odeint関数の戻り値は、時刻とその時刻の各状態の配列で、$ [[t_i], [x_i, y_i, z_i, vx_i, vy_i, vz_i]] $ となっています。
微分方程式の数値解のプロット¶
グラフに数値解のx,y平面のプロットをします。
plt.plot(sol[:, 0], sol[:, 1], 'r')
plt.grid() # 格子を描く
plt.gca().set_aspect('equal')
plt.show()
matplotlibのpyplotライブラリの関数plotに、x軸、y軸の数値列を指定、プロットの色と種類に'r' (赤色を指定、形状はデフォルトのsolid line）を指定しています。
横軸は、すべての時刻における状態の1つ目の値（sol[:, 0]）、縦軸は、すべての時刻における状態の2つ目の値（sol[:, 1]）です。
地球のプロット¶
グラフに、地球を円として描画します。
matplotlibには、大きく分けて2種類のインタフェース（API）があります。
- pyplotインタフェース
- axインタフェース
前者は、matplotlibライブラリの基となった MATLAB に近いインタフェースです。後者は、Pythonにおけるオブジェクト指向インタフェースです。
グラフ描画の細かな調整をするにはオブジェクト指向インタフェースが必要になるので、最初からaxインタフェースを理解するのがよいようです。
ですが、簡便のため今回は 前者のインタフェースで地球を描画します。
earth = plt.Circle((0, 0), 6400)
plt.gca().add_patch(earth)
plt.Circle関数で、中心を(0,0)、半径を6400とした円を作成します。plt.Circleは、実際にはmatplotlib.patches.Circleへの
plt.gca()関数で現在描画対象のAxisを取得し、add_patch関数でそこに図形を貼ります。
add_patch関数と同様の機能を持つadd_artist関数もあります。add_artist関数で貼った図形は、グラフの自動スケール対象外となります。