"""Checks the numbers quoted in ground-tracks.qmd.

The repeat-orbit table is computed in the figure script (fig-gt-repeat.py); this script checks everything
that follows from formulas.

Run from the sfm-slides directory:   python3 codes/ground-tracks-verify.py
"""
import numpy as np

RESULTS = []
MU, RE, J2 = 398600.4418, 6378.137, 1.08263e-3
DAY = 86400.0
wE = 2 * np.pi / 86164.0989                     # the Earth's rotation rate against the stars [rad/s]
Ysd = 365.256363
Ytr = 365.242189


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


print("ground tracks")
aI = RE + 420.0
TI = 2 * np.pi * np.sqrt(aI ** 3 / MU)
check("ISS period [min]", TI / 60, 92.97, 0.01)
check("ISS westward shift per orbit [deg]", np.degrees(wE * TI), 23.3, 0.05)
check("  at the equator [km]", wE * TI * RE, 2594, 2)
check("ISS orbits per day", DAY / TI, 15.49, 0.005)
a650 = RE + 650
T650 = 2 * np.pi * np.sqrt(a650 ** 3 / MU)
check("650 km orbit shift [deg]", np.degrees(wE * T650), 24.5, 0.05)
check("  at the equator [km]", wE * T650 * RE, 2727, 2)
nI = 2 * np.pi / TI; iI = np.radians(51.64)
check("ISS omega_E / n", wE / nI, 0.0647, 5e-5)
check("ISS heading at the equator [deg]", np.degrees(np.arctan2(nI * np.sin(iI), nI * np.cos(iI) - wE)), 54.7, 0.05)
for i_d, exact, approx in ((29, 3.84, 3.67), (10, 0.439, 0.436)):
    u = np.linspace(0, 2 * np.pi, 20001); i = np.radians(i_d)
    lam = np.arctan2(np.cos(i) * np.sin(u), np.cos(u)) - u
    lam = (lam + np.pi) % (2 * np.pi) - np.pi
    check(f"figure 8 at i = {i_d}: exact excursion [deg]", np.degrees(np.abs(lam).max()), exact, 0.006)
    check(f"  approximation i^2/4 [deg]", np.degrees(i ** 2 / 4), approx, 0.006)


def Omdot(a, e, i):
    n = np.sqrt(MU / a ** 3); return -1.5 * n * J2 * (RE / (a * (1 - e * e))) ** 2 * np.cos(i)


sun = np.radians(0.98565) / DAY
for h, istated in ((800, 98.60), (700, 98.19)):
    a = RE + h
    i_ss = np.degrees(np.arccos(-sun / (1.5 * np.sqrt(MU / a ** 3) * J2 * (RE / a) ** 2)))
    check(f"sun-synchronous inclination at {h} km [deg]", i_ss, istated, 0.005)
check("highest latitude of an 800 km SSO [deg]", 180 - 98.60, 81.4, 0.05)
for i_d, stated in ((90.0, -1.97), (51.6, -11.9)):
    rel = (np.degrees(Omdot(RE + (800 if i_d == 90 else 420), 0, np.radians(i_d))) * DAY - 0.98565) * 30 / 15
    check(f"local time of the node drift, i = {i_d} [h per 30 days]", rel, stated, 0.05)
check("ISS node drifts through 24 h in [days]", 360 / (4.95 + 0.9856), 61, 0.5)
aP = RE + 800; TP = 2 * np.pi * np.sqrt(aP ** 3 / MU)
check("polar 800 km: period [min]", TP / 60, 100.9, 0.05)
check("  orbits per day", DAY / TP, 14.3, 0.05)
check("  track spacing [deg]", np.degrees(wE * TP), 25.3, 0.05)
check("  at the equator [km]", wE * TP * RE, 2800, 25)
aM = (MU * (86164.0989 / 2 / (2 * np.pi)) ** 2) ** (1 / 3)
check("Molniya a [km]", aM, 26562, 1)
check("  perigee altitude (e = 0.74) [km]", aM * 0.26 - RE, 528, 1)
check("  apogee altitude [km]", aM * 1.74 - RE, 39840, 2)
Mg = np.linspace(0, 2 * np.pi, 200001); Eg = Mg.copy()
for _ in range(60):
    Eg = Eg - (Eg - 0.74 * np.sin(Eg) - Mg) / (1 - 0.74 * np.cos(Eg))
fg = 2 * np.arctan2(np.sqrt(1.74) * np.sin(Eg / 2), np.sqrt(0.26) * np.cos(Eg / 2))
latg = np.degrees(np.arcsin(np.sin(np.radians(63.4)) * np.sin(np.radians(270) + fg)))
check("  hours above 55 N per 12 h orbit (omega = 270 deg)", (latg > 55).mean() * 86164.0989 / 2 / 3600, 6.1, 0.05)
check("osculating a amplitude 1.5 J2 R^2/a for a 400 km orbit [km]", 1.5 * J2 * RE ** 2 / (RE + 400), 10, 0.5)
check("critical inclination [deg]", np.degrees(np.arccos(1 / np.sqrt(5))), 63.4, 0.05)
check("Landsat grid: 360/233 [deg]", 360 / 233, 1.545, 5e-4)
check("  at the equator [km]", np.radians(360 / 233) * RE, 172.0, 0.1)
check("  sidelap with a 185 km swath [%]", (185 - 172.0) / 185 * 100, 7, 0.1)
check("  orbits per day (233/16)", 233 / 16, 14.56, 0.005)

print(f"\n{sum(RESULTS)} PASS, {len(RESULTS) - sum(RESULTS)} FAIL")