"""Checks every number of lunar-trajectories-bates.qmd.

Part 1: the Earth-Moon constants and the simple (no lunar gravity) trajectories, against a numerical integration.
Part 2: the patched-conic procedure: reproduces the worked example of patched-conics.qmd, then Bate's example, and compares with a numerical
        integration in which the Moon pulls on the vehicle throughout (restricted problem, Moon on a circle).
Part 3: the non-coplanar geometry (spherical trigonometry, checked with explicit vectors) and the exercises.

Run from the sfm-slides directory:   python3 codes/lunar-verify.py
"""
import sys
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
sys.path.insert(0, "codes")
from _lunar import *

D2R = np.pi / 180
RESULTS = []


def check(name, value, stated, tol):
    ok = abs(value - stated) <= tol
    RESULTS.append(ok)
    print(f"{'PASS' if ok else 'FAIL'}  {name:<76} {value:16.8g}   note: {stated}")


def truth(name, cond):
    RESULTS.append(bool(cond))
    print(f"{'PASS' if cond else 'FAIL'}  {name}")


# ------------------------------------------------------------------ Part 1
print("Part 1: the Earth-Moon model and the simple trajectories")
check("mean motion n = sqrt((mu_E + mu_M)/D^3) [1e-6 rad/s]", N_MOON * 1e6, 2.6653, 5e-5)
check("  period [days]", 2 * np.pi / N_MOON / 86400, 27.285, 5e-4)
check("  Moon's speed n D [km/s]", V_MOON, 1.0245, 5e-5)
check("  with the Earth alone, sqrt(mu_E/D) [km/s]  (Bate's 1.018)", np.sqrt(MU_E / D), 1.0183, 5e-5)
check("Laplace radius D (mu_M/mu_E)^(2/5) [km]", R_SOI, 66183, 0.5)
check("  as a fraction of D", R_SOI / D, 0.172, 5e-4)
check("barycentre: D mu_M/(mu_E + mu_M) from the Earth's centre [km]", D * MU_M / (MU_E + MU_M), 4671, 1.0)
r0 = RE + 320.0
check("r0 [km]", r0, 6698.137, 5e-4)
vmin = v_min_to(r0, D)
vesc = np.sqrt(2 * MU_E / r0)
check("least injection speed (apogee at the Moon's distance) [km/s]", vmin, 10.8157, 5e-5)
check("escape speed at r0 [km/s]", vesc, 10.9096, 5e-5)
fm = flight(r0, vmin, 0.0, D)
check("minimum-energy trajectory: flight time [h]", fm["t"] / HOUR, 119.53, 0.005)
check("  eccentricity", fm["e"], (D - r0) / (D + r0), 1e-9); check("  eccentricity (value)", fm["e"], 0.9657, 5e-5)
check("  arrival speed [km/s]", fm["v"], 0.1885, 5e-5); check("  arrival speed = v_min r0/D [km/s]", fm["v"], vmin * r0 / D, 1e-9)
check("  sweep angle [deg]", fm["psi"] / D2R, 180.0, 1e-5)
check("  semi-major axis, a = (r0 + D)/2 [km]", (r0 + D) / 2, 195549.1, 0.05)
check("  half the period, pi sqrt(a^3/mu) [h]", np.pi * np.sqrt(((r0 + D) / 2) ** 3 / MU_E) / HOUR, fm["t"] / HOUR, 1e-6)


def integ(r0_, v0_, g0_, rT):
    y0 = [r0_, 0, v0_ * np.sin(g0_), v0_ * np.cos(g0_)]
    f = lambda t, y: [y[2], y[3], -MU_E * y[0] / np.hypot(y[0], y[1]) ** 3, -MU_E * y[1] / np.hypot(y[0], y[1]) ** 3]
    ev = lambda t, y: np.hypot(y[0], y[1]) - rT
    ev.terminal = True; ev.direction = 1
    s = solve_ivp(f, [0, 1e7], y0, events=ev, rtol=1e-12, atol=1e-9)
    y = s.y_events[0][0]
    return s.t_events[0][0], np.arctan2(y[1], y[0])


worst_t = worst_p = 0
for v_ in (10.82, 10.85, 10.95, 11.0, 11.2):
    t_i, p_i = integ(r0, v_, 0.0, D); fl = flight(r0, v_, 0.0, D)
    worst_t = max(worst_t, abs(t_i - fl["t"])); worst_p = max(worst_p, abs(p_i - fl["psi"]))
check("flight time: formula minus integration, worst of 5 speeds [s]", worst_t, 0.0, 1e-3)
check("sweep angle: formula minus integration, worst of 5 speeds [rad]", worst_p, 0.0, 1e-8)
print("      speed, flight time [h], sweep [deg], arrival speed [km/s], arrival flight-path angle [deg]")
rows = {}
for v_ in (10.82, 10.85, 10.90, 10.95, 11.00, 11.10, 11.20):
    fl = flight(r0, v_, 0.0, D); rows[v_] = fl
    print(f"      {v_:6.2f}  {fl['t']/HOUR:8.2f}  {fl['psi']/D2R:8.2f}  {fl['v']:8.4f}  {fl['g']/D2R:8.2f}   e = {fl['e']:.4f}")
check("v0 = 10.85: flight time [h]", rows[10.85]["t"] / HOUR, 66.28, 0.005); check("v0 = 11.00: flight time [h]", rows[11.00]["t"] / HOUR, 40.71, 0.005)
check("v0 = 11.00: sweep [deg]", rows[11.00]["psi"] / D2R, 158.98, 0.005); check("v0 = 11.20: flight time [h]", rows[11.20]["t"] / HOUR, 30.95, 0.005)
v72 = brentq(lambda v: flight(r0, v, 0.0, D)["t"] / HOUR - 72, vmin + 1e-4, 11.3)
check("injection speed for a 72 h flight [km/s]", v72, 10.8392, 5e-5)
check("  the excess over the minimum [m/s]", (v72 - vmin) * 1000, 23.5, 0.1)
check("  time saved against the minimum-energy flight [h]", fm["t"] / HOUR - 72, 47.5, 0.05)
check("  at escape speed the flight time [h]", flight(r0, vesc + 1e-6, 0.0, D)["t"] / HOUR, 50.71, 0.005)

print("\nMiss at the Moon's orbit caused by an injection error (direct launch; the Moon's travel during the changed flight time included)")


def S(v0_, rr=r0, g0_=0.0):
    fl = flight(rr, v0_, g0_, D); return fl["psi"] - N_MOON * fl["t"], fl["psi"], fl["t"]


h = 1e-5
vc_speed = brentq(lambda v: S(v + h)[0] - S(v - h)[0], 10.84, 11.3)
vc_alt = brentq(lambda v: S(v, r0 + 1e-2)[0] - S(v, r0 - 1e-2)[0], 10.84, 11.3)
check("speed errors cancel at injection speed [km/s]", vc_speed, 10.9345, 5e-4); check("  sweep angle there [deg]", S(vc_speed)[1] / D2R, 162.99, 0.01)
check("altitude errors cancel at [km/s]", vc_alt, 10.9296, 5e-4)
d11 = (S(11.0 + h)[0] - S(11.0 - h)[0]) / (2 * h) * D / 1000
check("at 11.0 km/s: net miss per m/s of speed [km]", d11, -80.7, 0.05)
check("  sweep-angle part [km per m/s]", (S(11.0 + h)[1] - S(11.0 - h)[1]) / (2 * h) * D / 1000, -364.7, 0.05)
check("  Moon's-travel part [km per m/s]", -(S(11.0 + h)[2] - S(11.0 - h)[2]) / (2 * h) * N_MOON * D / 1000, 284.0, 0.05)
dg = (S(11.0, r0, 1e-6)[0] - S(11.0, r0, -1e-6)[0]) / 2e-6 * D * D2R
check("at 11.0 km/s: miss per degree of flight-path angle [km]", dg, -13183, 1.0)


def miss_integrated(v0_, r_, g_):
    t_, p_ = integ(r_, v0_, g_, D); return p_ - N_MOON * t_


for dgd in (0.1, 1.0):
    d_int = (miss_integrated(11.0, r0, dgd * D2R) - miss_integrated(11.0, r0, -dgd * D2R)) / 2 * D
    check(f"  same, by integrating a +-{dgd} deg change [km per {dgd} deg]", d_int, dg * dgd, abs(dg * dgd) * 3e-3)
check("  speed error by integration, +-0.1 m/s [km per m/s]", (miss_integrated(11.0 + 1e-4, r0, 0.0) - miss_integrated(11.0 - 1e-4, r0, 0.0)) / 2e-4 * D / 1000, d11, 0.5)
check("  miss of 1 deg [km] (the book prints about 1300)", abs(dg), 13183, 1.0)

# ------------------------------------------------------------------ Part 2
print("\nPart 2: patched conics")
res = patch(6678.137, 10.84, 0.0, 60 * D2R)
check("patched-conics.qmd example: r1 [km]", res["r1"], 355953, 1.0); check("  eta1 [deg]", res["eta1"] / D2R, 9.27, 0.005)
check("  v1 [km/s]", res["v1"], 0.6087, 5e-5); check("  gamma1 [deg]", res["g1"] / D2R, 70.48, 0.005); check("  t1 [h]", res["t1"] / HOUR, 70.75, 0.005)
check("  v2 [km/s]", res["v2"], 0.9053, 5e-5); check("  |h2| [km^2/s]", abs(res["h2"]), 6371, 1.0); check("  e2", res["e2"], 1.4607, 5e-5)
check("  perilune radius [km]", res["rp"], 3364, 1.0); check("  phase angle at departure [deg]", res["gamma_dep"] / D2R, 125.66, 0.005)
check("  h2 is negative (retrograde about the Moon)", float(res["h2"] < 0), 1.0, 0)

r0b, v0b, lam = r0, 10.846, 30 * D2R
R = patch(r0b, v0b, 0.0, lam)
print("Bate's example: injection at r0 = 6698.137 km, v0 = 10.846 km/s, gamma0 = 0, lambda1 = 30 deg")
eps0 = v0b ** 2 / 2 - MU_E / r0b
check("  energy eps0 [km^2/s^2]", eps0, -0.6913, 5e-5); check("  h0 [km^2/s]", R["h0"], 72648.0, 0.5)
check("  e0", R["e0"], 0.9768, 5e-5); check("  a0 [km]", R["a0"], 288301, 1.0); check("  p0 [km]", R["p0"], 13240.7, 0.05)
check("  r1 = sqrt(D^2 + Rs^2 - 2 D Rs cos lambda1) [km]", R["r1"], 328754, 1.0)
check("  eta1 = asin(Rs sin(lambda1)/r1) [deg]", np.degrees(np.arcsin(R_SOI * np.sin(lam) / R["r1"])), 5.777, 0.0005); check("  eta1 from the patch point [deg]", R["eta1"] / D2R, 5.777, 0.0005)
check("  v1 [km/s]", R["v1"], 1.0209, 5e-5); check("  gamma1 [deg]", R["g1"] / D2R, 77.50, 0.005)
check("  true anomaly at the patch point f1 [deg]", R["f1"] / D2R, 169.28, 0.01)
check("  flight time to the patch point [h]", R["t1"] / HOUR, 50.97, 0.005)
check("  phase angle at departure [deg]", R["gamma_dep"] / D2R, 135.49, 0.005)
check("  v2 [km/s]  (Bate: 1.198)", R["v2"], 1.1983, 5e-5); check("  eps2 [km^2/s^2]  (Bate: 0.6437)", R["eps2"], 0.6439, 5e-5)
check("  h2 [km^2/s]", R["h2"], 8290.9, 0.1); check("  p2 [km]", R["p2"], 14020.5, 0.1); check("  e2", R["e2"], 2.1639, 5e-5)
check("  perilune radius [km]", R["rp"], 4431.3, 0.1); check("  perilune altitude [km]", R["rp"] - RM, 2693.9, 0.1)
check("  aiming angle zeta2 [deg]", R["zeta2"] / D2R, 6.001, 0.0005)
check("  h2 = r_SOI v2 sin(zeta2)", R_SOI * R["v2"] * np.sin(R["zeta2"]), R["h2"], 1e-6)
check("  v_infinity = sqrt(2 eps2) [km/s]", np.sqrt(2 * R["eps2"]), 1.1348, 5e-5)
check("  |a2| = mu_M/(2 eps2) [km]", MU_M / (2 * R["eps2"]), 3807.2, 0.05)
vcirc = np.sqrt(MU_M / R["rp"]); vp = np.sqrt(2 * (R["eps2"] + MU_M / R["rp"]))
check("  speed at perilune [km/s]", vp, 1.8710, 5e-5); check("  circular speed there [km/s]", vcirc, 1.0519, 5e-5)
check("  capture burn to a circular lunar orbit [km/s]", vp - vcirc, 0.8191, 5e-5)
lam_hit = brentq(lambda l: patch(r0b, v0b, 0.0, l)["h2"], 25 * D2R, 50 * D2R)
check("  lambda1 for a dead-centre hit [deg]", lam_hit / D2R, 36.16, 0.005)
lam_lo = brentq(lambda l: patch(r0b, v0b, 0.0, l)["rp"] - RM, 25 * D2R, lam_hit); lam_hi = brentq(lambda l: patch(r0b, v0b, 0.0, l)["rp"] - RM, lam_hit, 50 * D2R)
check("  hits the Moon for lambda1 from [deg]", lam_lo / D2R, 32.72, 0.005); check("  to [deg]", lam_hi / D2R, 39.72, 0.005)
dlam = 0.01 * D2R
check("  d(perilune altitude)/d lambda1 at 30 deg [km per deg]", (patch(r0b, v0b, 0.0, lam + dlam)["rp"] - patch(r0b, v0b, 0.0, lam - dlam)["rp"]) / (2 * dlam) * D2R, -1076, 1.0)

# restricted three-body check: the Moon pulls on the vehicle all the way
def restricted(r0_, v0_, gamma_dep_, tmax=200 * HOUR):
    def f(t, y):
        rm = D * np.array([np.cos(gamma_dep_ + N_MOON * t), np.sin(gamma_dep_ + N_MOON * t)])
        r = y[:2]; d = rm - r
        a = -MU_E * r / np.linalg.norm(r) ** 3 + MU_M * (d / np.linalg.norm(d) ** 3 - rm / np.linalg.norm(rm) ** 3)
        return [y[2], y[3], a[0], a[1]]
    def moon_dist(t, y):
        rm = D * np.array([np.cos(gamma_dep_ + N_MOON * t), np.sin(gamma_dep_ + N_MOON * t)])
        return np.dot(y[:2] - rm, y[2:] - D * N_MOON * np.array([-np.sin(gamma_dep_ + N_MOON * t), np.cos(gamma_dep_ + N_MOON * t)]))
    moon_dist.terminal = True; moon_dist.direction = 1
    s = solve_ivp(f, [0, tmax], [r0_, 0, 0, v0_], events=moon_dist, rtol=1e-11, atol=1e-6, max_step=300)
    t = s.t_events[0][0]; y = s.y_events[0][0]
    rm = D * np.array([np.cos(gamma_dep_ + N_MOON * t), np.sin(gamma_dep_ + N_MOON * t)])
    return t, np.linalg.norm(y[:2] - rm)


t_true, d_true = restricted(r0b, v0b, R["gamma_dep"])
check("restricted three-body: time of closest approach [h]", t_true / HOUR, 64.66, 0.005); check("  closest approach, altitude [km]", d_true - RM, 3192, 1.0)
res60 = patch(6678.137, 10.84, 0.0, 60 * D2R)
t_true60, d_true60 = restricted(6678.137, 10.84, res60["gamma_dep"])
print(f"      for patched-conics.qmd's example: closest approach {d_true60:.0f} km (altitude {d_true60 - RM:.0f}) at {t_true60/HOUR:.2f} h (that note: about 1030 km at 87.0 h)")
check("  patched-conics.qmd example: true perilune altitude [km]", d_true60 - RM, 1030, 60.0)

# ------------------------------------------------------------------ Part 3
print("\nPart 3: non-coplanar geometry")
rng = np.random.default_rng(5)
worst = 0.0
for _ in range(200):
    d0 = rng.uniform(-60, 60) * D2R; a0 = rng.uniform(0, 360) * D2R; b0 = rng.uniform(0, 180) * D2R; pt = rng.uniform(10, 350) * D2R
    P0 = np.array([np.cos(d0) * np.cos(a0), np.cos(d0) * np.sin(a0), np.sin(d0)])
    north = np.array([-np.sin(d0) * np.cos(a0), -np.sin(d0) * np.sin(a0), np.cos(d0)]); east = np.array([-np.sin(a0), np.cos(a0), 0.0])
    tdir = np.cos(b0) * north + np.sin(b0) * east
    P1 = np.cos(pt) * P0 + np.sin(pt) * tdir
    d1 = np.arcsin(P1[2]); a1 = np.arctan2(P1[1], P1[0])
    worst = max(worst, abs(d1 - intercept_declination(d0, pt, b0)))
    da = (a1 - a0 + np.pi) % (2 * np.pi) - np.pi
    if abs(np.cos(d1)) > 1e-3 and abs(np.cos(d0)) > 1e-3:
        worst = max(worst, abs(abs(da) - delta_alpha(d0, d1, pt)) if pt < np.pi else 0)
        bb = launch_azimuth_for(d0, d1, pt)
        if bb is not None and 0.05 < b0 < np.pi - 0.05:
            worst = max(worst, abs(bb - b0))
check("sin d1, cos(da) and cos(beta0) of Bate (7.5-2..4) against explicit vectors, 200 random cases, worst [rad]", worst, 0.0, 1e-9)
d0c = 28.5 * D2R
check("psi_t = 180: d1 = -d0 for every azimuth [deg]", max(abs(intercept_declination(d0c, np.pi, b * D2R) / D2R + 28.5) for b in (10, 40, 90, 115, 170)), 0.0, 1e-9)
check("Bate exercise 7.8: cos(beta0) for d1 = 15, psi_t = 150 (a number above 1: no solution)", (np.sin(15 * D2R) - np.sin(d0c) * np.cos(150 * D2R)) / (np.cos(d0c) * np.sin(150 * D2R)), 1.529, 5e-4)
truth("  so there is no launch azimuth at all", launch_azimuth_for(d0c, 15 * D2R, 150 * D2R) is None)
dmin, dmax = intercept_declination(d0c, 150 * D2R, 0.0) / D2R, intercept_declination(d0c, 150 * D2R, np.pi) / D2R
print("      declination range after a 150 deg sweep from 28.5 N, azimuth 0..180:", dmin, dmax)
dlo = intercept_declination(d0c, 150 * D2R, 40 * D2R) / D2R; dhi = intercept_declination(d0c, 150 * D2R, 115 * D2R) / D2R
check("  with the 40..115 window: declination from [deg]", dhi, -36.79, 0.005); check("  up to [deg]", dlo, -4.39, 0.005)
check("  with all azimuths 0..180: from [deg]", dmax, -58.5, 1e-9); check("  up to [deg]", dmin, 1.5, 1e-9)
# worked example: psi_t = 190, d1 = -25, d0 = 28.5
pt = 190 * D2R; d1 = -25 * D2R
b_ex = launch_azimuth_for(d0c, d1, pt); da_ex = delta_alpha(d0c, d1, pt)
check("worked example: launch azimuth [deg]", b_ex / D2R, 108.1, 0.05); check("  right-ascension difference [deg]", da_ex / D2R, 169.5, 0.05)
check("  azimuth is inside the window 40..115", float(40 <= b_ex / D2R <= 115), 1.0, 0)
# the flight times of that example
fl_ex = flight(r0, 10.9345, 0.0, D)
check("  free-flight sweep at 10.9345 km/s [deg]", fl_ex["psi"] / D2R, 162.99, 0.005); check("  free-flight time [h]", fl_ex["t"] / HOUR, 47.13, 0.005)
n_park = np.sqrt(MU_E / r0 ** 3)
check("  burn-plus-coast arc psi_c = 27 deg: time in the 320 km parking orbit [min]", 27 * D2R / n_park / 60, 6.82, 0.005)
check("  psi_t = psi_c + psi_ff [deg]", 27 + fl_ex["psi"] / D2R, 190.0, 0.015)
# Sriharikota-like site: reachable declinations
d0s = 13.7 * D2R
lo_s = intercept_declination(d0s, 190 * D2R, np.pi) / D2R; hi_s = intercept_declination(d0s, 190 * D2R, 0.0) / D2R
check("13.7 N, psi_t = 190: lowest declination reachable [deg]", min(lo_s, hi_s), -23.7, 0.05); check("  highest [deg]", max(lo_s, hi_s), -3.7, 0.05)
truth("  a Moon at -25 deg cannot be reached at any azimuth", launch_azimuth_for(d0s, -25 * D2R, 190 * D2R) is None)
# exercises
vp_h = np.sqrt(2 * MU_E * D / (RE * (RE + D))); va_h = vp_h * RE / D; vcm = np.sqrt(MU_E / D)
check("Ex 1: Hohmann perigee speed at r = R_E [km/s]", vp_h, 11.0883, 5e-5)
check("Ex 1: apogee speed [km/s]", va_h, 0.1840, 5e-5); check("Ex 1: the extra speed to match the Moon's orbit circle [km/s]", vcm - va_h, 0.8343, 5e-5)
e_soi = ((D - R_SOI) - RE) / ((D - R_SOI) + RE)
check("Ex 2: eccentricity of the orbit that just reaches the SOI from the surface", e_soi, 0.9607, 5e-5)
rp5 = RM + 262.0; v2_5 = 0.5; eps5 = v2_5 ** 2 / 2 - MU_M / R_SOI
vp5 = np.sqrt(2 * (eps5 + MU_M / rp5)); h5 = rp5 * vp5
check("Ex 5: energy at the SOI [km^2/s^2]", eps5, 0.0509, 5e-5); check("  perilune speed [km/s]", vp5, 2.2376, 5e-4)
check("  aiming angle zeta2 [deg]", np.degrees(np.arcsin(h5 / (R_SOI * v2_5))), 7.77, 0.01)
check("  circular speed at perilune [km/s]", np.sqrt(MU_M / rp5), 1.5658, 5e-4); check("  braking burn for the circular orbit [km/s]", vp5 - np.sqrt(MU_M / rp5), 0.672, 1e-3)
# Ex 7: arrival at lambda1 = 0 ...
P7 = np.array([D - R_SOI, 0.0]); v1v = 0.2 * np.array([np.cos(80 * D2R), np.sin(80 * D2R)])[::-1] * np.array([1, 1])
v1v = 0.2 * (np.cos(80 * D2R) * np.array([0.0, 1.0]) + np.sin(80 * D2R) * np.array([1.0, 0.0]))
v2v = v1v - np.array([0.0, V_MOON]); rho7 = P7 - np.array([D, 0.0]); h7 = rho7[0] * v2v[1] - rho7[1] * v2v[0]
check("Ex 7: speed relative to the Moon [km/s]", np.hypot(*v2v), 1.0092, 5e-4); check("  h2 [km^2/s] (positive: direct about the Moon)", h7, 65508, 5.0)
check("  aiming angle: angle between v2 and the line to the Moon's centre [deg]", np.degrees(np.arccos(v2v @ np.array([1.0, 0.0]) / np.hypot(*v2v))), 78.7, 0.05)
check("  energy w.r.t. the Moon [km^2/s^2]", np.hypot(*v2v) ** 2 / 2 - MU_M / R_SOI, 0.4352, 5e-4)
check("Ex 9: speed of a vertical shot to reach D from R_E [km/s]", np.sqrt(2 * MU_E * (1 / RE - 1 / D)), 11.0867, 5e-5)
check("  its speed relative to the Moon at apogee = the Moon's speed [km/s]", V_MOON, 1.0245, 5e-5)

# ---- extra checks quoted in the note
print("\nExtra checks")
h2_closed = R_SOI * (V_MOON * np.cos(lam) - R["v1"] * np.cos(lam + R["eta1"] - R["g1"]))
check("h2 = r_SOI [v_Moon cos(lambda1) - v1 cos(lambda1 + eta1 - gamma1)]  (Bate 7.4-20)", h2_closed, R["h2"], 1e-6)
v2_cos = np.sqrt(R["v1"] ** 2 + V_MOON ** 2 - 2 * R["v1"] * V_MOON * np.cos(R["g1"] - R["eta1"]))
check("v2 by the law of cosines with the angle gamma1 - eta1 [km/s]", v2_cos, R["v2"], 1e-9)
f_in = np.arccos((R["p2"] / R_SOI - 1) / R["e2"])
t_in = time_between(-f_in, 0.0, R["p2"], R["e2"], mu=MU_M)
check("time from the patch point to perilune on the hyperbola [h]", t_in / HOUR, 14.4, 0.05)
check("total flight time of the patched conic [h]", (R["t1"] + t_in) / HOUR, 65.34, 0.005)
check("  against the restricted three-body integration [h]", t_true / HOUR, 64.66, 0.005)
check("patched conic minus true perilune altitude [km]", (R["rp"] - RM) - (d_true - RM), -498, 2.0)
dr_alt = (S(11.0, r0 + 1e-2)[0] - S(11.0, r0 - 1e-2)[0]) / 2e-2 * D
check("at 11.0 km/s: miss per km of injection altitude [km]", dr_alt, -76.4, 0.05)
# first-order estimate of the flight-path-angle effect: f0 ~ 2 gamma0 for e ~ 1
e11 = flight(r0, 11.0, 0.0, D)["e"]
check("flight-path angle: sweep change per degree if f0 = 2 gamma0 [km]", 2 * D2R * D, 13418.1, 0.05)
# eastward right-ascension increase for a sweep beyond 180 degrees
d0_, a0_ = 28.5 * D2R, 0.3
P0 = np.array([np.cos(d0_) * np.cos(a0_), np.cos(d0_) * np.sin(a0_), np.sin(d0_)])
north = np.array([-np.sin(d0_) * np.cos(a0_), -np.sin(d0_) * np.sin(a0_), np.cos(d0_)]); east = np.array([-np.sin(a0_), np.cos(a0_), 0.0])
b_ = b_ex; tdir = np.cos(b_) * north + np.sin(b_) * east
sweep = np.linspace(0, 190 * D2R, 2000)
ras = np.unwrap([np.arctan2(*(np.cos(t) * P0 + np.sin(t) * tdir)[1::-1]) for t in sweep])
check("worked example: eastward RA increase after the 190 deg sweep [deg]", (ras[-1] - ras[0]) / D2R, 190.50, 0.005)
check("  = 360 deg minus the arccos value 169.50", 360 - da_ex / D2R, (ras[-1] - ras[0]) / D2R, 1e-3)
check("  declination reached [deg]", np.degrees(np.arcsin((np.cos(sweep[-1]) * P0 + np.sin(sweep[-1]) * tdir)[2])), -25.0, 1e-6)
check("Ex 7: 13.7 N, psi_t = 150: lowest declination for any azimuth [deg]", np.degrees(np.arcsin(np.sin(d0s) * np.cos(150 * D2R) - np.cos(d0s) * np.sin(150 * D2R))), -43.7, 0.05)
check("  highest [deg]", np.degrees(np.arcsin(np.sin(d0s) * np.cos(150 * D2R) + np.cos(d0s) * np.sin(150 * D2R))), 16.3, 0.05)
check("Ex 9: flight time at v_min + 50 m/s [h]", flight(r0, vmin + 0.05, 0.0, D)["t"] / HOUR, 60.45, 0.005)
check("  time saved [h]", fm["t"] / HOUR - flight(r0, vmin + 0.05, 0.0, D)["t"] / HOUR, 59.07, 0.005)
hh = 0.01 * D2R
check("d gamma_dep/d lambda1 at 30 deg", (patch(r0b, v0b, 0.0, lam + hh)["gamma_dep"] - patch(r0b, v0b, 0.0, lam - hh)["gamma_dep"]) / (2 * hh), -0.2407, 5e-4)
check("escape speed from the Moon at the edge of the sphere of influence [km/s]", np.sqrt(2 * MU_M / R_SOI), 0.3849, 5e-5)
check("flight-path-angle estimate over the exact value (13418/13182)", 13418.1 / 13182.1, 1.018, 5e-4)
check("Moon's distance range from e = 0.0549: perigee [km]", D * (1 - 0.0549), 363296, 1.0); check("  apogee [km]", D * (1 + 0.0549), 405504, 1.0)
check("closing speed of the Moon on the minimum-energy probe [km/s]", V_MOON - fm["v"], 0.836, 5e-4)
check("n_Moon t1 [rad]  (Example 1)", N_MOON * R["t1"], 0.4890, 5e-5); check("  in degrees", N_MOON * R["t1"] / D2R, 28.02, 0.005)
check("  phase angle from the rounded terms 169.28 - 5.78 - 28.02 [deg]", 169.28 - 5.78 - 28.02, 135.48, 1e-9)
check("  the two terms of h2: v_Moon cos(lambda1)", V_MOON * np.cos(lam), 0.8873, 5e-5)
check("  v1 cos(lambda1 + eta1 - gamma1)", R["v1"] * np.cos(lam + R["eta1"] - R["g1"]), 0.7620, 5e-5)
check("  sin(zeta2) = h2/(r_SOI v2)", R["h2"] / (R_SOI * R["v2"]), 0.1045, 5e-5)
check("Try it: 2 mu D [1e11]", 2 * MU_E * D / 1e11, 3.0644, 5e-5); check("  r0 (r0 + D) [1e9]", r0 * (r0 + D) / 1e9, 2.6196, 5e-5)
check("  sqrt(a0^3/mu) [s]", np.sqrt(((r0 + D) / 2) ** 3 / MU_E), 136970, 5.0); check("  pi sqrt(a0^3/mu) [s]", np.pi * np.sqrt(((r0 + D) / 2) ** 3 / MU_E), 430290, 5.0)
print("\n%d checks, %d failed" % (len(RESULTS), RESULTS.count(False)))