"""Verification of every identity and number quoted in the two-body note (two-body-updated.qmd).

Part A is symbolic (sympy); Part B checks the numbers in the text, the worked examples and the
tasks; Part C integrates the two-body problem and compares with the closed-form results.
Run:  python3 codes/two-body-verify.py        (a few seconds)
"""
import numpy as np
import sympy as sp
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
from scipy.special import jv

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

# ======================================================================================= A. symbolic
print("A. symbolic")
mu, h, p, a, r, t = sp.symbols("mu h p a r t", positive=True)
e, f, E, H, D = sp.symbols("e f E H D", real=True)

# A1. relative equation from the two N-body equations
G, m1, m2 = sp.symbols("G m1 m2", positive=True)
rv = sp.Matrix(sp.symbols("x y z", real=True))
R = sp.sqrt(rv.dot(rv))
acc1 = G*m2*rv/R**3            # body 1 pulled towards body 2 (r = R2 - R1)
acc2 = -G*m1*rv/R**3
ok("r'' = R2'' - R1'' = -G(m1+m2) r/r^3", sp.simplify(acc2 - acc1 + G*(m1 + m2)*rv/R**3) == sp.zeros(3, 1))
A = sp.Matrix([[m1/(m1 + m2), m2/(m1 + m2)], [-1, 1]])
ok("det A = 1 and A^-1 = [[1, -m2/M], [1, m1/M]]",
   sp.simplify(A.det()) == 1 and sp.simplify(A.inv() - sp.Matrix([[1, -m2/(m1 + m2)], [1, m1/(m1 + m2)]])) == sp.zeros(2, 2))

# A2. integrals with a generic state vector
X = sp.Matrix(sp.symbols("X1:4", real=True))
Vv = sp.Matrix(sp.symbols("V1:4", real=True))
rn = sp.sqrt(X.dot(X))
hv = X.cross(Vv)
ev = Vv.cross(hv)/mu - X/rn
eps = Vv.dot(Vv)/2 - mu/rn
ok("Casimir 1: h . e = 0", sp.simplify(hv.dot(ev)) == 0)
ok("Casimir 2: mu^2 (e^2 - 1) = 2 eps h^2", sp.simplify(mu**2*(ev.dot(ev) - 1) - 2*eps*hv.dot(hv)) == 0)
ok("e . v = -rdot (so e is not perpendicular to v)", sp.simplify(ev.dot(Vv) + X.dot(Vv)/rn) == 0)
ok("r . (v x h) = h^2", sp.simplify(X.dot(Vv.cross(hv)) - hv.dot(hv)) == 0)
# de/dt = 0 along the flow
accv = -mu*X/rn**3
de = sp.Matrix([sum(sp.diff(ev[k], X[j])*Vv[j] + sp.diff(ev[k], Vv[j])*accv[j] for j in range(3)) for k in range(3)])
ok("de/dt = 0 along r'' = -mu r/r^3", sp.simplify(de) == sp.zeros(3, 1))
dh = sp.Matrix([sum(sp.diff(hv[k], X[j])*Vv[j] + sp.diff(hv[k], Vv[j])*accv[j] for j in range(3)) for k in range(3)])
ok("dh/dt = 0", sp.simplify(dh) == sp.zeros(3, 1))
deps = sum(sp.diff(eps, X[j])*Vv[j] + sp.diff(eps, Vv[j])*accv[j] for j in range(3))
ok("d eps/dt = 0", sp.simplify(deps) == 0)
# hodograph form v = (mu/h) hhat x (rhat + e)
hn = sp.sqrt(hv.dot(hv))
vh = (mu/hn**2)*hv.cross(X/rn + ev)
ok("v = (mu/h) hhat x (rhat + e)", sp.simplify(vh - Vv) == sp.zeros(3, 1))

# A3. orbit equation, velocity, flight path angle
rf = p/(1 + e*sp.cos(f))
vr = mu/h*e*sp.sin(f)
vt = mu/h*(1 + e*sp.cos(f))
fdot = h/rf**2
ok("v_r = dr/dt = (mu/h) e sin f  (with p = h^2/mu)", sp.simplify((sp.diff(rf, f)*fdot - vr).subs(p, h**2/mu)) == 0)
ok("v_perp = r fdot = (mu/h)(1 + e cos f)", sp.simplify((rf*fdot - vt).subs(p, h**2/mu)) == 0)
ok("v^2 = (mu/p)(1 + 2e cos f + e^2)", sp.simplify((vr**2 + vt**2 - mu/p*(1 + 2*e*sp.cos(f) + e**2)).subs(p, h**2/mu)) == 0)
ok("vis-viva v^2 = mu(2/r - 1/a), a = p/(1-e^2)",
   sp.simplify((mu/p*(1 + 2*e*sp.cos(f) + e**2)) - mu*(2/rf - (1 - e**2)/p)) == 0)
tg = e*sp.sin(f)/(1 + e*sp.cos(f))
ok("d(tan gamma)/df = e(cos f + e)/(1 + e cos f)^2", sp.simplify(sp.diff(tg, f) - e*(sp.cos(f) + e)/(1 + e*sp.cos(f))**2) == 0)
ok("e = 1: tan gamma = tan(f/2)", sp.simplify((tg.subs(e, 1) - sp.tan(f/2)).rewrite(sp.tan)) == 0)
# hodograph components
vP = -mu/h*sp.sin(f)
vQ = mu/h*(e + sp.cos(f))
ok("perifocal components = rotation of (v_r, v_perp)",
   sp.simplify(vP - (vr*sp.cos(f) - vt*sp.sin(f))) == 0 and sp.simplify(vQ - (vr*sp.sin(f) + vt*sp.cos(f))) == 0)
ok("hodograph circle v_P^2 + (v_Q - mu e/h)^2 = (mu/h)^2", sp.simplify(vP**2 + (vQ - mu*e/h)**2 - (mu/h)**2) == 0)

# A4. energy at periapsis
ok("eps = h^2/(2 r_p^2) - mu/r_p = -mu/(2a)",
   sp.simplify((mu*a*(1 - e**2)/(2*a**2*(1 - e)**2) - mu/(a*(1 - e))) + mu/(2*a)) == 0)

# A5. eccentric anomaly
xE = a*(sp.cos(E) - e)
yE = a*sp.sqrt(1 - e**2)*sp.sin(E)
ok("r^2 = a^2 (1 - e cos E)^2", sp.simplify(xE**2 + yE**2 - a**2*(1 - e*sp.cos(E))**2) == 0)
ok("x dy - y dx = a^2 sqrt(1-e^2)(1 - e cos E) dE",
   sp.simplify(xE*sp.diff(yE, E) - yE*sp.diff(xE, E) - a**2*sp.sqrt(1 - e**2)*(1 - e*sp.cos(E))) == 0)
ok("1 + cos f = (1-e)(1+cos E)/(1 - e cos E)",
   sp.simplify(1 + (sp.cos(E) - e)/(1 - e*sp.cos(E)) - (1 - e)*(1 + sp.cos(E))/(1 - e*sp.cos(E))) == 0)
# numerical check of half-angle and df/dE = b/r
for ee in (0.1, 0.6, 0.95):
    EE = np.linspace(0.1, 6.1, 13)
    ff = 2*np.arctan2(np.sqrt(1 + ee)*np.sin(EE/2), np.sqrt(1 - ee)*np.cos(EE/2))
    cf = (np.cos(EE) - ee)/(1 - ee*np.cos(EE))
    sf = np.sqrt(1 - ee**2)*np.sin(EE)/(1 - ee*np.cos(EE))
    ok(f"half-angle E -> f matches cos f, sin f (e = {ee})", np.allclose(np.cos(ff), cf) and np.allclose(np.sin(ff), sf))
dfdE = sp.sqrt(1 - e**2)/(1 - e*sp.cos(E))
ftrue = 2*sp.atan(sp.sqrt((1 + e)/(1 - e))*sp.tan(E/2))
_fd = sp.lambdify((E, e), sp.diff(ftrue, E) - dfdE)
ok("df/dE = sqrt(1-e^2)/(1 - e cos E) = b/r", all(abs(_fd(EE, ee)) < 1e-12 for EE in (0.3, 1.2, 2.9, 4.0) for ee in (0.1, 0.7)))
# small-e expansions
ok("f = E + e sin E + O(e^2) (numerical)", all(abs(2*np.arctan(np.sqrt((1 + ee)/(1 - ee))*np.tan(EE/2)) - EE - ee*np.sin(EE)) < 3*ee**2
                                          for ee in (1e-3, 2e-3) for EE in (0.4, 1.3, 2.5)))
M_ = sp.Symbol("M")
for ee in (0.001, 0.002):
    for MM in (0.4, 1.3, 2.5):
        EE = brentq(lambda x: x - ee*np.sin(x) - MM, 0, np.pi)
        fft = 2*np.arctan(np.sqrt((1 + ee)/(1 - ee))*np.tan(EE/2))
        assert abs(fft - (MM + 2*ee*np.sin(MM))) < 5*ee**2
ok("f = M + 2e sin M + O(e^2) (numerical)", True)

# A6. hyperbolic anomaly
xH = a*(e - sp.cosh(H))
yH = a*sp.sqrt(e**2 - 1)*sp.sinh(H)
ok("r^2 = |a|^2 (e cosh H - 1)^2", sp.simplify((xH**2 + yH**2 - a**2*(e*sp.cosh(H) - 1)**2).rewrite(sp.exp)) == 0)
ok("x dy - y dx = |a|^2 sqrt(e^2-1)(e cosh H - 1) dH",
   sp.simplify((xH*sp.diff(yH, H) - yH*sp.diff(xH, H) - a**2*sp.sqrt(e**2 - 1)*(e*sp.cosh(H) - 1)).rewrite(sp.exp)) == 0)
for ee in (1.2, 2.0):
    HH = np.linspace(-3, 3, 13)
    rr = (e_ := ee)*np.cosh(HH) - 1
    cf = (ee - np.cosh(HH))/rr
    ff = 2*np.arctan(np.sqrt((ee + 1)/(ee - 1))*np.tanh(HH/2))
    ok(f"tanh(H/2) = sqrt((e-1)/(e+1)) tan(f/2) (e = {ee})", np.allclose(np.cos(ff), cf))

# A7. Barker and Cardano
Mp = sp.Symbol("M_p", real=True)
s_ = sp.sqrt(9*Mp**2/4 + 1)
Dc = sp.cbrt(3*Mp/2 + s_) - sp.cbrt(s_ - 3*Mp/2)
for val in (0.3, 1.0, 4.0):
    dval = float(Dc.subs(Mp, val))
    ok(f"Cardano root solves D^3 + 3D - 3M_p = 0 (M_p = {val})", abs(dval**3 + 3*dval - 3*val) < 1e-12)
ok("sec^4(f/2) df = 2(1 + D^2) dD", sp.simplify(sp.sec(f/2)**4 - 2*(1 + sp.tan(f/2)**2)*sp.diff(sp.tan(f/2), f)) == 0)

# A8. polar acceleration and Binet
ff_ = sp.Function("f")(t)
r_ = sp.Function("r")(t)
pos = r_*sp.Matrix([sp.cos(ff_), sp.sin(ff_)])
acc = pos.diff(t, 2)
rh = sp.Matrix([sp.cos(ff_), sp.sin(ff_)])
th = sp.Matrix([-sp.sin(ff_), sp.cos(ff_)])
ok("radial acceleration r'' - r f'^2", sp.simplify(acc.dot(rh) - (r_.diff(t, 2) - r_*ff_.diff(t)**2)) == 0)
ok("transverse acceleration r f'' + 2 r' f' = (1/r) d(r^2 f')/dt",
   sp.simplify(acc.dot(th) - (r_*ff_.diff(t, 2) + 2*r_.diff(t)*ff_.diff(t))) == 0 and
   sp.simplify(acc.dot(th) - (r_**2*ff_.diff(t)).diff(t)/r_) == 0)
u = sp.Function("u")(f)
ok("Binet: u = (1 + e cos f)/p solves u'' + u = mu/h^2 with p = h^2/mu",
   sp.simplify(((1 + e*sp.cos(f))/p).diff(f, 2) + (1 + e*sp.cos(f))/p - 1/p) == 0)

# A9. universal Kepler equation reduces to the elliptic one
for (aa, ee) in ((1.0, 0.3), (2.0, 0.8)):
    rp = aa*(1 - ee)
    for EE in (0.3, 1.5, 3.0):
        chi = np.sqrt(aa)*EE
        z = chi**2/aa
        S = (np.sqrt(z) - np.sin(np.sqrt(z)))/z**1.5
        lhs = rp*chi + (1 - rp/aa)*chi**3*S
        rhs = np.sqrt(aa**3)*(EE - ee*np.sin(EE))     # sqrt(mu)(t - tp) with mu = 1
        assert abs(lhs - rhs) < 1e-12
ok("universal: sqrt(mu)(t-tp) = r_p chi + (1 - r_p/a) chi^3 S(z) equals the elliptic Kepler equation", True)

# ======================================================================================= B. numbers
print("B. numbers in the text")
MU_E, MU_S, MU_J, MU_M = 398600.4418, 1.32712440018e11, 1.26686534e8, 4902.800
R_E, R_M, AU = 6371.0, 1737.4, 1.495978707e8

# mu of the Sun from the Earth's orbit
T_yr = 365.25636*86400
mu_s = 4*np.pi**2*AU**3/T_yr**2
ok(f"mu_sun from a = 1 AU, T = 1 sidereal year: {mu_s:.4e} km^3/s^2", close(mu_s/MU_S, 1, 1e-4))

# Pluto-Charon barycentre
q = 1.586e21/1.303e22
dPC = 19596*q/(1 + q)
ok(f"Pluto-Charon barycentre {dPC:.0f} km from Pluto's centre (> R_Pluto = 1188 km)", close(dPC, 2126, 3) and dPC > 1188)
ok("Sun-Jupiter: barycentre 1.07 R_sun", close(7.785e8/1048.35/6.957e5, 1.07, 0.01))

# circular and escape speeds
r_iss = R_E + 420
vc_iss = np.sqrt(MU_E/r_iss)
ok(f"ISS: v_c = {vc_iss:.2f} km/s (7.66)", close(vc_iss, 7.66, 0.005))
ok(f"ISS: v_esc = {np.sqrt(2)*vc_iss:.2f} km/s (10.83)", close(np.sqrt(2)*vc_iss, 10.83, 0.005))
T_iss = 2*np.pi*np.sqrt(r_iss**3/MU_E)
ok(f"ISS period {T_iss/60:.1f} min (92.8)", close(T_iss/60, 92.8, 0.05))
ok(f"escape at Earth's surface {np.sqrt(2*MU_E/R_E):.2f} km/s (11.19)", close(np.sqrt(2*MU_E/R_E), 11.19, 0.005))
ok(f"escape at r = 6700 km {np.sqrt(2*MU_E/6700):.2f} km/s (10.91)", close(np.sqrt(2*MU_E/6700), 10.91, 0.005))
ok(f"escape at GEO {np.sqrt(2*MU_E/42164):.2f} km/s (4.35)", close(np.sqrt(2*MU_E/42164), 4.35, 0.005))
ok(f"escape from the Sun at 1 AU {np.sqrt(2*MU_S/AU):.1f} km/s (42.1)", close(np.sqrt(2*MU_S/AU), 42.1, 0.05))
ok(f"Earth's orbital speed {np.sqrt(MU_S/AU):.1f} km/s (29.8)", close(np.sqrt(MU_S/AU), 29.8, 0.05))
ok(f"escape from the Moon's surface {np.sqrt(2*MU_M/R_M):.2f} km/s (2.38)", close(np.sqrt(2*MU_M/R_M), 2.38, 0.005))
r_geo = (MU_E*(86164.0905/(2*np.pi))**2)**(1/3)
ok(f"GEO radius from the sidereal day: {r_geo:.0f} km (42164)", close(r_geo, 42164, 1))

# precision of mu vs G
ok("GM_earth known to ~2e-9, G to ~2.2e-5 (CODATA 2018)", 0.0008/398600.4418 < 3e-9 and 0.00015/6.67430 > 2e-5)

# equal-a question
for ee in (0.1, 0.9):
    print(f"     e = {ee}: v_p/sqrt(mu/a) = {np.sqrt((1 + ee)/(1 - ee)):.3f}, h/sqrt(mu a) = {np.sqrt(1 - ee**2):.3f}")
ok("v_p ratio 4.36/1.11 = 3.9", close(np.sqrt(1.9/0.1)/np.sqrt(1.1/0.9), 3.94, 0.01))

# Example 2 (two equal masses, G = 1)
m = np.array([0.5, 0.5])
R1, R2 = np.array([0.5, 0, 0]), np.array([-0.5, 0, 0])
V1, V2 = np.array([0, 0.5, 0]), np.array([0, 0, 0.5])
rr0, vv0 = R2 - R1, V2 - V1
hh = np.cross(rr0, vv0)
en = vv0 @ vv0/2 - 1/np.linalg.norm(rr0)
aa = -1/(2*en)
ee_ = np.sqrt(1 + 2*en*(hh @ hh))
TT = 2*np.pi*np.sqrt(aa**3)
ok(f"Example: h = {hh} (= (0, 1/2, 1/2)), |h| = {np.linalg.norm(hh):.4f}", np.allclose(hh, [0, 0.5, 0.5]))
ok(f"Example: eps = {en}, a = {aa:.4f}, e = {ee_:.4f}, T = {TT:.4f}",
   close(en, -0.75, 1e-12) and close(aa, 2/3, 1e-12) and close(ee_, 0.5, 1e-12) and close(TT, 3.4201, 1e-4))
ok("Example: r(0) = 1 = a(1+e): apoapsis", close(aa*(1 + ee_), 1.0, 1e-12))
Vcm = (m[0]*V1 + m[1]*V2)/m.sum()
ok("Example: V_cm = (0, 1/4, 1/4) parallel to h", np.allclose(Vcm, [0, 0.25, 0.25]) and np.linalg.norm(np.cross(Vcm, hh)) < 1e-15)

# periapsis burn: d r_a / d v = 4 a^2 v_p / mu for a 400 km circular orbit
a4 = R_E + 407
coef = 4*a4**2*np.sqrt(MU_E/a4)/MU_E
ok(f"tangential burn on a ~400 km circular orbit: d r_a/dv = {coef:.0f} km per km/s (~3500)", close(coef, 3500, 50))

# GTO
rp, ra = 6678.0, 42164.0
a_g, e_g = (rp + ra)/2, (ra - rp)/(ra + rp)
vp, va = np.sqrt(MU_E/a_g*(1 + e_g)/(1 - e_g)), np.sqrt(MU_E/a_g*(1 - e_g)/(1 + e_g))
print(f"     GTO a = {a_g:.0f}, e = {e_g:.4f}, vp = {vp:.3f}, va = {va:.3f}, ratio {vp/va:.3f}")
ok("GTO numbers", close(a_g, 24421, 0.5) and close(e_g, 0.7265, 1e-4) and close(vp, 10.15, 0.005)
   and close(va, 1.608, 0.001) and close(vp/va, 6.31, 0.005))
ok(f"GTO steepest climb at f = {np.degrees(np.arccos(-e_g)):.1f} deg, gamma = {np.degrees(np.arcsin(e_g)):.1f} deg, v = {np.sqrt(MU_E/a_g):.2f}",
   close(np.degrees(np.arccos(-e_g)), 136.6, 0.05) and close(np.degrees(np.arcsin(e_g)), 46.6, 0.05) and close(np.sqrt(MU_E/a_g), 4.04, 0.005))
ok("GTO: v_p v_a = mu/a", close(vp*va, MU_E/a_g, 1e-9))

# gamma_max values
for ee, g in ((0.0167, 0.96), (0.2, 11.54), (0.8, 53.13), (0.6, 36.87)):
    ok(f"gamma_max(e = {ee}) = {np.degrees(np.arcsin(ee)):.2f} deg", close(np.degrees(np.arcsin(ee)), g, 0.01))
ok("e = 0.6: steepest at f = 126.87 deg", close(np.degrees(np.arccos(-0.6)), 126.87, 0.005))
gam = lambda ee, fd: np.degrees(np.arctan2(ee*np.sin(fd*d2r), 1 + ee*np.cos(fd*d2r)))
ok(f"e = 0.6: gamma(45) = {gam(0.6, 45):.1f}, gamma(135) = {gam(0.6, 135):.1f}, gamma(225) = {gam(0.6, 225):.1f}",
   close(gam(0.6, 45), 16.6, 0.05) and close(gam(0.6, 135), 36.4, 0.05) and close(gam(0.6, 225), -36.4, 0.05))

# hyperbola: gamma -> 90 deg at the asymptote
for ee in (1.2, 2.0, 5.0):
    finf = np.arccos(-1/ee)
    ok(f"e = {ee}: gamma -> {gam(ee, np.degrees(finf) - 1e-6):.4f} deg as f -> f_inf", gam(ee, np.degrees(finf) - 1e-6) > 89.99)
ok("e = 2: gamma(119 deg) ~ 89 deg", close(gam(2.0, 119), 89.0, 0.1))

# f - E lead
for ee, lead in ((0.2, 11.5), (0.95, 92.8)):
    EE = np.linspace(0, np.pi, 200001)
    ff = 2*np.arctan(np.sqrt((1 + ee)/(1 - ee))*np.tan(EE/2))
    k = np.argmax(ff - EE)
    rk = 1 - ee*np.cos(EE[k])
    ok(f"max lead f - E at e = {ee}: {np.degrees((ff - EE).max()):.1f} deg, at r/a = {rk:.4f} = b/a = {np.sqrt(1 - ee**2):.4f}",
       close(np.degrees((ff - EE).max()), lead, 0.1) and close(rk, np.sqrt(1 - ee**2), 1e-4))

# Bessel series vs power series
def kepler(M, ee):
    Ek = M if ee < 0.8 else np.pi
    for _ in range(60):
        Ek -= (Ek - ee*np.sin(Ek) - M)/(1 - ee*np.cos(Ek))
    return Ek
ee, MM = 0.9, 1.0
Eb = MM + sum(2/k*jv(k, k*ee)*np.sin(k*MM) for k in range(1, 400))
ok(f"Bessel series converges at e = 0.9 (> Laplace limit): {Eb:.10f} vs Newton {kepler(MM, ee):.10f}", close(Eb, kepler(MM, ee), 1e-8))
lap = brentq(lambda x: x*np.exp(np.sqrt(1 + x*x))/(1 + np.sqrt(1 + x*x)) - 1, 0.5, 0.8)
ok(f"Laplace limit = {lap:.6f} (0.662743)", close(lap, 0.662743, 1e-6))

# Newton task: e = 0.5, M = 1 rad
Es = [1.0]
for _ in range(4):
    Es.append(Es[-1] - (Es[-1] - 0.5*np.sin(Es[-1]) - 1.0)/(1 - 0.5*np.cos(Es[-1])))
print("     Newton e = 0.5, M = 1:", [f"{x:.9f}" for x in Es])
ok("Newton task converges to E = 1.498701", close(Es[-1], 1.498701, 1e-6))
f_task = 2*np.arctan(np.sqrt(3)*np.tan(Es[-1]/2))
ok(f"  ... f = {np.degrees(f_task):.2f} deg, r/a = {1 - 0.5*np.cos(Es[-1]):.4f}", close(np.degrees(f_task), 116.36, 0.01) and close(1 - 0.5*np.cos(Es[-1]), 0.9640, 1e-4))

# Tundra example
a_t, e_t = 42164.0, 0.25
n_t = np.sqrt(MU_E/a_t**3)
ok(f"Tundra: r_p {a_t*(1 - e_t):.0f}, r_a {a_t*(1 + e_t):.0f}, b {a_t*np.sqrt(1 - e_t**2):.0f}",
   close(a_t*(1 - e_t), 31623, 0.5) and close(a_t*(1 + e_t), 52705, 0.5) and close(a_t*np.sqrt(1 - e_t**2), 40825, 0.5))
ok(f"Tundra: n = {n_t:.5e}, T = {2*np.pi/n_t:.1f} s", close(n_t, 7.29216e-5, 2e-10) and close(2*np.pi/n_t, 86163.6, 0.5))
E0 = 2*np.arctan(np.sqrt((1 - e_t)/(1 + e_t))*np.tan(15*d2r))
M0 = E0 - e_t*np.sin(E0)
ok(f"Tundra: E0 = {np.degrees(E0):.4f} deg = {E0:.6f}, M0 = {M0:.6f} = {np.degrees(M0):.4f} deg",
   close(np.degrees(E0), 23.4508, 1e-4) and close(M0, 0.309804, 1e-6))
ok(f"Tundra: t0 - tp = {M0/n_t:.1f} s = {M0/n_t/3600:.3f} h, r0 = {a_t*(1 - e_t*np.cos(E0)):.0f}",
   close(M0/n_t, 4248.5, 0.2) and close(a_t*(1 - e_t*np.cos(E0)), 32494, 0.5))
M1 = M0 + n_t*21600
ok(f"Tundra: M1 = {M1:.6f} = {np.degrees(M1):.3f} deg", close(M1, 1.884910, 1e-6))
Ek = M1
rows = []
for k in range(4):
    g = Ek - e_t*np.sin(Ek) - M1
    dE = -g/(1 - e_t*np.cos(Ek))
    rows.append((Ek, g, dE))
    Ek += dE
for k, (Ek_, g, dE) in enumerate(rows):
    print(f"     k = {k}: E = {Ek_:.6f}, g = {g:+.3e}, dE = {dE:+.3e}")
ok("Tundra Newton table", close(rows[1][0], 2.105629, 1e-6) and close(rows[2][0], 2.100635, 1e-6) and close(rows[3][0], 2.100633, 1e-6)
   and close(rows[0][1], -2.378e-1, 1e-4) and close(rows[1][1], 5.630e-3, 1e-6) and close(rows[2][1], 2.685e-6, 1e-9))
E1 = Ek
f1 = 2*np.arctan(np.sqrt((1 + e_t)/(1 - e_t))*np.tan(E1/2))
r1 = a_t*(1 - e_t*np.cos(E1))
p_t = a_t*(1 - e_t**2)
ok(f"Tundra: E1 = {np.degrees(E1):.4f}, f1 = {np.degrees(f1):.3f}, r1 = {r1:.0f}, p = {p_t:.0f}",
   close(np.degrees(E1), 120.3574, 1e-4) and close(np.degrees(f1), 132.118, 1e-3) and close(r1, 47491, 0.5) and close(p_t, 39529, 1))
ok("Tundra: r from f agrees", close(p_t/(1 + e_t*np.cos(f1)), r1, 1e-6))
v1 = np.sqrt(MU_E*(2/r1 - 1/a_t))
ok(f"Tundra: v1 = {v1:.3f} km/s, gamma1 = {gam(e_t, np.degrees(f1)):.2f} deg, 1 - e cos E1 = {1 - e_t*np.cos(E1):.3f}",
   close(v1, 2.708, 0.001) and close(gam(e_t, np.degrees(f1)), 12.56, 0.01) and close(1 - e_t*np.cos(E1), 1.126, 0.001))

# parabola tasks (Sun, r_p = 0.5 AU)
rp_c = 0.5*AU
k_c = np.sqrt(MU_S/(2*rp_c**3))
t_1au = (1 + 1/3)/k_c
ok(f"parabolic comet, r_p = 0.5 AU: perihelion to 1 AU in {t_1au/86400:.1f} days", close(t_1au/86400, 38.8, 0.05))
Mp20 = k_c*20*86400
s20 = np.sqrt(9*Mp20**2/4 + 1)
D20 = np.cbrt(1.5*Mp20 + s20) + np.cbrt(1.5*Mp20 - s20)
f20 = 2*np.arctan(D20)
r20 = 2*rp_c/(1 + np.cos(f20))
print(f"     20 days: M_p = {Mp20:.4f}, D = {D20:.4f}, f = {np.degrees(f20):.2f} deg, r = {r20/AU:.4f} AU")
ok("comet 20 days after perihelion", close(Mp20, 0.6881, 1e-4) and close(D20, 0.6118, 1e-4) and close(np.degrees(f20), 62.91, 0.01) and close(r20/AU, 0.6871, 1e-4))
ok("comet speed at 1 AU = 42.1 km/s (local escape)", close(np.sqrt(2*MU_S/AU), 42.1, 0.05))

# Jupiter flyby question
e_j = 1 + 1e6*36/1.267e8
ok(f"Jupiter flyby: e = {e_j:.3f}, delta = {2*np.degrees(np.arcsin(1/e_j)):.1f} deg",
   close(e_j, 1.284, 0.001) and close(2*np.degrees(np.arcsin(1/e_j)), 102.3, 0.1))
ok(f"Jupiter's orbital speed {np.sqrt(MU_S/(5.2026*AU)):.1f} km/s (13.1)", close(np.sqrt(MU_S/(5.2026*AU)), 13.1, 0.05))

# hyperbolic departure task: v_inf = 3 km/s from 300 km altitude to 925 000 km
rp_h, vinf = R_E + 300, 3.0
a_h = MU_E/vinf**2
e_h = 1 + rp_h*vinf**2/MU_E
Hh = np.arccosh((925000/a_h + 1)/e_h)
Mh = e_h*np.sinh(Hh) - Hh
nh = np.sqrt(MU_E/a_h**3)
print(f"     departure: |a| = {a_h:.0f}, e = {e_h:.4f}, H = {Hh:.4f}, M_h = {Mh:.3f}, t = {Mh/nh:.0f} s = {Mh/nh/86400:.2f} d")
ok("departure hyperbola to 925 000 km takes ~3.1 days",
   close(a_h, 44289, 1) and close(e_h, 1.1506, 1e-4) and close(Mh/nh/86400, 3.11, 0.01))
vp_h = np.sqrt(vinf**2 + 2*MU_E/rp_h)
ok(f"departure periapsis speed {vp_h:.3f} km/s, delta = {2*np.degrees(np.arcsin(1/e_h)):.1f} deg",
   close(vp_h, 11.33, 0.01) and close(2*np.degrees(np.arcsin(1/e_h)), 120.7, 0.1))

# Kepler III: Jupiter's constant smaller by 1/(1 + m_J/M_sun)
ok("T^2/a^3 for Jupiter = 1/(1 + 1/1047) = 0.99905 of the massless value", close(1/(1 + 1/1047.35), 0.99905, 1e-5))

# ======================================================================================= C. integration
print("C. two-body integration vs closed form (mu = 1)")
def rhs(t_, y):
    rr_ = y[:3]
    return np.concatenate([y[3:], -rr_/np.linalg.norm(rr_)**3])
for (aa, ee) in ((1.0, 0.0), (1.0, 0.6), (-1.0, 1.6)):
    rp_ = aa*(1 - ee) if aa > 0 else -aa*(ee - 1)
    vp_ = np.sqrt(2/rp_ - 1/aa)
    y0 = np.array([rp_, 0, 0, 0, vp_, 0])
    tend = 2.0
    sol = solve_ivp(rhs, [0, tend], y0, method="DOP853", rtol=1e-12, atol=1e-13)
    rN = sol.y[:3, -1]
    if aa > 0:
        Ek = kepler(np.sqrt(1/aa**3)*tend, ee) if ee > 0 else np.sqrt(1/aa**3)*tend
        xp, yp = aa*(np.cos(Ek) - ee), aa*np.sqrt(1 - ee**2)*np.sin(Ek)
    else:
        Mh_ = np.sqrt(1/(-aa)**3)*tend
        Hk = np.arcsinh(Mh_/ee)
        for _ in range(60):
            Hk -= (ee*np.sinh(Hk) - Hk - Mh_)/(ee*np.cosh(Hk) - 1)
        xp, yp = -aa*(ee - np.cosh(Hk)), -aa*np.sqrt(ee**2 - 1)*np.sinh(Hk)
    ok(f"a = {aa}, e = {ee}: Kepler's equation reproduces the integration (|dr| = {np.hypot(rN[0]-xp, rN[1]-yp):.1e})",
       np.hypot(rN[0] - xp, rN[1] - yp) < 1e-9)

print("D. further numbers in the tasks")
ok(f"bound check: eps(7000 km, 8 km/s) = {8**2/2 - MU_E/7000:.1f} km^2/s^2 < 0", close(8**2/2 - MU_E/7000, -24.9, 0.05))
ve = np.sqrt(2*MU_E/6678)
ok(f"300 km orbit: v_esc = {ve:.3f}, v_inf at 11.0 km/s = {np.sqrt(11**2 - ve**2):.2f} km/s",
   close(ve, 10.926, 0.0005) and close(np.sqrt(11**2 - ve**2), 1.27, 0.005) and close(11 - ve, 0.074, 0.0005))
ok(f"ISS: v_esc - v_c = {np.sqrt(2)*vc_iss - vc_iss:.2f} km/s (3.2); {86400/T_iss:.1f} orbits/day",
   close(np.sqrt(2)*vc_iss - vc_iss, 3.17, 0.01) and close(86400/T_iss, 15.5, 0.05))
ok(f"Earth: max radial speed e sqrt(mu/p) = {0.0167*np.sqrt(MU_S/(AU*(1 - 0.0167**2))):.3f} km/s",
   close(0.0167*np.sqrt(MU_S/(AU*(1 - 0.0167**2))), 0.50, 0.005))
ok(f"GEO altitude {r_geo - 6378:.0f} km (35786), speed {np.sqrt(MU_E/r_geo):.3f} km/s",
   close(r_geo - 6378.137, 35786, 1) and close(np.sqrt(MU_E/r_geo), 3.07, 0.005))
Es_ = np.array(Es)
errs = np.abs(Es_ - 1.498701134)
print("     Newton errors:", [f"{x:.1e}" for x in errs])
ok("Newton errors 5e-1, 8e-2, 2e-3, 6e-7", close(errs[0], 0.5, 0.01) and close(errs[1], 0.078, 0.001)
   and close(errs[2], 1.5e-3, 1e-4) and close(errs[3], 5.9e-7, 0.2e-7))
ok(f"Earth: 2e = {np.degrees(2*0.0167):.2f} deg = {np.degrees(2*0.0167)/0.9856:.2f} days of solar motion",
   close(np.degrees(2*0.0167), 1.91, 0.01))
ok(f"Jupiter flyby: |dv| = 2 v_inf sin(delta/2) = {2*6/e_j:.2f} km/s", close(2*6/e_j, 9.35, 0.01))
ok(f"comet: sqrt(2 r_p^3/mu) = {1/k_c:.4g} s", close(1/k_c, 2.512e6, 0.001e6))
ok(f"departure: cosh H = {(925000/a_h + 1)/e_h:.2f}, n_h = {nh:.3e}", close((925000/a_h + 1)/e_h, 19.02, 0.01) and close(nh, 6.77e-5, 0.01e-5))
ok(f"Pluto-Charon with q = 0.122: {19596*0.122/1.122:.0f} km", close(19596*0.122/1.122, 2130, 5))
ok(f"Earth-Moon barycentre {384400/82.3:.0f} km", close(384400/82.3, 4671, 1))
ok(f"wedge half-angle = turning half-angle (e = 1.6): asin(1/e) = {np.degrees(np.arcsin(1/1.6)):.2f} = f_inf - 90 = {np.degrees(np.arccos(-1/1.6)) - 90:.2f}",
   close(np.degrees(np.arcsin(1/1.6)), np.degrees(np.arccos(-1/1.6)) - 90, 1e-9))
ok(f"Ex 5: fall into the Sun from 1 AU in {0.5**1.5/2*365.25:.1f} days", close(0.5**1.5/2*365.25, 64.6, 0.05))