"""Checks objective-questions.qmd.

(i)  Structure: 200 questions, four options each, 20/15/15 per level in every topic, and the
     letter in each answer agrees with the answer key.
(ii) Numbers: every numerical answer is recomputed here, and the recomputed value must appear,
     as printed, in the option the page marks as correct.
Run from the sfm-slides directory:   python3 codes/objective-questions-verify.py
"""
import re
import numpy as np
from scipy.optimize import brentq

RESULTS = []


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


# ------------------------------------------------------------------ data (as in the page's data box)
mu, Re, g0 = 398600.0, 6378.0, 9.80665
G, Ms, Me, Mm = 6.674e-11, 1.989e30, 5.972e24, 7.342e22
mus, AU, MJr, Rsun = 1.327e11, 1.496e8, 1047.6, 6.96e5
J2 = 1.08263e-3

# ------------------------------------------------------------------ parse the page
src = open("objective-questions.qmd").read()
blocks = re.findall(r":::\{\.question\}\n(.*?)\n:::\{\.answer\}\n(.*?)\n:::\n:::", src, re.S)
Q = {}
for qtext, atext in blocks:
    m = re.search(r"\*\*(\d+)\.\*\*", qtext)
    if not m:
        continue
    n = int(m.group(1))
    opts = dict(re.findall(r"(?m)^\(([ABCD])\) (.*)$", qtext))
    ans = re.match(r"\*\*\(([ABCD])\)\*\*", atext).group(1)
    Q[n] = dict(opts=opts, ans=ans)
check("200 questions, numbered 1-200", sorted(Q) == list(range(1, 201)))
check("every question has options A-D", all(sorted(q["opts"]) == list("ABCD") for q in Q.values()))
levels = re.findall(r"## Level (\d) .*?\(Q(\d+)--(\d+)\)", src)
sizes = [(int(l), int(b) - int(a) + 1) for l, a, b in levels]
check("each topic has 20 / 15 / 15 questions at levels 1 / 2 / 3",
      sizes == [(1, 20), (2, 15), (3, 15)] * 4)
keytab = dict((int(n), c) for n, c in re.findall(r"(\d+) ([ABCD])", src))
check("answer key agrees with every answer", keytab == {n: q["ans"] for n, q in Q.items()})
counts = {c: sum(q["ans"] == c for q in Q.values()) for c in "ABCD"}
print("     answer letters:", counts)
check("answer letters balanced (each 40-60)", all(40 <= v <= 60 for v in counts.values()))


def right(n):
    return Q[n]["opts"][Q[n]["ans"]]


def num(n, value, shown, x=None, tol=None):
    """value: recomputed here.  shown: the text that must appear in the correct option.
    x: the number that text stands for (parsed from it when it is a plain number).
    tol: allowed |value - x| (default: half a unit in the last printed digit)."""
    plain = shown.replace("\\,", "").replace("^\\circ", "")
    if x is None:
        x = float(plain)
    if tol is None:
        tol = 0.5 * 10.0 ** (-(len(plain.split(".")[1]) if "." in plain else 0))
    ok = shown in right(n) and abs(value - x) <= tol * 1.0001
    check(f"Q{n:<3} {value:.6g}  shown {shown}", ok)


# ------------------------------------------------------------------ n-body
num(9, 9 * 8 / 2, "36")
num(12, Rsun / AU, "4.7\\times10^{-3}", 4.7e-3, 0.05e-3)
num(14, 4 ** 1.5, "8")
num(19, G * Ms * Me / (AU * 1e3) ** 2, "3.5\\times10^{22}", 3.5e22, 0.05e22)
aJ = 5.203 * AU
num(21, aJ / (1 + MJr) / Rsun, "1.07")
num(22, 384400 / (1 + 81.3), "4670", tol=5)
num(23, mus / MJr / aJ ** 2 * 1e3, "2\\times10^{-7}", 2e-7, 0.15e-7)
num(25, np.hypot(1, 2), "2.24")
num(26, 1000 * 999 / 2, "499\\,500")
num(27, 3.33e5 * (384400 / AU) ** 2, "2.2")
num(32, 1 / MJr / 4.203 ** 2, "5\\times10^{-5}", 5e-5, 0.5e-5)
num(36, 1.898e27 / 5.683e26 * np.sqrt(5.203 / 9.537), "2.5")
num(37, 18 - 10, "8")
num(39, 4 + 2 / np.sqrt(2), "4+\\sqrt2", 4 + np.sqrt(2), 1e-12)
vJ = np.sqrt(mus / aJ)
check("Q40  Jupiter's speed 13.06 km/s", abs(vJ - 13.06) < 0.005)
num(40, vJ / MJr * 1e3, "12.5")
num(45, (Mm / Ms) * (AU / 384400) ** 3, "2.2")
d6 = np.radians(6)
num(46, np.degrees(np.arctan2(0.5 * np.sin(d6), 1 + 0.5 * np.cos(d6))), "2.0")
num(48, 2 ** 1.5, "2\\sqrt2", 2 * np.sqrt(2), 1e-12)

# ------------------------------------------------------------------ two-body
num(58, 2 * np.pi * np.sqrt(6678 ** 3 / mu) / 60, "90.5")
num(59, np.sqrt(mu / 7000), "7.55")
num(66, (mu * (23.934 * 3600) ** 2 / (4 * np.pi ** 2)) ** (1 / 3), "42\\,164", tol=1)
check("Q66  distractor 42 241 km is the 24 h radius",
      abs((mu * (24 * 3600) ** 2 / (4 * np.pi ** 2)) ** (1 / 3) - 42241) < 1)
num(69, 1.5 / 0.5, "3")
num(72, np.sqrt(2 * mu * 14000 / (7000 * 21000)), "8.71")
check("Q73  gamma at L is arctan(0.5) = 26.6 deg", abs(np.degrees(np.arctan(0.5)) - 26.57) < 0.01)
num(73, np.degrees(np.arcsin(0.5)), "30^\\circ", tol=1e-9)
num(75, np.sqrt(144 - 2 * mu / 7000), "5.49")
num(76, np.degrees(2 * np.arcsin(0.5)), "60^\\circ", tol=1e-9)
num(77, (1 / np.sqrt(1 + 1 / 81.3) - 1) * 100, "-0.6")
num(78, np.sqrt(1 + 2 * (-20) * 60000 ** 2 / mu ** 2), "0.31")
num(80, 1 / (2 / 8000 - 64 / mu), "11\\,180", tol=5)
num(82, 1 + 7000 * 9 / mu, "1.16")
e = 0.5
E = 2 * np.arctan(np.sqrt((1 - e) / (1 + e)))
num(86, (E - e * np.sin(E)) / (2 * np.pi), "0.098")
num(87, 0.5 - e / np.pi, "0.34")
nu = 7000 * 8.266 ** 2 / mu
check("Q88  r v^2/mu = 1.2", abs(nu - 1.2) < 1e-3)
num(88, np.sqrt(1 - 1.2 * 0.8 * np.cos(np.radians(10)) ** 2), "0.26")
num(90, 1 + 0.6 ** 2 / 2, "1.18")
E1 = 1 + 0.1 * np.sin(1) / (1 - 0.1 * np.cos(1))
num(91, E1, "1.0890")
check("Q91  converged root 1.0886", abs(brentq(lambda x: x - 0.1 * np.sin(x) - 1, 0, 3) - 1.0886) < 5e-5)
num(93, 27.32 / (4 * np.sqrt(2)), "4.8")
num(94, np.degrees(np.arccos(-1 / np.sqrt(2))), "135^\\circ", tol=1e-9)
num(95, np.degrees(np.arccos(-0.8)), "143.1")
num(96, 7000 * np.sqrt(1 + 2 * mu / (7000 * 9)), "25\\,900", tol=50)
num(98, (2 / 3) * 14000 ** 1.5 / np.sqrt(mu) / 60, "29.2")

# ------------------------------------------------------------------ orbital elements
num(116, np.degrees(np.arccos(0)), "90^\\circ", tol=1e-9)
num(121, np.degrees(np.arccos(1 / 2)), "60^\\circ", tol=1e-9)
num(122, np.degrees(np.arctan2(1, -1)), "135^\\circ", tol=1e-9)
num(125, 180 - 98, "82^\\circ", tol=1e-9)
nu = 7000 * 64 / mu
num(129, nu - 1, "0.12")
num(130, 7000 / (2 - nu), "7990", tol=5)
h = np.cross([0, 7000, 0], [0, 0, 7.546])
nv = np.cross([0, 0, 1], h)
check("Q133 h along +I, node along +J", h[0] > 0 and abs(h[1]) + abs(h[2]) < 1e-9 and nv[1] > 0 and abs(nv[0]) < 1e-9)
num(134, np.degrees(np.arcsin(np.sin(np.radians(51.6)) * 0.5)), "23.1")
num(135, 360 - 60, "300^\\circ", tol=1e-9)
num(137, np.degrees(np.arccos(np.cos(np.radians(13.72)) * np.sin(np.radians(135)))), "46.6")
num(140, np.degrees(1e-5 / 1e-4), "5.7")
e = 0.74
num(141, 100 * (1 - (np.arccos(e) - e * np.sqrt(1 - e * e)) / np.pi), "92")
a = Re + 700
n = np.sqrt(mu / a ** 3)
Od = 2 * np.pi / (365.2422 * 86400)
num(142, np.degrees(np.arccos(-2 * Od / (3 * n * J2 * (Re / a) ** 2))), "98.2")
r, v = np.array([7000.0, 0, 0]), np.array([1.0, np.sqrt(mu / 7000), 0])
hv = np.cross(r, v)
ev = np.cross(v, hv) / mu - r / 7000
f = np.degrees(np.arctan2(np.dot(np.cross(ev, r), hv) / np.linalg.norm(hv), np.dot(ev, r))) % 360
check(f"Q144 e = {np.linalg.norm(ev):.4f} (0.13), f = {f:.2f} deg (90)",
      abs(np.linalg.norm(ev) - 0.13) < 0.005 and abs(f - 90) < 0.01 and "0.13" in right(144) and "90" in right(144))
num(149, 50.3 * 100 / 3600, "1.4")

# ------------------------------------------------------------------ orbital maneuvers
num(152, 300 * g0 * np.log(2) / 1e3, "2.04")
r1, r2 = 6678.0, 42164.0
vc1, vc2 = np.sqrt(mu / r1), np.sqrt(mu / r2)
at = (r1 + r2) / 2
vt1, vt2 = np.sqrt(mu * (2 / r1 - 1 / at)), np.sqrt(mu * (2 / r2 - 1 / at))
num(155, vt1 - vc1 + vc2 - vt2, "3.89")
num(156, np.pi * np.sqrt(at ** 3 / mu) / 3600, "5.3")
num(158, 2 * 7.5 * np.sin(np.radians(30)), "7.50")
num(163, (np.sqrt(2) - 1) * vc1, "3.20")
num(171, np.sqrt(2 * mu * 14000 / (7000 * 21000)) - np.sqrt(mu / 7000), "1.17")
num(172, 23.934 * (1 - 45 / 360), "20.94")
check("Q173 GEO speeds 1.608, 3.075", abs(vt2 - 1.608) < 5e-4 and abs(vc2 - 3.075) < 5e-4)
num(173, np.sqrt(vt2 ** 2 + vc2 ** 2 - 2 * vt2 * vc2 * np.cos(np.radians(28.5))), "1.83")
check("Q173 separate burns 2.981", abs(vc2 - vt2 + 2 * vc2 * np.sin(np.radians(14.25)) - 2.981) < 5e-4)
num(174, np.sqrt(9 + 2 * vc1 ** 2) - vc1, "3.60")
num(176, 0.5 * ((1 + 1.524) / 2) ** 1.5 * 365.25, "259")
num(177, 1 / (1 / 365.25 - 1 / 686.98), "780")
num(181, 2 * np.sqrt(mu / 10000) * 0.3 * np.sin(np.radians(30)), "1.89")
num(182, 2 * 7.5 * np.sin(np.radians(30)) * np.sin(np.radians(30)), "3.75")
num(183, 100 * (1 - np.exp(-3000 / (450 * g0))), "49")
p1, p2 = 2 * 7000 * (11 / 12) ** (2 / 3) - 7000, 2 * 7000 * (23 / 24) ** (2 / 3) - 7000
check(f"Q186 perigee k=1 {p1:.0f} < 6378 < k=2 {p2:.0f}", p1 < Re < p2 and right(186).strip() == "$2$")
R = 20.0
g = np.sqrt(2 * R / (1 + R)) - 1 + 1 / np.sqrt(R) - np.sqrt(2 / (R * (1 + R)))
bp = (np.sqrt(2) - 1) * (1 + 1 / np.sqrt(R))
check(f"Q187 Hohmann {g:.3f}, bi-parabolic {bp:.3f}, saving {(g - bp) / g:.1%}",
      abs(g - 0.535) < 5e-4 and bp < g and abs((g - bp) / g - 0.05) < 0.005 and "bi-parabolic" in right(187))
roots = [x.real for x in np.roots([1, -15, -9, -1]) if abs(x.imag) < 1e-12 and x.real > 1]
check("Q188 single root of R^3-15R^2-9R-1 above 1 is 15.58", len(roots) == 1 and abs(roots[0] - 15.5817) < 1e-4)
s = np.sin(np.radians(25))
num(189, s / (1 - 2 * s), "2.73")
check("Q189 optimal 0.794 vc vs direct 0.845 vc",
      abs(2 * (2 * np.sqrt(2 * s * (1 - s)) - 1) - 0.794) < 5e-4 and abs(2 * s - 0.845) < 5e-4)
sm = np.sqrt(0.6)
num(192, 2 * np.degrees(np.arcsin(sm / (1 + sm))), "51.8")
num(193, (2 / 3) * 7000 * np.radians(10) / 86400 * 1e3, "9.4")
vcg = np.sqrt(mu / 42164)
num(197, vt1 - vc1 + vc2 - vt2 + np.sqrt(9 + 2 * vcg ** 2) - vcg, "6.10")
ra, rp = Re + 500, Re + 80
num(198, (np.sqrt(mu / ra) - np.sqrt(2 * mu * rp / (ra * (ra + rp)))) * 1e3, "121")
c = np.hypot(7000, 14000)
num(200, (7000 + 14000 + c) / 4, "9163")

print()
nf = sum(not ok for _, ok in RESULTS)
print(f"{len(RESULTS) - nf} PASS, {nf} FAIL")