火星に探査機を送りたいとき、「今日出発して8ヶ月後に火星に着くには、どんな軌道を飛べばよいか？」という問いに答える必要があります。出発地点（地球の位置）と到着地点（8ヶ月後の火星の位置）、そして飛行時間は決まっています。この条件から軌道を一意に決定する問題がランベルト問題（Lambert’s Problem）です。
ランベルト問題は、1761年にスイスの数学者ヨハン・ハインリヒ・ランベルトが定理を発表して以来、天体力学の中核的な問題として260年以上にわたって研究されてきました。この問題は、惑星間航行の軌道設計だけでなく、軌道決定やランデブー問題など、宇宙工学の幅広い場面で登場します。
ランベルト問題を理解すると、以下のような応用が見えてきます。
- 惑星間遷移軌道の設計: 地球から火星、木星などへの最適な打ち上げ時期と飛行経路の決定
- ポークチョッププロット: 出発日と到着日のあらゆる組み合わせに対するデルタVの可視化
- ランデブー・ドッキング: 宇宙ステーションや他の衛星に接近するための軌道計算
本記事の内容
- ランベルト問題の定式化
- ランベルトの定理と幾何学的意味
- ユニバーサル変数による定式化
- ニュートン法による反復解法
- Pythonでのランベルトソルバーの実装
- 火星遷移軌道の計算例
- ポークチョッププロットの可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ランベルト問題の定式化
問題の定義
ランベルト問題を正確に述べると、次のようになります。
与えられた条件: – 位置ベクトル $\bm{r}_1$（出発点） – 位置ベクトル $\bm{r}_2$（到着点） – 飛行時間 $\Delta t = t_2 – t_1$ – 中心天体の重力パラメータ $\mu$ – 軌道の方向（順行/逆行）
求めるもの: – 出発点での速度ベクトル $\bm{v}_1$ – 到着点での速度ベクトル $\bm{v}_2$
$\bm{v}_1$ と $\bm{v}_2$ が求まれば、軌道の6要素が全て決定されます。二体問題の解は位置と速度の組（6自由度）で一意に定まるからです。
なぜランベルト問題が重要か
ホーマン遷移は「2つの円軌道間の遷移」という特殊な場合の解です。それに対してランベルト問題は、任意の2点間の遷移を任意の飛行時間で扱えるため、はるかに一般的です。
惑星間航行では、出発惑星と到着惑星の位置は時刻によって変わるため、打ち上げ日と到着日のあらゆる組み合わせについてランベルト問題を解き、最もデルタVが少ない（つまり推進剤が少なくて済む）組み合わせを見つけます。この計算結果を2次元マップにしたものがポークチョッププロットです。
ランベルト問題の物理的意味が理解できたところで、次にこの問題の解の存在と性質を保証するランベルトの定理を紹介しましょう。
ランベルトの定理
定理の内容
ランベルトの定理は、二体問題における飛行時間が3つの幾何学的量だけで決まることを主張します。
ランベルトの定理: ケプラー軌道上の2点間の飛行時間 $\Delta t$ は、以下の3つの量のみに依存する。
- $r_1 + r_2$: 2つの動径の和
- $c = |\bm{r}_2 – \bm{r}_1|$: 弦の長さ（2点間の距離）
- $a$: 軌道の半長軸
この定理の重要性は、飛行時間が軌道の「形」と「大きさ」の情報のみで決まり、軌道の「向き」（傾斜角や昇交点赤経）に依存しないことを示している点です。
幾何学的パラメータ
2つの位置ベクトル $\bm{r}_1$ と $\bm{r}_2$ から、以下の幾何学的量を定義します。
$$ r_1 = |\bm{r}_1|, \quad r_2 = |\bm{r}_2|, \quad c = |\bm{r}_2 – \bm{r}_1| $$
半周長（semi-perimeter）を
$$ s = \frac{r_1 + r_2 + c}{2} $$
と定義します。これは2つの位置ベクトルと弦で作る三角形の半周長です。
さらに、ランベルトの定理を表す際に便利な補助変数 $\alpha$ と $\beta$ を次のように定義します。
$$ \sin^2\frac{\alpha}{2} = \frac{s}{2a}, \quad \sin^2\frac{\beta}{2} = \frac{s-c}{2a} $$
飛行時間の公式
ランベルトの定理による飛行時間の公式は、楕円軌道の場合
$$ \boxed{\sqrt{\mu}\,\Delta t = a^{3/2}\left[(\alpha – \sin\alpha) – (\beta – \sin\beta)\right]} $$
この式の形はケプラー方程式に似ています。実際、$\alpha$ と $\beta$ はそれぞれ出発点と到着点に対応する離心近点角の関数であり、ケプラー方程式の一般化と見なせます。
放物線（$a \to \infty$）の場合は
$$ \sqrt{\mu}\,\Delta t = \frac{1}{3}\left[s^{3/2} – (s-c)^{3/2}\right] $$
双曲線（$a < 0$）の場合は三角関数が双曲線関数に置き換わります。
$$ \sqrt{\mu}\,\Delta t = |a|^{3/2}\left[(\sinh\alpha_h – \alpha_h) – (\sinh\beta_h – \beta_h)\right] $$
ここで $\sinh^2(\alpha_h/2) = s/(2|a|)$, $\sinh^2(\beta_h/2) = (s-c)/(2|a|)$ です。
ランベルトの定理は飛行時間と半長軸の関係を与えますが、実際に解を求めるには半長軸 $a$ を未知数として反復法で解く必要があります。次に、より一般的で数値的に安定なユニバーサル変数による定式化を紹介します。
ユニバーサル変数による定式化
ユニバーサル変数の導入
楕円・放物線・双曲線の3つの場合を統一的に扱うために、ユニバーサル変数 $z$ を導入します。$z$ は半長軸 $a$ と関連する無次元パラメータで、次のように定義されます。
$$ z = \frac{\alpha^2}{a} \quad \text{（楕円）}, \quad z = -\frac{\alpha_h^2}{|a|} \quad \text{（双曲線）} $$
しかし実用上は、$z$ を直接的に扱う代わりに、Stumpff関数と呼ばれる補助関数を導入する方が便利です。
Stumpff関数
Stumpff関数 $c_2(z)$ と $c_3(z)$ は次のように定義されます。
$$ c_2(z) = \begin{cases} \frac{1 – \cos\sqrt{z}}{z} & (z > 0) \\[4pt] \frac{\cosh\sqrt{-z} – 1}{-z} & (z < 0) \\[4pt] \frac{1}{2} & (z = 0) \end{cases} $$
$$ c_3(z) = \begin{cases} \frac{\sqrt{z} – \sin\sqrt{z}}{z^{3/2}} & (z > 0) \\[4pt] \frac{\sinh\sqrt{-z} – \sqrt{-z}}{(-z)^{3/2}} & (z < 0) \\[4pt] \frac{1}{6} & (z = 0) \end{cases} $$
$z > 0$ は楕円軌道、$z = 0$ は放物線軌道、$z < 0$ は双曲線軌道に対応します。$z = 0$ のとき $c_2 = 1/2$, $c_3 = 1/6$ となるのは、テイラー展開の最低次の項と一致します。
ランベルト問題のユニバーサル変数形式
2点間の移動角 $\Delta\theta$ は、位置ベクトルから計算します。
$$ \cos\Delta\theta = \frac{\bm{r}_1 \cdot \bm{r}_2}{r_1 r_2} $$
軌道の方向（順行/逆行）により $\Delta\theta$ の範囲が決まります。
$A$ パラメータを次のように定義します。
$$ A = \sin\Delta\theta \sqrt{\frac{r_1 r_2}{1 – \cos\Delta\theta}} $$
$\Delta\theta = 180°$（対向する2点）の場合、$A = 0$ となり問題は退化します（無限個の解が存在）。
飛行時間 $\Delta t$ は $z$ の関数として次のように表されます。
$$ \sqrt{\mu}\,\Delta t = \left[\frac{y(z)}{c_2(z)}\right]^{3/2} c_3(z) + A\sqrt{y(z)} $$
ここで
$$ y(z) = r_1 + r_2 + A\frac{z\,c_3(z) – 1}{\sqrt{c_2(z)}} $$
この方程式は $z$ についての超越方程式であり、ニュートン法などの反復法で解きます。
ニュートン法による解法
$F(z) = 0$ の形に整理します。
$$ F(z) = \left[\frac{y(z)}{c_2(z)}\right]^{3/2} c_3(z) + A\sqrt{y(z)} – \sqrt{\mu}\,\Delta t $$
ニュートン法の更新則は
$$ z_{n+1} = z_n – \frac{F(z_n)}{F'(z_n)} $$
$F'(z)$ の計算は解析的に行えますが、やや複雑なため数値微分で代替することもあります。初期値は $z_0 = 0$（放物線に対応）から始めるのが安定した選択です。
$z$ が収束したら、ラグランジュ係数 $f, g, \dot{g}$ を計算して速度ベクトルを求めます。
$$ f = 1 – \frac{y}{r_1}, \quad g = A\sqrt{\frac{y}{\mu}}, \quad \dot{g} = 1 – \frac{y}{r_2} $$
速度ベクトルは
$$ \bm{v}_1 = \frac{\bm{r}_2 – f\bm{r}_1}{g}, \quad \bm{v}_2 = \frac{\dot{g}\bm{r}_2 – \bm{r}_1}{g} $$
理論的な準備が整いましたので、次にPythonでランベルトソルバーを実装しましょう。
Pythonによるランベルトソルバーの実装
Stumpff関数とソルバー本体
import numpy as np
def stumpff_c2(z):
"""Stumpff関数 c2(z)"""
if isinstance(z, np.ndarray):
result = np.zeros_like(z)
pos = z > 1e-10
neg = z < -1e-10
zero = ~pos & ~neg
result[pos] = (1 - np.cos(np.sqrt(z[pos]))) / z[pos]
result[neg] = (np.cosh(np.sqrt(-z[neg])) - 1) / (-z[neg])
result[zero] = 0.5
return result
else:
if z > 1e-10:
return (1 - np.cos(np.sqrt(z))) / z
elif z < -1e-10:
return (np.cosh(np.sqrt(-z)) - 1) / (-z)
else:
return 0.5
def stumpff_c3(z):
"""Stumpff関数 c3(z)"""
if isinstance(z, np.ndarray):
result = np.zeros_like(z)
pos = z > 1e-10
neg = z < -1e-10
zero = ~pos & ~neg
sz_pos = np.sqrt(z[pos])
result[pos] = (sz_pos - np.sin(sz_pos)) / (sz_pos**3)
sz_neg = np.sqrt(-z[neg])
result[neg] = (np.sinh(sz_neg) - sz_neg) / (sz_neg**3)
result[zero] = 1.0 / 6.0
return result
else:
if z > 1e-10:
sz = np.sqrt(z)
return (sz - np.sin(sz)) / (sz**3)
elif z < -1e-10:
sz = np.sqrt(-z)
return (np.sinh(sz) - sz) / (sz**3)
else:
return 1.0 / 6.0
def lambert_solver(r1_vec, r2_vec, dt, mu, prograde=True):
"""
ランベルト問題のソルバー（ユニバーサル変数法）
Parameters
----------
r1_vec : 出発位置ベクトル [m]
r2_vec : 到着位置ベクトル [m]
dt : 飛行時間 [s]
mu : 重力パラメータ [m^3/s^2]
prograde : 順行軌道ならTrue
Returns
-------
v1_vec : 出発速度ベクトル [m/s]
v2_vec : 到着速度ベクトル [m/s]
"""
r1 = np.linalg.norm(r1_vec)
r2 = np.linalg.norm(r2_vec)
# 移動角の計算
cos_dtheta = np.dot(r1_vec, r2_vec) / (r1 * r2)
cos_dtheta = np.clip(cos_dtheta, -1, 1)
# 軌道方向の判定（z成分の外積で判定）
cross = np.cross(r1_vec, r2_vec)
if prograde:
if cross[2] < 0:
dtheta = 2 * np.pi - np.arccos(cos_dtheta)
else:
dtheta = np.arccos(cos_dtheta)
else:
if cross[2] >= 0:
dtheta = 2 * np.pi - np.arccos(cos_dtheta)
else:
dtheta = np.arccos(cos_dtheta)
# Aパラメータ
A = np.sin(dtheta) * np.sqrt(r1 * r2 / (1 - cos_dtheta))
if abs(A) < 1e-14:
raise ValueError("対向する2点（Δθ=180°）: 解が退化します")
# ニュートン法
z = 0.0 # 初期値（放物線）
max_iter = 100
tol = 1e-10
for iteration in range(max_iter):
c2 = stumpff_c2(z)
c3 = stumpff_c3(z)
y = r1 + r2 + A * (z * c3 - 1) / np.sqrt(c2)
if y < 0:
# yが負になる場合、zを調整
z = z + 0.1
continue
F = (y / c2)**1.5 * c3 + A * np.sqrt(y) - np.sqrt(mu) * dt
# 数値微分
dz = 1e-7 * max(abs(z), 1)
c2p = stumpff_c2(z + dz)
c3p = stumpff_c3(z + dz)
yp = r1 + r2 + A * ((z+dz) * c3p - 1) / np.sqrt(c2p)
if yp < 0:
yp = abs(yp)
Fp = (yp / c2p)**1.5 * c3p + A * np.sqrt(yp) - np.sqrt(mu) * dt
dFdz = (Fp - F) / dz
if abs(dFdz) < 1e-20:
break
z_new = z - F / dFdz
if abs(z_new - z) < tol:
z = z_new
break
z = z_new
# ラグランジュ係数
c2 = stumpff_c2(z)
c3 = stumpff_c3(z)
y = r1 + r2 + A * (z * c3 - 1) / np.sqrt(c2)
f = 1 - y / r1
g = A * np.sqrt(y / mu)
g_dot = 1 - y / r2
# 速度ベクトル
v1_vec = (r2_vec - f * r1_vec) / g
v2_vec = (g_dot * r2_vec - r1_vec) / g
return v1_vec, v2_vec
検証: 既知の軌道での確認
実装したソルバーを検証するため、既知の楕円軌道上の2点を取り、ソルバーで得られた速度が元の速度と一致するか確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# 上のstumpff_c2, stumpff_c3, lambert_solver関数を使用
mu = 3.986e14 # 地球の重力パラメータ
R_earth = 6378e3
# テスト軌道: 楕円軌道
a_test = R_earth + 5000e3 # 半長軸
e_test = 0.3 # 離心率
p_test = a_test * (1 - e_test**2)
# 出発点: 真近点角 30°
nu1 = np.radians(30)
r1 = p_test / (1 + e_test * np.cos(nu1))
r1_vec = np.array([r1 * np.cos(nu1), r1 * np.sin(nu1), 0])
# 速度（ビザビバ方程式 + 軌道力学）
h = np.sqrt(mu * p_test)
vr1 = mu / h * e_test * np.sin(nu1)
vt1 = mu / h * (1 + e_test * np.cos(nu1))
v1_true = np.array([
vr1 * np.cos(nu1) - vt1 * np.sin(nu1),
vr1 * np.sin(nu1) + vt1 * np.cos(nu1),
0
])
# 到着点: 真近点角 150°（数値積分で飛行時間を計算）
nu2 = np.radians(150)
r2 = p_test / (1 + e_test * np.cos(nu2))
r2_vec = np.array([r2 * np.cos(nu2), r2 * np.sin(nu2), 0])
# 飛行時間の計算（ケプラー方程式）
def true_to_eccentric(nu, e):
return 2 * np.arctan(np.sqrt((1 - e) / (1 + e)) * np.tan(nu / 2))
E1 = true_to_eccentric(nu1, e_test)
E2 = true_to_eccentric(nu2, e_test)
M1 = E1 - e_test * np.sin(E1)
M2 = E2 - e_test * np.sin(E2)
n = np.sqrt(mu / a_test**3)
dt = (M2 - M1) / n
print(f"=== テスト条件 ===")
print(f"半長軸 a = {a_test/1000:.0f} km")
print(f"離心率 e = {e_test}")
print(f"出発真近点角: {np.degrees(nu1):.0f}°")
print(f"到着真近点角: {np.degrees(nu2):.0f}°")
print(f"飛行時間: {dt:.1f} s ({dt/60:.1f} min)")
# ランベルト問題を解く
v1_lambert, v2_lambert = lambert_solver(r1_vec, r2_vec, dt, mu)
print(f"\n=== 速度の比較 ===")
print(f"真の v1 = {v1_true/1000} km/s")
print(f"Lambert v1 = {v1_lambert/1000} km/s")
print(f"速度誤差 = {np.linalg.norm(v1_lambert - v1_true):.6f} m/s")
# 軌道の可視化
theta_full = np.linspace(0, 2*np.pi, 500)
r_full = p_test / (1 + e_test * np.cos(theta_full))
x_full = r_full * np.cos(theta_full) / 1000
y_full = r_full * np.sin(theta_full) / 1000
# ランベルト解の軌道を数値積分で確認
def eom(t, state):
rv = state[:3]
vv = state[3:]
rm = np.linalg.norm(rv)
return np.concatenate([vv, -mu / rm**3 * rv])
sol = solve_ivp(eom, [0, dt], np.concatenate([r1_vec, v1_lambert]),
method='RK45', rtol=1e-12, atol=1e-14,
t_eval=np.linspace(0, dt, 300))
fig, ax = plt.subplots(figsize=(10, 10))
# 地球
earth = plt.Circle((0, 0), R_earth/1000, color='#4FC3F7', alpha=0.4)
ax.add_patch(earth)
# 元の楕円軌道
ax.plot(x_full, y_full, '--', color='gray', linewidth=1, alpha=0.5,
label='元の楕円軌道')
# ランベルト解の軌道
ax.plot(sol.y[0]/1000, sol.y[1]/1000, '-', color='#F44336',
linewidth=2.5, label='ランベルト解')
# 出発点・到着点
ax.plot(r1_vec[0]/1000, r1_vec[1]/1000, 'o', color='#4CAF50',
markersize=12, markeredgecolor='black', markeredgewidth=1,
label=f'出発点 (ν={np.degrees(nu1):.0f}°)')
ax.plot(r2_vec[0]/1000, r2_vec[1]/1000, 's', color='#2196F3',
markersize=12, markeredgecolor='black', markeredgewidth=1,
label=f'到着点 (ν={np.degrees(nu2):.0f}°)')
# 焦点
ax.plot(0, 0, 'o', color='#FFC107', markersize=12,
markeredgecolor='black', markeredgewidth=1, label='焦点')
ax.set_xlabel('x [km]', fontsize=12)
ax.set_ylabel('y [km]', fontsize=12)
ax.set_title('ランベルト問題の解の検証', fontsize=14)
ax.legend(fontsize=11, loc='upper right')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
この検証結果から、ソルバーの精度を確認できます。
- ランベルト解の速度と真の速度の誤差は $10^{-6}$ m/s以下であり、数値的に十分な精度で解が得られています。これはニュートン法の収束精度（$10^{-10}$）を反映しています。
- ランベルト解の軌道（赤い実線）が元の楕円軌道（灰色の破線）と完全に重なっていることが視覚的に確認できます。ソルバーは正しく動作しています。
- ランベルト問題は既知の軌道を「再発見」することが確認されました。つまり、出発点と到着点の位置、および飛行時間の情報だけから、元の軌道が一意に復元されるのです。
飛行時間を変えた場合の軌道の変化
同じ2点間でも飛行時間を変えると、異なる軌道が得られます。これを可視化してみましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# stumpff_c2, stumpff_c3, lambert_solver関数を使用（上で定義済み）
mu = 3.986e14
R_earth = 6378e3
# 固定の出発点と到着点
r1_vec = np.array([8000e3, 0, 0])
r2_vec = np.array([-4000e3, 7000e3, 0])
r1 = np.linalg.norm(r1_vec)
r2 = np.linalg.norm(r2_vec)
# 飛行時間を変えてランベルト問題を解く
flight_times = [2000, 3000, 4000, 5000, 7000, 10000] # 秒
colors = ['#F44336', '#FF9800', '#FFEB3B', '#4CAF50', '#2196F3', '#9C27B0']
fig, ax = plt.subplots(figsize=(10, 10))
# 地球
earth = plt.Circle((0, 0), R_earth/1000, color='#4FC3F7', alpha=0.4)
ax.add_patch(earth)
for dt, color in zip(flight_times, colors):
try:
v1, v2 = lambert_solver(r1_vec, r2_vec, dt, mu)
# 数値積分で軌道を描画
def eom(t, state):
rv = state[:3]
vv = state[3:]
rm = np.linalg.norm(rv)
return np.concatenate([vv, -mu / rm**3 * rv])
sol = solve_ivp(eom, [0, dt],
np.concatenate([r1_vec, v1]),
method='RK45', rtol=1e-10, atol=1e-12,
t_eval=np.linspace(0, dt, 300))
# デルタV（円軌道からの出発を仮定）
v_circ1 = np.sqrt(mu / r1)
dv1 = np.linalg.norm(v1) - v_circ1
ax.plot(sol.y[0]/1000, sol.y[1]/1000, '-', color=color,
linewidth=2, label=f'Δt = {dt} s ({dt/60:.0f} min)')
except Exception as e:
print(f"dt = {dt} s: {e}")
# 出発点・到着点
ax.plot(r1_vec[0]/1000, r1_vec[1]/1000, 'o', color='#4CAF50',
markersize=15, markeredgecolor='black', markeredgewidth=1.5,
zorder=5, label='出発点')
ax.plot(r2_vec[0]/1000, r2_vec[1]/1000, 's', color='#E91E63',
markersize=15, markeredgecolor='black', markeredgewidth=1.5,
zorder=5, label='到着点')
ax.set_xlabel('x [km]', fontsize=12)
ax.set_ylabel('y [km]', fontsize=12)
ax.set_title('同じ2点間でも飛行時間によって軌道が変わる', fontsize=14)
ax.legend(fontsize=10, loc='upper left')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
この図から、ランベルト問題の本質的な性質が視覚的に理解できます。
- 飛行時間が短い軌道ほど「直線的」に見え、速度（エネルギー）が大きくなります。時間が短いほど速く飛ぶ必要があるため、必要なデルタVも大きくなります。
- 飛行時間が長い軌道は大きく迂回する弧を描きます。ゆっくり飛ぶことでデルタVは減りますが、ミッション期間が長くなります。
- ある飛行時間に対して唯一の軌道が決まる（同じ周回数の場合）ことが確認できます。これがランベルト問題の一意性です。ただし、1周以上する「マルチレボリューション」の場合は複数解が存在し得ます。
ランベルトソルバーが正しく動作することを確認できました。次に、このソルバーを使って実際の惑星間航行の問題を解いてみましょう。
惑星間遷移軌道の計算
地球-火星遷移の設定
地球と火星の間の遷移軌道をランベルト問題として解きます。ここでは簡略化のため、地球と火星の軌道を同一平面上の円軌道と仮定します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# stumpff_c2, stumpff_c3, lambert_solver関数を使用
# 太陽の重力パラメータ
mu_sun = 1.327e20 # m^3/s^2
# 惑星の軌道半径（円軌道近似）
r_earth = 1.496e11 # 1 AU [m]
r_mars = 2.279e11 # 1.524 AU [m]
# 地球の公転周期と角速度
T_earth = 2 * np.pi * np.sqrt(r_earth**3 / mu_sun)
omega_earth = 2 * np.pi / T_earth
# 火星の公転周期と角速度
T_mars = 2 * np.pi * np.sqrt(r_mars**3 / mu_sun)
omega_mars = 2 * np.pi / T_mars
# 地球出発時の位置
theta_earth_0 = 0 # 出発時の地球の角度
r1_vec = r_earth * np.array([np.cos(theta_earth_0),
np.sin(theta_earth_0), 0])
# ホーマン遷移との比較のため、ホーマン遷移時間を計算
a_hohmann = (r_earth + r_mars) / 2
dt_hohmann = np.pi * np.sqrt(a_hohmann**3 / mu_sun)
print(f"ホーマン遷移時間: {dt_hohmann/86400:.1f} 日")
# ホーマン遷移時に火星がいるべき位置
theta_mars_arrival_hohmann = omega_mars * dt_hohmann
# 出発時の火星の角度（到着時にdelta_thetaだけ進む）
theta_mars_0_hohmann = np.pi - theta_mars_arrival_hohmann
# 異なる飛行時間でランベルト問題を解く
flight_days = [180, 210, 259, 300, 360] # 日
colors = ['#F44336', '#FF9800', '#4CAF50', '#2196F3', '#9C27B0']
fig, ax = plt.subplots(figsize=(10, 10))
# 地球と火星の軌道
theta_orbit = np.linspace(0, 2*np.pi, 500)
ax.plot(r_earth * np.cos(theta_orbit) / 1e11,
r_earth * np.sin(theta_orbit) / 1e11,
'--', color='#64B5F6', linewidth=1, label='地球軌道')
ax.plot(r_mars * np.cos(theta_orbit) / 1e11,
r_mars * np.sin(theta_orbit) / 1e11,
'--', color='#EF9A9A', linewidth=1, label='火星軌道')
for dt_days, color in zip(flight_days, colors):
dt = dt_days * 86400 # 秒に変換
# 到着時の火星の位置（出発時の火星角度はホーマン遷移用を使用）
theta_mars = theta_mars_0_hohmann + omega_mars * dt
r2_vec = r_mars * np.array([np.cos(theta_mars),
np.sin(theta_mars), 0])
try:
v1, v2 = lambert_solver(r1_vec, r2_vec, dt, mu_sun)
# 軌道を数値積分
def eom(t, state):
rv = state[:3]
vv = state[3:]
rm = np.linalg.norm(rv)
return np.concatenate([vv, -mu_sun / rm**3 * rv])
sol = solve_ivp(eom, [0, dt],
np.concatenate([r1_vec, v1]),
method='RK45', rtol=1e-10, atol=1e-12,
t_eval=np.linspace(0, dt, 500))
# デルタV計算
v_earth = np.sqrt(mu_sun / r_earth)
v_earth_vec = np.array([0, v_earth, 0]) # 地球の公転速度
dv1 = np.linalg.norm(v1 - v_earth_vec)
v_mars = np.sqrt(mu_sun / r_mars)
theta_mars_rad = theta_mars
v_mars_vec = v_mars * np.array([-np.sin(theta_mars_rad),
np.cos(theta_mars_rad), 0])
dv2 = np.linalg.norm(v2 - v_mars_vec)
ax.plot(sol.y[0]/1e11, sol.y[1]/1e11, '-', color=color,
linewidth=2,
label=f'{dt_days}日 (ΔV={dv1/1000:.1f}+{dv2/1000:.1f} km/s)')
# 到着点
ax.plot(r2_vec[0]/1e11, r2_vec[1]/1e11, 'o', color=color,
markersize=8, markeredgecolor='black', markeredgewidth=0.5)
except Exception as e:
print(f"dt = {dt_days} days: {e}")
# 太陽、地球出発点
ax.plot(0, 0, 'o', color='#FFC107', markersize=20,
markeredgecolor='black', markeredgewidth=1, label='太陽')
ax.plot(r1_vec[0]/1e11, r1_vec[1]/1e11, 'o', color='#2196F3',
markersize=12, markeredgecolor='black', markeredgewidth=1,
label='地球（出発）')
ax.set_xlabel('x [×10¹¹ m]', fontsize=12)
ax.set_ylabel('y [×10¹¹ m]', fontsize=12)
ax.set_title('地球→火星 遷移軌道（ランベルト問題）', fontsize=14)
ax.legend(fontsize=9, loc='upper left')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)
plt.tight_layout()
plt.show()
この図から、惑星間遷移軌道の特性が読み取れます。
- ホーマン遷移（約259日）に近い飛行時間が最もデルタVが小さいことがわかります。ホーマン遷移は出発軌道と到着軌道の両方に接する楕円であり、エネルギー効率が最適です。
- 飛行時間がホーマン遷移時間より短い場合、軌道は「内回り」になり、より高いエネルギー（大きなデルタV）が必要です。時間が短いほど速く飛ぶ必要があるからです。
- 飛行時間が長い場合は「外回り」の軌道となり、こちらもデルタVが増加します。遠回りするために余分なエネルギーが必要です。
- 到着時のデルタV（火星軌道速度との差）も考慮する必要があるため、合計デルタVの最適化には出発と到着の両方のコストを考えなければなりません。
ポークチョッププロットの生成
最後に、出発日と飛行時間を系統的に変えてデルタVを計算し、ポークチョッププロットを作成します。
import numpy as np
import matplotlib.pyplot as plt
# stumpff_c2, stumpff_c3, lambert_solver関数を使用
mu_sun = 1.327e20
r_earth = 1.496e11
r_mars = 2.279e11
T_earth = 2 * np.pi * np.sqrt(r_earth**3 / mu_sun)
T_mars = 2 * np.pi * np.sqrt(r_mars**3 / mu_sun)
omega_earth = 2 * np.pi / T_earth
omega_mars = 2 * np.pi / T_mars
# 出発日の範囲（地球の角度で表現）
n_dep = 60
n_tof = 60
departure_days = np.linspace(0, 600, n_dep) # 出発日 [day]
tof_days = np.linspace(120, 400, n_tof) # 飛行時間 [day]
dv_total = np.full((n_tof, n_dep), np.nan)
for i, dep_day in enumerate(departure_days):
# 出発時の地球位置
theta_e = omega_earth * dep_day * 86400
r1 = r_earth * np.array([np.cos(theta_e), np.sin(theta_e), 0])
v_earth = np.sqrt(mu_sun / r_earth)
v_earth_vec = v_earth * np.array([-np.sin(theta_e),
np.cos(theta_e), 0])
for j, tof in enumerate(tof_days):
dt = tof * 86400
# 到着時の火星位置
theta_m = omega_mars * (dep_day + tof) * 86400
r2 = r_mars * np.array([np.cos(theta_m), np.sin(theta_m), 0])
v_mars = np.sqrt(mu_sun / r_mars)
v_mars_vec = v_mars * np.array([-np.sin(theta_m),
np.cos(theta_m), 0])
try:
v1, v2 = lambert_solver(r1, r2, dt, mu_sun)
dv1 = np.linalg.norm(v1 - v_earth_vec)
dv2 = np.linalg.norm(v2 - v_mars_vec)
dv = (dv1 + dv2) / 1000 # km/sに変換
if dv < 50: # 現実的な範囲に制限
dv_total[j, i] = dv
except Exception:
pass
fig, ax = plt.subplots(figsize=(12, 8))
# 等高線プロット
levels = np.arange(5, 35, 1)
contour = ax.contourf(departure_days, tof_days, dv_total,
levels=levels, cmap='RdYlGn_r')
contour_lines = ax.contour(departure_days, tof_days, dv_total,
levels=[8, 10, 12, 15, 20, 25, 30],
colors='black', linewidths=0.5)
ax.clabel(contour_lines, fmt='%.0f', fontsize=9)
cbar = plt.colorbar(contour, ax=ax, label='合計ΔV [km/s]')
ax.set_xlabel('出発日 [day from epoch]', fontsize=12)
ax.set_ylabel('飛行時間 [day]', fontsize=12)
ax.set_title('地球→火星 ポークチョッププロット（C3 = ΔV₁ + ΔV₂）',
fontsize=14)
ax.grid(True, alpha=0.2)
plt.tight_layout()
plt.show()
# 最小ΔVの場所を表示
min_idx = np.nanargmin(dv_total)
min_j, min_i = np.unravel_index(min_idx, dv_total.shape)
print(f"\n=== 最適遷移 ===")
print(f"出発日: {departure_days[min_i]:.0f} day")
print(f"飛行時間: {tof_days[min_j]:.0f} day")
print(f"合計ΔV: {dv_total[min_j, min_i]:.2f} km/s")
ポークチョッププロットは惑星間ミッション設計の最も重要なツールの一つです。この図から読み取れることは非常に多いです。
- デルタVが最小の「谷」が特定の出発日付近に存在します。これは地球と火星の相対配置（シナジック周期約780日≈26ヶ月ごと）に対応しており、約26ヶ月に1回の打ち上げウィンドウがあることを意味します。実際の火星探査ミッション（Mars 2020, Hope, 天問1号など）は全てこのウィンドウに合わせて打ち上げられています。
- 最適な飛行時間は約250日前後であり、これはホーマン遷移時間（約259日）に近い値です。しかし、ポークチョッププロットの谷は広がりを持っているため、最適から多少外れた飛行時間でも合理的なデルタVで遷移可能です。
- 等高線の形状が「ポークチョップ（豚の骨付き肉）」に似ていることから、この名前が付けられました。2つのローブ（短時間遷移と長時間遷移）の間に最小値の谷があります。
- 出発日を最適から大きく外すと、必要なデルタVが急増します。打ち上げウィンドウを逃すと、次の機会まで約26ヶ月待たなければなりません。
まとめ
本記事では、ランベルト問題の理論から実装、そして惑星間航行への応用までを解説しました。
- ランベルト問題は「2点間を指定時間で移動する軌道を求める」問題であり、軌道力学の最も基本的な境界値問題の一つ
- ランベルトの定理により、飛行時間は $r_1+r_2$, 弦の長さ $c$, 半長軸 $a$ のみに依存する
- ユニバーサル変数法とStumpff関数により、楕円・放物線・双曲線を統一的に扱える
- ニュートン法で反復的に解くことで、高精度な解が得られる
- ポークチョッププロットは出発日と飛行時間の全組み合わせについてデルタVを計算した2次元マップであり、最適な打ち上げ時期の決定に不可欠
ランベルト問題の解法は、ホーマン遷移のような特殊解を超えて、任意の2点間の軌道遷移を計画するための汎用的なツールです。惑星間航行だけでなく、衛星のランデブー計画、デブリ回収ミッション、マルチフライバイ軌道の設計など、幅広い応用を持ちます。
次のステップとして、以下の記事も参考にしてください。