楕円軌道上の時刻ごとの位置を数値的に求めるための下ごしらえ

楕円軌道は $r(\phi)$ で表されるので,$\phi$ を $t$ の関数として数値的に求めるための下ごしらえ。

万有引力の2体問題の運動方程式を数値的に解く前の下ごしらえ」では,万有引力の2体問題(は結局1体問題に帰着するんだけど)の運動方程式を数値的に解く前の下ごしらえをしていた。実際に Maxima で解くのは「下ごしらえした万有引力の2体問題の運動方程式を Maxima で数値的に解く」に書いている。

連立2階常微分方程式としての(規格化された)運動方程式を数値的に解くのも教育的ではあるが,この問題は楕円軌道になることがわかっている。ただし,解が

$$r = \frac{a (1 – e^2)}{1 + e \cos\phi}$$

と書けることはわかっているのだが,$r$ も $\phi$ も時間 $t$ の関数としてはあからさまに与えられていないため,ある時刻に楕円軌道上のどこにいるかは,数値的に解いてみないとわからない。

時刻 $t$ の関数として楕円軌道上の天体の位置を知るためには,ケプラー方程式を解く必要があるが,このあたりについては,カテゴリー「ケプラー方程式」にいくつか書いた。

 

ここでは,楕円軌道は $r(\phi)$ で表されるので,$\phi$ を $t$ の関数として数値的にもとめ,$\phi(t)$ が数値的に求まれば,

$$x(t) = r(\phi(t))\,\cos\phi(t), \quad  y(t) = r(\phi(t))\, \sin\phi(t)$$

のように,時刻を与えれば楕円軌道上の位置がわかるようにしてみる。以下,「参考:ニュートン力学における万有引力の2体問題」から要点を抜粋して利用する。

万有引力の2体問題は1体問題に帰着する

万有引力の2体問題は1体問題に帰着し,運動方程式は

$$\frac{d^2\boldsymbol{r}}{dt^2}  = – \frac{GM \boldsymbol{r}}{r^3}$$

保存量(運動の定数)

運動方程式から,2つの保存量(運動の定数)が得られる。1つは,(単位質量あたりの)角運動量の $z$ 成分に相当する $\ell$:

$$\ell \equiv r^2 \frac{d\phi}{dt}$$

もう1つは(単位質量あたりの)力学的エネルギーに相当する $\varepsilon$:

$$\varepsilon \equiv \frac{1}{2} \left\{\left(\frac{dr}{dt}\right)^2 + \frac{\ell^2}{r^2}\right\} – \frac{GM}{r}$$

解:楕円軌道

$$r = \frac{a(1-e^2)}{1 + e \cos\phi}$$

ここで,軌道長半径 $a$,離心率 $e$ と保存量の関係は

$$\ell = \sqrt{GM a (1-e^2)}, \quad \varepsilon = – \frac{GM}{2a}$$

解くべき方程式

基本的に以下の1階常微分方程式を解くだけでよい。

$$\frac{d\phi}{dt} = \frac{\ell}{r^2}, \quad r = \frac{a(1-e^2)}{1 + e \cos\phi}$$

無次元化

ケプラーの第3法則

$$\frac{a^3}{P^2} = \frac{GM}{4\pi^2}, \quad \therefore\ \ P = \frac{2\pi a^2}{\sqrt{GM a}}$$

軌道周期 $P$ を使って無次元化された時間 $T$ を

$$T \equiv \frac{t}{P}$$

とすると,

$$\frac{d\phi}{dT} = 2 \pi \frac{(1+e \cos\phi)^2}{(1-e^2)^{3/2}} \tag{1}$$

これを数値的に解き,$\phi(T)$ がわかったら,

\begin{eqnarray}
R(T) &\equiv& \frac{r}{a} = \frac{1-e^2}{1+e\cos\phi(T)} \tag{2}\\
X(T) &\equiv& R(T) \cos\phi(T) \tag{3}\\
Y(T) &\equiv& R(T) \sin\phi(T) \tag{4}
\end{eqnarray}

として,規格化された $X$ 座標,$Y$ 座標を求めればよい。