軌道遷移とデルタV — ホーマン遷移から二重インパルスまで

import numpy as np
import matplotlib.pyplot as plt

def hohmann_transfer(r1, r2, mu):
    """ホーマン遷移のデルタVと遷移時間を計算する"""
    # 遷移楕円の半長軸
    a_t = (r1 + r2) / 2.0

    # 各速度の計算
    v_c1 = np.sqrt(mu / r1)         # 内側円軌道速度
    v_c2 = np.sqrt(mu / r2)         # 外側円軌道速度
    v_tp = np.sqrt(mu * (2/r1 - 1/a_t))  # 遷移楕円の近点速度
    v_ta = np.sqrt(mu * (2/r2 - 1/a_t))  # 遷移楕円の遠点速度

    # デルタV
    dv1 = v_tp - v_c1
    dv2 = v_c2 - v_ta
    dv_total = abs(dv1) + abs(dv2)

    # 遷移時間（半周期）
    t_transfer = np.pi * np.sqrt(a_t**3 / mu)

    return {
        'dv1': dv1, 'dv2': dv2, 'dv_total': dv_total,
        't_transfer': t_transfer, 'a_t': a_t,
        'v_c1': v_c1, 'v_c2': v_c2, 'v_tp': v_tp, 'v_ta': v_ta
    }

# 地球の重力パラメータ
mu_earth = 3.986e14  # m^3/s^2
R_earth = 6378e3     # m

# LEO（高度200 km）からの遷移
r1 = R_earth + 200e3  # LEO半径

# 軌道半径比を変えてデルタVを計算
ratios = np.linspace(1.1, 15, 500)
r2_array = r1 * ratios

dv1_array = []
dv2_array = []
dv_total_array = []
t_array = []

for r2 in r2_array:
    result = hohmann_transfer(r1, r2, mu_earth)
    dv1_array.append(result['dv1'] / 1000)      # km/sに変換
    dv2_array.append(result['dv2'] / 1000)
    dv_total_array.append(result['dv_total'] / 1000)
    t_array.append(result['t_transfer'] / 3600)  # 時間に変換

fig, axes = plt.subplots(2, 1, figsize=(10, 8))

# デルタVのプロット
axes[0].plot(ratios, dv1_array, label='$\\Delta V_1$（第1噴射）',
             linewidth=2, color='#2196F3')
axes[0].plot(ratios, dv2_array, label='$\\Delta V_2$（第2噴射）',
             linewidth=2, color='#4CAF50')
axes[0].plot(ratios, dv_total_array, label='$\\Delta V_{total}$',
             linewidth=2.5, color='#F44336')
axes[0].set_xlabel('軌道半径比 $r_2/r_1$', fontsize=12)
axes[0].set_ylabel('$\\Delta V$ [km/s]', fontsize=12)
axes[0].set_title('ホーマン遷移のデルタV（LEO起点）', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(1, 15)

# 遷移時間のプロット
axes[1].plot(ratios, t_array, linewidth=2, color='#9C27B0')
axes[1].set_xlabel('軌道半径比 $r_2/r_1$', fontsize=12)
axes[1].set_ylabel('遷移時間 [時間]', fontsize=12)
axes[1].set_title('ホーマン遷移の所要時間', fontsize=14)
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(1, 15)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

mu = 3.986e14  # m^3/s^2
R_earth = 6378e3  # m

# LEO（高度200 km）→ GEO（高度35786 km）
r_leo = R_earth + 200e3
r_geo = R_earth + 35786e3

# ホーマン遷移の計算
result = hohmann_transfer(r_leo, r_geo, mu)
a_t = result['a_t']

print("=== LEO → GEO ホーマン遷移 ===")
print(f"LEO半径: {r_leo/1000:.0f} km")
print(f"GEO半径: {r_geo/1000:.0f} km")
print(f"遷移楕円半長軸: {a_t/1000:.0f} km")
print(f"ΔV1: {result['dv1']/1000:.3f} km/s")
print(f"ΔV2: {result['dv2']/1000:.3f} km/s")
print(f"ΔV合計: {result['dv_total']/1000:.3f} km/s")
print(f"遷移時間: {result['t_transfer']/3600:.2f} 時間")

# 可視化
theta = np.linspace(0, 2 * np.pi, 500)

# 地球
theta_earth = np.linspace(0, 2 * np.pi, 100)
x_earth = R_earth * np.cos(theta_earth) / 1000
y_earth = R_earth * np.sin(theta_earth) / 1000

# LEO
x_leo = r_leo * np.cos(theta) / 1000
y_leo = r_leo * np.sin(theta) / 1000

# GEO
x_geo = r_geo * np.cos(theta) / 1000
y_geo = r_geo * np.sin(theta) / 1000

# 遷移楕円（半分だけ描画）
e_t = 1 - r_leo / a_t
p_t = a_t * (1 - e_t**2)
theta_transfer = np.linspace(0, np.pi, 300)
r_transfer = p_t / (1 + e_t * np.cos(theta_transfer))
x_transfer = r_transfer * np.cos(theta_transfer) / 1000
y_transfer = r_transfer * np.sin(theta_transfer) / 1000

fig, ax = plt.subplots(figsize=(10, 10))

# 地球
earth_circle = plt.Circle((0, 0), R_earth / 1000,
                           color='#4FC3F7', alpha=0.4, label='地球')
ax.add_patch(earth_circle)

# 軌道
ax.plot(x_leo, y_leo, '--', color='#2196F3', linewidth=1.5,
        label=f'LEO (h=200 km)')
ax.plot(x_geo, y_geo, '--', color='#4CAF50', linewidth=1.5,
        label=f'GEO (h=35,786 km)')
ax.plot(x_transfer, y_transfer, '-', color='#F44336', linewidth=2.5,
        label='遷移楕円')

# 噴射点をマーク
ax.plot(r_leo / 1000, 0, 'v', color='#F44336', markersize=15,
        markeredgecolor='black', markeredgewidth=1)
ax.annotate(f'第1噴射\n$\\Delta V_1$ = {result["dv1"]/1000:.3f} km/s',
            xy=(r_leo / 1000, 0), fontsize=10,
            textcoords="offset points", xytext=(15, -30))

ax.plot(-r_geo / 1000, 0, '^', color='#F44336', markersize=15,
        markeredgecolor='black', markeredgewidth=1)
ax.annotate(f'第2噴射\n$\\Delta V_2$ = {result["dv2"]/1000:.3f} km/s',
            xy=(-r_geo / 1000, 0), fontsize=10,
            textcoords="offset points", xytext=(15, 10))

ax.set_xlim(-55000, 55000)
ax.set_ylim(-55000, 55000)
ax.set_xlabel('x [km]', fontsize=12)
ax.set_ylabel('y [km]', fontsize=12)
ax.set_title('LEO → GEO ホーマン遷移軌道', fontsize=14)
ax.legend(fontsize=11, loc='upper right')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def bielliptic_transfer(r1, r2, rb, mu):
    """バイエリプティック遷移のデルタVを計算"""
    a1 = (r1 + rb) / 2
    a2 = (r2 + rb) / 2

    v_c1 = np.sqrt(mu / r1)
    v_c2 = np.sqrt(mu / r2)

    # 楕円1の近点速度
    v1_p = np.sqrt(mu * (2/r1 - 1/a1))
    # 楕円1の遠点速度
    v1_a = np.sqrt(mu * (2/rb - 1/a1))
    # 楕円2の遠点速度
    v2_a = np.sqrt(mu * (2/rb - 1/a2))
    # 楕円2の近点速度
    v2_p = np.sqrt(mu * (2/r2 - 1/a2))

    dv1 = abs(v1_p - v_c1)
    dv2 = abs(v2_a - v1_a)
    dv3 = abs(v_c2 - v2_p)

    return dv1 + dv2 + dv3

mu = 1.0
r1 = 1.0

ratios = np.linspace(2, 30, 500)
dv_hohmann = []
dv_bielliptic = []

for R in ratios:
    r2 = r1 * R

    # ホーマン遷移
    res = hohmann_transfer(r1, r2, mu)
    dv_hohmann.append(res['dv_total'])

    # バイエリプティック（中間半径 rb を最適化）
    rb_values = np.linspace(r2 * 1.1, r2 * 50, 200)
    dv_bi_min = float('inf')
    for rb in rb_values:
        dv_bi = bielliptic_transfer(r1, r2, rb, mu)
        if dv_bi < dv_bi_min:
            dv_bi_min = dv_bi
    dv_bielliptic.append(dv_bi_min)

# 円軌道速度で正規化
v_c1 = np.sqrt(mu / r1)
dv_hohmann_norm = np.array(dv_hohmann) / v_c1
dv_bielliptic_norm = np.array(dv_bielliptic) / v_c1

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(ratios, dv_hohmann_norm, linewidth=2.5, color='#2196F3',
        label='ホーマン遷移（2インパルス）')
ax.plot(ratios, dv_bielliptic_norm, linewidth=2.5, color='#F44336',
        label='バイエリプティック遷移（3インパルス, 最適rb）')

# 交差点の近似表示
ax.axvline(x=11.94, color='gray', linestyle=':', linewidth=1.5,
           label='$r_2/r_1 = 11.94$（交差点）')

ax.set_xlabel('軌道半径比 $r_2/r_1$', fontsize=12)
ax.set_ylabel('$\\Delta V / v_{c1}$ [正規化]', fontsize=12)
ax.set_title('ホーマン遷移 vs バイエリプティック遷移のデルタV比較',
             fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_xlim(2, 30)
plt.tight_layout()
plt.show()

import numpy as np

mu = 3.986e14  # m^3/s^2
R_earth = 6378e3  # m

# 軌道パラメータ
r_leo = R_earth + 200e3   # LEO高度200 km
r_geo = R_earth + 35786e3  # GEO高度35786 km
i_launch = 28.5  # 打ち上げ緯度 [deg]
delta_i = np.radians(i_launch)  # 傾斜角変更量 [rad]

# === 方法1: ホーマン遷移 + LEOで傾斜角変更 ===
res = hohmann_transfer(r_leo, r_geo, mu)
v_leo = np.sqrt(mu / r_leo)
dv_plane_leo = 2 * v_leo * np.sin(delta_i / 2)
dv_method1 = res['dv_total'] + dv_plane_leo

# === 方法2: ホーマン遷移 + GEOで傾斜角変更 ===
v_geo = np.sqrt(mu / r_geo)
dv_plane_geo = 2 * v_geo * np.sin(delta_i / 2)
dv_method2 = res['dv_total'] + dv_plane_geo

# === 方法3: ホーマン遷移 + GTO遠点で複合マヌーバ ===
# 第1噴射はLEOで通常通り
dv1 = res['dv1']

# 第2噴射（GTO遠点）で高度変更と傾斜角変更を同時実行
v_ta = res['v_ta']  # GTO遠点速度
v_c2 = res['v_c2']  # GEO円軌道速度
# 複合マヌーバのデルタV（余弦定理）
dv2_combined = np.sqrt(v_ta**2 + v_c2**2 - 2 * v_ta * v_c2 * np.cos(delta_i))
dv_method3 = abs(dv1) + dv2_combined

print("=== LEO(28.5°) → GEO(0°) 遷移比較 ===")
print(f"\n方法1: ホーマン遷移 + LEOで傾斜角変更")
print(f"  ΔV(Hohmann) = {res['dv_total']/1000:.3f} km/s")
print(f"  ΔV(plane@LEO) = {dv_plane_leo/1000:.3f} km/s")
print(f"  ΔV(合計) = {dv_method1/1000:.3f} km/s")

print(f"\n方法2: ホーマン遷移 + GEOで傾斜角変更")
print(f"  ΔV(Hohmann) = {res['dv_total']/1000:.3f} km/s")
print(f"  ΔV(plane@GEO) = {dv_plane_geo/1000:.3f} km/s")
print(f"  ΔV(合計) = {dv_method2/1000:.3f} km/s")

print(f"\n方法3: GTO遠点で複合マヌーバ")
print(f"  ΔV1(LEO→GTO) = {abs(dv1)/1000:.3f} km/s")
print(f"  ΔV2(複合) = {dv2_combined/1000:.3f} km/s")
print(f"  ΔV(合計) = {dv_method3/1000:.3f} km/s")

print(f"\n方法3は方法1より {(dv_method1-dv_method3)/1000:.3f} km/s 節約")
print(f"方法3は方法2より {(dv_method2-dv_method3)/1000:.3f} km/s 節約")