前回、ある時刻での位置と速度から天体の軌道要素を具体的に計算する方法について説明しました。今回は、軌道要素が与えられたときに、ある時刻での位置と速度を求める方法について説明します。
このような問題は逆問題と呼ばれ、結論を示すと以下のような関係となります。
軌道要素 $a,e,t_{\pi},i,\Omega,\omega$ が与えられているとき、軌道面に設定した直交座標系の単位ベクトル $\B{p},\B{q},\B{w}$ は以下のように与えられる。なお、$\B{i},\B{j},\B{k}$ を基準面に設定した直交座標系の単位ベクトルとする。
$$
\left\{
\begin{split}
&\B{p}= (\cos\omega\cos\Omega-\sin\omega\cos i \sin\Omega)\B{i}-(\cos\omega\sin\Omega+\sin\omega\cos i\cos\Omega)\B{j}+\sin\omega\sin i\,\B{k}\EE
&\B{q}=(\sin\omega\cos\Omega+\cos\omega\cos i\sin\Omega)\B{i}-(\sin\omega\sin\Omega-\cos\omega\cos i\cos\Omega)\B{j}-\cos\omega\sin i\,\B{k} \EE
&\B{w}=\sin i\sin\Omega\B{i}+\sin i\cos\Omega\B{j}+\cos i\,\B{k} \EE
\end{split}
\right.
$$
また、ある時刻 $t$ での軌道上の位置ベクトル $\B{r}$ と速度ベクトル $\B{v}$ は以下のように与えられる。
$$
\left\{
\begin{eqnarray}
&&\B{r}=a(\cos E-e)\B{p}+a\sqrt{1-e^2}\sin E\B{q}\EE
&&\B{v}=-\sqrt{\ff{\mu}{a}}\ff{\sin E}{1-e\cos E}\,\B{p}+\sqrt{\ff{\mu(1-e^2)}{a}}\ff{\cos E}{1-e\cos E}\B{q}
\end{eqnarray}
\right.
$$
ただし、$\mu$ を重力定数、$E$ を離心近点離角として、$E$ は以下のケプラー方程式を満たす。
\begin{split}
\sqrt{\ff{\mu}{a^3}}(t-t_{\pi})=E-e\sin E\\
\,
\end{split}
上の表式から分かるように、位置と速度を求める際は離心近点離角 $E$ が重要な働きをします。上の関係を導く準備として、まずは離心近点離角と軌道要素との関係を導いていきます。
離心近点離角と真近点離角の関係
最初の準備として、離心近点離角 $E$ と真近点離角 $\nu$（ニュー）との関係を導くことを考えます。
まず、軌道上の天体が $\B{r}$ の位置にある状況について考えましょう。このとき、軌道面に設定した直交座標系の単位ベクトル $\B{p},\B{q}$ を用いて、$\B{r}$ が以下のように表示できます。
\begin{eqnarray}
\B{r}=r\cos\nu\,\B{p}+r\sin\nu\,\B{q}\tag{1}
\end{eqnarray}
さて、半長軸を $a$、離心率を $e$ とすると動径 $r$ は軌道方程式により以下のように表示できます。
\begin{eqnarray}
r=\ff{a(1-e^2)}{1+e\cos\nu}\tag{2}
\end{eqnarray}
ところで、ケプラー方程式の導出過程で現れた以下の式がありましたが、（※ $P$ はラプラスベクトルの大きさ）
\begin{split}
r=-\left(\ff{\mu}{2\eps}-\ff{P}{2\eps}\cos E \right)=-\ff{\mu}{2\eps}\left(1-\ff{P}{\mu}\cos E \right)
\end{split}
今、$a=\DL{-\ff{\mu}{2\eps}},\,e=\ff{P}{\mu}$ の関係が成り立つので、
\begin{eqnarray}
r=a(1-e\cos E)\tag{3}
\end{eqnarray}
という動径についての表示が得られます。これらの動径についての式$(2)$と$(3)$を等値すると、真近点離角と離心近点離角の関係、
\begin{split}
\cos\nu=\ff{\cos E-e}{1-e\cos E}
\end{split}
が得られます。$\sin\nu$ の関係については次節にて導出を行います。
離心近点離角による動径の表示
ここでの目的は、離心近点離角 $E$ を用いて動径の極座標表示を得ることです。そのために、$\sin\nu$ を求める必要があります。これを得る準備として式$(3)$の両辺を時間微分し、
\begin{split}
\ff{\diff r}{\diff t}=ae\sin E\ff{\diff E}{\diff t}
\end{split}
次に、ケプラー方程式の両辺を時間時微分すると、
\begin{eqnarray}
&&\sqrt{\ff{\mu}{a^3}}\diff t=(1-e\cos E)\diff E \EE
&&\therefore\,\ff{\diff E}{\diff t}=\sqrt{\ff{\mu}{a^3}}\ff{1}{1-e\cos E}\tag{4}
\end{eqnarray}
これを最初の $\DL{\ff{\diff r}{\diff t}}$ の式に適用すると、
\begin{split}
\ff{\diff r}{\diff t}=\ff{e\sqrt{\mu a}\sin E}{1-e\cos E}
\end{split}
が得られます。もう一度軌道方程式について考え、これの両辺を時間微分すると、
\begin{split}
\ff{\diff r}{\diff t}&=\ff{ep\sin\nu}{(1+e\cos \nu)^2}\ff{\diff \nu}{\diff t}\EE
&=\ff{e\sin\nu}{p}\cdot r^2\ff{\diff \nu}{\diff t}
\end{split}
これにケプラーの第二法則の証明過程で得られた角運動量 $h$ に関する関係、$r^2\dot{\nu}=h$ を適用すると、
\begin{split}
\ff{\diff r}{\diff t}&=\ff{e\sin\nu}{p}h=\sqrt{\ff{\mu}{p}}e\sin\nu
\end{split}
なお、$p=\DL{\ff{h^2}{\mu}}$ を用いています。最後にここまでで得られた $\DL{\ff{\diff r}{\diff t}}$ の関係を等値すると、
\begin{split}
\sin\nu=\ff{\sin E}{1-e\cos E}\sqrt{\ff{p}{a}}
\end{split}
今、$p=a(1-e^2)$ の関係にあるため、
\begin{split}
\sin\nu=\ff{\sqrt{1-e^2}\sin E}{1-e\cos E}
\end{split}
が成立します。以上より、動径ベクトルの離心近点離角を用いた極座標表示が以下のようにできることが分かります。
\begin{split}
\B{r}=r\left( \ff{\cos E-e}{1-e\cos E} \right)\B{p}+r\left( \ff{\sqrt{1-e^2}\sin E}{1-e\cos E} \right)\B{q}
\end{split}
これに再び式$(3)$を適用すると、
\begin{eqnarray}
\B{r}=a(\cos E-e)\B{p}+a\sqrt{1-e^2}\sin E\B{q}\tag{5}
\end{eqnarray}
が得られます。
軌道要素からの位置と速度の推算
いよいよ今回の本題である、軌道要素 $a,e,t_{\pi},i,\Omega,\omega$ が与えられているとき、ある時刻での軌道上の位置と速度を求めていきます。
今、軌道面に対して直交座標系を設定して、これの単位ベクトルを $\B{p},\B{q},\B{w}$ とします。このとき、軌道面の各方向の単位ベクトルの関係は、回転行列を用いることでこちらのように計算できることを用いて、
$$
\left\{
\begin{split}
&\B{p}= (\cos\omega\cos\Omega-\sin\omega\cos i \sin\Omega)\B{i}-(\cos\omega\sin\Omega+\sin\omega\cos i\cos\Omega)\B{j}+\sin\omega\sin i\,\B{k}\EE
&\B{q}=(\sin\omega\cos\Omega+\cos\omega\cos i\sin\Omega)\B{i}-(\sin\omega\sin\Omega-\cos\omega\cos i\cos\Omega)\B{j}-\cos\omega\sin i\,\B{k} \EE
&\B{w}=\sin i\sin\Omega\B{i}+\sin i\cos\Omega\B{j}+\cos i\,\B{k} \EE
\end{split}
\right.
$$
さらに $a,e,t_{\pi}$ および、時刻 $t$ が与えられたなら、この時刻での位置と速度を求めることができます。まず、離心近点離角 $E$ はケプラー方程式を用いて
\begin{split}
\sqrt{\ff{\mu}{a^3}}(t-t_{\pi})=E-e\sin E
\end{split}
の関係にあり、この式より $E$ の値が求められます。（ニュートンラプソン法などの数値計算により $E$ は求められます。）
さて、得られた $E$ と式$(5)$より、動径ベクトル $\B{r}$ が
\begin{eqnarray}
\B{r}=a(\cos E-e)\B{p}+a\sqrt{1-e^2}\sin E\B{q}
\end{eqnarray}
と求められます。次に、速度ベクトル $\B{v}$ については両辺を時間微分すれば得られ、
\begin{eqnarray}
\B{v}=\ff{\diff \B{v}}{\diff t}=-a\sin E\cdot \ff{\diff E}{\diff t}\B{p}+a\sqrt{1-e^2}\cos E\cdot\ff{\diff E}{\diff t}\B{q}
\end{eqnarray}
これに、式$(4)$の結果を適用すると、
\begin{eqnarray}
\B{v}=-\sqrt{\ff{\mu}{a}}\ff{\sin E}{1-e\cos E}\,\B{p}+\sqrt{\ff{\mu(1-e^2)}{a}}\ff{\cos E}{1-e\cos E}\B{q}
\end{eqnarray}
が得られます。以上より、軌道要素が与えられたときに、ある時刻 $t$ での位置と速度を以下のように表せます。
軌道要素 $a,e,t_{\pi},i,\Omega,\omega$ が与えられているとき、軌道面に設定した直交座標系の単位ベクトル $\B{p},\B{q},\B{w}$ は以下のように与えられる。なお、$\B{i},\B{j},\B{k}$ を基準面に設定した直交座標系の単位ベクトルとする。
$$
\left\{
\begin{split}
&\B{p}= (\cos\omega\cos\Omega-\sin\omega\cos i \sin\Omega)\B{i}-(\cos\omega\sin\Omega+\sin\omega\cos i\cos\Omega)\B{j}+\sin\omega\sin i\,\B{k}\EE
&\B{q}=(\sin\omega\cos\Omega+\cos\omega\cos i\sin\Omega)\B{i}-(\sin\omega\sin\Omega-\cos\omega\cos i\cos\Omega)\B{j}-\cos\omega\sin i\,\B{k} \EE
&\B{w}=\sin i\sin\Omega\B{i}+\sin i\cos\Omega\B{j}+\cos i\,\B{k} \EE
\end{split}
\right.
$$
また、ある時刻 $t$ での軌道上の位置ベクトル $\B{r}$ と速度ベクトル $\B{v}$ は以下のように与えられる。
$$
\left\{
\begin{eqnarray}
&&\B{r}=a(\cos E-e)\B{p}+a\sqrt{1-e^2}\sin E\B{q}\EE
&&\B{v}=-\sqrt{\ff{\mu}{a}}\ff{\sin E}{1-e\cos E}\,\B{p}+\sqrt{\ff{\mu(1-e^2)}{a}}\ff{\cos E}{1-e\cos E}\B{q}
\end{eqnarray}
\right.
$$
ただし、$\mu$ を重力定数、$E$ を離心近点離角として、$E$ は以下のケプラー方程式を満たす。
\begin{split}
\sqrt{\ff{\mu}{a^3}}(t-t_{\pi})=E-e\sin E\\
\,
\end{split}
位置と速度の導出
では、実際に人工衛星のある時刻での位置と速度を計算してみましょう。
ここでは、$e=0.1,i=40^{\circ},\Omega=50^{\circ},\omega=45^{\circ},a=42,000\,\RM{km}$ であったとします。
ここでは簡単のため、近心点を通過した時刻 $t_{\pi}$ を $0$ とします。なお、重力定数を $\mu=4.0\times 10^5\,\RM{km^3/s}$ とします。
今回は、近心点を通過して $2$ 時間後の位置と速度を求めてみます。
まず、離心近点離角 $E$ の導出を行います。今、$t=7,200\,\RM{s}$ と置けるため、ケプラー方程式は下のようになります。
\begin{split}
\sqrt{\ff{4.0\times 10^5}{42,000^3}}(7,200-0)=E-0.1\sin E
\end{split}
これを満たす $E$ をニュートン・ラプソン法により数値計算すると、$E\NEQ33.5^{\circ}$ と求められます。今、軌道面に置かれた単位ベクトルは次のように計算できるので、
$$
\left\{
\begin{split}
&\B{p}=0.04\B{i}-0.89\B{j}+0.45\B{k}\EE
&\B{q}=0.87\B{i}-0.19\B{j}-0.45\B{k}\EE
&\B{w}=0.49\B{i}+0.49\B{j}+0.77\B{k}
\end{split}
\right.
$$
これより、$\B{r}$ と $\B{v}$ が以下のように求められます。
$$
\left\{
\begin{split}
&\B{r}=37550\B{i}-35510\B{j}-4980\B{k}\,\RM{km}\EE
&\B{v}=2.36\B{i}+1.11\B{j}-2.11\B{k}\,\RM{km/s}
\end{split}
\right.
$$
これより、$r=51930\,\RM{km},v=3.36\,\RM{km/s}$ が得られます。
