二体問題を解くと軌道方程式と呼ばれる数式を導出できます。
$p$ を定数、離心率を $e$ として軌道方程式は次のように記述される。
\begin{eqnarray}
r = \frac{p}{1+e\cos\nu} \\
\,
\end{eqnarray}
この軌道方程式に離心率を代入すると、軌道の形を決定できます。
応用上、未来のある時刻 $t$ での天体の位置を計算する方が重要です。
例えば、日食や月食の起きる時期や人工衛星のある時刻での位置を知りたい場合などです。
その際に利用される計算式がケプラー方程式です。
平均近点離角を$M$、離心近点離角を$E$、離心率を $e$ として
ケプラー方程式は次のように表される。
\begin{eqnarray}
M = E \,-\, e\sin E \\
\,
\end{eqnarray}
今回はケプラー方程式の導出と、ケプラー方程式を利用したある時刻での天体の位置を計算する手法について解説します。
ケプラー方程式を解くことで、未来のある時刻での天体の位置を予測できるようになります。
※ ベクトルや時間微分の表記についてはこちらで解説しています。
ケプラー方程式導出の準備
エネルギー積分の変形
エネルギー積分より、動径ベクトル $\B{r}$、重力定数 $\mu$（ミュー）、全エネルギー $\varepsilon$（エプシロン）を使って次のように表されます。（$\eps<0$であることに注意してください。）
\begin{eqnarray}
\ff{1}{2}(\dot{\B{r}}\cdot\dot{\B{r}})= \ff{\mu}{r} + \eps \tag{1} \\
\end{eqnarray}
次に $\dot{\B{r}}\cdot\dot{\B{r}}$ について考えます。
ところで、角運動量 $\B{h}$ は $\B{r}\times\dot{\B{r}}$ と表せます。そこで、角運動量の内積をとると、
\begin{eqnarray}
\B{h}\cdot\B{h} &=& h^2 \EE
&=& (\B{r}\times\dot{\B{r}})\cdot(\B{r}\times\dot{\B{r}}) \EE
&=& (\dot{\B{r}}\cdot\dot{\B{r}})(\B{r}\cdot\B{r})-(\dot{\B{r}}\cdot\B{r})(\B{r}\cdot\dot{\B{r}}) \\
\end{eqnarray}
となります。さらに、
\begin{eqnarray}
\ff{\diff}{\diff t}r^2 &=& \ff{\diff}{\diff t}(\B{r}\cdot\B{r}) \EE
2r\dot{r} &=& 2(\B{r}\cdot\dot{\B{r}}) \EE
\therefore\,\,\,r\dot{r} &=& (\B{r}\cdot\dot{\B{r}})
\end{eqnarray}
の関係があるので、
\begin{eqnarray}
\B{h}\cdot\B{h} &=& (\dot{\B{r}}\cdot\dot{\B{r}})(\B{r}\cdot\B{r})-(\dot{\B{r}}\cdot\B{r})(\B{r}\cdot\dot{\B{r}}) \EE
h^2 &=& (\dot{\B{r}}\cdot\dot{\B{r}})r^2-r^2\dot{r}^2 \EE
\therefore \dot{\B{r}}\cdot\dot{\B{r}} &=& \left(\ff{\diff r}{\diff t}\right)^2 + \ff{h^2}{r^2}
\end{eqnarray}
とできます。
※ 上付きドットは時間微分を表します。詳しくはこちらで解説しています。
この結果を式(1)に代入すると、
\begin{eqnarray}
\ff{1}{2}\left(\ff{\diff r}{\diff t}\right)^2 &=& \ff{\mu}{r} + \eps-\ff{h^2}{2r^2} \EE
\ff{\diff r}{\diff t} &=& \pm\sqrt{2\eps + \ff{2\mu}{r}-\ff{h^2}{r^2} } \tag{2}
\end{eqnarray}
と整理されます。
ここで、式(2)を $\diff t$ について変形すると、
\begin{eqnarray}
\diff t &=& \pm\ff{\diff r}{\sqrt{2\eps + \ff{2\mu}{r}\,-\, \ff{h^2}{r^2}}} \\[3pt]
\therefore\,\,\diff t &=& \pm\ff{r\,\diff r}{\sqrt{2\eps r^2 + 2r\mu \,-\, h^2}} \tag{3}
\end{eqnarray}
となり微分方程式が得られます。
微分方程式の変形
式(3)の微分方程式を計算したいのですが、このままでは積分が難しいため、
右辺の分母を簡単にすることを考えます。
分母を平方完成により次のように変形します。
\begin{eqnarray}
2\eps r^2 + 2r\mu \,-\, h^2 &=& 2\eps\left( r^2 + \ff{\mu}{\eps}r \,-\, \ff{h^2}{2\eps} \right) \\[3pt]
&=& 2\eps\left\{ \left( r+\ff{\mu}{2\eps} \right)^2 \,-\, \left(\ff{\mu}{2\eps}\right)^2 \,-\, \ff{h^2}{2\eps} \right\} \EE
\end{eqnarray}
ところで、角運動量とエネルギーの関係から、
\begin{split}
&\ff{h^2}{2\eps} = -\left(\ff{\mu}{2\eps} \right)^2\left( 1\,-\,\ff{P^2}{\mu^2} \right)^2 \EE
\therefore\,\,\,&\left(\ff{\mu}{2\eps}\right)^2-\ff{h^2}{2\eps} = \left( \ff{P}{2\eps} \right)^2
\end{split}
と変形できます。
※ この関係式についての導出の詳細はこちらで解説しています。
この結果を先ほどの式に適用すると、
\begin{split}
&-2\eps\left\{\left(\ff{\mu}{2\eps}\right)^2 + \ff{h^2}{2\eps}-\left( r+\ff{\mu}{2\eps} \right)^2 \right\} \EE
=&-2\eps\left\{\left( \ff{P}{2\eps} \right)^2 \,-\, \left( r+\ff{\mu}{2\eps} \right)^2 \right\} \EE
\end{split}
とできます。
この結果を式(3)に代入して整理すると、
\begin{eqnarray}
\pm\sqrt{-2\epsilon}\,\diff t &=& \ff{r\,\diff r}{\sqrt{ \left( \ff{P}{2\epsilon} \right)^2 \,-\, \left( r+\ff{\mu}{2\epsilon} \right)^2}} \tag{4} \EE
\end{eqnarray}
が得られます。式(4)について積分を行いますが、このままでは厳しいため、変数変換を行い三角関数の形に持ちこみます。
ケプラー方程式の導出
式(4)の見通しを良くするため、変数変換を二回実施します。
まず、$\DL{z = r+\ff{\mu}{2\eps}}$ とすると、$\diff r = \diff z$ とできるため式(4)は、
\begin{eqnarray}
\pm\sqrt{-2\eps}\,\diff t &=& \ff{z-\ff{\mu}{2\eps}}{\sqrt{ \left( \ff{P}{2\eps} \right)^2 \,-\, z^2}} \diff z \EE
&=& \ff{2\eps}{P}\ff{z-\ff{\mu}{2\eps}}{\sqrt{ 1\,-\,\left( \ff{2\eps z}{P} \right)^2}} \diff z \tag{5}
\end{eqnarray}
と変形できます。
次に、$ \DL{\ff{2\eps z}{P} = \cos E} $ とすると、$\DL{\diff z = -\ff{P}{2\eps}\sin E \,\diff E} $ とできるので、式(5)は、
\begin{eqnarray}
\pm\sqrt{-2\eps}\,\diff t &=& \ff{2\eps}{P}\ff{z-\ff{\mu}{2\eps}}{\sqrt{ 1\,-\,\left( \ff{2\eps z}{P} \right)^2}} \diff z \\[5pt]
&=& -\ff{2\eps}{P}\ff{\ff{P}{2\eps}\cos E \,-\, \ff{\mu}{2\eps}}{\sqrt{1-\cos^2 E}}\cdot \ff{P}{2\eps}\sin E \,\diff E \\[5pt]
&=& \left(\ff{\mu}{2\eps} \,-\, \ff{P}{2\eps}\cos E \right) \diff E
\end{eqnarray}
と変数変換できます。
見通しの良い形に整理できたので、微分方程式を解いて行きましょう。
離心近点離角とは？
変数変換の過程で突然、角度 $E$ なるものが登場しました。
$E$ を図示すると以下のようになります。
楕円軌道に外接する円を描き、楕円軌道と動径ベクトルが交わる点$\RM{P}$から$x$軸に下した垂線と外接円との交点を$\RM{Q}$とします。
この線分 $\RM{OQ}$ と $x$ 軸の成す角が $E$ となります。
天体力学では $E$ を離心近点離角と呼びます。
動径ベクトル $\B{r}$ と、 $x$軸 の成す角 $\nu$（真近点離角）の方が物理的な意味を理解しやすいですが、
計算上は不便なため、離心近点離角 $E$ を使う場面があります。
では、先述の微分方程式を積分して、$E$ の表式を得ましょう。
さて、動径ベクトルが近点から$\RM{P}$まで動く間に、離心近点離角が $E’$ まで増加したとします。
また、時刻 $t_0$ にて近点にいて、時刻 $t’$ で$\RM{P}$に来たとします。
このとき、近点から$\RM{P}$に至る時間は式(5)を積分することにより計算できます。
具体的には、次のように計算できます。
\begin{eqnarray}
\pm\int_{t_0}^{t’} \sqrt{-2\eps}\,\diff t &=& \int_{0}^{E’} \left(\ff{\mu}{2\eps} \,-\, \ff{P}{2\eps}\cos E \right) \diff E \EE
\pm\sqrt{-2\eps}\,(t’ \,-\, t_0) &=& \ff{\mu}{2\eps}E’-\ff{P}{2\eps}\sin E’ \EE
\end{eqnarray}
両辺を $\DL{\ff{\mu}{2\eps}}$ で割ると、
\begin{eqnarray}
\pm\ff{2\eps}{\mu}\sqrt{-2\eps}\,(t’ \,-\, t_0) &=& E’ \,-\, \ff{P}{\mu}\sin E’ \EE
\end{eqnarray}
となります。
さらに、$\DL{a=-\ff{\mu}{2\eps}}$, $\DL{e=\ff{P}{\mu}}$ であることを利用し、符号が正の場合を採用すると
\begin{eqnarray}
\sqrt{\ff{\mu}{a^3}}\,(t’ \,-\, t_0) &=& E’ \,-\, e\sin E’ \EE
\end{eqnarray}
とできます。以上より、$\DL{M=\sqrt{\ff{\mu}{a^3}}(t’ \,-\, t_0)}$として$E=E’$とすると、
\begin{eqnarray}
M &=& E-e\sin E
\end{eqnarray}
と整理できます。この方程式をケプラー方程式と呼びます。
平均近点離角を$M$、離心近点離角を$E$、離心率を $e$ として
ケプラー方程式は次のように表される。
\begin{eqnarray}
M &=& E \,-\, e\sin E \EE
\,
\end{eqnarray}
以上、ケプラー方程式が得られました。
ところで、軌道を一回転したとき $E=2\pi$ となります。この周期を $T=(t’-t_0)$ とすると、
\begin{split}
&\sqrt{\ff{\mu}{a^3}}\,T = 2\pi \EE
&T = 2\pi\sqrt{\ff{a^3}{\mu}} \EE
\therefore\,\,\,&\ff{T^2}{a^3} = \ff{2\pi}{\mu} = const.
\end{split}
となり、確かにケプラーの第三法則が導出できることが分かります。
平均近点離角とは？
ケプラー方程式の左辺に現れる$M$の意味について考えます。
まず、$\DL{n=\sqrt{\ff{\mu}{a^3}}}$ とすると軌道を一周する時間（＝周期）は $T$ とできるので、
\begin{eqnarray}
nT &=& 2\pi \EE
\therefore n&=& \ff{2\pi}{T}
\end{eqnarray}
が成り立ちます。
$n$は角速度に相当するもので、天体力学では平均運動と呼ばれます。
天体が軌道を一周する時間を $T$ として、平均運動を次のように定義する。
\begin{eqnarray}
n&=& \ff{2\pi}{T} \\
\,
\end{eqnarray}
改めて $t=t’-t_0$ とすると、$M=nt$ とおけて、これを図示すると次のようになります。
なお、ケプラー方程式は平均近点離角$M$と離心近点離角$E$との関係を与える方程式といえます。
したがって、楕円軌道の離心率 $e$ と離心近点角が分かれば、平均近点離角が分かると言えます。
実用上、将来の時刻における天体の位置を計算したいという要請があります。
そのためには、求めたい時刻での離心近点離角 $E$ を知る必要があります。
しかし、観測データから平均運動$n$と離心率$e$は計算できる一方、将来の時点の$E$は観測からは算出できません。
そのため、$E$ については、ケプラー方程式方程式を解いて算出する必要があります。
$M$ と $e$ が与えられたとき、果たしてケプラー方程式は解けるでしょうか？
ケプラー方程式の解法
最古の超越方程式
結論からいうと、ケプラー方程式は解析的には解けません。（このような方程式を超越方程式と呼びます。）
と言っても、話が先に進まないので近似的に $E$ を計算する方法について考えましょう。
次のようにケプラー方程式を変形します。
\begin{eqnarray}
E \,-\, M = e\sin E \tag{6} \\
\end{eqnarray}
ケプラー方程式の右辺は近点（$E=0$）で$M=0$、遠点（$E=\pi$）で$M=\pi$となります。
また、式(6)の右辺が三角関数であるため、$E-M$は周期$2\pi$の周期関数と見なすことができます。
\begin{array}{c|ccccc}
E & 0 & \cdots & \pi & \cdots & 2\pi \\\hline
M & 0 & \cdots & \pi & \cdots & 2\pi \\\hline
E\,-M & 0 & \cdots & 0 & \cdots & 0
\end{array}
このとき、次のような関数を考えます。
\begin{eqnarray}
\sum_{l=1}^{\infty}A_l \sin lM = A_1 \sin M + A_2 \sin 2M + \cdots
\end{eqnarray}
この関数は、$M=0$ と $M=\pi$ で $0$ になる周期関数ですので、式(6)の右辺として適用できます。（フーリエ級数の考え方です。）
すると、
\begin{eqnarray}
E \,-\, M = \sum_{l=1}^{\infty}A_l \sin lM \tag{7} \\
\end{eqnarray}
とケプラー方程式を変形でき、$E$ を $M$ で表現できそうな見通しが立ちました。
このままでは、係数の $A_l$ が分からないので $A_l$ の具体的な表示を考えていきます。
ベッセル関数表示
式(7)の右辺の無限級数を解消して係数の計算を行います。
少々、トリッキーな変形を行いますが三角関数の性質（直交性）を巧みに利用して計算を行います。
まず、式(7)の両辺を $M$ に関して微分すると、
\begin{eqnarray}
\ff{\diff E}{\diff M} \,-\, 1 = \sum_{l=1}^{\infty}lA_l \cos lM
\end{eqnarray}
となり、両辺に $\cos kM$ を掛けると、
\begin{eqnarray}
\left(\ff{\diff E}{\diff M} \,-\, 1\right) \cos kM = \sum_{l=1}^{\infty}lA_l \cos lM \cos kM \EE
\end{eqnarray}
となります。
そして、$0$ から $\pi$ までの区間で積分を行うと、
\begin{eqnarray}
\int_{0}^{\pi}\left(\ff{\diff E}{\diff M} \,-\, 1\right) \cos kM \,\diff M &=& \int_{0}^{\pi} \sum_{l=1}^{\infty}lA_l \cos lM \cos kM \, \diff M \EE
&=& \sum_{l=1}^{\infty}lA_l \int_{0}^{\pi} \cos lM \cos kM \, \diff M \tag{8} \EE
\end{eqnarray}
となります。
右辺は三角関数の直交性より、
\begin{eqnarray}
\int_{0}^{\pi} \cos lM \cos kM \, \diff M =
\begin{cases}
0 & (l \neq k) \\
\ff{\pi}{2} & (l = k)
\end{cases}
\end{eqnarray}
とできて、式(8)は、
\begin{eqnarray}
\ff{\pi k}{2}A_k &=& \int_{0}^{\pi}\left(\ff{\diff E}{\diff M} \,-\, 1\right) \cos kM \,\diff M \EE
&=& \int_0^{\pi} \ff{\diff E}{\diff M} \cos kM \,\diff M \,-\, \int_0^{\pi} \cos kM \,\diff M \EE
&=& \int_0^{\pi} \cos kM \,\diff E \EE
&=& \int_0^{\pi} \cos k(E \,-\, e\sin E) \,\diff E \EE
\therefore A_k &=& \ff{2}{\pi k} \int_{0}^{\pi} \cos (kE \,-\, ek \sin E) \diff E
\end{eqnarray}
となります。
ところで、数学の世界にはベッセル関数と呼ばれる特殊関数があります。
ベッセル関数は次のように表されます。（$k$が自然数の場合）
$k$ を自然数として、ベッセル関数は次のように定義される。
\begin{eqnarray}
J_k(x) &=& \ff{1}{\pi} \int_0^{\pi} \cos (kE \,-\, x \sin E) \diff E \\[5pt]
&=& \sum_{s=0}^{\infty} \ff{(-1)^s}{s!(k+s)!}\left( \ff{x}{2} \right)^{k+2s} \\
\,
\end{eqnarray}
※ ベッセル関数についての詳しい解説はこちらでしています。
ベッセル関数を使うと $A_k$ は、
\begin{eqnarray}
A_k &=& \ff{2}{k}J_k(ek) \,\,\,\, (x=ek)
\end{eqnarray}
となり、離心近点離角 $E$ を次のように級数表示できます。
\begin{eqnarray}
E &=& M + 2\sum_{k=1}^{\infty} \ff{1}{k} J_k(ek) \sin kM \\
\,
\end{eqnarray}
なかなかに派手な式ですが、上式を計算することで理論上は離心近点離角を計算できます。
しかし、この式を直接計算するのは現実的ではありません。実用上はコンピュータを使い、ニュートン・ラプソン法により離心近点離角を計算します。
天体の未来の位置を計算するためには、相当な計算量が必要になることが実感できます。
→天体力学のまとめ記事はこちら
ベッセル関数やニュートン・ラプソン法については別の記事で詳細な解説を行います。
