今回は天体の軌道要素を具体的に計算する方法について説明します。まず、軌道要素は以下のような関係で結ばれています。
軌道要素 $a,e,t_{\pi},i,\Omega,\omega$ は以下のように表せる。
$$
\left\{
\begin{split}
&a=-\ff{\mu}{2\eps} \EE
&e=\ff{P}{\mu}=\sqrt{1-\ff{h^2}{\mu a}} \EE
&t_{\pi}=t-\sqrt{\ff{a^3}{\mu}}(E-e\sin E) \EE
&i=\cos^{-1}\left(\ff{\B{h}}{h}\cdot \B{k}\right) \EE
&\Omega=\tan^{-1}\left(\ff{\B{w}\cdot\B{i}}{\B{w}\cdot\B{j}} \right) \EE
&\omega=\cos^{-1}\left(\ff{\B{P}}{P}\cdot \B{n}\right)
\end{split}
\right.
$$
ただし、$\mu$ を重力定数、$\eps$ を軌道上の天体が持つ全力学的エネルギー、$P$ をラプラスベクトル(近心点ベクトル)の大きさ、$E$ を離心近点離角、$\B{h}$ を角運動量、$\B{i},\B{j},\B{k}$ を基準面に設定した直交座標系の単位ベクトル、$\B{p},\B{q},\B{w}$ を軌道面に設定した直交座標系の単位ベクトルとする。
まずは、上で示した物理量と軌道要素との関係を導出していきます。
軌道要素と位置と速度の関係
例えば、地球の観測局から人工衛星のある時刻での位置と速度が得られたとき、これらから機体の軌道を導くことは、産業や安全保障上、大きな意味があります。
さて、下図のように人工衛星が地球の上空を周回しているとします。
このとき、ある時刻 $t$ での人工衛星の動径ベクトルが $\B{r}$、速度ベクトルを $\B{v}$ と表せたとします。
始めに、人工衛星の軌道の長半径 $a$ については、こちらで計算した結果を用いることで以下のように表せます。
\begin{split}
a=-\ff{\mu}{2\eps}
\end{split}
ただし、$\eps$ を天体が持つ全力学的エネルギー、$\mu=G(m_1+m_2)$ を重力定数とします。
ところで、角運動量 $\B{h}$ は次のように計算でき、
\begin{split}
\B{h}=\B{r}\times\B{v}
\end{split}
角運動量 $\B{h}$ からは、ラプラスベクトル $\B{P}$(近心点ベクトル)と呼ばれるベクトルが導けます。
\begin{split}
\B{P}=-\ff{\mu}{r}\B{r}-\B{h}\times\B{v}
\end{split}
さらに、離心率 $e$ は軌道方程式の導出過程から分かるように $|\B{P}|$ を用いて次のように表せます。
\begin{split}
e=\ff{P}{\mu}=\sqrt{1-\ff{h^2}{\mu a}}
\end{split}
次に、近心点通過時刻 $t_{\pi}$ は、ある時刻 $t$ での離心近点離角 $E$ を用いると、ケプラー方程式から次のような関係を導けます。
\begin{split}
&\sqrt{\ff{\mu}{a^3}}(t-t_{\pi})=E-e\sin E\EE
\therefore&\,t_{\pi}=t-\sqrt{\ff{a^3}{\mu}}(E-e\sin E)
\end{split}
そして、軌道傾斜角 $i$ は軌道面に垂直な単位ベクトル $\B{w}$ と基準面の $z$ 軸方向の単位ベクトル $\B{k}$ の内積を用いて以下のように記述できます。
\begin{split}
\cos i&=\ff{\B{w}\cdot \B{k}}{|\B{w}||\B{k}|}=\B{w}\cdot \B{k}\EE
\therefore\,\,&i=\cos^{-1}(\B{w}\cdot \B{k})
\end{split}
なお、$\B{w}$ は $\B{h}$ と平行なため、$\DL{\B{w}=\ff{\B{h}}{h}}$ と表せます。したがって、
\begin{split}
i=\cos^{-1}\left(\ff{\B{h}}{h}\cdot \B{k}\right)
\end{split}
という関係となります。
昇交点経度 $\Omega$ については、軌道面の各方向の単位ベクトルの関係が回転行列を用いることでこちらのように表せることを用いて、以下のように記述できます。ただし、$\B{i},\B{j}$ を基準面の $x$ 軸、$y$ 軸方向の単位ベクトルとします。
$$
\left\{
\begin{split}
\B{w}\cdot\B{i}&=\sin\Omega\sin i\EE
\B{w}\cdot\B{j}&=\cos\Omega\sin i
\end{split}
\right.
$$
二式を割ると $i$ を消去でき、
\begin{split}
&\tan \Omega=\ff{\B{w}\cdot\B{i}}{\B{w}\cdot\B{j}}\EE
\therefore\,&\Omega=\tan^{-1}\left(\ff{\B{w}\cdot\B{i}}{\B{w}\cdot\B{j}} \right)
\end{split}
と $\Omega$ が求められます。最後に、近心点引数 $\omega$ ですが、これは $\B{P}$ と交線方向の単位ベクトル $\B{n}$ の内積として
\begin{split}
\B{P}\cdot\B{n}=|\B{P}||\B{n}|\cos\omega
\end{split}
とできることを用います。今、$|\B{n}|=1$ ことも用いると、
\begin{split}
\omega=\cos^{-1}\left(\ff{\B{P}}{P}\cdot \B{n}\right)
\end{split}
と求められます。以上の結果をまとめると、ある時刻での位置と速度の観測結果から、軌道要素が次のように記述できることが言えます。
軌道要素 $a,e,t_{\pi},i,\Omega,\omega$ は以下のように表せる。
$$
\left\{
\begin{split}
&a=-\ff{\mu}{2\eps} \EE
&e=\ff{P}{\mu}=\sqrt{1-\ff{h^2}{\mu a}} \EE
&t_{\pi}=t-\sqrt{\ff{a^3}{\mu}}(E-e\sin E) \EE
&i=\cos^{-1}\left(\ff{\B{h}}{h}\cdot \B{k}\right) \EE
&\Omega=\tan^{-1}\left(\ff{\B{w}\cdot\B{i}}{\B{w}\cdot\B{j}} \right) \EE
&\omega=\cos^{-1}\left(\ff{\B{P}}{P}\cdot \B{n}\right)
\end{split}
\right.
$$
ただし、$\mu$ を重力定数、$\eps$ を軌道上の天体が持つ全力学的エネルギー、$P$ をラプラスベクトル(近心点ベクトル)の大きさ、$E$ を離心近点離角、$\B{h}$ を角運動量、$\B{i},\B{j},\B{k}$ を基準面に設定した直交座標系の単位ベクトル、$\B{p},\B{q},\B{w}$ を軌道面に設定した直交座標系の単位ベクトルとする。
軌道要素の計算例題
それでは、上の位置と速度の関係から軌道要素の導出を行いましょう。このとき、人工衛星の位置ベクトルを $\B{r}$、速度ベクトル $\B{v}$ が次のように観測されたとします。
$$
\B{r}=
\begin{pmatrix}
x \\
y \\
z
\end{pmatrix}
=
\begin{pmatrix}
-3000\,\RM{km} \\
-6000\,\RM{km} \\
0
\end{pmatrix}
$$
$$
\B{v}=
\begin{pmatrix}
v_x \\
v_y \\
v_z
\end{pmatrix}
=
\begin{pmatrix}
0 \\
0 \\
10\,\RM{km/s}
\end{pmatrix}
$$
なお、重力定数を $\mu=4.0\times 10^5\,\RM{km^3/s}$ とします。
まず、$r$ と人工衛星の持つ全力学的エネルギー $\eps$ については次のように求められます。
\begin{split}
r&=\sqrt{x^2+y^2+z^2}\EE
&=\sqrt{3300^2+6600^2+0^2}\EE
&\NEQ 6700\,\RM{km} \EE
\eps&=\ff{1}{2}(\B{v}\cdot\B{v})-\ff{\mu}{r}\EE
&=\ff{1}{2}\times 10^2-\ff{4.0\times 10^5}{7380} \EE
&\NEQ -9.6\,\RM{km^2/s^2}
\end{split}
これより、長半径 $a$ が求められます。すなわち、
\begin{split}
a&=-\ff{\mu}{2\eps}=-\ff{4.0\times 10^5}{2\times(-4.2)}\EE
&=20770\,\RM{km}
\end{split}
次に、角運動量 $\B{h}$ は以下のように計算できます。
\begin{split}
\B{h}&=\B{r}\times\B{v}\\
&=
\begin{vmatrix}
\B{i} & \B{j} & \B{k} \\
x & y & 0 \\
0 & 0 & v_z
\end{vmatrix}\\
&=-60000\B{i}+30000\B{j}\,\,\RM{km^2/s}\EE
|\B{h}|&=\sqrt{60000^2+30000^2+0^2}\\
&\NEQ 67080\,\,\RM{km^2/s}
\end{split}
これより、ラプラスベクトル $\B{P}$ が求められます。
\begin{split}
\B{P}&=-\ff{\mu}{r}\B{r}-\B{h}\times\B{v}\EE
&=-121120\,\B{i}-242230\,\B{j}\,\,\RM{km^3/s^2}
\end{split}
また、離心率 $e$ が以下のように求められます。
\begin{split}
e&=\sqrt{1-\ff{h^2}{\mu a}}\EE
&=\sqrt{1-\ff{67080^2}{4.0\times 10^5\times47530 }}\EE
&\NEQ 0.677
\end{split}
近心点通過時刻 $t_{\pi}$ は後ほど計算するとして、$i,\Omega,\omega$ の計算を先に行います。
まず、軌道傾斜角 $i$ については、$\B{h}\cdot\B{k}=0$ であることより、$i=\DL{\ff{\pi}{2}}$ が得られます。そして、昇交点経度 $\Omega$ については $\B{w}$ が
\begin{split}
\B{w}=\ff{\B{h}}{h}=-0.894\B{i}+0.447\B{j}
\end{split}
となるので、
\begin{split}
\Omega&=\tan^{-1}\left(\ff{\B{w}\cdot\B{i}}{\B{w}\cdot\B{j}} \right)\EE
&=\tan^{-1}\ff{-0.894}{0.447}\NEQ -63.4^{\circ}
\end{split}
となり、$\omega$ については、$\B{n}$ が
\begin{split}
\B{n}&=\ff{\B{h}}{h}\times\B{k}\EE
&=-0.707\B{i}+0.707\B{j}
\end{split}
となることより、
\begin{split}
\omega&=\cos^{-1}\left(\ff{\B{P}}{P}\cdot \B{n}\right)\EE
&\NEQ 108^{\circ}
\end{split}
が得られます。最後の $t_{\pi}$ については、
$$
\left\{
\begin{split}
&e\cos E=1-\ff{r}{a}=0.677\EE
&e\sin E=\ff{\B{r}\cdot\B{v}}{\sqrt{\mu a}}=0 \EE
&\tan E=\ff{e\sin E}{e\cos E}=0
\end{split}
\right.
$$
これより、$E=0$ と言え、これを近点通過時刻 $t_{\pi}$ の計算式に適用すると、
\begin{split}
t_{\pi}&=t-\sqrt{\ff{a^3}{\mu}}(0-0.677\times \sin 0)\EE
\therefore\,\,t_{\pi}&=t
\end{split}
となるので、観測時刻がそのまま近点通過時刻となることが分かります。