2023年11月24日金曜日

ケプラー方程式

楕円軌道からの続き

軌道の形ではなく,時間発展を考える。
出発点は,$\bm{r}(t) = (x(t),y(t)) = ( a(\cos\theta(t) -e), a\sqrt{1-e^2} \sin \theta(t))$と,面積速度が一定であるということだ。
長半径$a$,短半径$b$,離心率$e$の 楕円の面積$S$は $S=\pi a b = \pi a^2 \sqrt{1-e^2}$なので,周期を$T$とすると,面積速度は,$\dfrac{dS}{dt}= \dfrac{S}{T} = \dfrac{\pi a^2 \sqrt{1-e^2}}{T}$である。

次に,楕円上の位置ベクトル$\bm{r}(t)$から面積速度を計算する。
$\dfrac{dS}{dt}= \frac{1}{2}(\bm{r} \times \dot{\bm{r}})_z = \frac{1}{2}(x \dot{y} - \dot{x} y) $
$= \frac{1}{2}  \{ a(\cos\theta-e) \cdot a\sqrt{1-e^2} \cos\theta \dot{\theta} - (-a \sin \theta \dot{\theta} ) \cdot a\sqrt{1-e^2}\sin \theta \}$
$= \dfrac{a^2 \dot{\theta}\sqrt{1-e^2}}{2} (\cos^2 \theta -e \cos\theta + \sin^2 \theta) =  \dfrac{a^2 \dot{\theta}\sqrt{1-e^2}}{2} (1 -e \cos\theta ) $

この2つの式が等しいので,$\dfrac{2\pi}{T} = \dot{\theta} ( 1- e \cos\theta)$となる。この両辺を時間$t$で積分して,$t=0$で$\theta=0$とすれば,次のケプラー方程式が得られる。
$\dfrac{2\pi}{T} t = \theta -\sin \theta \quad ( 0 \le t \le T \ \ \rightarrow\ \  0 \le \theta \le 2\pi) $
この解$\theta(t)$ を$\bm{r}(t) = (x(t),y(t)) = ( a(\cos\theta(t) -e), a\sqrt{1-e^2} \sin \theta(t))$に代入すれば,位置ベクトルが時間の関数として表される。

Mathematicaで計算してみた。
f[t_, e_] := FindRoot[u - e Sin[u] == 2 Pi t, {u, 0}]
g1[a_, e_] := 
  Table[{a (Cos[u] - e), a Sqrt[1 - e^2] Sin[u]} /. f[k/52., e], {k, 1, 52}];
gp1 = Graphics[{PointSize -> Large, Red, Point[g1[1, 0.2]]}];
g2[a_, e_] := 
 Plot[{a Sqrt[1 - e^2] Sqrt[1 - (x/a + e)^2], -a Sqrt[1 - e^2] Sqrt[
     1 - (x/a + e)^2]}, {x, -a (1 + e), a (1 - e)}, 
  AspectRatio -> Automatic, PlotStyle -> Blue]
gp2 = g2[1, 0.2];
Show[gp2, gp1]

図:ケプラー軌道の計算例(a=1, b=0.98,  e=0.2, r_ap=1.2, r_pe=0.8)

0 件のコメント: