"""Checks every number and identity quoted in phasing.qmd. Prints FAIL on a mismatch.
Run from anywhere:   python3 codes/phasing-verify.py"""
import numpy as np

fails = 0


def ok(name, x, y, tol):
    global fails
    good = abs(float(x) - float(y)) <= tol
    fails += not good
    print("ok  " if good else "FAIL", name, round(float(x), 6), y)


mu, Re = 398600.4418, 6378.137
vis = lambda r, a: np.sqrt(mu*(2/r - 1/a))
period = lambda a: 2*np.pi*np.sqrt(a**3/mu)

r = 6778.0; dth = np.radians(10.0)
T = period(r); vc = np.sqrt(mu/r)
ok("altitude", r - Re, 400, 0.5); ok("T s", T, 5553.46, 5e-3); ok("T min", T/60, 92.558, 5e-4)
Tp = T*(1 - dth/(2*np.pi)); ap = r*(Tp/T)**(2/3)
ok("T_ph", Tp, 5399.19, 5e-3); ok("a_ph", ap, 6651.89, 5e-3); ok("a_ph fraction", 6778*(35/36)**(2/3), ap, 1e-8); ok("period of the ellipse", period(ap), Tp, 1e-6)
ok("perigee radius", 2*ap - r, 6525.79, 5e-3); ok("perigee alt", 2*ap - r - Re, 147.6, 0.05)
ok("v_c", vc, 7.668636, 5e-7); ok("v_ph", vis(r, ap), 7.595597, 5e-7); total = 2*(vc - vis(r, ap)); ok("total", total, 0.146078, 5e-6); ok("m/s", total*1e3, 146.1, 0.05)

# the two spacecraft really meet, by following them
ang_chaser_end = 0.0                                   # the chaser is back at the burn point after one period of its ellipse
ang_target_end = dth + 2*np.pi*Tp/T
ok("target at the burn point after T_ph", np.degrees((ang_target_end - ang_chaser_end) % (2*np.pi)) % 360, 0, 1e-9)

# the approximate law
ok("small-gap formula k = 1, m/s", 1e3*vc*dth/(3*np.pi), 142, 0.5); ok("within 3 percent", 100*(1 - vc*dth/(3*np.pi)/total), 3, 0.5)
for k_ in (10, 50):
    ap_ = r*(1 - dth/(2*np.pi*k_))**(2/3)
    ok(f"law within 0.5 percent at k={k_}", 2*(vc - vis(r, ap_))/(vc*dth/(3*np.pi*k_)), 1, 5e-3)

# the table
tab = ((1, 6651.89, 147.6, 146.1, 1.50), (2, 6715.09, 274.1, 72.0, 3.04), (3, 6736.10, 316.1, 47.8, 4.59), (5, 6752.87, 349.6, 28.6, 7.67), (10, 6765.44, 374.7, 14.2, 15.4))
res = {}
for k, A, P, D, H in tab:
    ph = T*(1 - dth/(2*np.pi*k)); a_ = r*(ph/T)**(2/3); d = 2*(vc - vis(r, a_))*1e3; res[k] = (d, k*ph/3600)
    ok(f"k={k} a", a_, A, 6e-3); ok(f"k={k} perigee", 2*a_ - r - Re, P, 0.06); ok(f"k={k} dv", d, D, 0.06); ok(f"k={k} h", k*ph/3600, H, 0.06)
ok("saved going to ten", res[1][0] - res[10][0], 132, 0.5); ok("extra hours", res[10][1] - res[1][1], 14, 0.5)

# the 200 km limit used in the text and the figure: k = 1 fails and k = 2 passes for the 10 deg gap
ok("k = 1 below the limit", (2*r*(1 - dth/(2*np.pi))**(2/3) - r - Re) < 200, 1, 0)
ok("k = 2 above the limit", (2*r*(1 - dth/(4*np.pi))**(2/3) - r - Re) > 200, 1, 0)

# the chaser ahead of the target
Ta = T*(1 + dth/(2*np.pi)); aa = r*(Ta/T)**(2/3)
ok("a_ph ahead", aa, 6902.94, 6e-3); ok("fraction", 6778*(37/36)**(2/3), aa, 1e-8); ok("apoapsis", 2*aa - r, 7027.89, 6e-3)
ok("speed", vis(r, aa), 7.73773, 6e-6); ok("v_c, 5 places", vc, 7.66864, 6e-6); ok("total m/s", 2e3*(vis(r, aa) - vc), 138.2, 0.05)

# phasing on an ellipse
a, e, dt = 10000.0, 0.3, 600.0
T = period(a); ok("T", T, 9952.01, 5e-3); aph = a*(1 + dt/T)**(2/3); ok("a_ph", aph, 10397.99, 5e-3)
vals = []
for r_, b, c, dv, tot in ((a*(1 - e), 8.603825, 8.692036, 0.088211, 0.17642), (a*(1 + e), 4.632829, 4.794662, 0.161834, 0.32367)):
    ok("before", vis(r_, a), b, 6e-7); ok("after", vis(r_, aph), c, 6e-7); ok("each", vis(r_, aph) - vis(r_, a), dv, 6e-7); ok("total", 2*(vis(r_, aph) - vis(r_, a)), tot, 6e-6)
    vals.append(2*(vis(r_, aph) - vis(r_, a)))
ok("radii", a*(1 - e), 7000, 1e-9); ok("radii", a*(1 + e), 13000, 1e-9)
ok("periapsis less, percent", 100*(1 - vals[0]/vals[1]), 45, 0.5); ok("ratio", vals[1]/vals[0], 1.835, 5e-4); ok("prediction", (1 + e)/(1 - e), 1.857, 5e-4)
ok("v_p/v_a", vis(a*(1 - e), a)/vis(a*(1 + e), a), (1 + e)/(1 - e), 1e-12); ok("half a period", T/2, 4976, 0.5); ok("saving m/s", 1e3*(vals[1] - vals[0]), 147, 0.5)
ok("the differential at small da", (vis(7000, a + 1e-3) - vis(7000, a))/1e-3, mu/(2*vis(7000, a)*a*a), 1e-8)

# the exercise
Tg, rg = 86164.0989, 42164.0                   # one sidereal day, s
# period set by 2 deg per revolution east over 5 revolutions
ok("sidereal day h", Tg/3600, 23.934, 5e-4)
Tn = Tg*(1 - 2/360); an = (mu*(Tn/(2*np.pi))**2)**(1/3)     # Kepler III with the new period
ok("period h", Tn/3600, 23.80, 5e-3); ok("a", an, 42008, 0.5); ok("perigee radius", 2*an - rg, 41852, 0.5)
d = np.sqrt(mu/rg) - vis(rg, an); ok("each m/s", d*1e3, 5.72, 5e-3); ok("total m/s", 2e3*d, 11.4, 0.05)

print("\nfailures:", fails)
raise SystemExit(1 if fails else 0)