import numpy as np
import matplotlib.pyplot as plt

# 定数
Re = 6371.0  # 地球半径 [km]

def walker_constellation(T, P, F, h, inc):
    """Walkerデルタパターンのコンステレーションを生成

    Parameters
    ----------
    T : int - 総衛星数
    P : int - 軌道面数
    F : int - 位相パラメータ
    h : float - 軌道高度 [km]
    inc : float - 軌道傾斜角 [deg]

    Returns
    -------
    sats : list of dict - 各衛星の軌道要素
    """
    S = T // P  # 1面あたりの衛星数
    inc_rad = np.radians(inc)
    sats = []

    for p in range(P):
        raan = 360.0 * p / P  # 昇交点赤経 [deg]
        for s in range(S):
            M = 360.0 * s / S + 360.0 * F * p / T  # 平均近点角 [deg]
            M = M % 360.0
            sats.append({
                'plane': p,
                'index': s,
                'raan': raan,
                'M': M,
                'inc': inc,
                'h': h,
            })
    return sats

def sat_subpoint(sat, t=0):
    """衛星の直下点（緯度・経度）を計算（簡易モデル）"""
    h = sat['h']
    a = Re + h
    mu = 398600.4418
    T_orb = 2 * np.pi * np.sqrt(a**3 / mu)
    n = 2 * np.pi / T_orb

    inc_rad = np.radians(sat['inc'])
    raan_rad = np.radians(sat['raan'])
    M_rad = np.radians(sat['M']) + n * t

    # 衛星の位置（軌道面内）
    lat = np.degrees(np.arcsin(np.sin(inc_rad) * np.sin(M_rad)))
    lon_inertial = np.degrees(np.arctan2(
        np.sin(M_rad) * np.cos(inc_rad), np.cos(M_rad)))
    # 地球自転を考慮
    omega_e = 2 * np.pi / 86164.1
    lon = (lon_inertial + np.degrees(raan_rad) - np.degrees(omega_e * t)) % 360
    if lon > 180:
        lon -= 360

    return lat, lon

def coverage_angle(h, eps_min_deg):
    """カバー範囲の半頂角 [rad]"""
    eps_min = np.radians(eps_min_deg)
    theta = np.arccos(Re / (Re + h) * np.cos(eps_min)) - eps_min
    return theta

# Starlinkライクなコンステレーション
T = 72   # 総衛星数（簡易版）
P = 6    # 軌道面数
F = 1    # 位相パラメータ
h = 550  # 軌道高度 [km]
inc = 53 # 軌道傾斜角 [deg]

sats = walker_constellation(T, P, F, h, inc)

# 被覆率の計算
eps_min = 25  # 最低仰角 [deg]
theta_cover = coverage_angle(h, eps_min)
A_cover = 2 * np.pi * Re**2 * (1 - np.cos(theta_cover))
A_earth = 4 * np.pi * Re**2
N_min = 2 / (1 - np.cos(theta_cover))

print(f"=== コンステレーション設計 ===")
print(f"Walker {inc}:{T}/{P}/{F}")
print(f"軌道高度: {h} km")
print(f"最低仰角: {eps_min} deg")
print(f"カバー半頂角: {np.degrees(theta_cover):.1f} deg")
print(f"1衛星のカバー面積: {A_cover:.0f} km^2")
print(f"地球全表面積: {A_earth:.0f} km^2")
print(f"必要最小衛星数（下限）: {N_min:.0f}")

# 可視化
fig, axes = plt.subplots(1, 2, figsize=(16, 6))

# 衛星の配置（t=0）
lats = []
lons = []
colors = []
for sat in sats:
    lat, lon = sat_subpoint(sat, t=0)
    lats.append(lat)
    lons.append(lon)
    colors.append(sat['plane'])

scatter = axes[0].scatter(lons, lats, c=colors, cmap='tab10', s=30, zorder=5)
axes[0].set_xlim(-180, 180)
axes[0].set_ylim(-90, 90)
axes[0].set_xlabel('Longitude [deg]')
axes[0].set_ylabel('Latitude [deg]')
axes[0].set_title(f'Walker {inc}:{T}/{P}/{F} at h={h}km')
axes[0].grid(True, alpha=0.3)

# 各衛星のカバー範囲を円で描画
for lat, lon in zip(lats, lons):
    theta_deg = np.degrees(theta_cover)
    circle_lats = lat + theta_deg * np.cos(np.linspace(0, 2*np.pi, 50))
    circle_lons = lon + theta_deg / np.cos(np.radians(lat)) * np.sin(np.linspace(0, 2*np.pi, 50))
    circle_lats = np.clip(circle_lats, -90, 90)
    axes[0].plot(circle_lons, circle_lats, 'b-', alpha=0.1, linewidth=0.5)

plt.colorbar(scatter, ax=axes[0], label='Orbital Plane')

# 高度と必要衛星数の関係
altitudes = np.linspace(300, 2000, 100)
for eps in [10, 25, 40]:
    N_mins = []
    for alt in altitudes:
        theta = coverage_angle(alt, eps)
        N_min_val = 2 / (1 - np.cos(theta))
        N_mins.append(N_min_val)
    axes[1].semilogy(altitudes, N_mins, label=f'eps_min = {eps} deg')

axes[1].set_xlabel('Altitude [km]')
axes[1].set_ylabel('Minimum Number of Satellites')
axes[1].set_title('Minimum Satellites for Global Coverage')
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.show()