人工衛星を低軌道(LEO)から静止軌道(GEO)へ移したいとき、どうすればよいでしょうか。直感的には「上に向かって噴射すればいい」と思うかもしれません。しかし、軌道力学の世界では 「前に加速する」ことで「上に行く」 という、地上の感覚からすると意外な事実があります。車を運転していてアクセルを踏んだら高度が上がるようなものです。この直感に反する振る舞いこそが軌道力学の面白さであり、その核心に位置するのが ホーマン遷移軌道(Hohmann transfer orbit)です。
ホーマン遷移を理解すると、以下のような幅広い場面で役立ちます。
- 衛星の軌道投入設計: 低軌道から静止軌道への遷移(LEO → GEO)に必要な燃料を見積もることができます。気象衛星ひまわりや通信衛星は、まさにこの手法で静止軌道へ投入されています
- 惑星間航行の基本計画: 地球から火星や金星への探査機の軌道を概算する際の出発点となります。はやぶさ2やパーサヴィアランスの軌道計画も、最初の概算はホーマン遷移がベースです
- 宇宙ミッションの燃料予算(ΔVバジェット)の策定: ロケットの打ち上げ能力と衛星の搭載燃料を決定する際に、必要なΔV(速度増分)の下限値を与えるのがホーマン遷移です
- 軌道力学全般の基礎理解: より高度な軌道遷移手法(バイエリプティック遷移、ランベルト問題、低推力遷移など)を学ぶための前提知識となります
本記事の内容
- ホーマン遷移の直感的な理解と原理
- ΔVの丁寧な導出(ビザビバ方程式を使用)
- 遷移時間の計算(ケプラーの第3法則との関係)
- 複数の具体例(LEO → GEO、LEO → MEO、地球 → 火星)
- バイエリプティック遷移との比較 — いつホーマン遷移が最適でなくなるか
- Pythonでの軌道描画とΔV解析
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 万有引力とケプラーの法則 — 円軌道・楕円軌道の基本
- ビザビバ方程式 — 軌道上の速度計算の道具
- 軌道要素の定義 — 軌道を記述する6つのパラメータ
ホーマン遷移とは — 高速道路のインターチェンジで車線変更するイメージ
日常のアナロジー: 車線変更と軌道変更
高速道路で内側車線(低軌道)から外側車線(高軌道)に車線変更する場面を想像してみてください。いきなり真横に飛び出すのではなく、まず アクセルを踏んで速度を上げ、斜めに外側へ向かう軌跡を描き、外側車線に合流したところで 速度を合わせる のが安全で効率的です。
ホーマン遷移もまったく同じ発想です。内側の円軌道上のある点で 進行方向に加速(第1噴射)して楕円軌道に乗り、その楕円軌道を半周した後、外側の円軌道に到達した時点で再び 進行方向に加速(第2噴射)して目標の円軌道に乗ります。
ただし地上の車線変更と違う点が一つあります。地上では加速してもハンドルを切らなければ真っ直ぐ走り続けますが、宇宙では加速すると 重力によって自然に軌道が膨らみ、高度が上がっていきます。これは二体問題における楕円軌道の性質に起因しています。つまり宇宙では「アクセルを踏む = ハンドルを切る」のような効果を持つのです。
ホーマン遷移の3ステップ
この「高速道路の車線変更」を正式に定式化すると、ホーマン遷移は以下の3ステップで構成されます。
- 第1噴射(近地点でのΔV₁): 初期の円軌道上のある点で、進行方向に加速します。これにより衛星は円軌道から離脱し、遷移楕円軌道に乗ります
- 遷移楕円軌道上を半周(180°飛行): 推力なしの慣性飛行(コースティング)で楕円軌道を半周します。この間、衛星は徐々に高度を上げていきます
- 第2噴射(遠地点でのΔV₂): 遷移楕円軌道の最高点(遠地点)に到達したとき、再び進行方向に加速して、目標の円軌道速度に合わせます
この手法は1925年にドイツの工学者ヴァルター・ホーマンが著書 Die Erreichbarkeit der Himmelskörper(天体への到達可能性)で提案したもので、以来ほぼ100年にわたって宇宙ミッションの基本設計に使われ続けています。
遷移楕円軌道の幾何学的特徴
遷移楕円軌道は、初期軌道と目標軌道にちょうど 接する 楕円です。具体的には次の関係が成り立ちます。
- 近地点半径 = 初期軌道半径 $r_1$(遷移楕円が内側の円軌道に内接)
- 遠地点半径 = 目標軌道半径 $r_2$(遷移楕円が外側の円軌道に外接)
この「接する」という性質が重要です。接点では楕円軌道の速度方向と円軌道の速度方向が一致するため、噴射は純粋に速度の大きさを変えるだけで済みます。方向を変える必要がないので、エネルギーの無駄がありません。
ここまでホーマン遷移の概要を直感的に理解しました。しかし、なぜこの方法が「最も燃料効率が良い」と言えるのでしょうか。次のセクションでは、そのエネルギー的な根拠を掘り下げます。
なぜホーマン遷移は最小ΔVを与えるのか
エネルギー的な直感
「最小のΔVで軌道を変える」ということは、「最小の燃料で目的の軌道に到達する」ということです。では、なぜたった2回の噴射、しかもこの特定のタイミングで噴射するのが最適なのでしょうか。
鍵となるのは オベルト効果(Oberth effect)と呼ばれる原理です。ロケットが生み出す軌道エネルギーの増分は、噴射時の速度が大きいほど大きくなります。これは運動エネルギーが速度の2乗に比例するためです。速度 $v$ で飛んでいる物体に速度増分 $\Delta v$ を加えると、運動エネルギーの増分は次のようになります。
$$ \Delta E_k = \frac{1}{2}m(v + \Delta v)^2 – \frac{1}{2}mv^2 = mv\Delta v + \frac{1}{2}m(\Delta v)^2 $$
この式の第1項 $mv\Delta v$ に注目すると、同じ $\Delta v$ でも $v$ が大きいほど得られるエネルギー増分が大きいことがわかります。つまり、ロケットは 速く飛んでいるときに噴射するほど効率が良い のです。
楕円軌道において、速度が最大になるのは 近地点、最小になるのは 遠地点 です。エネルギー保存則と角運動量の保存からこれは証明できます。ホーマン遷移の第1噴射は近地点(内側軌道上の点)で、第2噴射は遠地点(外側軌道上の点)で行われますが、それぞれの噴射は その時点で取りうる最も効率的な場所 で行われているのです。
ΔVの最適性の条件
数学的に厳密な証明はここでは省略しますが、結論を述べます。同一平面上の2つの円軌道間を 2回のインパルス噴射 で遷移する場合、ホーマン遷移は次の条件下で最小のΔVを与えます。
$$ \frac{r_2}{r_1} < 11.94 $$
ここで $r_1$ は内側の円軌道半径、$r_2$ は外側の円軌道半径です。この閾値 $11.94$ は、ホーマン遷移のΔVとバイエリプティック遷移(3回噴射)のΔVが等しくなる軌道半径比です。
地球周りの典型的なミッションでは $r_2/r_1$ がこの値を超えることは稀です。たとえば LEO(高度200 km)から GEO(高度35,786 km)への遷移では $r_2/r_1 \approx 6.4$ であり、ホーマン遷移が最適です。一方、惑星間遷移では軌道半径比が大きくなるため、バイエリプティック遷移が有利になる場合もあります。この比較は記事の後半で詳しく取り上げます。
ホーマン遷移が最小ΔVを与える条件がわかったところで、次はその具体的なΔVの値を数式で求めていきましょう。ビザビバ方程式が威力を発揮する場面です。
ΔVの導出 — ビザビバ方程式による体系的アプローチ
準備: ビザビバ方程式と円軌道速度
ΔVを導出するために、2つの基本的な道具を用意します。
まず、ビザビバ方程式(Vis-viva equation)です。これは二体問題のエネルギー保存則から導かれる式で、楕円軌道上の任意の点における速度 $v$ を、その点の距離 $r$ と軌道の半長径 $a$ で表します。
$$ \begin{equation} v^2 = \mu \left(\frac{2}{r} – \frac{1}{a}\right) \end{equation} $$
ここで $\mu = GM$ は中心天体の重力パラメータです(地球の場合 $\mu = 3.986 \times 10^{14}$ m³/s²)。この式は「軌道のどの位置にいるか($r$)」と「軌道全体の大きさ($a$)」さえわかれば、速度が即座に求まるという強力な道具です。直感的には、距離 $r$ が小さいほど速く、半長径 $a$ が大きいほどエネルギーが高いことを表しています。ビザビバ方程式の導出についてはビザビバ方程式の記事で詳しく解説しています。
次に、円軌道速度を確認しておきます。円軌道では $a = r$(半長径 = 軌道半径)なので、ビザビバ方程式に代入すると次のようになります。
$$ v_{\text{circ}}^2 = \mu \left(\frac{2}{r} – \frac{1}{r}\right) = \frac{\mu}{r} $$
したがって、
$$ \begin{equation} v_{\text{circ}} = \sqrt{\frac{\mu}{r}} \end{equation} $$
この式は高度が上がると軌道速度が下がることを示しています。静止軌道(高度約36,000 km)の衛星は約3.1 km/sで飛んでいますが、LEO(高度200 km)の衛星は約7.8 km/sで飛んでいます。高い軌道のほうが遅いというのは、山の上のほうが重力が弱く、遠心力とのバランスに必要な速度が小さくて済むからです。
遷移楕円軌道の半長径
遷移楕円軌道は、近地点が内側軌道に接し($r_p = r_1$)、遠地点が外側軌道に接する($r_a = r_2$)楕円です。軌道要素の定義から、楕円軌道の半長径は近地点距離と遠地点距離の平均です。
$$ \begin{equation} a_t = \frac{r_p + r_a}{2} = \frac{r_1 + r_2}{2} \end{equation} $$
この式は直感的にも理解できます。遷移楕円は内側の円と外側の円の「間」を通るので、その「平均的な大きさ」は2つの円軌道半径の平均になります。
第1噴射のΔV — 内側軌道から遷移楕円へ
第1噴射は $r = r_1$ の位置で行います。噴射前の速度は内側円軌道の速度 $v_1$、噴射後の速度は遷移楕円軌道の近地点速度 $v_{t1}$ です。
噴射前(内側円軌道上)の速度:
$$ v_1 = \sqrt{\frac{\mu}{r_1}} $$
噴射後(遷移楕円の近地点)の速度:
ビザビバ方程式に $r = r_1$、$a = a_t$ を代入します。
$$ v_{t1} = \sqrt{\mu\left(\frac{2}{r_1} – \frac{1}{a_t}\right)} $$
ここで $a_t = (r_1 + r_2)/2$ を代入すると、$1/a_t = 2/(r_1 + r_2)$ なので、
$$ v_{t1} = \sqrt{\mu\left(\frac{2}{r_1} – \frac{2}{r_1 + r_2}\right)} $$
括弧内を通分して整理します。
$$ \frac{2}{r_1} – \frac{2}{r_1 + r_2} = \frac{2(r_1 + r_2) – 2r_1}{r_1(r_1 + r_2)} = \frac{2r_2}{r_1(r_1 + r_2)} $$
したがって、
$$ \begin{equation} v_{t1} = \sqrt{\frac{2\mu r_2}{r_1(r_1 + r_2)}} \end{equation} $$
第1噴射のΔV:
$$ \Delta v_1 = v_{t1} – v_1 $$
ここで $v_1 = \sqrt{\mu/r_1}$ を共通因子としてくくり出すと、
$$ \Delta v_1 = \sqrt{\frac{\mu}{r_1}} \left(\sqrt{\frac{2r_2}{r_1 + r_2}} – 1\right) $$
$\sqrt{2r_2/(r_1+r_2)} > 1$(なぜなら $r_2 > r_1$ のとき $2r_2 > r_1 + r_2$)なので $\Delta v_1 > 0$ です。これは、遷移楕円に乗るためには 加速する 必要がある(円軌道速度より速くなる)ことを意味しています。
$$ \begin{equation} \Delta v_1 = \sqrt{\frac{\mu}{r_1}} \left(\sqrt{\frac{2r_2}{r_1 + r_2}} – 1\right) \end{equation} $$
第2噴射のΔV — 遷移楕円から外側軌道へ
第2噴射は $r = r_2$(遷移楕円の遠地点)の位置で行います。
噴射前(遷移楕円の遠地点)の速度:
ビザビバ方程式に $r = r_2$、$a = a_t$ を代入します。
$$ v_{t2} = \sqrt{\mu\left(\frac{2}{r_2} – \frac{1}{a_t}\right)} $$
$a_t = (r_1 + r_2)/2$ を代入すると、
$$ v_{t2} = \sqrt{\mu\left(\frac{2}{r_2} – \frac{2}{r_1 + r_2}\right)} $$
先ほどと同様に括弧内を通分します。
$$ \frac{2}{r_2} – \frac{2}{r_1 + r_2} = \frac{2(r_1 + r_2) – 2r_2}{r_2(r_1 + r_2)} = \frac{2r_1}{r_2(r_1 + r_2)} $$
したがって、
$$ \begin{equation} v_{t2} = \sqrt{\frac{2\mu r_1}{r_2(r_1 + r_2)}} \end{equation} $$
噴射後(外側円軌道上)の速度:
$$ v_2 = \sqrt{\frac{\mu}{r_2}} $$
第2噴射のΔV:
$$ \Delta v_2 = v_2 – v_{t2} $$
$v_2 = \sqrt{\mu/r_2}$ を共通因子としてくくり出すと、
$$ \begin{equation} \Delta v_2 = \sqrt{\frac{\mu}{r_2}} \left(1 – \sqrt{\frac{2r_1}{r_1 + r_2}}\right) \end{equation} $$
$\sqrt{2r_1/(r_1+r_2)} < 1$(なぜなら $r_1 < r_2$ のとき $2r_1 < r_1 + r_2$)なので $\Delta v_2 > 0$ です。遠地点では遷移楕円上の速度が外側円軌道の速度より遅いため、ここでも 加速する 必要があります。
合計ΔV
ホーマン遷移に必要な合計ΔVは、2回の噴射の和として計算されます。
$$ \begin{equation} \Delta v_{\text{total}} = \Delta v_1 + \Delta v_2 = \sqrt{\frac{\mu}{r_1}} \left(\sqrt{\frac{2r_2}{r_1 + r_2}} – 1\right) + \sqrt{\frac{\mu}{r_2}} \left(1 – \sqrt{\frac{2r_1}{r_1 + r_2}}\right) \end{equation} $$
この式は $\mu$、$r_1$、$r_2$ の3つのパラメータだけでΔVが完全に決まることを示しています。中心天体と2つの軌道半径が決まれば、それだけで必要な燃料(の指標)がわかるという、非常に強力な結果です。
ΔVの式が得られたところで、次は遷移にかかる 時間 を求めましょう。衛星が遷移楕円を半周するのにどれだけの時間がかかるかは、ケプラーの第3法則から求められます。
遷移時間 — ケプラーの第3法則との関係
遷移時間の直感
遷移時間は「遷移楕円を半周する時間」です。ケプラーの第3法則によると、軌道の周期 $T$ は半長径 $a$ の3/2乗に比例します。
$$ T = 2\pi\sqrt{\frac{a^3}{\mu}} $$
大きな軌道ほど一周するのに時間がかかるというのは直感的にも自然です。惑星で考えると、地球は太陽の周りを365日で一周しますが、火星は687日、木星は約12年かかります。遠い惑星ほど周期が長いのは、軌道が大きいためです。
遷移時間の公式
ホーマン遷移では楕円軌道を 半周 するので、遷移時間は周期の半分です。
$$ t_{\text{transfer}} = \frac{T_t}{2} = \pi\sqrt{\frac{a_t^3}{\mu}} $$
ここで $a_t = (r_1 + r_2)/2$ を代入すると、
$$ \begin{equation} t_{\text{transfer}} = \pi\sqrt{\frac{(r_1 + r_2)^3}{8\mu}} \end{equation} $$
この式からわかる重要なポイントは、遷移時間が $r_1$ と $r_2$ の 和の3/2乗 に比例するという点です。つまり、目標軌道が遠くなるほど遷移時間は急速に長くなります。LEO → GEO は約5.3時間ですが、地球 → 火星は約8.5か月もかかります。
遷移時間とΔVのトレードオフ
ここで一つの重要な視点を提示しておきます。ホーマン遷移は最小ΔVを与えますが、必ずしも最短時間の遷移ではありません。ΔVを増やしてもよいなら、より短時間で目標軌道に到達することも可能です(ただし楕円1本では到達できず、より複雑な軌道設計が必要になります)。実際のミッション設計では、燃料(ΔV)と時間のトレードオフ を考慮して軌道を選択します。
具体的な数値で実感を得るために、次はいくつかの代表的なミッションシナリオでΔVと遷移時間を計算してみましょう。
具体例1: LEO → GEO(気象衛星・通信衛星の軌道投入)
パラメータ設定
最も代表的なホーマン遷移の例として、低軌道から静止軌道への遷移を計算します。
| パラメータ | 値 |
|---|---|
| $\mu$(地球の重力パラメータ) | $3.986 \times 10^{14}$ m³/s² |
| LEO高度 | 200 km → $r_1 = 6{,}578$ km |
| GEO高度 | 35,786 km → $r_2 = 42{,}164$ km |
| 軌道半径比 | $r_2/r_1 \approx 6.41$ |
軌道半径比が $11.94$ より小さいので、ホーマン遷移が最適な2回噴射遷移です。
計算過程
遷移楕円の半長径:
$$ a_t = \frac{6{,}578 + 42{,}164}{2} = 24{,}371 \text{ km} $$
各速度の計算:
$$ v_1 = \sqrt{\frac{3.986 \times 10^{14}}{6.578 \times 10^6}} = 7{,}784 \text{ m/s} = 7.784 \text{ km/s} $$
$$ v_{t1} = \sqrt{3.986 \times 10^{14} \times \left(\frac{2}{6.578 \times 10^6} – \frac{1}{24.371 \times 10^6}\right)} = 10{,}252 \text{ m/s} = 10.252 \text{ km/s} $$
$$ v_{t2} = \sqrt{3.986 \times 10^{14} \times \left(\frac{2}{42.164 \times 10^6} – \frac{1}{24.371 \times 10^6}\right)} = 1{,}597 \text{ m/s} = 1.597 \text{ km/s} $$
$$ v_2 = \sqrt{\frac{3.986 \times 10^{14}}{42.164 \times 10^6}} = 3{,}075 \text{ m/s} = 3.075 \text{ km/s} $$
ΔVと遷移時間:
| 項目 | 値 |
|---|---|
| $\Delta v_1$(LEOでの加速) | 2.468 km/s |
| $\Delta v_2$(GEOでの加速) | 1.478 km/s |
| $\Delta v_{\text{total}}$ | 3.946 km/s |
| 遷移時間 | 5.26 時間 |
結果の解釈
注目すべき点がいくつかあります。まず、第1噴射(2.468 km/s)が第2噴射(1.478 km/s)より大きくなっています。これは内側の円軌道速度と遷移楕円の近地点速度の差が、外側の円軌道速度と遷移楕円の遠地点速度の差よりも大きいためです。
また、遷移時間は約5.3時間と、地球を半周する程度の時間です。衛星バスの設計ではこの間のバッテリー持続時間や熱設計を考慮する必要があります。
さらに、合計ΔVの3.946 km/sは衛星の推進系にとって大きな値です。実際の静止衛星は、打ち上げロケットが一部のΔVを担い、衛星自身のアポジモータ(遠地点噴射エンジン)が残りを担うという分担が一般的です。GTO(静止遷移軌道)に直接投入するロケット(H3、Falcon 9など)を使えば、衛星が担うΔVを大幅に減らせます。
LEO → GEO以外のミッションではどうなるでしょうか。次に、GPS軌道への遷移と惑星間遷移の例を見てみましょう。
具体例2: LEO → MEO(GPS衛星の軌道投入)
GPS衛星は高度約20,200 km の中軌道(MEO: Medium Earth Orbit)に配置されています。GPS衛星のホーマン遷移を計算してみましょう。
| パラメータ | 値 |
|---|---|
| LEO高度 | 200 km → $r_1 = 6{,}578$ km |
| MEO高度(GPS) | 20,200 km → $r_2 = 26{,}571$ km |
| 軌道半径比 | $r_2/r_1 \approx 4.04$ |
遷移楕円の半長径:
$$ a_t = \frac{6{,}578 + 26{,}571}{2} = 16{,}575 \text{ km} $$
計算結果をまとめます。
| 項目 | 値 |
|---|---|
| $\Delta v_1$ | 2.053 km/s |
| $\Delta v_2$ | 1.282 km/s |
| $\Delta v_{\text{total}}$ | 3.335 km/s |
| 遷移時間 | 3.73 時間 |
GEOへの遷移と比較すると、合計ΔVは0.6 km/sほど小さく、遷移時間も1.5時間ほど短くなっています。これは目標軌道が近いためであり、直感に合う結果です。
具体例3: 地球 → 火星(惑星間ホーマン遷移)
ホーマン遷移の考え方は惑星間航行にもそのまま拡張できます。太陽を中心天体とし、地球軌道と火星軌道をそれぞれ円軌道と近似して計算します。
| パラメータ | 値 |
|---|---|
| $\mu_{\text{sun}}$ | $1.327 \times 10^{20}$ m³/s² |
| 地球軌道半径 | $r_1 = 1.496 \times 10^{11}$ m(1 AU) |
| 火星軌道半径 | $r_2 = 2.279 \times 10^{11}$ m(1.524 AU) |
| 軌道半径比 | $r_2/r_1 \approx 1.52$ |
遷移楕円の半長径:
$$ a_t = \frac{1.496 + 2.279}{2} \times 10^{11} = 1.888 \times 10^{11} \text{ m} \approx 1.262 \text{ AU} $$
計算結果をまとめます。
| 項目 | 値 |
|---|---|
| 地球の公転速度 $v_1$ | 29.78 km/s |
| 遷移楕円の近日点速度 $v_{t1}$ | 32.73 km/s |
| $\Delta v_1$(地球出発時) | 2.95 km/s |
| 遷移楕円の遠日点速度 $v_{t2}$ | 21.48 km/s |
| 火星の公転速度 $v_2$ | 24.13 km/s |
| $\Delta v_2$(火星到着時) | 2.65 km/s |
| $\Delta v_{\text{total}}$ | 5.60 km/s |
| 遷移時間 | 約 8.5 か月(258.8 日) |
この遷移時間の約8.5か月は「打ち上げウィンドウ」の計算基礎になります。地球と火星の相対位置が遷移楕円の両端と一致するタイミング(約26か月ごとに訪れる)を狙って打ち上げる必要があります。
なお、ここで求めたΔVは 太陽中心系 でのΔVです。実際の探査機は地球の重力圏内からスタートするため、地球の重力を脱出するためのΔVも必要です(パッチドコニック近似)。それでも、ホーマン遷移のΔVは惑星間ミッションの燃料予算を概算する出発点として不可欠です。
3つの具体例を通して、ホーマン遷移のΔVと遷移時間の感覚が掴めました。ここからは、ホーマン遷移が最適でなくなる場合 — バイエリプティック遷移 — について考察します。
バイエリプティック遷移との比較 — いつホーマン遷移が最適でなくなるか
バイエリプティック遷移とは
バイエリプティック遷移(Bi-elliptic transfer)は、ホーマン遷移の拡張版とも言える手法です。3回の噴射 を使い、一度目標軌道よりもはるか遠くまで飛んでから戻ってくるという、一見非効率に見える方法です。
具体的には次の3ステップです。
- 第1噴射: 内側軌道から第1遷移楕円に投入(遠地点 $r_b$ は $r_2$ より十分大きい)
- 第2噴射: 遠地点 $r_b$ で第2遷移楕円に投入(近地点が $r_2$)
- 第3噴射: 近地点 $r_2$ で目標の円軌道に投入
直感的には「遠回りするほうが効率が良い」というのは不思議に思えます。しかし、これもオベルト効果で説明できます。$r_2/r_1$ が非常に大きい場合、ホーマン遷移の第2噴射は非常に遠い地点(速度が遅い点)で行われるため、大きなΔVが必要です。一方バイエリプティック遷移では、さらに遠い中間点 $r_b$ で噴射して方向を変え、目標軌道に「落ちてくる」ことで、合計のΔVを削減できる場合があるのです。
臨界軌道半径比
ホーマン遷移とバイエリプティック遷移のΔVが等しくなる軌道半径比は、中間軌道半径 $r_b \to \infty$ の極限で次の値になります。
$$ \frac{r_2}{r_1} = 11.94 $$
これを境にして、
- $r_2/r_1 < 11.94$: ホーマン遷移 のほうがΔVが小さい
- $r_2/r_1 > 11.94$: バイエリプティック遷移 のほうがΔVが小さくなりうる($r_b$ が十分大きい場合)
ただし、バイエリプティック遷移の欠点は 遷移時間が非常に長い ことです。$r_b$ を大きくするほどΔVは節約できますが、楕円を半周する時間が非常に長くなります。実用上は、時間制約とのトレードオフでホーマン遷移が選ばれることがほとんどです。
具体的な比較
$r_2/r_1 = 15$ の場合を考えてみましょう。$r_1 = 6{,}578$ km(LEO)とすると $r_2 = 98{,}670$ km です。この場合のΔVを比較すると、
- ホーマン遷移: $\Delta v_{\text{total}} \approx 5.52$ km/s
- バイエリプティック遷移($r_b = 200{,}000$ km): $\Delta v_{\text{total}} \approx 5.41$ km/s
ΔVの差は約0.11 km/sですが、バイエリプティック遷移の遷移時間はホーマン遷移の何倍にもなります。この微小なΔV削減が、長大な遷移時間に見合うかどうかはミッション要件次第です。
バイエリプティック遷移との比較で、ホーマン遷移の適用範囲が明確になりました。それでは、ここまでの理論をPythonで実装し、軌道を可視化してみましょう。
Pythonでの実装
ΔV計算関数と軌道描画
まず、ホーマン遷移のΔVと遷移時間を計算する関数を実装し、LEO → GEO の遷移軌道を可視化します。
import numpy as np
import matplotlib.pyplot as plt
# 定数
mu = 3.986e14 # 地球の重力パラメータ [m^3/s^2]
R_earth = 6371e3 # 地球半径 [m]
def hohmann_transfer(r1, r2, mu):
"""
ホーマン遷移のΔVと遷移時間を計算する。
Parameters
----------
r1 : float
初期円軌道の半径 [m]
r2 : float
目標円軌道の半径 [m]
mu : float
中心天体の重力パラメータ [m^3/s^2]
Returns
-------
dv1 : float
第1噴射のΔV [m/s]
dv2 : float
第2噴射のΔV [m/s]
a_t : float
遷移楕円の半長径 [m]
t_transfer : float
遷移時間(半周期)[s]
"""
# 円軌道速度
v1 = np.sqrt(mu / r1)
v2 = np.sqrt(mu / r2)
# 遷移楕円の半長径
a_t = (r1 + r2) / 2
# 遷移楕円上の速度(ビザビバ方程式)
vt1 = np.sqrt(mu * (2 / r1 - 1 / a_t))
vt2 = np.sqrt(mu * (2 / r2 - 1 / a_t))
# ΔV
dv1 = vt1 - v1
dv2 = v2 - vt2
# 遷移時間(楕円軌道の半周期)
t_transfer = np.pi * np.sqrt(a_t**3 / mu)
return dv1, dv2, a_t, t_transfer
# === LEO → GEO ===
r1 = R_earth + 200e3 # LEO: 高度 200 km
r2 = R_earth + 35786e3 # GEO: 高度 35,786 km
dv1, dv2, a_t, t_transfer = hohmann_transfer(r1, r2, mu)
print("=== ホーマン遷移: LEO → GEO ===")
print(f"r1 (LEO) = {r1/1e3:.0f} km")
print(f"r2 (GEO) = {r2/1e3:.0f} km")
print(f"軌道半径比 r2/r1 = {r2/r1:.2f}")
print(f"遷移楕円の半長径 a_t = {a_t/1e3:.0f} km")
print(f"ΔV1 = {dv1/1e3:.3f} km/s")
print(f"ΔV2 = {dv2/1e3:.3f} km/s")
print(f"ΔV合計 = {(dv1 + dv2)/1e3:.3f} km/s")
print(f"遷移時間 = {t_transfer/3600:.2f} 時間")
このコードを実行すると、先ほど手計算で求めた値(ΔV合計 ≈ 3.946 km/s、遷移時間 ≈ 5.26 時間)と一致する結果が得られます。手計算と数値計算が一致することで、導出した公式の正しさが確認できます。
続いて、軌道を図示します。
import numpy as np
import matplotlib.pyplot as plt
# 定数
mu = 3.986e14
R_earth = 6371e3
r1 = R_earth + 200e3
r2 = R_earth + 35786e3
a_t = (r1 + r2) / 2
fig, ax = plt.subplots(figsize=(10, 10))
theta = np.linspace(0, 2 * np.pi, 500)
# 地球
earth = plt.Circle((0, 0), R_earth / 1e6, color='skyblue', alpha=0.8,
label='Earth')
ax.add_patch(earth)
# LEO(初期軌道)— 青い実線
ax.plot(r1 * np.cos(theta) / 1e6, r1 * np.sin(theta) / 1e6,
'b-', linewidth=2, label=f'LEO (200 km)')
# GEO(目標軌道)— 赤い実線
ax.plot(r2 * np.cos(theta) / 1e6, r2 * np.sin(theta) / 1e6,
'r-', linewidth=2, label=f'GEO (35,786 km)')
# 遷移楕円(半周)— 緑の破線
e_t = (r2 - r1) / (r1 + r2) # 離心率
theta_transfer = np.linspace(0, np.pi, 300)
r_transfer = a_t * (1 - e_t**2) / (1 + e_t * np.cos(theta_transfer))
ax.plot(r_transfer * np.cos(theta_transfer) / 1e6,
r_transfer * np.sin(theta_transfer) / 1e6,
'g--', linewidth=3, label='Transfer orbit')
# 噴射点のマーカー
dv1_val = np.sqrt(mu / r1) * (np.sqrt(2 * r2 / (r1 + r2)) - 1)
dv2_val = np.sqrt(mu / r2) * (1 - np.sqrt(2 * r1 / (r1 + r2)))
ax.plot(r1 / 1e6, 0, 'g^', markersize=15,
label=f'ΔV₁ = {dv1_val/1e3:.2f} km/s')
ax.plot(-r2 / 1e6, 0, 'rv', markersize=15,
label=f'ΔV₂ = {dv2_val/1e3:.2f} km/s')
# 速度ベクトルの矢印(噴射方向を示す)
arrow_scale = 3.0
ax.annotate('', xy=(r1 / 1e6, arrow_scale), xytext=(r1 / 1e6, 0),
arrowprops=dict(arrowstyle='->', color='green', lw=2))
ax.annotate('', xy=(-r2 / 1e6, -arrow_scale), xytext=(-r2 / 1e6, 0),
arrowprops=dict(arrowstyle='->', color='red', lw=2))
ax.set_xlabel('x [×1000 km]', fontsize=12)
ax.set_ylabel('y [×1000 km]', fontsize=12)
ax.set_title('Hohmann Transfer Orbit: LEO → GEO', fontsize=14)
ax.set_aspect('equal')
ax.legend(fontsize=11, loc='upper right')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
この図では、青い円がLEO、赤い円がGEO、緑の破線が遷移楕円軌道を示しています。三角形のマーカーは噴射点で、上向き三角(ΔV₁)がLEO上の出発点、下向き三角(ΔV₂)がGEO到達点に対応します。遷移楕円がLEOの右端に内接し、GEOの左端に外接していることが視覚的に確認できます。また、緑と赤の矢印は噴射の方向(進行方向)を示しており、2回とも加速方向であることがわかります。
複数ミッションのΔV比較
LEO → GEO だけでなく、LEO → MEO(GPS軌道)や地球 → 火星のΔVも一度に計算してみましょう。
import numpy as np
mu_earth = 3.986e14
mu_sun = 1.327e20
R_earth = 6371e3
# ミッション定義
missions = {
'LEO → MEO (GPS)': {
'r1': R_earth + 200e3,
'r2': R_earth + 20200e3,
'mu': mu_earth
},
'LEO → GEO': {
'r1': R_earth + 200e3,
'r2': R_earth + 35786e3,
'mu': mu_earth
},
'LEO → GEO×2': {
'r1': R_earth + 200e3,
'r2': R_earth + 71572e3,
'mu': mu_earth
},
'Earth → Mars': {
'r1': 1.496e11,
'r2': 2.279e11,
'mu': mu_sun
},
'Earth → Jupiter': {
'r1': 1.496e11,
'r2': 7.783e11,
'mu': mu_sun
},
}
print(f"{'ミッション':<22} {'ΔV1 [km/s]':>12} {'ΔV2 [km/s]':>12} "
f"{'ΔV合計 [km/s]':>14} {'遷移時間':>14} {'r2/r1':>8}")
print("-" * 90)
for name, params in missions.items():
r1, r2, mu = params['r1'], params['r2'], params['mu']
v1 = np.sqrt(mu / r1)
v2 = np.sqrt(mu / r2)
a_t = (r1 + r2) / 2
vt1 = np.sqrt(mu * (2 / r1 - 1 / a_t))
vt2 = np.sqrt(mu * (2 / r2 - 1 / a_t))
dv1 = vt1 - v1
dv2 = v2 - vt2
t = np.pi * np.sqrt(a_t**3 / mu)
# 遷移時間を適切な単位で表示
if t < 86400:
t_str = f"{t/3600:.2f} 時間"
elif t < 86400 * 365:
t_str = f"{t/86400:.1f} 日"
else:
t_str = f"{t/86400/365.25:.2f} 年"
ratio = r2 / r1
print(f"{name:<22} {dv1/1e3:>12.3f} {dv2/1e3:>12.3f} "
f"{(dv1+dv2)/1e3:>14.3f} {t_str:>14} {ratio:>8.2f}")
このコードの出力を見ると、いくつかの興味深い傾向がわかります。第一に、軌道半径比 $r_2/r_1$ が大きくなるほど合計ΔVが増加しますが、その増加は線形ではありません。LEO → GEO($r_2/r_1 \approx 6.4$)の合計ΔVは約3.9 km/sですが、LEO → GEO×2($r_2/r_1 \approx 11.9$)でも約4.9 km/sであり、軌道半径が2倍になってもΔVは1 km/s程度しか増えません。これはビザビバ方程式の平方根の性質に起因しています。第二に、惑星間遷移の遷移時間が地球周回軌道遷移に比べて桁違いに長いことが確認できます。これは太陽系のスケールの巨大さを反映しています。
軌道高度とΔVの関係グラフ
ΔVが目標軌道高度に対してどのように変化するかをグラフ化します。
import numpy as np
import matplotlib.pyplot as plt
mu = 3.986e14
R_earth = 6371e3
r1 = R_earth + 200e3 # LEO
# 目標軌道高度の範囲
altitudes = np.linspace(400, 100000, 1000) # km
r2_range = R_earth + altitudes * 1e3
dv1_arr = []
dv2_arr = []
dv_total_arr = []
t_arr = []
for r2 in r2_range:
v1 = np.sqrt(mu / r1)
v2 = np.sqrt(mu / r2)
a_t = (r1 + r2) / 2
vt1 = np.sqrt(mu * (2 / r1 - 1 / a_t))
vt2 = np.sqrt(mu * (2 / r2 - 1 / a_t))
dv1 = vt1 - v1
dv2 = v2 - vt2
dv1_arr.append(dv1 / 1e3)
dv2_arr.append(dv2 / 1e3)
dv_total_arr.append((dv1 + dv2) / 1e3)
t_arr.append(np.pi * np.sqrt(a_t**3 / mu) / 3600)
fig, axes = plt.subplots(1, 2, figsize=(15, 6))
# 左パネル: ΔV vs 高度
axes[0].plot(altitudes / 1e3, dv1_arr, 'b--', linewidth=1.5,
label='ΔV₁ (1st burn)')
axes[0].plot(altitudes / 1e3, dv2_arr, 'r--', linewidth=1.5,
label='ΔV₂ (2nd burn)')
axes[0].plot(altitudes / 1e3, dv_total_arr, 'k-', linewidth=2.5,
label='ΔV_total')
# 代表的な軌道を縦線で表示
axes[0].axvline(x=20.2, color='orange', linestyle=':', alpha=0.7,
label='MEO (GPS)')
axes[0].axvline(x=35.786, color='green', linestyle=':', alpha=0.7,
label='GEO')
axes[0].set_xlabel('Target altitude [×1000 km]', fontsize=12)
axes[0].set_ylabel('ΔV [km/s]', fontsize=12)
axes[0].set_title('Hohmann Transfer ΔV from LEO (200 km)', fontsize=14)
axes[0].legend(fontsize=10, loc='right')
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0, 100)
# 右パネル: 遷移時間 vs 高度
axes[1].plot(altitudes / 1e3, t_arr, 'purple', linewidth=2)
axes[1].axvline(x=20.2, color='orange', linestyle=':', alpha=0.7,
label='MEO (GPS)')
axes[1].axvline(x=35.786, color='green', linestyle=':', alpha=0.7,
label='GEO')
axes[1].set_xlabel('Target altitude [×1000 km]', fontsize=12)
axes[1].set_ylabel('Transfer time [hours]', fontsize=12)
axes[1].set_title('Hohmann Transfer Time from LEO (200 km)', fontsize=14)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(0, 100)
plt.tight_layout()
plt.show()
左パネルのΔVグラフから、いくつかの特徴が読み取れます。まず、合計ΔV(黒実線)は目標高度に対して単調増加していますが、その増加率は高度が上がるにつれて鈍化しています。これは $\Delta v \propto \sqrt{r_2/(r_1+r_2)}$ のような平方根の性質に起因します。次に、ΔV₁(青破線)とΔV₂(赤破線)の内訳を見ると、ΔV₁は高度とともに増加する一方、ΔV₂はある高度を超えると減少に転じています。これは遠い軌道ほど円軌道速度が遅くなるため、遷移楕円の遠地点速度との差が縮小するためです。右パネルの遷移時間は高度の約1.5乗に比例して増加しており、GEO到達に約5.3時間、高度100,000 kmまでは約46時間かかることが確認できます。
ホーマン遷移とバイエリプティック遷移のΔV比較
最後に、ホーマン遷移とバイエリプティック遷移のΔVを軌道半径比の関数としてプロットし、$r_2/r_1 = 11.94$ 付近での交差を確認します。
import numpy as np
import matplotlib.pyplot as plt
mu = 3.986e14
R_earth = 6371e3
r1 = R_earth + 200e3
# 軌道半径比の範囲
ratios = np.linspace(2, 25, 500)
r2_vals = r1 * ratios
# ホーマン遷移のΔV
dv_hohmann = []
for r2 in r2_vals:
v1 = np.sqrt(mu / r1)
v2 = np.sqrt(mu / r2)
a_t = (r1 + r2) / 2
vt1 = np.sqrt(mu * (2 / r1 - 1 / a_t))
vt2 = np.sqrt(mu * (2 / r2 - 1 / a_t))
dv = (vt1 - v1) + (v2 - vt2)
dv_hohmann.append(dv / 1e3)
def bielliptic_dv(r1, r2, rb, mu):
"""バイエリプティック遷移のΔVを計算"""
# 第1遷移楕円: r1 → rb
a1 = (r1 + rb) / 2
vt1_peri = np.sqrt(mu * (2 / r1 - 1 / a1))
vt1_apo = np.sqrt(mu * (2 / rb - 1 / a1))
# 第2遷移楕円: rb → r2
a2 = (r2 + rb) / 2
vt2_apo = np.sqrt(mu * (2 / rb - 1 / a2))
vt2_peri = np.sqrt(mu * (2 / r2 - 1 / a2))
v1 = np.sqrt(mu / r1)
v2 = np.sqrt(mu / r2)
dv1 = vt1_peri - v1 # 第1噴射
dv2 = vt2_apo - vt1_apo # 第2噴射(中間点)
dv3 = v2 - vt2_peri # 第3噴射
return abs(dv1) + abs(dv2) + abs(dv3)
# 複数のrb値でバイエリプティック遷移のΔVを計算
rb_factors = [2, 5, 50] # rb = rb_factor × r2
dv_bielliptic = {f: [] for f in rb_factors}
for r2 in r2_vals:
for f in rb_factors:
rb = f * r2
dv = bielliptic_dv(r1, r2, rb, mu)
dv_bielliptic[f].append(dv / 1e3)
fig, ax = plt.subplots(figsize=(12, 7))
ax.plot(ratios, dv_hohmann, 'b-', linewidth=2.5, label='Hohmann')
colors = ['orange', 'green', 'red']
for (f, dvs), c in zip(dv_bielliptic.items(), colors):
ax.plot(ratios, dvs, '--', color=c, linewidth=1.5,
label=f'Bi-elliptic ($r_b = {f} r_2$)')
ax.axvline(x=11.94, color='gray', linestyle=':', alpha=0.7,
label='$r_2/r_1 = 11.94$ (critical ratio)')
ax.set_xlabel('Orbit radius ratio $r_2/r_1$', fontsize=12)
ax.set_ylabel('Total ΔV [km/s]', fontsize=12)
ax.set_title('Hohmann vs Bi-elliptic Transfer ΔV', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_xlim(2, 25)
plt.tight_layout()
plt.show()
このグラフから、ホーマン遷移と各バイエリプティック遷移のΔVの関係が視覚的に理解できます。青い実線(ホーマン遷移)は $r_2/r_1$ に対して滑らかに増加しています。一方、バイエリプティック遷移の破線は $r_b$ の値によって異なるカーブを描きます。$r_b = 2r_2$(橙の破線)は小さな軌道半径比ではホーマンより明らかにΔVが大きいですが、軌道半径比が10を超えるあたりから差が縮まり始めます。$r_b = 50r_2$(赤の破線)では $r_2/r_1 \approx 11.94$ 付近でホーマン遷移のΔVを下回り、大きな軌道半径比ではバイエリプティック遷移のほうが効率的になることが確認できます。ただし、$r_b$ を大きくするとΔVは下がるものの遷移時間は急増するため、実用的には限界があります。
遷移楕円のアニメーション的可視化
衛星が遷移楕円上をどのように動くかを、位置の時間変化として可視化してみましょう。ケプラー方程式を数値的に解くことで、等時間間隔での衛星位置をプロットします。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq
mu = 3.986e14
R_earth = 6371e3
r1 = R_earth + 200e3
r2 = R_earth + 35786e3
a_t = (r1 + r2) / 2
e_t = (r2 - r1) / (r1 + r2)
def solve_kepler(M, e, tol=1e-10):
"""ケプラー方程式 M = E - e*sin(E) を数値的に解く"""
f = lambda E: E - e * np.sin(E) - M
return brentq(f, 0, 2 * np.pi)
# 遷移楕円上の等時間間隔点を計算
n_points = 20
t_transfer = np.pi * np.sqrt(a_t**3 / mu)
times = np.linspace(0, t_transfer, n_points)
n_mean = np.pi / t_transfer # 平均運動
positions_x = []
positions_y = []
for t in times:
M = n_mean * t # 平均近点角
E = solve_kepler(M, e_t) # 離心近点角
# 真近点角
nu = 2 * np.arctan2(np.sqrt(1 + e_t) * np.sin(E / 2),
np.sqrt(1 - e_t) * np.cos(E / 2))
r = a_t * (1 - e_t * np.cos(E))
positions_x.append(r * np.cos(nu) / 1e6)
positions_y.append(r * np.sin(nu) / 1e6)
fig, ax = plt.subplots(figsize=(10, 10))
theta = np.linspace(0, 2 * np.pi, 500)
# 地球
earth = plt.Circle((0, 0), R_earth / 1e6, color='skyblue', alpha=0.8)
ax.add_patch(earth)
# LEO と GEO
ax.plot(r1 * np.cos(theta) / 1e6, r1 * np.sin(theta) / 1e6,
'b-', linewidth=1, alpha=0.5)
ax.plot(r2 * np.cos(theta) / 1e6, r2 * np.sin(theta) / 1e6,
'r-', linewidth=1, alpha=0.5)
# 遷移楕円(全体を薄く表示)
theta_t = np.linspace(0, np.pi, 300)
r_t = a_t * (1 - e_t**2) / (1 + e_t * np.cos(theta_t))
ax.plot(r_t * np.cos(theta_t) / 1e6, r_t * np.sin(theta_t) / 1e6,
'g-', linewidth=1, alpha=0.3)
# 等時間間隔点をプロット
scatter = ax.scatter(positions_x, positions_y, c=np.linspace(0, 1, n_points),
cmap='viridis', s=80, zorder=5, edgecolors='black',
linewidths=0.5)
cbar = plt.colorbar(scatter, ax=ax, shrink=0.6, label='Normalized time')
ax.set_xlabel('x [×1000 km]', fontsize=12)
ax.set_ylabel('y [×1000 km]', fontsize=12)
ax.set_title('Satellite positions at equal time intervals\n'
'(Hohmann transfer LEO → GEO)', fontsize=14)
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
この図は、ホーマン遷移楕円上の衛星位置を等時間間隔で打点したものです。色が暗い(近地点付近)ほど点が密集し、明るい(遠地点付近)ほど点が疎になっています。これはケプラーの第2法則(面積速度一定の法則)を反映しています。近地点では速度が大きいため短時間で多くの距離を移動し、遠地点では速度が小さくなるため同じ時間ではあまり移動しません。この速度差こそが、遷移楕円軌道の本質的な特徴です。
実用上の注意点と発展的話題
軌道面変更を伴う場合
ここまでの議論では、初期軌道と目標軌道が 同一平面上 にあることを仮定しました。しかし実際のミッションでは、軌道面の傾斜角を変更する必要がある場合があります。たとえば、赤道直上のロケット射場から打ち上げない限り、LEOの軌道傾斜角は射場の緯度以上になります。GEOは赤道面上(傾斜角0°)の軌道ですから、軌道面変更のためのΔVが追加で必要です。
軌道面変更のΔVは次の式で与えられます。
$$ \Delta v_{\text{plane}} = 2v\sin\frac{\Delta i}{2} $$
ここで $v$ は噴射時の速度、$\Delta i$ は傾斜角変更量です。この式から、速度が遅い場所(遠地点)で軌道面変更を行うほうが効率的であることがわかります。そのため実際のGEO投入では、遠地点噴射で 軌道面変更と円軌道化を同時に 行うことが多いです。
有限推力の影響
ホーマン遷移の理論は 瞬間的な噴射(インパルス近似)を仮定しています。しかし実際のロケットエンジンには有限の推力があり、噴射には時間がかかります。噴射時間が軌道周期に比べて長くなると、噴射中に衛星の位置が変わるため、理想的なインパルスと比べてΔVの損失(重力損失)が生じます。
電気推進のような低推力エンジンを使う場合、この損失が顕著になります。低推力遷移ではホーマン遷移の概念をそのまま適用することはできず、連続的な推力を考慮したスパイラル軌道遷移の設計が必要です。
J2摂動の影響
地球は完全な球ではなく、赤道方向にわずかに膨らんだ回転楕円体です。この扁平率が原因で生じるJ2摂動は、軌道の昇交点赤経や近地点引数の歳差を引き起こします。ホーマン遷移の計算では二体問題を仮定しているため、長時間の遷移ではJ2摂動による軌道のずれを補正する追加のΔVが必要になる場合があります。特に遷移時間が長い(遠方の軌道への遷移の)場合に影響が大きくなります。
まとめ
本記事では、ホーマン遷移軌道の理論と計算について、直感的な理解から数式の導出、Pythonでの実装まで体系的に解説しました。
- ホーマン遷移の本質: 2回のインパルス噴射で同一平面上の円軌道間を最小のΔVで遷移する手法です。「前に加速することで高度を上げる」という軌道力学ならではの直感に反する振る舞いが鍵です
- ΔVの導出: ビザビバ方程式を用いて、$\Delta v_1 = \sqrt{\mu/r_1}(\sqrt{2r_2/(r_1+r_2)} – 1)$、$\Delta v_2 = \sqrt{\mu/r_2}(1 – \sqrt{2r_1/(r_1+r_2)})$ を導きました
- 遷移楕円の半長径: $a_t = (r_1 + r_2)/2$。内側と外側の円軌道半径の単純平均で求まります
- 遷移時間: ケプラーの第3法則から $t = \pi\sqrt{a_t^3/\mu}$。LEO → GEO で約5.3時間、地球 → 火星で約8.5か月です
- 最適性の限界: 軌道半径比 $r_2/r_1 > 11.94$ ではバイエリプティック遷移のほうがΔVを削減できる場合がありますが、遷移時間とのトレードオフがあります
- 実用上の考慮: 軌道面変更、有限推力、J2摂動などの影響により、実際のミッションではホーマン遷移の理論値からの補正が必要になります
ホーマン遷移は軌道力学の最も基本的な概念の一つですが、その理解は惑星間航行やコンステレーション設計など、はるかに高度な軌道設計問題への入口となります。
次のステップとして、以下の記事も参考にしてください。