衛星通信の基礎 — GEOからLEOコンステレーションまで

import numpy as np
import matplotlib.pyplot as plt

# 物理定数
Re = 6371        # 地球半径 [km]
c = 3e5          # 光速 [km/s]
G = 6.674e-11    # 万有引力定数 [m^3/(kg*s^2)]
Me = 5.972e24    # 地球質量 [kg]
mu = G * Me      # 地球の標準重力パラメータ [m^3/s^2]

def orbital_velocity(h_km):
    """軌道速度 [km/s]"""
    r = (Re + h_km) * 1000  # m
    v = np.sqrt(mu / r) / 1000  # km/s
    return v

def orbital_period(h_km):
    """軌道周期 [分]"""
    r = (Re + h_km) * 1000  # m
    T = 2 * np.pi * np.sqrt(r**3 / mu) / 60  # 分
    return T

def one_way_delay(h_km, elev_deg=90):
    """片道伝搬遅延 [ms]（最短は天頂方向）"""
    elev = np.radians(elev_deg)
    # 衛星までのスラント距離
    r_sat = Re + h_km
    slant = -Re * np.sin(elev) + np.sqrt((Re * np.sin(elev))**2 +
             2 * Re * h_km + h_km**2)
    delay = slant / c * 1000  # ms
    return delay

def fspl_sat(h_km, f_ghz, elev_deg=90):
    """自由空間損失 [dB]"""
    elev = np.radians(elev_deg)
    r_sat = Re + h_km
    slant = -Re * np.sin(elev) + np.sqrt((Re * np.sin(elev))**2 +
             2 * Re * h_km + h_km**2)
    fspl = 20 * np.log10(slant * 1000) + 20 * np.log10(f_ghz * 1e9) + \
           20 * np.log10(4 * np.pi / (3e8))
    return fspl

def max_doppler(h_km, f_ghz):
    """最大ドップラーシフト [kHz]（低仰角通過時）"""
    v = orbital_velocity(h_km) * 1000  # m/s
    # 最大ドップラーは衛星が真横を通過するとき（放射方向速度が最大）
    f_doppler = v * f_ghz * 1e9 / (3e8) / 1000  # kHz
    return f_doppler

def visibility_time(h_km, min_elev_deg=10):
    """1回の可視時間 [分]（最低仰角以上の時間）"""
    elev = np.radians(min_elev_deg)
    # 可視半角
    r_sat = Re + h_km
    half_angle = np.arccos(Re * np.cos(elev) / r_sat) - elev
    # 可視弧の長さ
    T = orbital_period(h_km)
    vis = T * (2 * half_angle) / (2 * np.pi)
    return vis

# 高度範囲
h = np.linspace(200, 40000, 1000)

fig, axes = plt.subplots(2, 3, figsize=(16, 10))

# 1. 往復遅延
delay_90 = 2 * one_way_delay(h, 90)  # 天頂
delay_10 = 2 * one_way_delay(h, 10)  # 仰角10度
axes[0, 0].semilogy(h, delay_90, 'b-', linewidth=2, label='Elevation = 90°')
axes[0, 0].semilogy(h, delay_10, 'r--', linewidth=2, label='Elevation = 10°')
axes[0, 0].axhline(y=250, color='gray', linestyle=':', alpha=0.5)
axes[0, 0].text(5000, 280, 'GEO RTT ≈ 250 ms', fontsize=9, color='gray')
axes[0, 0].set_xlabel('Altitude [km]', fontsize=11)
axes[0, 0].set_ylabel('Round-Trip Delay [ms]', fontsize=11)
axes[0, 0].set_title('Propagation Delay', fontsize=12)
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)

# 2. 自由空間損失（Ku帯 12 GHz）
fspl_90 = fspl_sat(h, 12, 90)
fspl_10 = fspl_sat(h, 12, 10)
axes[0, 1].plot(h, fspl_90, 'b-', linewidth=2, label='Elevation = 90°')
axes[0, 1].plot(h, fspl_10, 'r--', linewidth=2, label='Elevation = 10°')
axes[0, 1].set_xlabel('Altitude [km]', fontsize=11)
axes[0, 1].set_ylabel('FSPL [dB]', fontsize=11)
axes[0, 1].set_title('Free Space Path Loss (12 GHz)', fontsize=12)
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)

# 3. ドップラーシフト
doppler_ku = max_doppler(h, 12)
doppler_ka = max_doppler(h, 30)
axes[0, 2].semilogy(h, doppler_ku, 'b-', linewidth=2, label='Ku (12 GHz)')
axes[0, 2].semilogy(h, doppler_ka, 'r-', linewidth=2, label='Ka (30 GHz)')
axes[0, 2].set_xlabel('Altitude [km]', fontsize=11)
axes[0, 2].set_ylabel('Max Doppler Shift [kHz]', fontsize=11)
axes[0, 2].set_title('Maximum Doppler Shift', fontsize=12)
axes[0, 2].legend(fontsize=9)
axes[0, 2].grid(True, alpha=0.3)

# 4. 可視時間
vis_10 = visibility_time(h, 10)
vis_20 = visibility_time(h, 20)
axes[1, 0].plot(h, vis_10, 'b-', linewidth=2, label='Min elev = 10°')
axes[1, 0].plot(h, vis_20, 'r--', linewidth=2, label='Min elev = 20°')
axes[1, 0].set_xlabel('Altitude [km]', fontsize=11)
axes[1, 0].set_ylabel('Visibility Time [min]', fontsize=11)
axes[1, 0].set_title('Single Pass Visibility Time', fontsize=12)
axes[1, 0].legend(fontsize=9)
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].set_xlim([200, 5000])

# 5. 軌道速度と周期
ax5a = axes[1, 1]
ax5b = ax5a.twinx()
ax5a.plot(h, orbital_velocity(h), 'b-', linewidth=2)
ax5b.plot(h, orbital_period(h) / 60, 'r-', linewidth=2)
ax5a.set_xlabel('Altitude [km]', fontsize=11)
ax5a.set_ylabel('Orbital Velocity [km/s]', fontsize=11, color='b')
ax5b.set_ylabel('Orbital Period [hours]', fontsize=11, color='r')
ax5a.set_title('Orbital Velocity & Period', fontsize=12)
ax5a.grid(True, alpha=0.3)

# 6. 代表的な衛星の位置をプロット
sat_systems = {
    'Starlink': 550,
    'OneWeb': 1200,
    'O3b': 8062,
    'GPS': 20200,
    'GEO': 35786,
}

categories = ['Delay\n[ms]', 'FSPL\n[dB]', 'Doppler\n[kHz]',
               'Visibility\n[min]', 'Velocity\n[km/s]']
values = {}
for name, alt in sat_systems.items():
    values[name] = [
        2 * one_way_delay(alt, 30),
        fspl_sat(alt, 12, 30),
        max_doppler(alt, 12),
        visibility_time(alt, 10) if alt < 36000 else np.inf,
        orbital_velocity(alt)
    ]

# テーブル形式で比較
cell_text = []
for name, alt in sat_systems.items():
    v = values[name]
    vis_str = f'{v[3]:.1f}' if v[3] < 1000 else 'N/A'
    cell_text.append([name, f'{alt}', f'{v[0]:.1f}', f'{v[1]:.1f}',
                       f'{v[2]:.1f}', vis_str, f'{v[4]:.2f}'])

axes[1, 2].axis('off')
table = axes[1, 2].table(
    cellText=cell_text,
    colLabels=['System', 'Alt [km]', 'RTD [ms]', 'FSPL [dB]',
                'Doppler\n[kHz]', 'Vis [min]', 'Vel [km/s]'],
    loc='center',
    cellLoc='center'
)
table.auto_set_font_size(False)
table.set_fontsize(9)
table.scale(1.2, 1.8)
axes[1, 2].set_title('Satellite System Comparison (elev=30°)', fontsize=12)

plt.tight_layout()
plt.savefig('satellite_orbit_parameters.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

Re = 6371  # km

def coverage_radius(h_km, min_elev_deg):
    """衛星のカバレッジ半径 [km]（地表面上の距離）"""
    elev = np.radians(min_elev_deg)
    r_sat = Re + h_km

    # 中心角
    central_angle = np.arccos(Re * np.cos(elev) / r_sat) - elev

    # 地表面上の距離
    ground_dist = Re * central_angle
    return ground_dist

def coverage_area(h_km, min_elev_deg):
    """カバレッジ面積 [km^2]"""
    central_angle = coverage_radius(h_km, min_elev_deg) / Re
    area = 2 * np.pi * Re**2 * (1 - np.cos(central_angle))
    return area

def min_satellites_for_coverage(h_km, min_elev_deg):
    """全球カバーに必要な最小衛星数（概算）"""
    earth_surface = 4 * np.pi * Re**2
    single_area = coverage_area(h_km, min_elev_deg)
    # 球面を重なりなくカバーするのは不可能なので、重複係数1.5を掛ける
    n_min = earth_surface / single_area * 1.5
    return n_min

# 高度範囲
h = np.linspace(200, 40000, 500)

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# カバレッジ半径
for elev in [5, 10, 20, 30]:
    r = coverage_radius(h, elev)
    axes[0].plot(h, r, linewidth=2, label=f'Min elev = {elev}°')

axes[0].set_xlabel('Altitude [km]', fontsize=12)
axes[0].set_ylabel('Coverage Radius [km]', fontsize=12)
axes[0].set_title('Ground Coverage Radius', fontsize=13)
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)

# カバレッジ面積（地球表面積の割合で）
earth_area = 4 * np.pi * Re**2  # km^2

for elev in [5, 10, 20]:
    area_frac = coverage_area(h, elev) / earth_area * 100
    axes[1].semilogy(h, area_frac, linewidth=2, label=f'Min elev = {elev}°')

axes[1].set_xlabel('Altitude [km]', fontsize=12)
axes[1].set_ylabel('Coverage Area [% of Earth]', fontsize=12)
axes[1].set_title('Coverage Area Fraction', fontsize=13)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)

# 必要衛星数
h_leo = np.linspace(300, 2000, 200)
for elev in [10, 20, 30]:
    n_sat = min_satellites_for_coverage(h_leo, elev)
    axes[2].plot(h_leo, n_sat, linewidth=2, label=f'Min elev = {elev}°')

axes[2].set_xlabel('Altitude [km]', fontsize=12)
axes[2].set_ylabel('Minimum Satellites for Global Coverage', fontsize=12)
axes[2].set_title('Required Satellites (LEO range)', fontsize=13)
axes[2].legend(fontsize=10)
axes[2].grid(True, alpha=0.3)
axes[2].set_ylim([0, 3000])

# Starlinkの位置をマーク
axes[2].axvline(x=550, color='gray', linestyle=':', alpha=0.5)
axes[2].annotate('Starlink\n(550 km)', xy=(550, 1200),
                  fontsize=10, ha='center')

plt.tight_layout()
plt.savefig('satellite_coverage.png', dpi=150, bbox_inches='tight')
plt.show()

# 具体的な数値
for name, alt in [('Starlink', 550), ('OneWeb', 1200), ('GEO', 35786)]:
    r = coverage_radius(alt, 10)
    area = coverage_area(alt, 10)
    frac = area / earth_area * 100
    print(f"{name} ({alt} km): カバレッジ半径 {r:.0f} km, "
          f"面積 {area:.0f} km², 地球の {frac:.2f}%")

import numpy as np
import matplotlib.pyplot as plt

# 各衛星システムの代表的パラメータ
systems = {
    'L-band\n(Iridium/GPS)': {'f': 1.5, 'alt': 780, 'color': '#4CAF50'},
    'S-band': {'f': 2.5, 'alt': 1200, 'color': '#8BC34A'},
    'C-band\n(VSAT)': {'f': 6, 'alt': 35786, 'color': '#2196F3'},
    'Ku-band\n(BSat)': {'f': 12, 'alt': 35786, 'color': '#FF9800'},
    'Ka-band\n(Starlink)': {'f': 20, 'alt': 550, 'color': '#F44336'},
    'Ka-band\n(GEO HTS)': {'f': 20, 'alt': 35786, 'color': '#9C27B0'},
}

Re = 6371

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 仰角による自由空間損失の変化（各システム）
elevations = np.linspace(5, 90, 100)

for name, params in systems.items():
    f = params['f']
    h = params['alt']
    fspl_list = []
    for elev in elevations:
        elev_rad = np.radians(elev)
        slant = -Re * np.sin(elev_rad) + np.sqrt(
            (Re * np.sin(elev_rad))**2 + 2 * Re * h + h**2)
        fspl_val = 20 * np.log10(slant * 1000) + 20 * np.log10(f * 1e9) + \
                   20 * np.log10(4 * np.pi / 3e8)
        fspl_list.append(fspl_val)

    axes[0].plot(elevations, fspl_list, linewidth=2,
                  color=params['color'], label=name)

axes[0].set_xlabel('Elevation Angle [deg]', fontsize=12)
axes[0].set_ylabel('Free Space Path Loss [dB]', fontsize=12)
axes[0].set_title('FSPL vs Elevation for Various Systems', fontsize=13)
axes[0].legend(fontsize=8, ncol=2)
axes[0].grid(True, alpha=0.3)
axes[0].invert_yaxis()

# 周波数帯ごとの降雨マージン必要量（概算）
# ITU-R P.618に基づく概算（GEO、仰角30度、東京、0.01%）
freq_bands = ['L\n(1.5 GHz)', 'S\n(2.5 GHz)', 'C\n(6 GHz)',
               'Ku\n(12 GHz)', 'Ka-up\n(30 GHz)', 'Ka-dn\n(20 GHz)',
               'V\n(50 GHz)']
rain_atten = [0.01, 0.05, 0.3, 3.5, 18, 8, 35]  # dB概算
atm_atten = [0.05, 0.05, 0.08, 0.15, 0.4, 0.3, 2.0]  # dB概算

x = np.arange(len(freq_bands))
width = 0.35

bars1 = axes[1].bar(x - width/2, rain_atten, width, label='Rain (0.01%)',
                      color='#2196F3', alpha=0.8)
bars2 = axes[1].bar(x + width/2, atm_atten, width, label='Atmospheric',
                      color='#FF9800', alpha=0.8)

axes[1].set_xlabel('Frequency Band', fontsize=12)
axes[1].set_ylabel('Attenuation [dB]', fontsize=12)
axes[1].set_title('Rain & Atmospheric Attenuation\n'
                    '(GEO, elev=30°, Tokyo, 0.01% exceedance)', fontsize=13)
axes[1].set_xticks(x)
axes[1].set_xticklabels(freq_bands, fontsize=10)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3, axis='y')
axes[1].set_yscale('log')
axes[1].set_ylim([0.005, 50])

plt.tight_layout()
plt.savefig('satellite_frequency_comparison.png', dpi=150, bbox_inches='tight')
plt.show()