"""Verification of every number and identity quoted in the orbital-manoeuvres note (orbital-maneuvers.qmd).

Run:  python3 codes/orbital-maneuvers-verify.py        (a few seconds)
"""
import numpy as np
import sympy as sp
from scipy.optimize import brentq, minimize_scalar

ok = lambda label, cond: print(f"  [{'ok' if cond else 'FAIL'}] {label}")
close = lambda x, y, tol: abs(x - y) <= tol
d2r = np.pi/180

MU = 398600.4418
RE = 6378.137
r1, r2 = 6678.137, 42164.0
g0 = 9.80665e-3
vis = lambda r, a: np.sqrt(MU*(2/r - 1/a))
vc = lambda r: np.sqrt(MU/r)


def hohmann(ra_, rb_):
    at = (ra_ + rb_)/2
    return abs(vis(ra_, at) - vc(ra_)), abs(vc(rb_) - vis(rb_, at)), np.pi*np.sqrt(at**3/MU)


# ------------------------------------------------------------------ Part 0
print("Part 0: impulsive model")
dvH1, dvH2, tH = hohmann(r1, r2)
dvH = dvH1 + dvH2
ve = 315*g0
frac = 1 - np.exp(-dvH/ve)
ok(f"rocket: v_e = {ve:.4f} km/s, propellant fraction {frac:.3f} (0.716)", close(ve, 3.0891, 1e-4) and close(frac, 0.716, 0.0005))
v, dv, beta = sp.symbols("v dv beta", positive=True)
ok("d(eps) = v dv cos(beta) + dv^2/2",
   sp.simplify(((v + dv*sp.cos(beta))**2 + (dv*sp.sin(beta))**2)/2 - v**2/2 - (v*dv*sp.cos(beta) + dv**2/2)) == 0)
a_o, e_o = 12000.0, 0.5
vp_o, va_o = vis(a_o*(1 - e_o), a_o), vis(a_o*(1 + e_o), a_o)
ok(f"Oberth figure: v_p/v_a = {vp_o/va_o:.3f} (3)", close(vp_o/va_o, 3, 1e-9))
ok(f"  100 m/s buys {vp_o*0.1 + 0.005:.3f} vs {va_o*0.1 + 0.005:.3f} km^2/s^2 (ratio {(vp_o*0.1 + 0.005)/(va_o*0.1 + 0.005):.2f})",
   close((vp_o*0.1 + 0.005)/(va_o*0.1 + 0.005), 2.98, 0.02))

# ------------------------------------------------------------------ Part I: Hohmann
print("Part I: Hohmann")
at = (r1 + r2)/2
ok(f"LEO-GEO: a_t = {at:.2f}, e_t = {(r2 - r1)/(r2 + r1):.5f}", close(at, 24421.07, 0.01) and close((r2 - r1)/(r2 + r1), 0.72654, 1e-5))
ok(f"  v_c1 {vc(r1):.5f}, v_t1 {vis(r1, at):.5f}, v_t2 {vis(r2, at):.5f}, v_c2 {vc(r2):.5f}",
   close(vc(r1), 7.72576, 1e-5) and close(vis(r1, at), 10.15149, 1e-5) and close(vis(r2, at), 1.60784, 1e-5) and close(vc(r2), 3.07467, 1e-5))
ok(f"  dv1 {dvH1:.5f}, dv2 {dvH2:.5f}, total {dvH:.5f}, t_H {tH:.0f} s = {tH/3600:.3f} h",
   close(dvH1, 2.42573, 1e-5) and close(dvH2, 1.46682, 1e-5) and close(dvH, 3.89255, 1e-5) and close(tH, 18990, 1) and close(tH/3600, 5.275, 0.001))
g = lambda R: np.sqrt(2*R/(1 + R)) - 1 + (1 - np.sqrt(2/(1 + R)))/np.sqrt(R)
res = minimize_scalar(lambda R: -g(R), bounds=(2, 50), method="bounded", options={"xatol": 1e-10})
ok(f"Hohmann maximum at R = {res.x:.6f}, g = {g(res.x):.6f}", close(res.x, 15.581719, 1e-5) and close(g(res.x), 0.536258, 1e-6))
ok(f"R(LEO-GEO) = {r2/r1:.3f}", close(r2/r1, 6.31, 0.005))


def apsidal(rA, rB, a1, a2):
    at_ = (rA + rB)/2
    return abs(vis(rA, at_) - vis(rA, a1)), abs(vis(rB, a2) - vis(rB, at_)), np.pi*np.sqrt(at_**3/MU)/3600, at_


a1, e1, a2, e2 = 12000.0, 0.30, 42000.0, 0.50
rp1, ra1, rp2, ra2 = a1*(1 - e1), a1*(1 + e1), a2*(1 - e2), a2*(1 + e2)
A1 = apsidal(rp1, ra2, a1, a2)
A2 = apsidal(ra1, rp2, a1, a2)
print(f"     A1 {A1}, A2 {A2}")
ok("Example 2 A1: 1.29675 + 0.55850 = 1.85525, 9.324 h",
   close(A1[0], 1.29675, 1e-5) and close(A1[1], 0.55850, 1e-5) and close(A1[0] + A1[1], 1.85525, 1e-5) and close(A1[2], 9.324, 0.001))
ok("Example 2 A2: 1.18573 + 1.31337 = 2.49909, 3.422 h",
   close(A2[0], 1.18573, 1e-5) and close(A2[1], 1.31337, 1e-5) and close(A2[0] + A2[1], 2.49909, 1e-5) and close(A2[2], 3.422, 0.001))
ok(f"  e_t: {(ra2 - rp1)/(ra2 + rp1):.4f}, {(rp2 - ra1)/(rp2 + ra1):.4f}", close((ra2 - rp1)/(ra2 + rp1), 0.7647, 1e-4) and close((rp2 - ra1)/(rp2 + ra1), 0.1475, 1e-4))
ok(f"  A1 cheaper by {1000*((A2[0] + A2[1]) - (A1[0] + A1[1])):.0f} m/s, A2 faster x{A1[2]/A2[2]:.2f}, A2 costs {100*((A2[0]+A2[1])/(A1[0]+A1[1]) - 1):.1f}% more",
   close((A2[0] + A2[1]) - (A1[0] + A1[1]), 0.644, 0.001) and close(A1[2]/A2[2], 2.72, 0.005) and close(100*((A2[0]+A2[1])/(A1[0]+A1[1]) - 1), 34.7, 0.05))
T2 = apsidal(rp1, rp2, a1, a2)   # anti-aligned: arrival at r_p2 (speed of orbit 2 at its periapsis)
T1 = apsidal(ra1, ra2, a1, a2)
print(f"     T2 {T2}, T1 {T1}")
ok("Example 3 T2: 0.37924 + 2.04250 = 2.42174, 2.464 h; T1: 2.17083 + 0.19386 = 2.36469, 10.769 h",
   close(T2[0], 0.37924, 1e-5) and close(T2[1], 2.04250, 1e-5) and close(T2[2], 2.464, 0.001) and
   close(T1[0], 2.17083, 1e-5) and close(T1[1], 0.19386, 1e-5) and close(T1[2], 10.769, 0.001))
ok(f"  T1 cheaper by {1000*((T2[0]+T2[1]) - (T1[0]+T1[1])):.0f} m/s, T2 faster x{T1[2]/T2[2]:.2f}, T2 +{100*((T2[0]+T2[1])/(T1[0]+T1[1]) - 1):.1f}%",
   close((T2[0]+T2[1]) - (T1[0]+T1[1]), 0.057, 0.001) and close(T1[2]/T2[2], 4.37, 0.005) and close(100*((T2[0]+T2[1])/(T1[0]+T1[1]) - 1), 2.4, 0.05))
ok("orientation penalty 2.36469/1.85525 = +27%", close(T1[0] + T1[1], 2.365, 0.001) and close((T1[0]+T1[1])/(A1[0]+A1[1]) - 1, 0.27, 0.005))
ok("c1 + c2 = 24600 = 39300 - 14700; c2 - c1 = 17400 = 35700 - 18300",
   close(a1*e1 + a2*e2, 24600, 1e-9) and close(T1[3] - T2[3], 24600, 1e-9) and close(a2*e2 - a1*e1, A1[3] - A2[3], 1e-9))
# e1 = 0.45 at comparable radii reverses the anti-aligned ordering
a1b = 12000.0
T2b = apsidal(a1b*(1 - 0.45), rp2, a1b, a2)
T1b = apsidal(a1b*(1 + 0.45), ra2, a1b, a2)
ok(f"e1 = 0.45: T2 {T2b[0]+T2b[1]:.4f} vs T1 {T1b[0]+T1b[1]:.4f} -> T2 cheaper and faster",
   (T2b[0] + T2b[1]) < (T1b[0] + T1b[1]) and T2b[2] < T1b[2])
# A12 theorem, random test
rng = np.random.default_rng(1)
bad = 0
for _ in range(20000):
    rp_1 = rng.uniform(6600, 40000); ra_1 = rp_1*rng.uniform(1, 6)
    rp_2 = rp_1*rng.uniform(1.0001, 5); ra_2 = max(ra_1*rng.uniform(1.0001, 5), rp_2)
    aa1, aa2 = (rp_1 + ra_1)/2, (rp_2 + ra_2)/2
    x = apsidal(rp_1, ra_2, aa1, aa2); y = apsidal(ra_1, rp_2, aa1, aa2)
    if x[0] + x[1] > y[0] + y[1] + 1e-12:
        bad += 1
ok(f"A12: aligned A1 never more expensive in 20000 random cases (violations: {bad})", bad == 0)

# Example 4: crossing ellipses
rpA, raA, rpB, raB = 8000.0, 30000.0, 12000.0, 24000.0
aA, eA = (rpA + raA)/2, (raA - rpA)/(raA + rpA)
aB, eB = (rpB + raB)/2, (raB - rpB)/(raB + rpB)
pA, pB = aA*(1 - eA**2), aB*(1 - eB**2)
cf = (pB - pA)/(pA*eB - pB*eA)
fs = np.arccos(cf)
rs = pA/(1 + eA*cf)
ok(f"Example 4: a1 {aA}, e1 {eA:.6f}, p1 {pA:.2f}; a2 {aB}, e2 {eB:.6f}, p2 {pB:.0f}",
   close(eA, 0.578947, 1e-6) and close(pA, 12631.58, 0.01) and close(pB, 16000, 1e-6))
ok(f"  cos f* = {cf:.6f}, f* = {np.degrees(fs):.4f}, r* = {rs:.2f}", close(cf, -2/3, 1e-6) and close(np.degrees(fs), 131.8103, 1e-4) and close(rs, 20571.43, 0.01))
hA, hB = np.sqrt(MU*pA), np.sqrt(MU*pB)
vrA, vtA = MU/hA*eA*np.sin(fs), hA/rs
vrB, vtB = MU/hB*eB*np.sin(fs), hB/rs
dvs = np.hypot(vrB - vrA, vtB - vtA)
print(f"     orbit1 vr {vrA:.5f} vt {vtA:.5f} v {np.hypot(vrA, vtA):.5f} g {np.degrees(np.arctan2(vrA, vtA)):.3f}; orbit2 vr {vrB:.5f} vt {vtB:.5f} v {np.hypot(vrB, vtB):.5f} g {np.degrees(np.arctan2(vrB, vtB)):.3f}")
ok(f"  single impulse {dvs:.5f} km/s (1.26058)", close(dvs, 1.26058, 1e-5) and close(vrA, 2.42406, 1e-5) and close(vtB, 3.88208, 1e-5))
X1 = apsidal(rpA, raB, aA, aB)
X2 = apsidal(raA, rpB, aA, aB)
print(f"     two-impulse: {X1}, {X2}")
ok("  two-impulse routes 0.67039 (2.797 h) and 0.62376 (4.206 h); ratio 2.02",
   close(X1[0] + X1[1], 0.67039, 1e-5) and close(X1[2], 2.797, 0.001) and close(X2[0] + X2[1], 0.62376, 1e-5)
   and close(X2[2], 4.206, 0.001) and close(dvs/(X2[0] + X2[1]), 2.02, 0.005))
gA, gB = np.arctan2(vrA, vtA), np.arctan2(vrB, vtB)
vbar = (np.hypot(vrA, vtA) + np.hypot(vrB, vtB))/2
ok(f"  rotation estimate 2 vbar sin(dgamma/2) = {2*vbar*np.sin((gA - gB)/2):.3f} ({100*2*vbar*np.sin((gA-gB)/2)/dvs:.0f}%)",
   close(2*vbar*np.sin((gA - gB)/2), 1.253, 0.001) and close(np.degrees(gA - gB), 17.38, 0.01))

# Hohmann optimality: eq-dv-eh and a brute-force check
eps = sp.symbols("epsilon", real=True)
h, r = sp.symbols("h r", positive=True)
mus = sp.symbols("mu", positive=True)
vr2_ = 2*(eps + mus/r) - h**2/r**2
expr = (vr2_ + (h/r - sp.sqrt(mus/r))**2) - (2*eps + 3*mus/r - 2*h*sp.sqrt(mus)/r**sp.Rational(3, 2))
ok("eq-dv-eh is exact", sp.simplify(expr) == 0)
def cost_h(hh):
    F = lambda rr: hh**2/(2*rr**2) - MU/rr
    ee = max(F(r1), F(r2))
    return sum(np.sqrt(max(2*ee + 3*MU/rr - 2*hh*np.sqrt(MU)/rr**1.5, 0)) for rr in (r1, r2))
grid = np.linspace(0.8*np.sqrt(MU*r1), 1.2*np.sqrt(MU*r2), 4001)
k0 = int(np.argmin([cost_h(x) for x in grid]))
rr_ = minimize_scalar(cost_h, bounds=(grid[max(k0-2, 0)], grid[min(k0+2, 4000)]), method="bounded", options={"xatol": 1e-9})
hH = np.sqrt(2*MU*r1*r2/(r1 + r2))
ok(f"brute force over h: min two-impulse cost {rr_.fun:.5f} at h = {rr_.x:.2f} = Hohmann {dvH:.5f} at h_H = {hH:.2f}",
   close(rr_.fun, dvH, 1e-6) and close(rr_.x, hH, 0.5))

# ------------------------------------------------------------------ Part II: bi-elliptic
print("Part II: bi-elliptic")


def bi(R, Rs):
    d1 = np.sqrt(2*Rs/(1 + Rs)) - 1
    d2 = np.sqrt(2/Rs)*(np.sqrt(R/(R + Rs)) - np.sqrt(1/(1 + Rs)))
    d3 = np.sqrt(2*Rs/(R*(R + Rs))) - 1/np.sqrt(R)
    return d1, d2, d3


ok("bi(R, R) = Hohmann", all(close(sum(bi(R, R)), g(R), 1e-12) for R in (2, 6.31, 15, 40)))
rb = 400000.0
R5 = 15.0
d = np.array(bi(R5, rb/r1))*vc(r1)
r2_5 = 15*r1
ok(f"Example 5: r2 = {r2_5:.2f}, a_t1 = {(r1 + rb)/2:.0f}, a_t2 = {(r2_5 + rb)/2:.0f}",
   close(r2_5, 100172.06, 0.01) and close((r1 + rb)/2, 203339, 1) and close((r2_5 + rb)/2, 250086, 1))
print(f"     burns {d}, total {d.sum():.5f}, Hohmann {g(R5)*vc(r1):.5f}")
ok(f"  burns 3.11004, 0.45088, 0.52801 -> {d.sum():.5f} vs {g(R5)*vc(r1):.5f}",
   close(d[0], 3.11004, 1e-5) and close(d[1], 0.45088, 1e-5) and close(d[2], 0.52801, 1e-5)
   and close(d.sum(), 4.08892, 1e-5) and close(g(R5)*vc(r1), 4.14269, 1e-5))
at1, at2 = (r1 + rb)/2, (r2_5 + rb)/2
ok(f"  speeds at r_b: {vis(rb, at1):.5f} -> {vis(rb, at2):.5f}; at r2: {vis(r2_5, at2):.5f} -> {vc(r2_5):.5f}; at r1 -> {vis(r1, at1):.5f}",
   close(vis(rb, at1), 0.18091, 1e-5) and close(vis(rb, at2), 0.63178, 1e-5) and close(vis(r2_5, at2), 2.52279, 1e-5)
   and close(vc(r2_5), 1.99478, 1e-5) and close(vis(r1, at1), 10.83580, 1e-5))
tbi = (np.pi*np.sqrt(at1**3/MU) + np.pi*np.sqrt(at2**3/MU))/3600
tH5 = np.pi*np.sqrt(((r1 + r2_5)/2)**3/MU)/3600
ok(f"  saving {d.sum() - g(R5)*vc(r1):+.4f} ({100*(1 - d.sum()/(g(R5)*vc(r1))):.1f}%), time {tH5:.1f} h -> {tbi:.1f} h = {tbi/24:.1f} d",
   close(g(R5)*vc(r1) - d.sum(), 0.0538, 0.0001) and close(tH5, 17.1, 0.05) and close(tbi, 299.6, 0.1))
ok(f"  periapsis raised by {r2_5 - r1:.0f} km", close(r2_5 - r1, 93494, 1))
biinf = lambda R: (np.sqrt(2) - 1)*(1 + 1/np.sqrt(R))
Rdag = brentq(lambda R: biinf(R) - g(R), 5, 15)
ok(f"R-dagger = {Rdag:.6f} (11.938765)", close(Rdag, 11.938765, 1e-5))
Rs_ = sp.symbols("Rs", positive=True)
Rsym = sp.symbols("R", positive=True)
d2s = sp.sqrt(2/Rs_)*(sp.sqrt(Rsym/(Rsym + Rs_)) - sp.sqrt(1/(1 + Rs_)))
d3s = sp.sqrt(2*Rs_/(Rsym*(Rsym + Rs_))) - 1/sp.sqrt(Rsym)
dd = sp.simplify((sp.diff(d2s + d3s, Rsym)).subs(Rs_, Rsym))
ok("d(dv_bi)/dR = 0 at R* = R", dd == 0)
ok("  the two partial derivatives are +-1/(4 R^1.5)",
   sp.simplify(sp.diff(sp.sqrt(2/Rs_)*sp.sqrt(Rsym/(Rsym + Rs_)), Rsym).subs(Rs_, Rsym) - 1/(4*Rsym**sp.Rational(3, 2))) == 0)
# bi-elliptic r_b = 25 r1 curve vs Hohmann between 14.10 and 25
Rg = np.linspace(12, 25, 130001)
diff_ = np.array([g(R) - sum(bi(R, 25.0)) for R in Rg])
lo = Rg[np.argmax(diff_ > 0)]
ok(f"r_b = 25 r1 beats Hohmann from R = {lo:.2f} to 25; max gap {diff_.max():.4f} v_c1 at R = {Rg[np.argmax(diff_)]:.1f} ({100*diff_.max()/g(Rg[np.argmax(diff_)]):.2f}%)",
   close(lo, 14.10, 0.01) and close(diff_.max(), 0.0030, 0.0001) and close(Rg[np.argmax(diff_)], 18.5, 0.1) and close(100*diff_.max()/g(Rg[np.argmax(diff_)]), 0.56, 0.01))
dg = (g(25 + 1e-6) - g(25 - 1e-6))/2e-6
ok(f"|dg/dR| at R = 25: {abs(dg):.1e} (7.6e-4)", close(abs(dg), 7.6e-4, 0.05e-4))

# stopover
R30, Rb = 30.0, 25.0
r2s, rbs = 30*r1, 25*r1
at1s, at2s = (r1 + rbs)/2, (rbs + r2s)/2
s1 = vis(r1, at1s) - vc(r1)
s2 = vis(rbs, at2s) - vis(rbs, at1s)
s3 = vc(r2s) - vis(r2s, at2s)
print(f"     stopover burns {s1:.5f} {s2:.5f} {s3:.5f} total {s1+s2+s3:.5f}; Hohmann {g(30)*vc(r1):.5f}")
ok("stopover R = 30, r_b = 25 r1: 2.98794 + 1.18531 + 0.06564 = 4.23889 vs 4.07470 (+4.03%)",
   close(s1, 2.98794, 1e-5) and close(s2, 1.18531, 1e-5) and close(s3, 0.06564, 1e-5) and close(s1 + s2 + s3, 4.23889, 1e-5)
   and close(g(30)*vc(r1), 4.07470, 1e-5) and close(100*((s1+s2+s3)/(g(30)*vc(r1)) - 1), 4.03, 0.01))
ok(f"  r2 = {r2s:.2f}, r_b = {rbs:.2f}, a_t1 = {at1s:.2f}, a_t2 = {at2s:.2f}",
   close(r2s, 200344.11, 0.01) and close(rbs, 166953.43, 0.01) and close(at1s, 86815.78, 0.01) and close(at2s, 183648.77, 0.01))
ok(f"  speeds: at r_b {vis(rbs, at1s):.5f} -> {vis(rbs, at2s):.5f}; at r2 {vis(r2s, at2s):.5f} -> {vc(r2s):.5f}",
   close(vis(rbs, at1s), 0.42855, 1e-5) and close(vis(rbs, at2s), 1.61386, 1e-5) and close(vis(r2s, at2s), 1.34488, 1e-5) and close(vc(r2s), 1.41052, 1e-5))
tH30 = np.pi*np.sqrt(((r1 + r2s)/2)**3/MU)/86400
tst = (np.pi*np.sqrt(at1s**3/MU) + np.pi*np.sqrt(at2s**3/MU))/86400
ok(f"  time {tH30:.2f} d -> {tst:.2f} d", close(tH30, 1.92, 0.005) and close(tst, 6.01, 0.005))
h1 = vis(r1, (r1 + r2s)/2) - vc(r1)
ok(f"  Hohmann first burn {h1:.5f}; saved {1000*(h1 - s1):.1f} m/s; climbed {r2s - rbs:.0f} km less", close(h1, 3.02245, 1e-5) and close(1000*(h1 - s1), 34.5, 0.05) and close(r2s - rbs, 33391, 1))
# 'pays 198.7 m/s': the extra burn-2 cost relative to a burn that only matches
ok(f"  speed ratio at r_b vs r_1 on ellipse 1 = {vis(r1, at1s)/vis(rbs, at1s):.4f} (= 25)", close(vis(r1, at1s)/vis(rbs, at1s), 25, 1e-9))
st = []
for Rs in np.linspace(1.0001, 29.999, 30000):
    rbx = Rs*r1
    a1x, a2x = (r1 + rbx)/2, (rbx + r2s)/2
    st.append((vis(r1, a1x) - vc(r1)) + (vis(rbx, a2x) - vis(rbx, a1x)) + (vc(r2s) - vis(r2s, a2x)))
st = np.array(st)
Rgrid = np.linspace(1.0001, 29.999, 30000)
ok(f"  worst stopover at R* = {Rgrid[np.argmax(st)]:.2f}, +{100*(st.max()/(g(30)*vc(r1)) - 1):.0f}%",
   close(Rgrid[np.argmax(st)], 4.02, 0.02) and close(100*(st.max()/(g(30)*vc(r1)) - 1), 34, 0.5))
b40 = np.array(bi(30, 40))*vc(r1)
ok(f"  r_b = 40 r1 bi-elliptic: {b40.sum():.5f} ({100*(1 - b40.sum()/(g(30)*vc(r1))):.1f}% less)", close(b40.sum(), 4.02458, 1e-5) and close(100*(1 - b40.sum()/(g(30)*vc(r1))), 1.2, 0.05))
# decomposition of the 198.7 m/s
hb2 = g(30)*vc(r1) - h1
ok(f"  burns 2+3 cost {1000*(s2 + s3 - hb2):.1f} m/s more than Hohmann's second burn ({1000*hb2:.1f})", close(1000*(s2 + s3 - hb2), 198.7, 0.05))
print(f"     burn 2 = {1000*s2:.1f} m/s, burn 3 = {1000*s3:.1f} m/s, net penalty {1000*(s1+s2+s3 - g(30)*vc(r1)):.1f} m/s")

# ------------------------------------------------------------------ Part III: phasing
print("Part III: phasing")
r6 = 6778.0
T6 = 2*np.pi*np.sqrt(r6**3/MU)
Tph = T6*(1 - 10/360)
aph = r6*(1 - 1/36)**(2/3)
vph = vis(r6, aph)
ok(f"Example 6: T {T6:.2f} s = {T6/60:.3f} min, T_ph {Tph:.2f} s = {Tph/60:.3f} min",
   close(T6, 5553.46, 0.01) and close(T6/60, 92.558, 0.001) and close(Tph, 5399.19, 0.01) and close(Tph/60, 89.987, 0.001))
ok(f"  a_ph {aph:.2f}, r_p {2*aph - r6:.2f} ({2*aph - r6 - RE:.1f} km alt), v_ph {vph:.6f}, v_c {vc(r6):.6f}, dv {2000*(vc(r6) - vph):.1f} m/s",
   close(aph, 6651.89, 0.01) and close(2*aph - r6, 6525.79, 0.01) and close(2*aph - r6 - RE, 147.6, 0.05)
   and close(vph, 7.595597, 1e-6) and close(vc(r6), 7.668636, 1e-6) and close(2000*(vc(r6) - vph), 146.1, 0.05))
aph3 = r6*(1 - 1/(36*3))**(2/3)
ok(f"  k = 3: a_ph {aph3:.2f}, perigee alt {2*aph3 - r6 - RE:.0f} km, dv {2000*(vc(r6) - vis(r6, aph3)):.1f} m/s",
   close(aph3, 6736.10, 0.01) and close(2*aph3 - r6 - RE, 316, 0.5) and close(2000*(vc(r6) - vis(r6, aph3)), 47.8, 0.05))
a7, e7 = 10000.0, 0.3
T7 = 2*np.pi*np.sqrt(a7**3/MU)
aph7 = a7*(1 + 600/T7)**(2/3)
rp7, ra7 = a7*(1 - e7), a7*(1 + e7)
dp = vis(rp7, aph7) - vis(rp7, a7)
da = vis(ra7, aph7) - vis(ra7, a7)
ok(f"Example 7: T {T7:.2f}, a_ph {aph7:.2f}", close(T7, 9952.01, 0.01) and close(aph7, 10397.99, 0.01))
ok(f"  periapsis {vis(rp7, a7):.6f} -> {vis(rp7, aph7):.6f} ({dp:.6f}, total {2*dp:.5f}); apoapsis {vis(ra7, a7):.6f} -> {vis(ra7, aph7):.6f} ({da:.6f}, total {2*da:.5f})",
   close(vis(rp7, a7), 8.603825, 1e-6) and close(vis(rp7, aph7), 8.692036, 1e-6) and close(2*dp, 0.17642, 1e-5)
   and close(vis(ra7, a7), 4.632829, 1e-6) and close(vis(ra7, aph7), 4.794662, 1e-6) and close(2*da, 0.32367, 1e-5))
ok(f"  saving {100*(1 - dp/da):.0f}% (v_a/v_p = {(1-e7)/(1+e7):.3f} -> {100*(1 - (1-e7)/(1+e7)):.0f}%); ratio {da/dp:.3f} vs (1+e)/(1-e) = {(1+e7)/(1-e7):.3f}",
   close(100*(1 - dp/da), 45, 0.5) and close(da/dp, 1.835, 0.001) and close((1 + e7)/(1 - e7), 1.857, 0.001))

# ------------------------------------------------------------------ Part IV: apse rotation
print("Part IV: apse-line rotation")
a8, e8 = 12000.0, 0.4
p8 = a8*(1 - e8**2)
s8 = np.sqrt(MU/p8)
f8 = 15*d2r
vr8, vt8 = s8*e8*np.sin(f8), s8*(1 + e8*np.cos(f8))
ok(f"Example 8: p {p8:.0f}, sqrt(mu/p) {s8:.6f}, r {p8/(1 + e8*np.cos(f8)):.2f}, v {np.hypot(vr8, vt8):.6f}",
   close(s8, 6.288378, 1e-6) and close(p8/(1 + e8*np.cos(f8)), 7270.79, 0.01) and close(np.hypot(vr8, vt8), 8.742294, 1e-6))
ok(f"  v_r {vr8:.6f}, v_theta {vt8:.6f}, dv {2*vr8:.5f}, gamma {np.degrees(np.arctan2(vr8, vt8)):.3f}",
   close(vr8, 0.651021, 1e-6) and close(vt8, 8.718020, 1e-6) and close(2*vr8, 1.30204, 1e-5) and close(np.degrees(np.arctan2(vr8, vt8)), 4.271, 0.001))
f8b = 195*d2r
vr8b, vt8b = s8*e8*np.sin(f8b), s8*(1 + e8*np.cos(f8b))
ok(f"  second crossing: v {np.hypot(vr8b, vt8b):.3f}, gamma {np.degrees(np.arctan2(vr8b, vt8b)):.3f}, dv {2*abs(vr8b):.5f}",
   close(np.hypot(vr8b, vt8b), 3.913, 0.001) and close(np.degrees(np.arctan2(vr8b, vt8b)), -9.576, 0.001) and close(2*abs(vr8b), 1.30204, 1e-5))
ra8 = a8*(1 + e8)
alt = 2*(vc(ra8) - vis(ra8, a8))
ok(f"  alternative 2({vc(ra8):.6f} - {vis(ra8, a8):.6f}) = {alt:.5f}", close(vc(ra8), 4.870956, 1e-6) and close(vis(ra8, a8), 3.773027, 1e-6) and close(alt, 2.19586, 1e-5))
eta_c = 2*np.arcsin(alt/(2*s8*e8))
ok(f"  crossover eta = {np.degrees(eta_c):.1f} deg; eta = 180: {2*s8*e8:.2f}", close(np.degrees(eta_c), 51.8, 0.05) and close(2*s8*e8, 5.03, 0.005))

# ------------------------------------------------------------------ Part VI: plane change
print("Part VI: plane change")
di = 28.5*d2r
pc = lambda vv, th: 2*vv*np.sin(th/2)
ok(f"Example 11: {pc(vc(r1), di):.5f}, {pc(vc(r2), di):.5f}, {pc(vis(r2, at), di):.5f}; ratio {pc(vc(r1), di)/pc(vis(r2, at), di):.1f}",
   close(pc(vc(r1), di), 3.80344, 1e-5) and close(pc(vc(r2), di), 1.51368, 1e-5) and close(pc(vis(r2, at), di), 0.79155, 1e-5)
   and close(pc(vc(r1), di)/pc(vis(r2, at), di), 4.8, 0.05))
lc = lambda x, y, th: np.sqrt(x**2 + y**2 - 2*x*y*np.cos(th))
c_geo = lc(vis(r2, at), vc(r2), di)
c_leo = lc(vc(r1), vis(r1, at), di)
ok(f"Example 12: separate {dvH + pc(vc(r2), di):.5f}; folded into LEO burn {c_leo:.5f} + {dvH2:.5f} = {c_leo + dvH2:.5f}; folded into GEO burn {dvH1:.5f} + {c_geo:.5f} = {dvH1 + c_geo:.5f}",
   close(dvH + pc(vc(r2), di), 5.40623, 1e-5) and close(c_leo, 4.98923, 1e-5) and close(c_leo + dvH2, 6.45606, 1e-5)
   and close(c_geo, 1.83023, 1e-5) and close(dvH1 + c_geo, 4.25596, 1e-5))
ok(f"  saving {dvH + pc(vc(r2), di) - (dvH1 + c_geo):.2f} km/s = {100*(1 - (dvH1 + c_geo)/(dvH + pc(vc(r2), di))):.0f}%",
   close(dvH + pc(vc(r2), di) - (dvH1 + c_geo), 1.15, 0.005) and close(100*(1 - (dvH1 + c_geo)/(dvH + pc(vc(r2), di))), 21, 0.5))
tot = lambda dl: lc(vc(r1), vis(r1, at), dl) + lc(vis(r2, at), vc(r2), di - dl)
rs_ = minimize_scalar(tot, bounds=(0, di), method="bounded", options={"xatol": 1e-12})
dl = rs_.x
L = vc(r1)*vis(r1, at)*np.sin(dl)/lc(vc(r1), vis(r1, at), dl)
Rr = vis(r2, at)*vc(r2)*np.sin(di - dl)/lc(vis(r2, at), vc(r2), di - dl)
ok(f"Example 13: delta* = {np.degrees(dl):.2f}, dv1 {lc(vc(r1), vis(r1, at), dl):.5f}, dv2 {lc(vis(r2, at), vc(r2), di - dl):.5f}, total {tot(dl):.5f}; sides {L:.5f} {Rr:.5f}",
   close(np.degrees(dl), 2.20, 0.005) and close(lc(vc(r1), vis(r1, at), dl), 2.44945, 1e-5) and close(lc(vis(r2, at), vc(r2), di - dl), 1.78186, 1e-5)
   and close(tot(dl), 4.23131, 1e-5) and close(L, 1.22925, 1e-4) and close(1000*(dvH1 + c_geo - tot(dl)), 24.7, 0.06))
ok("plane change of 60 deg costs v; 180 deg costs 2v", close(pc(1, 60*d2r), 1, 1e-12) and close(pc(1, np.pi), 2, 1e-12))
i14, dO = 51.6*d2r, 10*d2r
ca = np.cos(i14)**2 + np.sin(i14)**2*np.cos(dO)
al = np.arccos(ca)
ok(f"Example 14: cos alpha {ca:.6f}, alpha {np.degrees(al):.3f}, dv {pc(vc(r6), al):.4f}",
   close(ca, 0.990669, 1e-6) and close(np.degrees(al), 7.833, 0.001) and close(pc(vc(r6), al), 1.0476, 1e-4))
J2 = 1.08262668e-3
n6 = np.sqrt(MU/r6**3)
Od = -1.5*J2*(RE/r6)**2*n6*np.cos(i14)*86400/d2r
ok(f"J2 node rate {Od:.4f} deg/day; 10 deg in {10/abs(Od):.1f} days", close(Od, -5.0027, 0.0005) and close(10/abs(Od), 2.0, 0.05))

# ------------------------------------------------------------------ Part VII: escape
print("Part VII: escape")
ok("escape cost sqrt(2x) - sqrt(2x-1) decreasing in x", all(np.diff(np.sqrt(2*np.linspace(0.5, 50, 1000)) - np.sqrt(2*np.linspace(0.5, 50, 1000) - 1)) < 0))
rows = {0: (3.2001, 3.2001), 0.3: (2.1172, 3.2742), 0.5: (1.4638, 3.1540), 0.7: (0.8527, 2.8122), 0.9: (0.2766, 1.9461)}
for e_, (pp, aa) in rows.items():
    a_ = r1/(1 - e_)
    rap = r1*(1 + e_)/(1 - e_)
    cp = np.sqrt(2*MU/r1) - vis(r1, a_)
    ca_ = np.sqrt(2*MU/rap) - vis(rap, a_)
    ok(f"  e = {e_}: periapsis {cp:.4f}, apoapsis {ca_:.4f}, penalty +{100*(ca_/cp - 1):.0f}%", close(cp, pp, 1e-4) and close(ca_, aa, 1e-4))
ok("  e = 0.9: v_p/v_esc = {:.2f}".format(vis(r1, r1/0.1)/np.sqrt(2*MU/r1)), close(vis(r1, r1/0.1)/np.sqrt(2*MU/r1), 0.97, 0.005))
q, Rq = sp.symbols("q R", positive=True)
uu = sp.symbols("u", positive=True)
tq = sp.integrate((uu + q)/sp.sqrt(uu), (uu, 0, Rq - q))/sp.sqrt(2*sp.Symbol("mu", positive=True))
ok("parabolic time of flight t = (2/3) sqrt((R-q)/(2mu)) (R+2q)",
   sp.simplify(tq - sp.Rational(2, 3)*sp.sqrt((Rq - q)/(2*sp.Symbol("mu", positive=True)))*(Rq + 2*q)) == 0)
ok("dt/dq proportional to R - 2q", sp.simplify(sp.diff(sp.Rational(2, 3)*sp.sqrt(Rq - q)*(Rq + 2*q), q)*sp.sqrt(Rq - q) - (Rq - 2*q)) == 0)
RSOI = 1.495978707e8*3.0034896e-6**0.4
ok(f"R_SOI = {RSOI:.3e} km, R/2 = {RSOI/2:.3e}", close(RSOI, 9.25e5, 0.005e5) and close(RSOI/2, 4.62e5, 0.01e5))
# savings quoted at e = 0.9, 0.95
for e_, sv in ((0.9, 1.67), (0.95, 1.34)):
    a_ = r1/(1 - e_); rap = r1*(1 + e_)/(1 - e_)
    s_ = (np.sqrt(2*MU/rap) - vis(rap, a_)) - (np.sqrt(2*MU/r1) - vis(r1, a_))
    ok(f"  e = {e_}: periapsis burn saves {s_:.2f} km/s; half period {np.pi*np.sqrt(a_**3/MU)/86400:.2f} d", close(s_, sv, 0.005))

# ------------------------------------------------------------------ Part VIII
print("Part VIII: low thrust")
ed = np.sqrt(vc(r1)**2 + vc(r2)**2 - 2*vc(r1)*vc(r2)*np.cos(np.pi/2*di))
ok(f"Edelbaum {ed:.3f} km/s; propellant {100*(1 - np.exp(-ed/(1800*g0))):.1f}% at 1800 s vs {100*(1 - np.exp(-tot(dl)/(315*g0))):.1f}% at 315 s; +{100*(ed/tot(dl) - 1):.0f}% dv",
   close(ed, 5.951, 0.001) and close(100*(1 - np.exp(-ed/(1800*g0))), 28.6, 0.05) and close(100*(1 - np.exp(-tot(dl)/(315*g0))), 74.6, 0.05)
   and close(100*(ed/tot(dl) - 1), 41, 0.5))
ok("Edelbaum with di = 0 is |v1 - v2|", close(np.sqrt(vc(r1)**2 + vc(r2)**2 - 2*vc(r1)*vc(r2)), vc(r1) - vc(r2), 1e-12))

print("Task answers added in the note")
ok(f"Oberth task: v_p {vp_o:.3f}, v_a {va_o:.3f}, d(eps) {vp_o*0.1+0.005:.3f} vs {va_o*0.1+0.005:.3f}",
   close(vp_o, 9.982, 0.001) and close(va_o, 3.327, 0.001) and close(vp_o*0.1 + 0.005, 1.003, 0.001) and close(va_o*0.1 + 0.005, 0.338, 0.001))
ok(f"g(6.31) {g(r2/r1):.4f}, g(15) {g(15):.4f}, g(100) {g(100):.4f}", close(g(r2/r1), 0.5038, 1e-4) and close(g(15), 0.5362, 1e-4) and close(g(100), 0.4931, 1e-4))
f0, f1 = np.exp(-dvH/ve), np.exp(-(dvH - 0.05)/ve)
ok(f"50 m/s saving: final mass {100*f0:.2f}% -> {100*f1:.2f}%, +{3000*(f1 - f0):.0f} kg on 3 t, +{100*(f1/f0 - 1):.1f}%",
   close(100*f0, 28.36, 0.005) and close(100*f1, 28.83, 0.005) and close(3000*(f1 - f0), 14, 0.5) and close(100*(f1/f0 - 1), 1.6, 0.05))