"""宇宙科学の学習用。Python 3の標準機能のみ使用。
実際の飛行計画ではありません。太陽以外の重力・噴射などは省略します。
実行: python3 orbit_practice.py
"""
import math
AU = 149_597_870_700.0  # m
MU = 1.32712440018e20  # m^3/s^2（太陽の重力パラメータ）
DAY = 86400.0

def acceleration(x, y):
    r = math.hypot(x, y)
    return -MU*x/r**3, -MU*y/r**3

def orbit(speed, days, dt=3600.0):
    if not (20000 <= speed <= 50000 and 0 < days <= 1000 and 0 < dt <= DAY):
        raise ValueError('この教材の計算範囲外です')
    x, y, vx, vy = AU, 0., 0., speed
    ax, ay = acceleration(x, y)
    energy_start = speed**2/2 - MU/AU
    time, error = 0., 0.
    while time < days*DAY:
        h = min(dt, days*DAY-time)
        nx = x + vx*h + ax*h*h/2
        ny = y + vy*h + ay*h*h/2
        bx, by = acceleration(nx, ny)
        vx += (ax+bx)*h/2
        vy += (ay+by)*h/2
        x, y, ax, ay = nx, ny, bx, by
        time += h
        energy = (vx*vx+vy*vy)/2-MU/math.hypot(x,y)
        error = max(error, abs(energy-energy_start)/(MU/AU))
    return x/AU, y/AU, error

def main():
    v_circle = math.sqrt(MU/AU)
    r_mars = 1.524*AU
    a = (AU+r_mars)/2
    v_depart = math.sqrt(MU*(2/AU-1/a))
    v_arrive = math.sqrt(MU*(2/r_mars-1/a))
    time = math.pi*math.sqrt(a**3/MU)
    phase = 180-360*time/(2*math.pi*math.sqrt(r_mars**3/MU))
    print('円軌道の速さ:', round(v_circle/1000, 3), 'km/s')
    print('遷移の出発:', round(v_depart/1000, 3), 'km/s')
    print('遷移の到着:', round(v_arrive/1000, 3), 'km/s')
    print('太陽中心の円から遷移への速度差:', round((v_depart-v_circle)/1000, 3), 'km/s')
    print('遷移から外側の円への速度差:', round((math.sqrt(MU/r_mars)-v_arrive)/1000, 3), 'km/s')
    print('遷移時間:', round(time/DAY, 3), '日')
    print('出発時の火星の先行角:', round(phase, 3), '度')
    for step in [3600., 1800.]:
        x,y,error = orbit(v_depart, time/DAY, step)
        miss = math.hypot(x+1.524, y)
        print(f'刻み{step:.0f}秒: 終点({x:.8f}, {y:.8f}) AU, 理論終点との差{miss:.3e} AU, エネルギー検算値{error:.3e}')
    for angle in [45, 135]:
        rad = math.radians(angle)
        x, y = 30+10*math.cos(rad), 10*math.sin(rad)
        print(f'退出角{angle}度: 太陽基準の速さ{math.hypot(x,y):.3f} km/s')
    print('地表での弾道: v0=10m/s、45度、g=9.8m/s²、空気抵抗なし')
    vy = 10/math.sqrt(2)
    print('飛行時間:', round(2*vy/9.8, 3), '秒; 飛距離:', round(vy*2*vy/9.8,3), 'm')

if __name__ == '__main__':
    main()