ケプラー方程式とは
Wikipedia のケプラー方程式の項を読んでも，天文業界でない私には今一つ意味がわからないので，以下のように自分が納得できるように噛み砕いてみた。
ケプラー方程式とは，楕円軌道を媒介変数表示する際のパラメータである離心近点離角 $u$ と時間 $t$ を関係づける式であり，楕円軌道の離心率を $e$，周期を $T$ とすると，以下のような式のことである。
$$ u – e \sin u = \frac{2\pi}{T} t $$
楕円軌道の媒介変数表示
まず，ケプラーの第1法則から，惑星の軌道は楕円軌道である。（平面上にあるのでそれを $xy$ 平面とする。）
楕円の中心を原点としたデカルト座標 $X, Y$ では，長半径 $a$，離心率 $e$ （したがって短半径 $b = a \sqrt{1-e^2}$）の楕円の式は
$$\frac{X^2}{a^2} + \frac{Y^2}{a^2(1-e^2)} = 1$$
である。この楕円は，離心近点離角 $u$ を使って，以下のように媒介変数表示できる。
\begin{eqnarray}
X &=& a \cos u\\
Y &=& a\sqrt{1 – e^2} \sin u
\end{eqnarray}
楕円の中心から $X$ 軸上に $a e$ だけ離れた楕円の焦点を原点としたデカルト座標 $x, y$ で書くと
\begin{eqnarray}
x &=& X – a e = a (\cos u – e) \\
y &=& Y = a\sqrt{1 – e^2} \sin u\\ \ \\
r^2 &=& x^2 + y^2 \\
&=& a^2 (\cos^2 u – 2 e \cos u + e^2) + a^2 (1 – e^2) \sin^2 u \\
&=& a^2 (1 – e \cos u)^2 \\
\therefore\ \ r &=& a (1 – e\cos u)
\end{eqnarray}
ケプラーの第2法則（面積速度一定の法則）
ケプラーの第2法則とは，楕円の焦点に位置する太陽と楕円上を運動する惑星とを結ぶ線分が単位時間に描く面積，すなわち面積速度が一定であることを表す。楕円の面積 $S = \pi a^2 \sqrt{1-e^2}$ の時間変化率 $\displaystyle \frac{dS}{dt}$ が一定であるから
$$\frac{dS}{dt} = \mbox{const.} = \frac{S}{T} = \frac{\pi a^2 \sqrt{1-e^2}}{T} \tag{1}$$
一方で，面積速度 $\displaystyle \frac{dS}{dt}$ は，力学的な視点からは単位質量あたりの角運動量の面に垂直な成分の半分であるから
\begin{eqnarray}
\frac{dS}{dt} &=&\frac{1}{2} \left( \boldsymbol{r}\times \dot{\boldsymbol{r}}\right)_z \\
&=& \frac{1}{2}\left( x \dot{y} – \dot{x} y\right)\\
&=& \frac{1}{2} a^2 \sqrt{1 – e^2}\left((\cos u – e)\cdot \cos u\, \dot{u}
– (- \sin u \, \dot{u}) \cdot \sin u\right)\\
&=& \frac{1}{2}a^2 \sqrt{1 – e^2} (1 – e \cos u) \dot{u} \tag{2}
\end{eqnarray}
つまり，面積速度一定の法則とは，角運動量保存則のことであった。
ケプラー方程式の導出
ケプラー方程式とは，面積速度一定の法則の積分のことである。
$(1)$ および $(2)$ 式より
\begin{eqnarray}
\frac{1}{2}a^2 \sqrt{1 – e^2} (1 – e \cos u) \dot{u} &=& \frac{\pi a^2 \sqrt{1-e^2}}{T}\\
(1 – e \cos u) \frac{du}{dt} &=& \frac{2 \pi}{T}\\
\int (1 – e \cos u)du &=& \frac{2 \pi}{T} \int dt\\
\therefore\ \ u – e \sin u &=& \frac{2 \pi}{T} t
\end{eqnarray}
上記で積分定数は $t = 0$ で $u = 0$ という初期条件でゼロとしている。
これがケプラー方程式。
ケプラー運動の時間平均
ちなみに，
$$ r = a (1 – e \cos u)$$
であったから，
$$(1 – e \cos u) \frac{du}{dt} = \frac{2 \pi}{T}$$
を使うと以下の関係が導かれる。
$${\color{blue}\frac{1}{T} dt} = {\color{red}\frac{1}{2\pi} \frac{r}{a} du}$$
この式はケプラー運動の時間的な平均値を計算するときに使う。たとえば
\begin{eqnarray}
\left\langle\frac{1}{r} \right\rangle &\equiv& {\color{blue}\frac{1}{T}} \int_0^T \frac{1}{r} {\color{blue}dt} \\
&=& {\color{red}\frac{1}{2\pi}} \int_0^{2\pi} \frac{1}{r} {\color{red}\frac{r}{a} du} \\
&=& \frac{1}{a}
\end{eqnarray}
Maxima-Jupyter でケプラー方程式を数値的に解く
ケプラー方程式を Maxima で数値的に解く
楕円の焦点を原点としたデカルト座標 $x, y$ は媒介変数 $u$ （離心近点離角）によって，時間 $t$ と関係づけられているが，$t$ の陽関数としては表すことができない。
\begin{eqnarray}
x &=& r \cos\varphi = a(\cos u – e)\\
y &=& r \sin\varphi = a \sqrt{1-e^2} \sin u\\
\frac{2\pi}{T} t &=& u – e \sin u
\end{eqnarray}
上記の3行目の式が，離心近点離角 $u$ と時間 $t$ を関係付ける ケプラー方程式 である。
そこで，周期 $T$ を $N$ 等分し，
$$t_i = \frac{T}{N} \times i, \quad (i = 0, 1, \dots, N)$$
に対して，$u_i = u(t_i)$ を数値的に求め，等しい時間間隔 $\displaystyle \Delta t = \frac{T}{N}$ ごとの位置 $x(u_i), y(u_i)$ を求めてみる。
/* 軌道長半径 */
a: 5$
/* 離心率 */
e: 6/10$
/* 分割数 */
N: 36$
$$g(u_i) \equiv u_i – e \sin u_i = \frac{2\pi}{T} t_i = \frac{2\pi}{N} \times i$$を $u_i$ について find_root() 関数で数値的に解きます。
g(u):= u - e*sin(u);
for i:0 thru N do
u[i]: find_root(g(u) = 2*%pi/N * i, u, 0, 2*%pi)$
xu(u):= a*(cos(u) - e);
yu(u):= a*sqrt(1-e**2)*sin(u);
xy: makelist([xu(u[i]), yu(u[i])], i, 0, N)$
plot2d([discrete, xy], [style, [points, 1.5]],
[same_xy], [x, -9, 3], [y, -6, 6])$
楕円
\begin{eqnarray}
r &=& \frac{a (1-e^2)}{1 + e \cos\varphi}\\
x &=& r \cos\varphi\\
y &=& r \sin\varphi
\end{eqnarray}
も重ねて描いてみます。
r(phi,a,e):= a*(1-e**2)/(1 + e*cos(phi));
x(phi):= r(phi,a,e)*cos(phi);
y(phi):= r(phi,a,e)*sin(phi);
plot2d([[discrete, xy],
[parametric, x(phi), y(phi), [phi, 0, 2*%pi]]],
[xlabel, "x"], [ylabel, "y"], [legend, "", ""],
[same_xy], [x, -9, 3], [y, -6, 6],
[style, [points, 1.5], lines])$
上の図から，一定時間間隔ごとの位置のプロットなのに，原点（焦点）に近いときは間隔が広い，つまりすばやく動き，原点から離れているときは感覚が狭い，つまりゆっくり動いていることがわかる。