"""Verification of every number quoted in Tutorial 1 (tutorial-01.qmd).

Run:  python3 codes/tutorial-01-verify.py
"""
import numpy as np
from scipy.integrate import quad

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

muS, muE, muM = 1.3271e20, 3.9860e14, 4.905e12           # m^3/s^2
aSE, aEM = 1.496e11, 3.844e8                              # m

print("Problem 1")
aEM_ = muE/aEM**2
aSM = muS/aSE**2
ok(f"a_E->M = {aEM_:.3e}, a_S->M = {aSM:.3e}, ratio {aSM/aEM_:.2f}", close(aEM_, 2.70e-3, 0.005e-3) and close(aSM, 5.93e-3, 0.005e-3) and close(aSM/aEM_, 2.20, 0.005))
ok(f"(3.844e8)^2 = {aEM**2:.3e}, (1.496e11)^2 = {aSE**2:.3e}", close(aEM**2, 1.478e17, 0.001e17) and close(aSE**2, 2.238e22, 0.001e22))
tid = 2*muS*aEM/aSE**3
ok(f"tidal 2 mu_S r_EM / r_SE^3 = {tid:.2e} = {100*tid/aEM_:.1f}% of Earth's pull", close(tid, 3.0e-5, 0.05e-5) and close(100*tid/aEM_, 1.1, 0.05))
# exact difference of the Sun's pulls on Moon (at new moon, Moon closer to Sun) and Earth
exact = muS/(aSE - aEM)**2 - muS/aSE**2
ok(f"exact difference at new moon {exact:.3e} (first-order {tid:.3e})", close(exact/tid, 1, 0.01))

print("Problem 2")
aME = muM/aEM**2
ok(f"a_M->E = {aME:.3e}, a_S->E/a_M->E = {aSM/aME:.0f}", close(aME, 3.32e-5, 0.005e-5) and close(aSM/aME, 179, 0.5))
ok(f"ratio of ratios = {(aSM/aME)/(aSM/aEM_):.1f} = mu_E/mu_M = {muE/muM:.1f}", close(muE/muM, 81.3, 0.05))

print("Problem 3")
ok(f"a_SE/a_EM = {aSE/aEM:.0f}", close(aSE/aEM, 389, 0.5))
ok(f"D_S/D_M = {1.392e6/3475:.0f}", close(1.392e6/3475, 400, 1))
ok(f"angular sizes: Sun {np.degrees(1.392e6/1.496e8):.3f} deg, Moon {np.degrees(3475/384400):.3f} deg", close(np.degrees(1.392e6/1.496e8), 0.533, 0.001) and close(np.degrees(3475/384400), 0.518, 0.001))

print("Problem 5")
mu = 398600.44
r = np.array([10000.0, 0, 0]); v = np.array([0, 5.15494, 0])
h = np.cross(r, v); hn = np.linalg.norm(h)
eps = v @ v/2 - mu/10000
a = -mu/(2*eps)
e = np.sqrt(1 + 2*eps*hn**2/mu**2)
ok(f"h = {hn:.1f}, v^2 = {v @ v:.4f}, eps = {eps:.4f}, a = {a:.1f}, e = {e:.4f}", close(hn, 51549.4, 0.05) and close(v @ v, 26.5734, 1e-4) and close(eps, -26.5734, 1e-4) and close(a, 7500, 0.1) and close(e, 1/3, 1e-4))
rp = a*(1 - e); p = hn**2/mu
ok(f"r_p = {rp:.1f}, p = {p:.1f}", close(rp, 5000, 0.5) and close(p, 6666.7, 0.05))
RE5 = 6378.0
cf = (p/RE5 - 1)/e
fimp = np.degrees(np.arccos(cf))
vimp = np.sqrt(mu*(2/RE5 - 1/a))
gimp = np.degrees(np.arctan2(e*np.sin(-fimp*d2r), 1 + e*cf))
print(f"     impact: cos f = {cf:.4f}, f = +-{fimp:.1f} deg, v = {vimp:.2f} km/s, gamma = {gimp:.1f} deg")
ok("impact at f = 277.8 deg (i.e. -82.2), v = 8.48 km/s (the deck had 122 deg and 9.7 km/s)",
   close(cf, 0.1358, 1e-4) and close(fimp, 82.2, 0.05) and close(vimp, 8.48, 0.005))
ok(f"impact flight path angle {gimp:.1f} deg (descending)", close(gimp, -17.5, 0.05))
ok(f"dips {RE5 - rp:.0f} km below the surface", close(RE5 - rp, 1378, 0.5))

print("Problem 7 (Sun)")
muS_km, AU, RS = 1.3271e11, 1.496e8, 6.957e5
PhiS, PhiE = -muS_km/RS, -muS_km/AU
ok(f"Phi(R_S) = {PhiS:.4g}, Phi(1 AU) = {PhiE:.1f}, ratio {PhiS/PhiE:.0f}", close(PhiS, -1.908e5, 0.001e5) and close(PhiE, -887.1, 0.05) and close(PhiS/PhiE, 215, 0.5))
ok(f"fraction of the climb done at 1 AU: {100*(1 - PhiE/PhiS):.1f}%", close(100*(1 - PhiE/PhiS), 99.5, 0.05))
ok(f"v_esc surface {np.sqrt(-2*PhiS):.1f}, at 1 AU {np.sqrt(-2*PhiE):.1f}, v_E {np.sqrt(-PhiE):.1f}",
   close(np.sqrt(-2*PhiS), 617.7, 0.1) and close(np.sqrt(-2*PhiE), 42.1, 0.05) and close(np.sqrt(-PhiE), 29.8, 0.05))
vE = np.sqrt(-PhiE)
ok(f"(sqrt2 - 1) v_E = {(np.sqrt(2) - 1)*vE:.1f}", close((np.sqrt(2) - 1)*vE, 12.3, 0.05))
va = np.sqrt(2*muS_km*RS/(AU*(AU + RS)))
ok(f"graze: v_a = {va:.2f}, dv = {vE - va:.1f}; ratio to escape {(vE - va)/((np.sqrt(2) - 1)*vE):.2f}, radial fall {vE/((np.sqrt(2)-1)*vE):.1f}",
   close(va, 2.87, 0.005) and close(vE - va, 26.9, 0.05) and close(vE/((np.sqrt(2) - 1)*vE), 2.4, 0.05))
r7, v7 = AU, 10.0
eps7 = v7**2/2 - muS_km/r7
h7 = r7*v7
a7 = -muS_km/(2*eps7)
e7 = np.sqrt(1 + 2*eps7*h7**2/muS_km**2)
rp7 = a7*(1 - e7); p7 = h7**2/muS_km
ok(f"(h) eps {eps7:.1f}, h {h7:.3e}, a {a7:.3e} = {a7/AU:.2f} AU, e {e7:.3f}, r_p {rp7:.3e} = {rp7/AU:.3f} AU = {rp7/RS:.1f} R_S",
   close(eps7, -837.1, 0.05) and close(a7, 7.93e7, 0.005e7) and close(a7/AU, 0.53, 0.005) and close(e7, 0.887, 0.0005)
   and close(rp7, 8.94e6, 0.005e6) and close(rp7/AU, 0.060, 0.0005) and close(rp7/RS, 12.8, 0.05))
ok(f"  2 eps h^2/mu^2 = {2*eps7*h7**2/muS_km**2:.4f}, p = {p7:.3e}, p/(1+e) = {p7/(1 + e7):.3e}", close(2*eps7*h7**2/muS_km**2, -0.2128, 1e-4) and close(p7, 1.686e7, 0.001e7) and close(p7/(1 + e7), rp7, 1))
ok(f"  burn {vE - 10:.1f} km/s, two-thirds of v_E ({(vE - 10)/vE:.2f})", close(vE - 10, 19.8, 0.05) and close((vE - 10)/vE, 0.66, 0.01))

print("Problem 8")
mu8, RE8 = 398600.0, 6378.0
rp8, ra8 = RE8 + 400, RE8 + 4000
e8 = (ra8 - rp8)/(ra8 + rp8); a8 = (rp8 + ra8)/2; p8 = 2*rp8*ra8/(rp8 + ra8); h8 = np.sqrt(mu8*p8)
ok(f"e {e8:.5f}, p {p8:.2f}, h {h8:.1f}", close(e8, 0.20984, 1e-5) and close(p8, 8200.29, 0.01) and close(h8, 57172, 1))
vp8, va8 = h8/rp8, h8/ra8
ok(f"v_p {vp8:.3f}, v_a {va8:.3f}, v_p v_a {vp8*va8:.2f} = mu/a {mu8/a8:.2f}", close(vp8, 8.435, 0.001) and close(va8, 5.509, 0.001) and close(vp8*va8, 46.47, 0.01) and close(mu8/a8, 46.47, 0.01))
T8 = 2*np.pi*np.sqrt(a8**3/mu8)
ok(f"a {a8}, T {T8:.1f} s = {T8/60:.2f} min = {T8/3600:.3f} h", close(a8, 8578, 1e-9) and close(T8, 7906.6, 0.1) and close(T8/60, 131.78, 0.01) and close(T8/3600, 2.196, 0.001))
I, _ = quad(lambda f: 1/(1 + e8*np.cos(f)), 0, 2*np.pi)
rbar = p8*I/(2*np.pi)
ok(f"true-anomaly average r = {rbar:.2f} = b = {a8*np.sqrt(1 - e8**2):.2f} = sqrt(r_p r_a) = {np.sqrt(rp8*ra8):.2f}; alt {rbar - RE8:.0f}",
   close(rbar, 8387.02, 0.01) and close(a8*np.sqrt(1 - e8**2), rbar, 1e-6) and close(np.sqrt(rp8*ra8), rbar, 1e-6) and close(rbar - RE8, 2009, 0.5))
tbar = quad(lambda M: a8*(1 - e8*np.cos(__import__('scipy.optimize', fromlist=['brentq']).brentq(lambda E: E - e8*np.sin(E) - M, 0, 2*np.pi))), 0, 2*np.pi)[0]/(2*np.pi)
ok(f"time-averaged r = {tbar:.1f} = a(1 + e^2/2) = {a8*(1 + e8**2/2):.1f}", close(tbar, 8766.9, 0.1) and close(a8*(1 + e8**2/2), 8766.9, 0.1))
cf8 = (p8/rbar - 1)/e8
f8 = np.degrees(np.arccos(cf8))
ok(f"cos f = {cf8:.5f}, f = {f8:.2f}, {360 - f8:.2f}", close(cf8, -0.10610, 1e-5) and close(f8, 96.09, 0.01) and close(360 - f8, 263.91, 0.01))
v8 = np.sqrt(mu8*(2/rbar - 1/a8))
tg = e8*np.sin(f8*d2r)/(1 + e8*np.cos(f8*d2r))
ok(f"v {v8:.3f}, tan gamma {tg:.5f}, gamma {np.degrees(np.arctan(tg)):.2f}", close(v8, 6.970, 0.001) and close(tg, 0.21341, 1e-5) and close(np.degrees(np.arctan(tg)), 12.05, 0.005))
fg = np.degrees(np.arccos(-e8))
ok(f"gamma_max at f = {fg:.2f} ({360 - fg:.2f}), gamma_max = {np.degrees(np.arcsin(e8)):.2f}", close(fg, 102.11, 0.005) and close(360 - fg, 257.89, 0.005) and close(np.degrees(np.arcsin(e8)), 12.11, 0.005))