"""Verification of every number and identity quoted in the solutions of endsem-2025.qmd (End-Semester, 8 December 2025).

Q1  (p) the Sun's and the Earth's pulls at the Earth's sphere of influence; (h) two orbits of equal energy that still need a Delta v;
    (i) a plane change by rotation about the line of nodes (a, e and omega unchanged); (g) the energy bought by an impulse
Q2  the aiming radius of a hyperbola
Q3  the six elements of the state vector of Curtis's Example 4.3
Q4  the identity e^2 = 1 + sigma (sigma - 2) cos^2(gamma), sigma = r v^2 / mu
Q5  the bi-elliptic transfer 7000 km -> 105 000 km with an apogee of 210 000 km, against the Hohmann transfer
Q6  the free-flight range of the ballistic missile (canonical units, with a numerical integration)

Run:  python3 courses/sfm/exam/endsem-2025-verify.py
"""
import numpy as np
import sympy as sp
from scipy.integrate import solve_ivp

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

MU = 398600.0                      # km^3/s^2: the paper gives no value; the midsem data table has 3.986e5, Curtis 398 600
MU_E, R_E = 398600.4418, 6378.137


def elements(r, v, mu=MU):
    """The six elements from a state vector (the algorithm of the solution): h, i, Omega, e, omega, theta (degrees), a."""
    r, v = np.asarray(r, float), np.asarray(v, float)
    rn, vn = np.linalg.norm(r), np.linalg.norm(v)
    vr = r @ v/rn
    h = np.cross(r, v); hn = np.linalg.norm(h)
    i = np.arccos(h[2]/hn)
    N = np.cross([0, 0, 1.0], h); Nn = np.linalg.norm(N)
    Om = np.arccos(N[0]/Nn) if N[1] >= 0 else 2*np.pi - np.arccos(N[0]/Nn)
    e = ((vn**2 - mu/rn)*r - rn*vr*v)/mu; en = np.linalg.norm(e)
    om = np.arccos(N @ e/(Nn*en)) if e[2] >= 0 else 2*np.pi - np.arccos(N @ e/(Nn*en))
    th = np.arccos(e @ r/(en*rn)) if vr >= 0 else 2*np.pi - np.arccos(e @ r/(en*rn))
    a = hn**2/(mu*(1 - en**2))
    return dict(h=hn, i=i*r2d, Om=Om*r2d, e=en, om=om*r2d, th=th*r2d, a=a, hvec=h, N=N, evec=e, vr=vr, eps=vn**2/2 - mu/rn, rn=rn, vn=vn)


# =========================================================================================== Q1
print("Q1: some of the blanks")
# (p) the Sun's pull against the Earth's, at the Earth's sphere of influence
MU_S, AU = 1.32712440018e11, 1.495978707e8
rsoi = AU*(MU_E/MU_S)**0.4
g_sun, g_earth = MU_S/AU**2, MU_E/rsoi**2
r_eq = AU*np.sqrt(MU_E/MU_S)
ok(f"r_SOI = {rsoi:,.0f} km; at that distance the Sun pulls with {g_sun:.3e} km/s^2 and the Earth with {g_earth:.3e}: the Sun's pull is {g_sun/g_earth:.1f} times larger",
   close(rsoi, 924647, 5) and close(g_sun, 5.930e-6, 5e-10) and close(g_earth, 4.662e-7, 5e-11) and close(g_sun/g_earth, 12.72, 0.01))
ok(f"the two pulls are equal at r = {r_eq:,.0f} km = {r_eq/rsoi:.2f} r_SOI: from there to r_SOI the Sun pulls harder, inside the sphere", close(r_eq, 259_300, 50) and r_eq < rsoi)
# Laplace's criterion compares the *perturbation* with the main acceleration: tidal Sun on the Earth-centred motion, at r_SOI
q_tide = 2*MU_S*rsoi/AU**3/g_earth
ok(f"the Sun's tidal (differential) acceleration at r_SOI, 2 mu_S r/R^3 = {2*MU_S*rsoi/AU**3:.3e} km/s^2, is {100*q_tide:.1f}% of the Earth's pull there",
   close(q_tide, 0.1575, 5e-4))
# (g) an impulse bought energy v dv + dv^2/2
v_, dv_ = sp.symbols("v dv", positive=True)
ok("impulse along v: d(eps) = (v + dv)^2/2 - v^2/2 = v dv + dv^2/2", sp.simplify(((v_ + dv_)**2/2 - v_**2/2) - (v_*dv_ + dv_**2/2)) == 0)
# (h) a circle of radius a and an ellipse with the same semi-major axis a: same energy, they cross, equal speeds, different directions
a_, e_ = 10000.0, 0.5
p_ = a_*(1 - e_**2)
f_x = np.arccos(-e_)                               # where r = a on the ellipse
r_x = p_/(1 + e_*np.cos(f_x))
ve = np.sqrt(MU/p_)*np.array([e_*np.sin(f_x), 1 + e_*np.cos(f_x)])          # (v_r, v_theta) on the ellipse
vc = np.sqrt(MU/a_)*np.array([0.0, 1.0])
gam = np.arctan2(ve[0], ve[1])
dv_h = np.linalg.norm(ve - vc)
ok(f"e = 0.5: the crossing is at f = {f_x*r2d:.1f} deg, r = {r_x:.1f} km = a; speeds {np.linalg.norm(ve):.4f} and {np.linalg.norm(vc):.4f} km/s; "
   f"gamma = {gam*r2d:.2f} deg = arcsin e; dv = {dv_h:.4f} = 2 v_c sin(gamma/2) = {2*np.linalg.norm(vc)*np.sin(gam/2):.4f} km/s ({dv_h/np.linalg.norm(vc):.4f} v_c)",
   close(r_x, a_, 1e-9) and close(np.linalg.norm(ve), np.linalg.norm(vc), 1e-12) and close(gam*r2d, 30.0, 1e-9) and close(dv_h, 2*np.linalg.norm(vc)*np.sin(gam/2), 1e-12)
   and close(dv_h/np.linalg.norm(vc), 0.5176, 5e-5))
# (i) plane change of 45 deg: rotate the state about the line of nodes
def rotate_about(x, axis, ang):
    axis = np.asarray(axis, float)/np.linalg.norm(axis)
    return x*np.cos(ang) + np.cross(axis, x)*np.sin(ang) + axis*(axis @ x)*(1 - np.cos(ang))


r_n = np.array([9000.0, 0.0, 0.0])                  # a point on the line of nodes (the x axis), the ascending node
# a state with the node exactly on the +x axis: v = (v_r, v_t cos i, v_t sin i)
i0 = 30*d2r
vr_, vt_ = 1.1, 6.0
v_n = np.array([vr_, vt_*np.cos(i0), vt_*np.sin(i0)])
el0 = elements(r_n, v_n)
v_new = rotate_about(v_n, [1, 0, 0], 45*d2r)
dv_vec = v_new - v_n
el1 = elements(r_n, v_new)
ok(f"rotation of v about the line of nodes by 45 deg: i {el0['i']:.2f} -> {el1['i']:.2f} deg, a {el0['a']:.3f} -> {el1['a']:.3f}, e {el0['e']:.5f} -> {el1['e']:.5f}, omega {el0['om']:.3f} -> {el1['om']:.3f}, Omega {el0['Om']:.3f} -> {el1['Om']:.3f}",
   close(el1["i"] - el0["i"], 45.0, 1e-9) and close(el0["a"], el1["a"], 1e-9) and close(el0["e"], el1["e"], 1e-12) and close(el0["om"], el1["om"], 1e-9) and close(el0["Om"], el1["Om"], 1e-9))
vt_node = el0["h"]/np.linalg.norm(r_n)
ok(f"dv = {np.linalg.norm(dv_vec):.5f} km/s = 2 (h/r) sin 22.5 deg = {2*vt_node*np.sin(22.5*d2r):.5f} (the radial part {vr_} km/s is untouched; h/r = {vt_node:.3f} km/s); circular orbit: {2*np.sin(22.5*d2r):.4f} v",
   close(np.linalg.norm(dv_vec), 2*vt_node*np.sin(22.5*d2r), 1e-12) and close(2*np.sin(22.5*d2r), 0.7654, 5e-5))

# =========================================================================================== Q2
print("Q2: the aiming radius")
for a_h, e_h in ((20000.0, 1.5), (35000.0, 3.0), (8000.0, 1.1)):
    p_h = a_h*(e_h**2 - 1)
    th_inf = np.arccos(-1/e_h)
    # the perpendicular distance from the focus to the asymptote, from the geometry: centre at a e from the focus, the asymptote at beta = arccos(1/e)
    beta = np.arccos(1/e_h)
    Delta_geo = a_h*e_h*np.sin(beta)
    # the limit  Delta = lim r sin(theta_inf - theta)
    th = th_inf - 1e-7
    Delta_lim = p_h/(1 + e_h*np.cos(th))*np.sin(th_inf - th)
    h_h = np.sqrt(MU*p_h); vinf = np.sqrt(MU/a_h)
    ok(f"a = {a_h:,.0f}, e = {e_h}: a sqrt(e^2-1) = {a_h*np.sqrt(e_h**2 - 1):.4f}; geometry (a e sin beta) {Delta_geo:.4f}; limit of r sin(theta_inf - theta) {Delta_lim:.4f}; h/v_inf {h_h/vinf:.4f}; turn angle {2*np.arcsin(1/e_h)*r2d:.3f} deg, theta_inf - 90 = {(th_inf - np.pi/2)*r2d:.3f}",
       close(Delta_geo, a_h*np.sqrt(e_h**2 - 1), 1e-8) and close(Delta_lim, a_h*np.sqrt(e_h**2 - 1), 1e-2*a_h*1e-4) and close(h_h/vinf, a_h*np.sqrt(e_h**2 - 1), 1e-7)
       and close(2*np.arcsin(1/e_h), 2*(th_inf - np.pi/2), 1e-12))

# =========================================================================================== Q3
print("Q3: the six elements of r = (-6045, -3490, 2500) km, v = (-3.457, 6.618, 2.533) km/s")
r3, v3 = np.array([-6045.0, -3490.0, 2500.0]), np.array([-3.457, 6.618, 2.533])
E3 = elements(r3, v3)
print(f"      |r| = {E3['rn']:.1f} km, |v| = {E3['vn']:.4f} km/s, v_r = {E3['vr']:.4f} km/s, h vector = {E3['hvec']}, N = {E3['N']}, e vector = {E3['evec']}")
ok(f"h = {E3['h']:.0f} km^2/s, i = {E3['i']:.2f}, Omega = {E3['Om']:.2f}, e = {E3['e']:.4f}, omega = {E3['om']:.2f}, theta = {E3['th']:.2f} deg; a = {E3['a']:.0f} km (Curtis: 58 310, 153.2, 255.3, 0.1712, 20.07, 28.45)",
   close(E3["h"], 58310, 5) and close(E3["i"], 153.2, 0.05) and close(E3["Om"], 255.3, 0.05) and close(E3["e"], 0.1712, 5e-5) and close(E3["om"], 20.07, 0.01) and close(E3["th"], 28.45, 0.01))
ok(f"the same with mu = 3.986e5 and 398 600.4418: i, Omega, omega, theta unchanged to 0.01 deg, e = {elements(r3, v3, 3.986e5)['e']:.4f}, {elements(r3, v3, MU_E)['e']:.4f}",
   all(close(elements(r3, v3, m)[k], E3[k], 0.01) for m in (3.986e5, MU_E) for k in ("i", "Om", "om", "th")) and all(close(elements(r3, v3, m)["e"], E3["e"], 1e-4) for m in (3.986e5, MU_E)))
rp3, ra3 = E3["a"]*(1 - E3["e"]), E3["a"]*(1 + E3["e"])
T3 = 2*np.pi*np.sqrt(E3["a"]**3/MU)
ok(f"quadrants: N_y = {E3['N'][1]:.0f} < 0 so Omega in the third/fourth; e_z = {E3['evec'][2]:.4f} > 0 so omega < 180; v_r = {E3['vr']:.3f} > 0 so theta < 180. "
   f"Extras: eps = {E3['eps']:.3f} km^2/s^2, a = {E3['a']:.0f} km, r_p = {rp3:.0f} km, r_a = {ra3:.0f} km, T = {T3/3600:.3f} h; i > 90: retrograde",
   E3["N"][1] < 0 and E3["evec"][2] > 0 and E3["vr"] > 0 and close(E3["a"], 8788, 2) and E3["i"] > 90)
Nn3 = np.linalg.norm(E3["N"])
ok(f"the quoted intermediates: |r| = {E3['rn']:.2f} km, |v| = {E3['vn']:.4f} km/s, r.v = {r3 @ v3:.3f} km^2/s, v_r = {E3['vr']:.4f} km/s, v^2 - mu/r = {E3['vn']**2 - MU/E3['rn']:.4f} km^2/s^2; "
   f"h = ({E3['hvec'][0]:.2f}, {E3['hvec'][1]:.3f}, {E3['hvec'][2]:.2f}), N = ({E3['N'][0]:.3f}, {E3['N'][1]:.2f}, 0), |N| = {Nn3:.1f}; e = ({E3['evec'][0]:.5f}, {E3['evec'][1]:.5f}, {E3['evec'][2]:.5f}); "
   f"cos i = {E3['hvec'][2]/E3['h']:.4f}, cos Omega = {E3['N'][0]/Nn3:.4f}, cos omega = {E3['N'] @ E3['evec']/(Nn3*E3['e']):.4f}, cos theta = {E3['evec'] @ r3/(E3['e']*E3['rn']):.4f}",
   close(E3["rn"], 7414.32, 5e-3) and close(E3["vn"], 7.8845, 5e-5) and close(r3 @ v3, 4133.245, 5e-3) and close(E3["vr"], 0.5575, 5e-5) and close(E3["vn"]**2 - MU/E3["rn"], 8.4040, 5e-5)
   and close(E3["hvec"][0], -25385.17, 5e-3) and close(E3["hvec"][1], 6669.485, 5e-4) and close(E3["hvec"][2], -52070.74, 5e-3) and close(Nn3, 26246.7, 0.05)
   and close(E3["evec"][0], -0.09160, 5e-6) and close(E3["evec"][1], -0.14221, 5e-6) and close(E3["evec"][2], 0.02644, 5e-6)
   and close(E3["hvec"][2]/E3["h"], -0.8930, 5e-5) and close(E3["N"][0]/Nn3, -0.2541, 5e-5) and close(E3["N"] @ E3["evec"]/(Nn3*E3["e"]), 0.9393, 5e-5) and close(E3["evec"] @ r3/(E3["e"]*E3["rn"]), 0.8793, 5e-5))

# =========================================================================================== Q4
print("Q4: e = sqrt(1 + sigma (sigma - 2) cos^2 gamma) with sigma = r v^2 / mu")
r_s, v_s, g_s, mu_s = sp.symbols("r v gamma mu", positive=True)
h_s = r_s*v_s*sp.cos(g_s)
eps_s = v_s**2/2 - mu_s/r_s
e2_s = 1 + 2*eps_s*h_s**2/mu_s**2
sig = r_s*v_s**2/mu_s
ok("e^2 = 1 + 2 eps h^2/mu^2 with h = r v cos(gamma) equals 1 + sigma (sigma - 2) cos^2(gamma), sigma = r v^2/mu", sp.simplify(e2_s - (1 + sig*(sig - 2)*sp.cos(g_s)**2)) == 0)
rng = np.random.default_rng(1)
good = True
for _ in range(6):
    r_t = rng.uniform(7000, 40000)
    v_t = rng.uniform(2, 11)
    g_t = rng.uniform(-80, 80)*d2r
    rv = np.array([r_t, 0, 0]); vv = v_t*np.array([np.sin(g_t), np.cos(g_t), 0])
    e_vec = np.linalg.norm(((vv @ vv - MU/r_t)*rv - (rv @ vv)*vv)/MU)          # the eccentricity vector, directly
    sg = r_t*v_t**2/MU
    good = good and close(e_vec, np.sqrt(1 + sg*(sg - 2)*np.cos(g_t)**2), 1e-9)
ok("six random states: the eccentricity of the state vector equals the formula", good)
ok("checks: sigma = 1, gamma = 0 gives e = 0 (circle); sigma = 2 gives e = 1 for every gamma (parabola); gamma = 0 gives e = |sigma - 1| (apse)",
   sp.simplify(e2_s.subs({v_s: sp.sqrt(mu_s/r_s), g_s: 0})) == 0 and sp.simplify(sp.sqrt(e2_s.subs(v_s, sp.sqrt(2*mu_s/r_s))) - 1) == 0)

# =========================================================================================== Q5
print("Q5: bi-elliptic 7000 -> 105 000 km via 210 000 km, against Hohmann (mu = 398 600)")
rA, rC, rB = 7000.0, 105000.0, 210000.0
vis = lambda r, a: np.sqrt(MU*(2/r - 1/a))
vc = lambda r: np.sqrt(MU/r)
a1, a2, aH = (rA + rB)/2, (rB + rC)/2, (rA + rC)/2
dv1 = vis(rA, a1) - vc(rA)
dv2 = vis(rB, a2) - vis(rB, a1)
dv3 = vc(rC) - vis(rC, a2)                          # negative: a brake
T = lambda a: 2*np.pi*np.sqrt(a**3/MU)
t_bi = T(a1)/2 + T(a2)/2
dvH1 = vis(rA, aH) - vc(rA); dvH2 = vc(rC) - vis(rC, aH)
t_H = T(aH)/2
print(f"      a1 = {a1:,.0f}, a2 = {a2:,.0f}, a_H = {aH:,.0f} km; v_cA = {vc(rA):.4f}, v_A1 = {vis(rA, a1):.4f}, v_B1 = {vis(rB, a1):.4f}, v_B2 = {vis(rB, a2):.4f}, v_C2 = {vis(rC, a2):.4f}, v_cC = {vc(rC):.4f} km/s")
ok(f"bi-elliptic: dv_A = {dv1:.4f} (prograde), dv_B = {dv2:.4f} (prograde), dv_C = {dv3:.4f} (retrograde: v_C2 > v_c); total {dv1 + dv2 + abs(dv3):.4f} km/s",
   close(dv1, 2.9521, 5e-5) and close(dv2, 0.7750, 5e-5) and close(dv3, -0.3014, 5e-5) and close(dv1 + dv2 + abs(dv3), 4.0285, 5e-5))
ok(f"  flight times: half of ellipse 1 {T(a1)/2/3600:.2f} h + half of ellipse 2 {T(a2)/2/3600:.2f} h = {t_bi/3600:.2f} h = {t_bi/86400:.2f} d",
   close(T(a1)/2/3600, 49.40, 5e-3) and close(T(a2)/2/3600, 86.40, 5e-3) and close(t_bi/3600, 135.80, 5e-3))
ok(f"Hohmann: dv_A = {dvH1:.4f}, dv_C = {dvH2:.4f}; total {dvH1 + dvH2:.4f} km/s in {t_H/3600:.2f} h",
   close(dvH1, 2.7868, 5e-5) and close(dvH2, 1.2595, 5e-5) and close(dvH1 + dvH2, 4.0463, 5e-5) and close(t_H/3600, 18.32, 5e-3))
ok(f"saving {1000*(dvH1 + dvH2 - dv1 - dv2 - abs(dv3)):.1f} m/s ({100*(1 - (dv1 + dv2 + abs(dv3))/(dvH1 + dvH2)):.2f}%) for {t_bi/t_H:.2f} times the time; R = r_C/r_A = {rC/rA:.0f}, r_B/r_A = {rB/rA:.0f}",
   close(1000*(dvH1 + dvH2 - dv1 - dv2 - abs(dv3)), 17.8, 0.05) and close(t_bi/t_H, 7.41, 5e-3) and rB/rA > 15.58 and rC/rA > 11.94)
ok(f"Hohmann speeds: v_A = {vis(rA, aH):.4f} (circular {vc(rA):.4f}), v_D = {vis(rC, aH):.4f} (circular {vc(rC):.4f}) km/s", close(vis(rA, aH), 10.3329, 5e-5) and close(vis(rC, aH), 0.6889, 5e-5) and close(vc(rA), 7.5460, 5e-5) and close(vc(rC), 1.9484, 5e-5))
from scipy.optimize import brentq
tot_bi = lambda rb_: (vis(rA, (rA + rb_)/2) - vc(rA)) + (vis(rb_, (rb_ + rC)/2) - vis(rb_, (rA + rb_)/2)) + abs(vc(rC) - vis(rC, (rb_ + rC)/2))
rb_star = brentq(lambda rb_: tot_bi(rb_) - (dvH1 + dvH2), 106000, 600000)
ok(f"at R = 15 the bi-elliptic transfer beats the Hohmann one only if its apogee exceeds {rb_star:,.0f} km = {rb_star/rA:.2f} r_A (here 210 000 km = 30 r_A)", close(rb_star, 127332, 5) and rB > rb_star)
for lab, val in (("dv_A", dv1), ("dv_B", dv2), ("dv_C", dv3), ("Hohmann dv_A", dvH1), ("Hohmann dv_C", dvH2)):
    print(f"      {lab} = {val:.5f}")
print(f"      bi-elliptic total {dv1 + dv2 + abs(dv3):.5f}, Hohmann total {dvH1 + dvH2:.5f}, times {t_bi/3600:.4f} h and {t_H/3600:.4f} h; ratio {t_bi/t_H:.3f}")

# =========================================================================================== Q6
print("Q6: the free-flight range, canonical units (mu = 1 DU^3/TU^2)")
rbo, vbo, ra6 = 1.0 + 1/5, 2/3, 1.0 + 0.5
mu6 = 1.0
eps6 = vbo**2/2 - mu6/rbo
a6 = -mu6/(2*eps6)
e6 = ra6/a6 - 1
p6 = a6*(1 - e6**2)
cosf = (p6/rbo - 1)/e6
fbo = np.arccos(cosf)
psi = 2*(np.pi - fbo)
Q6 = rbo*vbo**2/mu6
cg2 = (1 - e6**2)/(Q6*(2 - Q6))*(-1)*(-1) if False else (e6**2 - 1)/(Q6*(Q6 - 2))
gam6 = np.arccos(np.sqrt(cg2))
tan_half = Q6*np.sin(gam6)*np.cos(gam6)/(1 - Q6*np.cos(gam6)**2)
print(f"      eps = {eps6:.6f}, a = {a6:.6f} DU, e = {e6:.6f}, p = {p6:.6f}, r_p = {a6*(1 - e6):.6f} DU, cos f_bo = {cosf:.4f}, f_bo = {fbo*r2d:.4f} deg")
ok(f"a = {a6:.5f} DU, e = {e6:.5f} = 5/6, p = {p6:.5f} = 1/4, r_p = {a6*(1 - e6):.5f} DU (< 1: below the surface), cos f_bo = {cosf:.4f}",
   close(a6, 9/11, 1e-12) and close(e6, 5/6, 1e-12) and close(p6, 0.25, 1e-12) and close(cosf, -0.95, 1e-12))
ok(f"f_bo = {fbo*r2d:.3f} deg, free-flight range angle psi = 2 (pi - f_bo) = {psi*r2d:.3f} deg = {psi:.5f} rad; range = {psi:.5f} DU = {psi*6378:.1f} km (6378 km/DU)",
   close(fbo*r2d, 161.805, 5e-4) and close(psi*r2d, 36.390, 5e-4) and close(psi, 0.63512, 5e-6) and close(psi*6378, 4050.8, 0.1))
ok(f"cross-check with the Q4 identity: Q = {Q6:.4f}, cos^2 gamma = {cg2:.4f} -> gamma_bo = {gam6*r2d:.2f} deg, and tan(psi/2) = Q sin gamma cos gamma/(1 - Q cos^2 gamma) = {tan_half:.5f} gives psi = {2*np.arctan(tan_half)*r2d:.3f} deg",
   close(gam6*r2d, 51.32, 5e-3) and close(2*np.arctan(tan_half), psi, 1e-12))
# numerical integration of the two-body motion from the burnout state until the radius returns to r_bo
def rhs(t, y):
    r = np.hypot(y[0], y[1])
    return [y[2], y[3], -mu6*y[0]/r**3, -mu6*y[1]/r**3]


y0 = [rbo, 0.0, vbo*np.sin(gam6), vbo*np.cos(gam6)]            # position on the x axis; gamma measured from the local horizontal
ev = lambda t, y: np.hypot(y[0], y[1]) - rbo
ev.terminal, ev.direction = True, -1
sol6 = solve_ivp(rhs, [0, 20], y0, events=ev, rtol=1e-12, atol=1e-12, dense_output=True)
ang_end = np.arctan2(sol6.y_events[0][0][1], sol6.y_events[0][0][0])
rmax = max(np.hypot(*sol6.sol(t)[:2]) for t in np.linspace(0, sol6.t_events[0][0], 4000))
ok(f"integration: the vehicle returns to r = {rbo} DU after {sol6.t_events[0][0]:.4f} TU = {sol6.t_events[0][0]*806.78/60:.2f} min, having swept {ang_end*r2d:.3f} deg; apogee {rmax:.5f} DU",
   close(ang_end*r2d, psi*r2d, 2e-4) and close(rmax, ra6, 1e-6))
# flight time from the Kepler equation, for the record
E_bo = 2*np.arctan2(np.sqrt(1 - e6)*np.sin(fbo/2), np.sqrt(1 + e6)*np.cos(fbo/2))
t_ff = 2*np.sqrt(a6**3/mu6)*((np.pi - e6*np.sin(np.pi)) - (E_bo - e6*np.sin(E_bo)))
ok(f"flight time from Kepler's equation: {t_ff:.4f} TU = {t_ff*806.78/60:.2f} min (TU = sqrt(R^3/mu) = {np.sqrt(6378.137**3/MU_E):.2f} s)",
   close(t_ff, sol6.t_events[0][0], 1e-8))
ok("1 DU/TU = sqrt(mu/R) = 7.905 km/s", close(np.sqrt(MU_E/6378.0), 7.905, 5e-4))

# =========================================================================================== extra checks of the claims made in the text
print("Extra: numbers that appear in the solutions and are not tied to one question's chain")
# the Sun's tidal acceleration at r_SOI, quoted as 7.3e-8 km/s^2
ok(f"tidal acceleration 2 mu_S r_SOI/R^3 = {2*MU_S*rsoi/AU**3:.3e} km/s^2 (quoted 7.3e-8)", close(2*MU_S*rsoi/AU**3, 7.3e-8, 0.05e-8))
# Q1(h), second example: two ellipses with the same a and e, apse lines turned through eta: dv = 2 sqrt(mu/p) e sin(eta/2)
for e_h, eta_deg in ((0.3, 40.0), (0.6, 75.0)):
    a_h = 12000.0; p_h = a_h*(1 - e_h**2); eta = eta_deg*d2r
    th_x = eta/2                                       # r1 = r2 where e cos(th) = e cos(th - eta): th = eta/2
    def vel(th, om):                                    # inertial velocity on the ellipse with apse line at angle om, at inertial polar angle th
        f = th - om; k = np.sqrt(MU/p_h)
        vr, vt = k*e_h*np.sin(f), k*(1 + e_h*np.cos(f))
        return np.array([vr*np.cos(th) - vt*np.sin(th), vr*np.sin(th) + vt*np.cos(th)])
    r1 = p_h/(1 + e_h*np.cos(th_x - 0)); r2 = p_h/(1 + e_h*np.cos(th_x - eta))
    dv_rot = np.linalg.norm(vel(th_x, eta) - vel(th_x, 0))
    ok(f"two ellipses a = {a_h:.0f} km, e = {e_h}, apse lines {eta_deg:.0f} deg apart: they meet at r = {r1:.1f} = {r2:.1f} km; dv = {dv_rot:.5f} km/s = 2 sqrt(mu/p) e sin(eta/2) = {2*np.sqrt(MU/p_h)*e_h*np.sin(eta/2):.5f}",
       close(r1, r2, 1e-9) and close(dv_rot, 2*np.sqrt(MU/p_h)*e_h*np.sin(eta/2), 1e-12))
# Q5: the radius ratios at which the bi-elliptic transfer can beat the Hohmann one (units: mu = 1, r_A = 1, R = r_C/r_A, beta = r_B/r_A)
visn = lambda r, a: np.sqrt(2/r - 1/a)
dvH_n = lambda R: (visn(1, (1 + R)/2) - 1) + (1/np.sqrt(R) - visn(R, (1 + R)/2))
dvbi_n = lambda R, b: (visn(1, (1 + b)/2) - 1) + (visn(b, (b + R)/2) - visn(b, (1 + b)/2)) + abs(1/np.sqrt(R) - visn(R, (b + R)/2))
gap_inf = lambda R: dvbi_n(R, 1e12) - dvH_n(R)                           # r_B -> infinity: the best the bi-elliptic transfer can do
R1 = brentq(gap_inf, 5, 14)
gap_near = lambda R: dvbi_n(R, R*(1 + 1e-6)) - dvH_n(R)                  # r_B just above r_C: the bi-elliptic transfer is almost the Hohmann one
R2 = brentq(lambda R: gap_near(R)/(1e-6*R), 13, 20)
ok(f"bi-elliptic can beat Hohmann for R > {R1:.2f} (r_B -> infinity); for R > {R2:.2f} it beats it for every r_B > r_C (quoted 11.94 and 15.58)", close(R1, 11.94, 0.005) and close(R2, 15.58, 0.005))
bs = np.logspace(np.log10(1.001), 3, 400)
ok("R = 11.5: Hohmann better for every r_B; R = 13: bi-elliptic better only for r_B above a limit; R = 16: bi-elliptic better for every r_B > r_C",
   all(dvbi_n(11.5, b*11.5) > dvH_n(11.5) for b in bs) and any(dvbi_n(13, b*13) < dvH_n(13) for b in bs) and any(dvbi_n(13, b*13) > dvH_n(13) for b in bs)
   and all(dvbi_n(16, b*16) < dvH_n(16) for b in bs))
# does the quoted precision depend on mu? Q5 and Q6 with mu = 398 600 (used), 3.986e5 (identical) and 398 600.4418
def q5_strings(mu):
    v_ = lambda r, a: np.sqrt(mu*(2/r - 1/a)); c_ = lambda r: np.sqrt(mu/r)
    d1, d2, d3 = v_(rA, a1) - c_(rA), v_(rB, a2) - v_(rB, a1), c_(rC) - v_(rC, a2)
    h1, h2 = v_(rA, aH) - c_(rA), c_(rC) - v_(rC, aH)
    t = lambda a: 2*np.pi*np.sqrt(a**3/mu)
    vals = [c_(rA), v_(rA, a1), v_(rB, a1), v_(rB, a2), v_(rC, a2), c_(rC), v_(rA, aH), v_(rC, aH), d1, d2, abs(d3), d1 + d2 + abs(d3), h1, h2, h1 + h2]
    return [f"{x:.4f}" for x in vals] + [f"{(t(a1) + t(a2))/7200:.2f}", f"{t(aH)/7200:.2f}", f"{1000*(h1 + h2 - d1 - d2 - abs(d3)):.1f}", f"{(t(a1) + t(a2))/t(aH):.2f}"]
base = q5_strings(398600.0)
alt = q5_strings(398600.4418)
ok(f"Q5 with mu = 3.986e5 is identical to 398 600; with 398 600.4418 every quoted digit is the same except v_c(7000 km): {base[0]} -> {alt[0]} (7.54605 is on a rounding boundary). Quoted: {base}",
   q5_strings(3.986e5) == base and alt[1:] == base[1:] and base[0] == "7.5460" and alt[0] == "7.5461")
tu = lambda mu: f"{np.sqrt(6378.0**3/mu):.1f}"
ok(f"Q6: 1 TU = {tu(398600.0)} s (6378 km) for all three values of mu: {tu(3.986e5)}, {tu(398600.4418)}", tu(398600.0) == tu(3.986e5) == tu(398600.4418) == "806.8")
print("done")