夜空を見上げて、惑星がどのように動いているのか不思議に思ったことはないでしょうか。古代から人類は惑星の軌道を理解しようとしてきましたが、その運動を正確に記述する法則を見出したのは、17世紀のヨハネス・ケプラーでした。そしてその法則がなぜ成り立つのかを万有引力の理論から証明したのが、アイザック・ニュートンです。
ケプラーの法則を理解することは、単なる歴史の学習にとどまりません。人工衛星の軌道設計では、衛星をどの高度に投入すればどのくらいの周期で地球を周回するかを計算する基盤がケプラーの第3法則です。また、惑星探査ミッションでは、地球から火星へ最小エネルギーで到達するホーマン遷移軌道の設計に楕円軌道の理論が不可欠です。さらに近年注目される系外惑星の検出でも、恒星のふらつき（ドップラー法）から惑星の軌道パラメータを推定する際に、ケプラーの法則が基礎理論として活躍しています。
本記事では、万有引力の法則とニュートンの運動法則を出発点として、ケプラーの3法則を1つずつ数学的に導出します。導出の各ステップでは物理的な意味を丁寧に解説し、Pythonによるシミュレーションと可視化で理論を実感できるようにしました。
本記事の内容
- ケプラーの法則の歴史的背景とティコ・ブラーエの観測
- 万有引力の法則と運動方程式の準備
- 第2法則（面積速度一定）の証明 ― 角運動量保存から
- 第1法則（楕円軌道）の導出 ― 極座標の運動方程式とビネの方程式
- 第3法則（周期と半長軸の関係）の証明 ― 円軌道からの導出と楕円への一般化
- Pythonによる楕円軌道の描画、面積速度一定の可視化、第3法則の検証
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ケプラーの法則の歴史的背景
ティコ・ブラーエの精密観測
ケプラーの法則は、いきなり天才のひらめきから生まれたわけではありません。その背景には、デンマークの天文学者ティコ・ブラーエ（1546–1601）による20年以上にわたる惑星の精密観測がありました。
ティコは望遠鏡が発明される以前の時代に、肉眼観測用の巨大な分度器のような器具を自ら設計し、惑星の位置を角度にして約1分角（1/60度）の精度で記録し続けました。当時の他の天文学者の観測精度が10分角程度だったことを考えると、ティコの精度は一桁違う驚異的なものでした。特に火星の軌道データが膨大に蓄積されており、このデータが後にケプラーの発見を可能にします。
ケプラーの苦闘と3法則の発見
1600年にティコの助手となったケプラー（1571–1630）は、師の死後にその観測データを引き継ぎます。ケプラーは当初、古代ギリシャ以来の伝統にしたがい、火星の軌道を「正円」の組み合わせで説明しようと試みました。しかし、何年もの計算の末、どうしても8分角の誤差が残ってしまいます。
ここでケプラーは決定的な判断を下しました。「ティコの観測精度は8分角の誤差を許さない。軌道が正円であるという仮定が間違っているのだ」と結論づけたのです。円を捨て、楕円を試したところ、火星の観測データと完全に一致しました。こうして1609年に第1法則（楕円軌道）と第2法則（面積速度一定）が、1619年に第3法則（周期と半長軸の関係）が発表されました。
ケプラーの3法則を簡潔にまとめると次のとおりです。
- 第1法則（楕円軌道の法則）: 惑星は太陽を1つの焦点とする楕円軌道を描く
- 第2法則（面積速度一定の法則）: 太陽と惑星を結ぶ線分が単位時間に掃く面積は一定である
- 第3法則（調和の法則）: 惑星の公転周期の2乗は、軌道の半長軸の3乗に比例する
これらはあくまで観測データから帰納的に見出された「経験法則」です。ケプラー自身はなぜこれらの法則が成り立つのかを説明できませんでした。その理由を明らかにしたのが、ケプラーの死後に生まれたニュートン（1643–1727）です。ニュートンは万有引力の法則と運動の法則を定式化し、そこからケプラーの3法則すべてを数学的に演繹できることを示しました。
これから、ニュートンの理論的枠組みを使って、ケプラーの3法則を1つずつ証明していきましょう。まずは出発点となる万有引力の法則を確認します。
万有引力の法則と運動方程式
万有引力の直感的理解
2つの物体は、質量を持つだけで互いに引き合います。日常生活ではこの力を感じることはありませんが、地球とリンゴのように一方の質量が巨大になると、私たちが「重力」と呼ぶ力として実感できるようになります。ニュートンの偉大な洞察は、リンゴを地面に落とす力と、月を軌道に留めている力が同じ種類の力であると見抜いたことでした。
万有引力の大きさは、2つの物体の質量の積に比例し、距離の2乗に反比例します。これをベクトルの形で書くと次のようになります。
$$ \begin{equation} \bm{F} = -\frac{GMm}{r^2}\hat{\bm{r}} \end{equation} $$
ここで $G$ は万有引力定数（$6.674 \times 10^{-11} \, \mathrm{N \cdot m^2/kg^2}$）、$M$ は中心天体（太陽や地球）の質量、$m$ は軌道上の物体（惑星や衛星）の質量、$r$ は2物体間の距離、$\hat{\bm{r}}$ は中心天体から物体へ向かう単位ベクトルです。マイナス符号は、力の向きが $\hat{\bm{r}}$ と逆、すなわち引力であることを表しています。
運動方程式への変換
ニュートンの運動法則 $\bm{F} = m\bm{a}$ を適用すると、質量 $m$ の物体の加速度が得られます。
$$ m\ddot{\bm{r}} = -\frac{GMm}{r^2}\hat{\bm{r}} $$
両辺を $m$ で割ると、軌道上の物体の質量が消えます。
$$ \ddot{\bm{r}} = -\frac{GM}{r^2}\hat{\bm{r}} $$
この結果は非常に重要です。惑星の軌道は惑星自身の質量によらず、中心天体の質量 $M$ だけで決まるということを意味しています。実際、地球の周りを回る人工衛星は、質量が100kgでも10,000kgでも同じ軌道を描きます（厳密には二体問題として扱うべきですが、$m \ll M$ の場合には非常に良い近似です）。
ここで 重力パラメータ $\mu = GM$ を導入すると、運動方程式はさらに簡潔になります。
$$ \begin{equation} \ddot{\bm{r}} = -\frac{\mu}{r^2}\hat{\bm{r}} = -\frac{\mu}{r^3}\bm{r} \end{equation} $$
最後の等式では $\hat{\bm{r}} = \bm{r}/r$ を代入しました。地球の場合、$\mu_{\oplus} = 3.986 \times 10^{14} \, \mathrm{m^3/s^2}$、太陽の場合、$\mu_{\odot} = 1.327 \times 10^{20} \, \mathrm{m^3/s^2}$ です。
この運動方程式が、これからケプラーの3法則すべてを導出する出発点です。万有引力の重要な特徴は、力が常に中心天体の方向を向いている中心力であるということです。この性質から、まず第2法則（面積速度一定）を証明できます。
第2法則: 面積速度一定の法則
なぜ第2法則から始めるのか
ケプラーの法則は第1、第2、第3の順番で番号が振られていますが、証明の論理的順序としては第2法則を先に示すのが自然です。第2法則（面積速度一定）は角運動量保存の直接的な帰結であり、その結果が第1法則（楕円軌道）の導出で使われるからです。
面積速度の直感的イメージ
面積速度とは、太陽と惑星を結ぶ「動径ベクトル」が単位時間に掃く面積の速さです。日常のアナロジーで考えると、回転ドアの扉を想像してください。扉（動径）の先端が素早く動くと扉が掃く面積は大きくなり、ゆっくり動くと小さくなります。ケプラーの第2法則は、この掃く面積の速さが軌道上のどの位置でも一定であると主張しています。
つまり、惑星が太陽に近いとき（近日点付近）は速く動き、遠いとき（遠日点付近）は遅く動くのですが、その速さの変化がちょうど距離の変化を補って、面積速度が一定に保たれるのです。この性質は、万有引力が中心力であること、すなわち角運動量が保存されることから自然に導かれます。
角運動量保存の証明
中心力のもとで運動する物体の角運動量が保存されることを示します。角運動量は $\bm{L} = m\bm{r} \times \bm{v}$ で定義されますが、ここでは単位質量あたりの比角運動量 $\bm{h} = \bm{r} \times \bm{v}$ を使います。
$\bm{h}$ の時間微分を計算しましょう。
$$ \dot{\bm{h}} = \frac{d}{dt}(\bm{r} \times \bm{v}) = \dot{\bm{r}} \times \bm{v} + \bm{r} \times \dot{\bm{v}} $$
右辺の第1項を見ると、$\dot{\bm{r}} = \bm{v}$ なので $\bm{v} \times \bm{v} = \bm{0}$ です（任意のベクトルと自分自身の外積はゼロ）。
第2項では $\dot{\bm{v}} = \bm{a}$ が加速度です。万有引力による加速度は $\bm{a} = -(\mu/r^3)\bm{r}$ であり、$\bm{r}$ と平行です。平行なベクトル同士の外積はゼロなので、$\bm{r} \times \bm{a} = \bm{0}$ となります。
したがって、
$$ \dot{\bm{h}} = \bm{0} + \bm{0} = \bm{0} $$
$$ \begin{equation} \bm{h} = \bm{r} \times \bm{v} = \text{const} \end{equation} $$
比角運動量ベクトル $\bm{h}$ が時間によらず一定であることが示されました。この結果には2つの重要な意味があります。
- $\bm{h}$ の大きさが一定 → 面積速度が一定（第2法則）
- $\bm{h}$ の向きが一定 → 軌道は $\bm{h}$ に垂直な平面内に限られる（軌道面が一定）
2番目の帰結は、惑星の運動が3次元空間ではなく2次元平面内で起こることを意味しています。これにより、以降の解析を極座標 $(r, \theta)$ で行うことが正当化されます。
面積速度の計算
極座標において、位置ベクトルの成分は $\bm{r} = r\hat{\bm{r}}$、速度ベクトルは $\bm{v} = \dot{r}\hat{\bm{r}} + r\dot{\theta}\hat{\bm{\theta}}$ です。外積を計算すると、
$$ \bm{h} = \bm{r} \times \bm{v} = r\hat{\bm{r}} \times (\dot{r}\hat{\bm{r}} + r\dot{\theta}\hat{\bm{\theta}}) = r^2\dot{\theta}\,(\hat{\bm{r}} \times \hat{\bm{\theta}}) $$
$\hat{\bm{r}} \times \hat{\bm{\theta}}$ は軌道面に垂直な単位ベクトルなので、比角運動量の大きさは、
$$ h = |\bm{h}| = r^2\dot{\theta} $$
微小時間 $dt$ の間に動径が掃く微小面積は、底辺 $r$、高さ $r\,d\theta$ の三角形として近似でき、$dA = \frac{1}{2}r^2\,d\theta$ です。したがって面積速度は、
$$ \begin{equation} \frac{dA}{dt} = \frac{1}{2}r^2\dot{\theta} = \frac{h}{2} = \text{const} \end{equation} $$
これがケプラーの第2法則の数学的表現です。面積速度 $dA/dt$ は軌道上のどの位置でも $h/2$ に等しく、一定です。 $\square$
ここで得られた $h = r^2\dot{\theta} = \text{const}$ という関係式は、第1法則の導出で繰り返し使われます。角運動量保存という1つの原理が、軌道の形を決定する鍵となるのです。次に、この関係式を使って惑星の軌道が楕円であることを証明しましょう。
第1法則: 楕円軌道の導出
楕円軌道の直感的理解
なぜ惑星の軌道が楕円になるのでしょうか。まず「もし重力が存在しなかったら」を考えると、惑星は初速度の方向に直線運動をします。一方「もし惑星が静止していたら」、重力によってまっすぐ太陽に向かって落下します。実際にはこの2つの効果が同時に働いています。惑星は「落ちながら進む」ことで、直線でも直接落下でもない、曲がった軌道を描くのです。
初速度と方向がちょうど適切であれば正円軌道になりますが、一般的には楕円になります。速度が大きすぎれば放物線や双曲線になって太陽から逃げてしまいます。これらはすべて円錐曲線と呼ばれる曲線の仲間であり、逆2乗則の力（万有引力）のもとでの軌道が必ず円錐曲線になることを、これから数学的に示します。
極座標での運動方程式
軌道面が一定であることが示されたので、その平面内で極座標 $(r, \theta)$ を用います。極座標における加速度の $r$ 成分と $\theta$ 成分はそれぞれ次のように表されます。
$r$ 方向（動径方向）の運動方程式:
$$ \ddot{r} – r\dot{\theta}^2 = -\frac{\mu}{r^2} $$
左辺の $-r\dot{\theta}^2$ は遠心加速度の項です。右辺は万有引力による加速度です。
$\theta$ 方向（接線方向）の運動方程式:
$$ r\ddot{\theta} + 2\dot{r}\dot{\theta} = 0 $$
この式は $\frac{1}{r}\frac{d}{dt}(r^2\dot{\theta}) = 0$ と変形でき、$r^2\dot{\theta} = h = \text{const}$ を意味します。これは先ほど証明した角運動量保存そのものです。
$r$ 方向の方程式に含まれる $\dot{\theta}$ を消去するため、角運動量保存の式 $\dot{\theta} = h/r^2$ を利用します。これにより、$r$ 方向の方程式は $r$ と時間 $t$ だけの微分方程式になります。しかし、このままでは非線形で解きにくいため、巧妙な変数変換を導入します。
ビネの変数変換
ここで $u = 1/r$ とおくビネの変数変換を導入します。この変換の意図は、独立変数を時間 $t$ から角度 $\theta$ に切り替え、微分方程式を線形化することです。
まず $r = 1/u$ から $\dot{r}$ を $\theta$ の関数として表します。$h = r^2\dot{\theta}$ より $\dot{\theta} = hu^2$ なので、
$$ \dot{r} = \frac{dr}{dt} = \frac{dr}{d\theta}\cdot\dot{\theta} $$
$r = 1/u$ を $\theta$ で微分すると $\frac{dr}{d\theta} = -\frac{1}{u^2}\frac{du}{d\theta}$ です。したがって、
$$ \dot{r} = -\frac{1}{u^2}\frac{du}{d\theta}\cdot hu^2 = -h\frac{du}{d\theta} $$
次に $\ddot{r}$ を求めます。$\dot{r} = -h\frac{du}{d\theta}$ をさらに時間で微分します。
$$ \ddot{r} = \frac{d\dot{r}}{dt} = \frac{d\dot{r}}{d\theta}\cdot\dot{\theta} = -h\frac{d^2u}{d\theta^2}\cdot hu^2 = -h^2u^2\frac{d^2u}{d\theta^2} $$
ビネの方程式の導出
得られた $\ddot{r}$ と $\dot{\theta} = hu^2$ を $r$ 方向の運動方程式に代入します。
$r$ 方向の方程式 $\ddot{r} – r\dot{\theta}^2 = -\mu/r^2$ において、$r = 1/u$ を使うと右辺は $-\mu u^2$ です。左辺に代入すると、
$$ -h^2u^2\frac{d^2u}{d\theta^2} – \frac{1}{u}(hu^2)^2 = -\mu u^2 $$
左辺第2項を整理すると $-h^2u^3$ となります。
$$ -h^2u^2\frac{d^2u}{d\theta^2} – h^2u^3 = -\mu u^2 $$
両辺を $-h^2u^2$（これは $r > 0$ より常にゼロでない）で割ると、
$$ \begin{equation} \frac{d^2u}{d\theta^2} + u = \frac{\mu}{h^2} \end{equation} $$
これがビネの方程式です。右辺 $\mu/h^2$ は定数なので、これは $u$ についての非斉次の単振動方程式（調和振動子方程式）です。驚くべきことに、元の非線形微分方程式が、変数変換によって線形の定係数微分方程式に変換されました。
ビネの方程式の解
単振動方程式の一般解は容易に書き下せます。斉次解（$\cos$ または $\sin$）に特殊解（定数 $\mu/h^2$）を加えたものです。
$$ u(\theta) = \frac{\mu}{h^2} + C\cos(\theta – \theta_0) $$
ここで $C$ は初期条件で決まる定数、$\theta_0$ は近点方向を定める角度です。座標を適切に選んで $\theta_0 = 0$ とすると、
$$ u = \frac{\mu}{h^2}(1 + e\cos\theta) $$
ここで $e = Ch^2/\mu$ は離心率と呼ばれる無次元のパラメータです。$r = 1/u$ に戻すと、
$$ \begin{equation} r = \frac{h^2/\mu}{1 + e\cos\theta} = \frac{p}{1 + e\cos\nu} \end{equation} $$
ここで $p = h^2/\mu$ は半通径（セミラタスレクタム）、$\nu = \theta – \theta_0$ は真近点角（近点からの角度）です。
円錐曲線の分類
式(6)は円錐曲線の極座標表現です。離心率 $e$ の値によって軌道の形が決まります。
| 離心率 $e$ の範囲 | 軌道の形 | 力学的エネルギー |
|---|---|---|
| $e = 0$ | 円 | $E < 0$ |
| $0 < e < 1$ | 楕円 | $E < 0$ |
| $e = 1$ | 放物線 | $E = 0$ |
| $e > 1$ | 双曲線 | $E > 0$ |
惑星は太陽の重力に束縛されており、力学的エネルギーが負なので $e < 1$ です。したがって惑星の軌道は楕円であり、これがケプラーの第1法則です。 $\square$
この導出の本質は、「逆2乗則の中心力」という性質がビネの方程式を通じて「単振動方程式」に帰着し、その解が円錐曲線を与えるという点にあります。もし力の法則が逆2乗則でなければ、ビネの方程式の形が変わり、軌道は一般に閉じた曲線にならないことが知られています（ベルトランの定理）。
楕円の幾何学的パラメータ
導出された軌道方程式 $r = p/(1 + e\cos\nu)$ から、楕円の幾何学的パラメータを読み取ることができます。
近点（太陽に最も近い点）では $\nu = 0$ なので $r_{\min} = p/(1+e)$、遠点（太陽から最も遠い点）では $\nu = \pi$ なので $r_{\max} = p/(1-e)$ です。
楕円の半長軸 $a$ は近点距離と遠点距離の平均です。
$$ a = \frac{r_{\min} + r_{\max}}{2} = \frac{p}{2}\left(\frac{1}{1+e} + \frac{1}{1-e}\right) = \frac{p}{1 – e^2} $$
最初の等式では近点と遠点の距離を代入し、通分して整理すると $a = p/(1-e^2)$ が得られます。
逆に解くと $p = a(1-e^2)$ であり、これを使って各パラメータを半長軸 $a$ と離心率 $e$ で表せます。
$$ \begin{equation} a = \frac{p}{1 – e^2}, \quad b = a\sqrt{1 – e^2}, \quad r_{\min} = a(1-e), \quad r_{\max} = a(1+e) \end{equation} $$
ここで $b$ は半短軸で、$b^2 = a^2(1-e^2) = ap$ の関係があります。半短軸は楕円の「幅」に対応し、離心率が大きいほど楕円は細長くなります。$e = 0$ のとき $a = b$ となり、楕円は正円に退化します。
これらの幾何学的パラメータは軌道要素の定義として体系的に整理されています。
ここまでで第1法則（楕円軌道）と第2法則（面積速度一定）を証明しました。残る第3法則は、面積速度と楕円の面積を組み合わせることで導くことができます。
第3法則: 周期と半長軸の関係
円軌道での簡単な導出
第3法則の本質をつかむために、まず最も簡単なケースである円軌道で考えてみましょう。
半径 $r$ の円軌道を一定の速さ $v$ で周回する場合、遠心加速度 $v^2/r$ と重力加速度 $\mu/r^2$ がつり合います。
$$ \frac{v^2}{r} = \frac{\mu}{r^2} $$
これを $v$ について解くと、円軌道速度は $v = \sqrt{\mu/r}$ です。この結果は直感に合います。太陽から遠いほど重力が弱いため、つり合いに必要な軌道速度が小さくなるのです。この関係はビザビバ方程式で $a = r$ とした特殊ケースに対応します。
周期 $T$ は円周 $2\pi r$ を速さ $v$ で割ったものです。
$$ T = \frac{2\pi r}{v} = \frac{2\pi r}{\sqrt{\mu/r}} = 2\pi\sqrt{\frac{r^3}{\mu}} $$
2乗すると、
$$ T^2 = \frac{4\pi^2}{\mu}r^3 $$
円軌道では $r = a$（半長軸＝半径）なので、$T^2 \propto a^3$ が得られました。これがケプラーの第3法則の円軌道版です。
ここで注目すべきは、比例定数 $4\pi^2/\mu$ が軌道上の物体の質量に依存せず、中心天体の質量 $M$（$\mu = GM$）だけで決まるという点です。したがって、同じ中心天体を周回するすべての物体について $T^2/a^3$ は同一の値を持ちます。
楕円軌道への一般化
円軌道での導出は直感的でしたが、楕円軌道でも同じ関係が成り立つことを示す必要があります。ここで第2法則の結果を活用します。
楕円の面積は $A = \pi ab$ です。面積速度が $dA/dt = h/2$ で一定なので、1周期 $T$ の間に掃く面積が楕円の全面積に等しいことから、
$$ A = \int_0^T \frac{dA}{dt}\,dt = \frac{h}{2}\cdot T $$
これを $T$ について解くと、
$$ T = \frac{2A}{h} = \frac{2\pi ab}{h} $$
ここから $T^2$ を計算するために、$b$ と $h$ を $a$ と $e$ で表します。
半短軸は $b = a\sqrt{1-e^2}$ です。半通径の定義から $p = h^2/\mu = a(1-e^2)$ なので、$h^2 = \mu a(1-e^2)$ です。
$T = 2\pi ab/h$ を2乗すると、
$$ T^2 = \frac{4\pi^2 a^2 b^2}{h^2} $$
分子に $b^2 = a^2(1-e^2)$ を代入すると $4\pi^2 a^2 \cdot a^2(1-e^2) = 4\pi^2 a^4(1-e^2)$ です。
分母に $h^2 = \mu a(1-e^2)$ を代入すると、
$$ T^2 = \frac{4\pi^2 a^4(1-e^2)}{\mu a(1-e^2)} $$
分子と分母の $(1-e^2)$ が約分され、$a^4/a = a^3$ が残ります。
$$ \begin{equation} T^2 = \frac{4\pi^2}{\mu}a^3 \end{equation} $$
離心率 $e$ がきれいに消えて、周期は半長軸 $a$ のみで決まることが示されました。これがケプラーの第3法則の一般的な証明です。 $\square$
この結果の意味を具体的な数値で確認してみましょう。地球の重力パラメータ $\mu = 3.986 \times 10^{14} \, \mathrm{m^3/s^2}$ を使うと、静止軌道（$T = 24$ 時間）の半長軸が計算できます。
$$ a = \left(\frac{\mu T^2}{4\pi^2}\right)^{1/3} = \left(\frac{3.986 \times 10^{14} \times (86400)^2}{4\pi^2}\right)^{1/3} \approx 42{,}164 \, \mathrm{km} $$
地球の半径が約6,371 kmなので、静止軌道は地表から約35,786 kmの高度に位置することがわかります。これは実際の静止衛星の軌道高度と一致しており、第3法則が実用的な計算に直結していることを実感できます。
なお、実際の衛星軌道では、地球が完全な球ではないことによるJ2摂動などの効果が加わり、ケプラー軌道からのずれが生じます。しかし、ケプラーの法則は最も基本的な近似として、軌道設計の出発点であり続けています。
ここまでで3つの法則すべてを数学的に証明しました。次に、Pythonを使ってこれらの法則を視覚的に確認していきましょう。
Pythonでの実装: 楕円軌道の描画
離心率を変えた楕円軌道
まず、異なる離心率を持つ楕円軌道を描画し、離心率が軌道の形状にどのように影響するかを確認します。
import numpy as np
import matplotlib.pyplot as plt
# 重力パラメータ（地球）
mu = 3.986e14 # [m^3/s^2]
# 半長軸を固定し、離心率を変えた軌道を描画
a = 20000e3 # 半長軸 20,000 km [m]
eccentricities = [0.0, 0.3, 0.6, 0.8]
fig, ax = plt.subplots(figsize=(9, 9))
for e in eccentricities:
p = a * (1 - e**2) # 半通径
nu = np.linspace(0, 2 * np.pi, 1000)
r = p / (1 + e * np.cos(nu))
x = r * np.cos(nu)
y = r * np.sin(nu)
ax.plot(x / 1e6, y / 1e6, linewidth=2, label=f'e = {e}')
# 焦点（原点）を表示
ax.plot(0, 0, 'o', color='orange', markersize=12, label='焦点（中心天体）')
ax.set_xlabel('x [×1000 km]', fontsize=13)
ax.set_ylabel('y [×1000 km]', fontsize=13)
ax.set_title('離心率による楕円軌道の変化（半長軸 a = 20,000 km 固定）', fontsize=14)
ax.set_aspect('equal')
ax.legend(fontsize=11, loc='upper right')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このグラフから、いくつかの重要な点が読み取れます。
- $e = 0$ のとき、軌道は完全な正円になり、焦点は円の中心と一致します。これは楕円の特殊ケースです。
- 離心率が大きくなるほど軌道は細長くなり、焦点（中心天体）は楕円の中心から大きくずれます。$e = 0.8$ の軌道では、近点と遠点の距離が $a(1-e) = 4{,}000$ km と $a(1+e) = 36{,}000$ km で、9倍もの差があります。
- すべての軌道が同じ焦点を共有しています。これはケプラーの第1法則の「太陽（中心天体）が楕円の焦点に位置する」という主張を視覚的に確認しています。
楕円軌道の詳細表示
次に、1つの楕円軌道について、近点・遠点・半長軸・半短軸などの幾何学的パラメータを可視化します。
import numpy as np
import matplotlib.pyplot as plt
mu = 3.986e14 # [m^3/s^2]
a = 20000e3 # 半長軸 [m]
e = 0.6 # 離心率
p = a * (1 - e**2) # 半通径
b = a * np.sqrt(1 - e**2) # 半短軸
c = a * e # 焦点距離
# 軌道
nu = np.linspace(0, 2 * np.pi, 1000)
r = p / (1 + e * np.cos(nu))
x = r * np.cos(nu)
y = r * np.sin(nu)
# 楕円の中心座標（焦点から -c の位置）
cx = -c
fig, ax = plt.subplots(figsize=(12, 7))
# 軌道描画
ax.plot(x / 1e6, y / 1e6, 'b-', linewidth=2, label='楕円軌道')
# 焦点
ax.plot(0, 0, 'o', color='orange', markersize=14, label='焦点（地球）', zorder=5)
# 近点・遠点
r_min = a * (1 - e)
r_max = a * (1 + e)
ax.plot(r_min / 1e6, 0, 'r^', markersize=12, label=f'近点 r_min = {r_min/1e3:.0f} km')
ax.plot(-r_max / 1e6, 0, 'gv', markersize=12, label=f'遠点 r_max = {r_max/1e3:.0f} km')
# 半長軸（楕円の中心を通る水平線）
ax.plot([cx / 1e6 - a / 1e6, cx / 1e6 + a / 1e6], [0, 0], 'k--', alpha=0.5)
ax.annotate('2a', xy=(cx / 1e6, -1.5), fontsize=12, ha='center', color='black')
# 半短軸
ax.plot([cx / 1e6, cx / 1e6], [-b / 1e6, b / 1e6], 'k--', alpha=0.5)
ax.annotate('2b', xy=(cx / 1e6 + 1.0, 0), fontsize=12, ha='left', color='black')
ax.set_xlabel('x [×1000 km]', fontsize=13)
ax.set_ylabel('y [×1000 km]', fontsize=13)
ax.set_title(f'楕円軌道の幾何学 (a={a/1e3:.0f} km, e={e}, b={b/1e3:.0f} km)', fontsize=14)
ax.set_aspect('equal')
ax.legend(fontsize=10, loc='upper right')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 軌道パラメータの数値確認
T = 2 * np.pi * np.sqrt(a**3 / mu)
print(f"半長軸 a = {a/1e3:.0f} km")
print(f"半短軸 b = {b/1e3:.0f} km")
print(f"離心率 e = {e}")
print(f"近点距離 r_min = {r_min/1e3:.0f} km")
print(f"遠点距離 r_max = {r_max/1e3:.0f} km")
print(f"周期 T = {T/3600:.2f} 時間")
この図は楕円軌道の全体像を示しています。焦点（地球）が楕円の中心からずれていることが一目でわかります。$e = 0.6$ の場合、近点距離 8,000 km と遠点距離 32,000 km の差が大きく、焦点の偏心が顕著です。半短軸 $b$ は半長軸 $a$ よりも短く、$b = a\sqrt{1-e^2} \approx 16{,}000$ km です。
出力される周期は約 8.93 時間であり、第3法則 $T = 2\pi\sqrt{a^3/\mu}$ の計算結果と一致することが確認できます。
次に、第2法則（面積速度一定）を可視化して確認しましょう。
Pythonでの実装: 面積速度一定の可視化
ケプラー方程式による時間発展
面積速度が一定であることを可視化するには、「等しい時間間隔で惑星がどの位置にいるか」を計算する必要があります。真近点角 $\nu$ は時間に対して一様には変化しない（近点付近で速く、遠点付近で遅い）ため、ケプラー方程式を用いて時間と位置の対応を求めます。
ケプラー方程式は、平均近点角 $M$（時間に比例して一様に増加する角度）と離心近点角 $E$ の関係を与えます。
$$ M = E – e\sin E $$
この方程式は $E$ について解析的に解けないため、数値的に解く必要があります。ここではニュートン法を使います。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
def solve_kepler(M, e, tol=1e-12):
"""ケプラー方程式 M = E - e*sin(E) をニュートン法で解く"""
E = M.copy() # 初期値
for _ in range(100):
dE = (E - e * np.sin(E) - M) / (1 - e * np.cos(E))
E -= dE
if np.max(np.abs(dE)) < tol:
break
return E
def eccentric_to_true(E, e):
"""離心近点角 E → 真近点角 nu"""
return 2 * np.arctan2(np.sqrt(1 + e) * np.sin(E / 2),
np.sqrt(1 - e) * np.cos(E / 2))
# 軌道パラメータ
mu = 3.986e14
a = 20000e3
e = 0.6
p = a * (1 - e**2)
T = 2 * np.pi * np.sqrt(a**3 / mu)
# 等時間間隔で12分割
N_sectors = 12
M_boundaries = np.linspace(0, 2 * np.pi, N_sectors + 1)
E_boundaries = solve_kepler(M_boundaries, e)
nu_boundaries = eccentric_to_true(E_boundaries, e)
# 軌道全体の描画データ
nu_fine = np.linspace(0, 2 * np.pi, 1000)
r_fine = p / (1 + e * np.cos(nu_fine))
x_fine = r_fine * np.cos(nu_fine)
y_fine = r_fine * np.sin(nu_fine)
# 各セクターの面積を色分けで表示
fig, ax = plt.subplots(figsize=(10, 10))
colors = plt.cm.tab20(np.linspace(0, 1, N_sectors))
areas = []
for i in range(N_sectors):
nu_sec = np.linspace(nu_boundaries[i], nu_boundaries[i + 1], 200)
r_sec = p / (1 + e * np.cos(nu_sec))
x_sec = r_sec * np.cos(nu_sec)
y_sec = r_sec * np.sin(nu_sec)
# 扇形のポリゴン（焦点→軌道→焦点）
verts = [(0, 0)]
verts += list(zip(x_sec / 1e6, y_sec / 1e6))
verts.append((0, 0))
polygon = Polygon(verts, alpha=0.4, facecolor=colors[i], edgecolor='gray')
ax.add_patch(polygon)
# 面積の数値計算（扇形の面積 = 0.5 * ∫ r^2 dnu）
area = 0.5 * np.trapz(r_sec**2, nu_sec)
areas.append(area)
# 軌道の輪郭
ax.plot(x_fine / 1e6, y_fine / 1e6, 'k-', linewidth=2)
ax.plot(0, 0, 'o', color='orange', markersize=14, zorder=5, label='焦点')
ax.set_xlabel('x [×1000 km]', fontsize=13)
ax.set_ylabel('y [×1000 km]', fontsize=13)
ax.set_title(f'面積速度一定の可視化（e = {e}, 等時間間隔 {N_sectors} 分割）', fontsize=14)
ax.set_aspect('equal')
ax.legend(fontsize=12)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 面積の比較
areas = np.array(areas)
print("各セクターの面積 [km^2]:")
for i, area in enumerate(areas):
print(f" セクター {i+1:2d}: {area/1e6:.2f} ×10^6 km^2")
print(f"\n面積の最大/最小比: {areas.max()/areas.min():.6f}")
print(f"面積の標準偏差/平均: {areas.std()/areas.mean()*100:.4f} %")
このコードは、1周期を12等分した時間間隔ごとに、動径が掃く扇形の領域を色分けして表示しています。
グラフを見ると、近点付近の扇形は幅が狭く遠くまで伸びず、遠点付近の扇形は幅が広く遠くまで伸びていることがわかります。しかし面積を数値的に計算すると、すべてのセクターの面積はほぼ等しいことが確認できます。面積の最大/最小比はほぼ1.0であり、数値誤差の範囲で完全に一致しています。これがケプラーの第2法則の視覚的な確認です。
物理的に解釈すると、惑星は近点付近で速く移動するため狭い角度を短時間で掃きますが、距離が近い分だけ面積への寄与が小さくなります。遠点付近では逆に遅く移動しますが、距離が遠いため同じ面積を掃くことができるのです。この「速さと距離のトレードオフ」が角運動量保存の帰結であることを、先ほどの数学的証明で確認しました。
次に、第3法則を太陽系の惑星データを使って検証してみましょう。
Pythonでの実装: 第3法則の検証
太陽系惑星での検証
ケプラーの第3法則 $T^2 = (4\pi^2/\mu)a^3$ が太陽系の惑星データで成り立つかを検証します。理論的には、全惑星について $T^2/a^3$ が同一の値（$4\pi^2/\mu_{\odot}$）になるはずです。
import numpy as np
import matplotlib.pyplot as plt
# 太陽の重力パラメータ
mu_sun = 1.327e20 # [m^3/s^2]
# 太陽系惑星データ（半長軸 [AU], 公転周期 [年]）
planets = {
'水星': (0.387, 0.241),
'金星': (0.723, 0.615),
'地球': (1.000, 1.000),
'火星': (1.524, 1.881),
'木星': (5.203, 11.86),
'土星': (9.537, 29.46),
'天王星': (19.19, 84.01),
'海王星': (30.07, 164.8),
}
names = list(planets.keys())
a_au = np.array([planets[n][0] for n in names])
T_yr = np.array([planets[n][1] for n in names])
# 第3法則の検証: T^2 vs a^3
a_cubed = a_au**3
T_squared = T_yr**2
fig, axes = plt.subplots(1, 2, figsize=(15, 6))
# 左: T^2 vs a^3 の直線プロット
axes[0].scatter(a_cubed, T_squared, s=100, c='royalblue', zorder=5)
for i, name in enumerate(names):
axes[0].annotate(name, (a_cubed[i], T_squared[i]),
textcoords="offset points", xytext=(8, 5), fontsize=10)
# 理論直線（AU, 年の単位系では T^2 = a^3）
a_theory = np.linspace(0, 30000, 100)
axes[0].plot(a_theory, a_theory, 'r--', linewidth=1.5,
label=r'理論: $T^2 = a^3$（AU, 年）')
axes[0].set_xlabel(r'$a^3$ [AU$^3$]', fontsize=13)
axes[0].set_ylabel(r'$T^2$ [年$^2$]', fontsize=13)
axes[0].set_title(r'ケプラーの第3法則: $T^2$ vs $a^3$', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
# 右: 対数プロット（log T vs log a）
axes[1].scatter(np.log10(a_au), np.log10(T_yr), s=100, c='royalblue', zorder=5)
for i, name in enumerate(names):
axes[1].annotate(name, (np.log10(a_au[i]), np.log10(T_yr[i])),
textcoords="offset points", xytext=(8, 5), fontsize=10)
# 理論直線: log T = (3/2) log a
log_a = np.linspace(-0.5, 1.6, 100)
axes[1].plot(log_a, 1.5 * log_a, 'r--', linewidth=1.5,
label=r'理論: 傾き 3/2')
axes[1].set_xlabel(r'$\log_{10}(a$ / AU$)$', fontsize=13)
axes[1].set_ylabel(r'$\log_{10}(T$ / 年$)$', fontsize=13)
axes[1].set_title(r'対数プロット: $\log T = \frac{3}{2}\log a$', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# T^2/a^3 の値を表で出力
print(f"{'惑星':>6s} {'a [AU]':>8s} {'T [年]':>8s} {'T^2/a^3':>10s}")
print("-" * 40)
for i, name in enumerate(names):
ratio = T_yr[i]**2 / a_au[i]**3
print(f"{name:>6s} {a_au[i]:8.3f} {T_yr[i]:8.3f} {ratio:10.4f}")
このコードは、太陽系の8惑星について第3法則を2つの方法で可視化しています。
左のグラフは $T^2$ と $a^3$ の関係を直接プロットしたもので、すべての惑星が理論直線 $T^2 = a^3$（AU・年の単位系）の上に正確に乗っていることがわかります。右のグラフは対数プロットで、$\log T$ と $\log a$ の関係が傾き $3/2$ の直線になっています。$T \propto a^{3/2}$ という関係が全惑星にわたって精密に成り立っていることを視覚的に確認できます。
表に出力される $T^2/a^3$ の値は、すべての惑星でほぼ 1.00（AU・年の単位系）になっています。水星から海王星まで、半長軸が約78倍、周期が約684倍も異なるにもかかわらず、$T^2/a^3$ がすべて同じ値になるという事実は、ケプラーの第3法則の正確さを印象的に示しています。
地球周回衛星での検証
同じ法則が地球を周回する人工衛星にも適用できることを確認します。ここでは、低軌道からの半長軸を変化させたときの周期の変化を連続的にプロットし、いくつかの代表的な衛星軌道を重ねて表示します。
import numpy as np
import matplotlib.pyplot as plt
mu_earth = 3.986e14 # [m^3/s^2]
R_earth = 6371e3 # 地球半径 [m]
# 連続的な T vs a の曲線
a_range = np.linspace(R_earth + 200e3, 50000e3, 500) # LEOから高軌道まで
T_range = 2 * np.pi * np.sqrt(a_range**3 / mu_earth)
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(a_range / 1e3, T_range / 3600, 'b-', linewidth=2,
label=r'$T = 2\pi\sqrt{a^3/\mu}$')
# 代表的な衛星軌道
satellites = {
'ISS (LEO)': 6771e3, # 高度 ~400 km
'GPS': 26560e3, # 半長軸 ~26,560 km
'静止軌道 (GEO)': 42164e3, # 半長軸 ~42,164 km
}
for name, a_sat in satellites.items():
T_sat = 2 * np.pi * np.sqrt(a_sat**3 / mu_earth)
ax.plot(a_sat / 1e3, T_sat / 3600, 'o', markersize=10, zorder=5)
ax.annotate(f'{name}\n(T={T_sat/3600:.1f}h)',
xy=(a_sat / 1e3, T_sat / 3600),
textcoords="offset points", xytext=(10, 5), fontsize=10)
ax.axhline(y=24, color='gray', linestyle=':', alpha=0.5, label='T = 24 時間')
ax.set_xlabel('半長軸 a [km]', fontsize=13)
ax.set_ylabel('周期 T [時間]', fontsize=13)
ax.set_title('地球周回衛星のケプラー第3法則', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このグラフから、半長軸の増加に伴い周期が $a^{3/2}$ に比例して増大する様子が読み取れます。ISS（国際宇宙ステーション）は地球に近いため約1.5時間で1周しますが、GPS衛星は約12時間、静止衛星はちょうど24時間です。静止衛星の周期が地球の自転周期と一致するため、地上から見ると静止しているように見えます。この「ちょうど良い高度」を求めるのに第3法則を逆に使えることが、軌道設計における第3法則の実用的価値です。
各法則の物理的意味の深掘り
ここまでの数学的導出と可視化を踏まえて、3つの法則それぞれが物理的に何を語っているのかを改めて整理しましょう。
第1法則: なぜ楕円になるのか
導出の過程を振り返ると、軌道が楕円になる理由は「力が逆2乗則に従うこと」に帰着します。ビネの方程式において、右辺が定数（$\mu/h^2$）になったのは、力が $1/r^2$ に比例するからでした。もし力が $1/r^3$ や $1/r$ に比例していたら、ビネの方程式は異なる形になり、軌道は一般に閉じた曲線にはなりません。
実際、ベルトランの定理によれば、すべての束縛軌道が閉じるのは、逆2乗則の力と調和振動子的な力（$r$ に比例する力）の2つの場合に限られます。太陽系で惑星が閉じた楕円軌道を描くのは、万有引力がちょうど逆2乗則に従っているおかげなのです。
離心率 $e$ が軌道のエネルギーに対応していることも重要です。エネルギー保存則を適用すると、軌道の力学的エネルギー $E$ と離心率の間に次の関係が成り立ちます。
$$ e = \sqrt{1 + \frac{2Eh^2}{\mu^2}} $$
エネルギーが負（束縛状態）のとき $e < 1$ で楕円、ゼロのとき $e = 1$ で放物線（ぎりぎり脱出）、正のとき $e > 1$ で双曲線（脱出軌道）になります。衛星を軌道に投入する際の速度をわずかに変えるだけで、束縛軌道から脱出軌道に変わりうるのです。
第2法則: 角運動量保存の幾何学的表現
第2法則は角運動量保存をわかりやすく言い換えたものです。角運動量 $h = r^2\dot{\theta}$ が一定であることから、$r$ が大きい（太陽から遠い）ときには $\dot{\theta}$ が小さくなり（角速度が遅くなり）、$r$ が小さいときには $\dot{\theta}$ が大きくなります。
このことは、フィギュアスケーターのスピンに似ています。スケーターが腕を縮めると回転が速くなり、腕を広げると回転が遅くなります。これは角運動量 $L = I\omega$（慣性モーメント × 角速度）が保存されるからです。惑星の場合、太陽に近づくと「腕を縮める」のと同じ効果で角速度が増し、遠ざかると角速度が減ります。
また、第2法則は万有引力だけでなく、あらゆる中心力のもとで成り立つことに注意してください。力が距離の関数であっても（$1/r^2$ でなくても）、方向が常に中心を向いている限り、角運動量は保存され面積速度は一定です。したがって、第2法則は第1法則や第3法則よりも一般的な法則と言えます。
第3法則: 普遍的なスケーリング則
第3法則 $T^2 = (4\pi^2/\mu)a^3$ の最も重要な含意は、比例定数が中心天体の質量だけで決まるということです。これにより、第3法則は天体の質量を測定するための強力なツールになります。
例えば、地球の周りを回る月の公転周期と軌道半径から地球の質量を、木星の衛星の運動から木星の質量を推定できます。歴史的にも、ケプラーの第3法則と万有引力定数 $G$ の精密測定（キャベンディッシュの実験、1798年）を組み合わせて、太陽や地球の質量が初めて正確に求められました。
近年では、系外惑星の検出においてもこの法則が活躍しています。恒星のドップラーシフトから惑星の公転周期 $T$ を測定し、恒星の質量 $M$ を推定すれば、第3法則から惑星の軌道半径 $a$ を求めることができます。ケプラー宇宙望遠鏡が数千個の系外惑星を発見する際にも、この基本法則が使われました。
数値シミュレーション: 軌道の時間発展
運動方程式の数値積分
最後に、万有引力の運動方程式を数値的に解いて軌道を直接シミュレーションし、結果がケプラーの法則と一致することを確認します。初期条件から運動方程式を時間積分すれば、理論に頼らずに軌道が得られます。これにより、導出で使った解析的な方法と独立に法則を検証できます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
mu = 3.986e14 # 地球の重力パラメータ [m^3/s^2]
def equations_of_motion(t, state):
"""二体問題の運動方程式"""
x, y, vx, vy = state
r = np.sqrt(x**2 + y**2)
ax = -mu * x / r**3
ay = -mu * y / r**3
return [vx, vy, ax, ay]
# 初期条件: 近点から出発
a = 20000e3 # 半長軸 [m]
e = 0.6 # 離心率
r_p = a * (1 - e) # 近点距離
v_p = np.sqrt(mu * (2 / r_p - 1 / a)) # 近点速度（ビザビバ方程式）
state0 = [r_p, 0, 0, v_p]
T = 2 * np.pi * np.sqrt(a**3 / mu)
# 3周期分を積分
sol = solve_ivp(equations_of_motion, [0, 3 * T], state0,
method='DOP853', rtol=1e-12, atol=1e-14,
dense_output=True)
# 描画
t_plot = np.linspace(0, 3 * T, 5000)
state_plot = sol.sol(t_plot)
x_num = state_plot[0]
y_num = state_plot[1]
# 解析解との比較
nu_fine = np.linspace(0, 2 * np.pi, 1000)
p = a * (1 - e**2)
r_analytical = p / (1 + e * np.cos(nu_fine))
x_analytical = r_analytical * np.cos(nu_fine)
y_analytical = r_analytical * np.sin(nu_fine)
fig, axes = plt.subplots(1, 2, figsize=(15, 7))
# 左: 数値解の軌道
axes[0].plot(x_num / 1e6, y_num / 1e6, 'b-', linewidth=1, alpha=0.7,
label='数値解（3周回）')
axes[0].plot(x_analytical / 1e6, y_analytical / 1e6, 'r--', linewidth=2,
label='解析解（楕円）')
axes[0].plot(0, 0, 'o', color='orange', markersize=12, label='地球')
axes[0].set_xlabel('x [×1000 km]', fontsize=13)
axes[0].set_ylabel('y [×1000 km]', fontsize=13)
axes[0].set_title('数値解 vs 解析解', fontsize=14)
axes[0].set_aspect('equal')
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)
# 右: 角運動量の時間変化
vx_plot = state_plot[2]
vy_plot = state_plot[3]
h_num = x_num * vy_plot - y_num * vx_plot # h = r × v のz成分
h_theory = np.sqrt(mu * p)
axes[1].plot(t_plot / T, h_num / h_theory, 'b-', linewidth=1.5)
axes[1].set_xlabel('時間 [周期]', fontsize=13)
axes[1].set_ylabel(r'$h / h_{\mathrm{theory}}$', fontsize=13)
axes[1].set_title('角運動量の保存（数値解）', fontsize=14)
axes[1].set_ylim(0.999999, 1.000001)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"角運動量の相対誤差: {np.abs(h_num/h_theory - 1).max():.2e}")
左のグラフでは、数値積分による軌道（青の実線）と解析解の楕円（赤の破線）が完全に重なっていることがわかります。3周回分を積分しても軌道はぴったり閉じており、数値的にもケプラーの第1法則が確認できます。
右のグラフは角運動量の時間変化を示しています。理論値で規格化した $h/h_{\text{theory}}$ がほぼ1.0で一定に保たれており、角運動量保存（第2法則の基礎）が数値積分でも精度よく維持されていることがわかります。角運動量の相対誤差は $10^{-9}$ 程度であり、高精度の数値積分法（DOP853）の性能を反映しています。
このように、解析的な導出と数値シミュレーションの両面から、ケプラーの法則が万有引力から自然に導かれることを確認できました。
ケプラーの法則の適用限界
ここまでケプラーの法則の導出を見てきましたが、これらの法則が正確に成り立つのは「2つの質点が万有引力のみで相互作用する」という理想的な状況に限られます。実際の太陽系や衛星の運動では、いくつかの効果がケプラー軌道からのずれ（摂動）を引き起こします。
多体効果: 太陽系には複数の惑星が存在し、互いに重力を及ぼし合います。木星のような大質量惑星の重力は、他の惑星の軌道に測定可能な影響を与えます。このため、惑星の軌道は厳密には閉じた楕円にはならず、わずかに歳差運動します。
天体の形状: 地球は完全な球ではなく、赤道方向にわずかに膨らんだ扁平楕円体です。この非球形性はJ2摂動と呼ばれ、人工衛星の軌道面の歳差や近点の回転を引き起こします。低軌道衛星では数度/日の軌道面変化が生じるため、ミッション設計で必ず考慮されます。
一般相対性理論: ニュートン力学による予測と実際の観測のわずかなずれは、アインシュタインの一般相対性理論によって説明されます。最も有名な例は水星の近日点移動で、ニュートン力学では説明できない43秒角/世紀のずれが相対論的効果として正確に説明されます。
大気抵抗・太陽輻射圧: 低軌道衛星では大気抵抗が、高軌道衛星では太陽光の輻射圧が軌道に影響を与えます。これらは非保存力であるため、軌道エネルギーが変化し、長期的に軌道が変化します。
これらの摂動を考慮した高精度な軌道計算では、ケプラー軌道を出発点として摂動の効果を補正する手法が一般的です。そのため、ケプラーの法則は「近似」でありながらも、軌道力学の最も重要な基盤として位置づけられています。
まとめ
本記事では、万有引力の法則とニュートンの運動法則を出発点として、ケプラーの3法則を数学的に導出しました。
- 第2法則（面積速度一定）: 万有引力が中心力であるため角運動量が保存され、面積速度 $dA/dt = h/2$ が一定となることを証明しました。この法則はすべての中心力に共通する性質です。
- 第1法則（楕円軌道）: 極座標の運動方程式にビネの変数変換を施すと単振動方程式に帰着し、その解が円錐曲線 $r = p/(1+e\cos\nu)$ を与えることを示しました。束縛軌道（$e < 1$）が楕円になるのは、万有引力が逆2乗則に従うことの直接的な帰結です。
- 第3法則（$T^2 \propto a^3$）: 楕円の面積と面積速度の関係から $T^2 = (4\pi^2/\mu)a^3$ を導きました。離心率に依存しない普遍的な関係であり、天体の質量推定や衛星の軌道設計に直結する実用的な法則です。
Pythonによる可視化では、楕円軌道の形状、面積速度の一定性、太陽系惑星での第3法則の精密な成立を確認しました。数値シミュレーションによる独立検証でも、解析解と完全に一致する結果が得られました。
ケプラーの法則は、17世紀に発見された法則でありながら、現代の宇宙工学においても不可欠な基礎理論です。衛星軌道の設計、惑星間航行の計画、系外惑星の検出に至るまで、その応用範囲は広がり続けています。
次のステップとして、以下の記事も参考にしてください。