"""The Sync Test: ISS position and an imitated signal on four maps.

  ../.venv/bin/python build.py   -> data/sync.json, site/index.html

Orbit: data/iss_25544_20261028.tle, fetched 2026-10-29 from
https://celestrak.org/NORAD/elements/gp.php?CATNR=25544&FORMAT=TLE
(epoch 2026-10-28 14:18 UTC), propagated with SGP4 (python sgp4).

For every second of 16 orbits:
  station   SGP4 position -> Earth-fixed (GMST rotation) -> geodetic lat/lon
  map speed speed of the point beneath the station across each flat map
            (models.py: North-pole AE, dual-pole AE on 0N 0E, Mercator with its
            east and west edges joined)
  signal    starts level with the station at each northbound equator
            crossing and moves along the map's track at the measured speed:
            7.07 km/s, the straight-line (level-flight) reading of the SatNOGS
            recordings of 3 Oct 2026 (satnogs.py), which over a plane is also
            the speed of the point beneath. lag = time since the station was
            where the signal is.
On the globe the same recordings fit the published orbit (7.36 km/s at
altitude, x R/r on the ground); check_globe() confirms the signal stays level
there, as it must, since both come from the orbit.
"""
import json, math, os, sys
import numpy as np
from sgp4.api import Satrec, jday

HERE = os.path.dirname(os.path.abspath(__file__))
SRC = os.path.join(HERE, "..", "src")
sys.path.insert(0, SRC)
from models import R  # 6371.0 km sphere shared with the cable survey
from build import land  # Natural Earth 110m, rings of [lon, lat]

TLE = os.path.join(HERE, "data", "iss_25544_20261028.tle")
START = (2026, 10, 28, 15, 0, 0)   # first northbound equator crossing after this
N_ORBITS = 16
MARGIN_S = 2 * 3600                # track after the last orbit, for a signal running ahead
STEP = 10                          # s between samples sent to the page
G0 = 9.80665e-3                    # km/s^2
OMEGA = 7.2921150e-5               # rad/s
WGS_A, WGS_F = 6378.137, 1 / 298.257223563
MIN_POP = 1_000_000
DEDUP_KM = 30                      # merge districts of one city
SATNOGS = os.path.join(HERE, "data", "satnogs_20261003.json")   # satnogs.py
F_UHF = 437.55e6                   # Hz, the recorded downlink
C_KMS = 299_792.458


def gmst(jd):
    t = (jd - 2451545.0) / 36525.0
    s = 67310.54841 + (876600 * 3600 + 8640184.812866) * t + 0.093104 * t**2 - 6.2e-6 * t**3
    return np.radians((s % 86400.0) / 240.0)


def propagate(sat, t0jd, t0fr, secs):
    jd = np.full(secs.shape, t0jd)
    fr = t0fr + secs / 86400.0
    e, r, v = sat.sgp4_array(jd, fr)
    assert not e.any(), "sgp4 error"
    g = gmst(jd + fr)
    c, s = np.cos(g), np.sin(g)
    x = c * r[:, 0] + s * r[:, 1]
    y = -s * r[:, 0] + c * r[:, 1]
    z = r[:, 2]
    vx = c * v[:, 0] + s * v[:, 1] + OMEGA * y
    vy = -s * v[:, 0] + c * v[:, 1] - OMEGA * x
    vz = v[:, 2]
    return np.stack([x, y, z], 1), np.stack([vx, vy, vz], 1)


def geodetic(p):
    x, y, z = p.T
    e2 = WGS_F * (2 - WGS_F)
    lon = np.degrees(np.arctan2(y, x))
    rho = np.hypot(x, y)
    lat = np.arctan2(z, rho * (1 - e2))
    for _ in range(6):
        n = WGS_A / np.sqrt(1 - e2 * np.sin(lat) ** 2)
        h = rho / np.cos(lat) - n
        lat = np.arctan2(z, rho * (1 - e2 * n / (n + h)))
    return np.degrees(lat), lon, h


def unit(lat, lon):
    p, l = np.radians(lat), np.radians(lon)
    return np.stack([np.cos(p) * np.cos(l), np.cos(p) * np.sin(l), np.sin(p)], -1)


def to_latlon(u):
    u = u / np.linalg.norm(u, axis=-1, keepdims=True)
    return np.degrees(np.arcsin(np.clip(u[..., 2], -1, 1))), np.degrees(np.arctan2(u[..., 1], u[..., 0]))


# flat maps, vectorised forms of models.ae_xy / dualpole_xy / mercator_xy (km)
def ae_xy(lat, lon):
    r = R * np.radians(90 - lat)
    t = np.radians(lon)
    return r * np.sin(t), -r * np.cos(t)


def dualpole_xy(lat, lon):
    p, l = np.radians(lat), np.radians(lon)
    c = np.arccos(np.clip(np.cos(p) * np.cos(l), -1, 1))
    k = np.where(c < 1e-9, 1.0, c / np.where(c < 1e-9, 1.0, np.sin(c)))
    return R * k * np.cos(p) * np.sin(l), R * k * np.sin(p)


def mercator_xy(lat, lon):
    return R * np.radians(lon), R * np.log(np.tan(np.pi / 4 + np.radians(np.clip(lat, -85, 85)) / 2))


MAPS = {"ae": ae_xy, "dualpole": dualpole_xy, "mercator": mercator_xy}
WRAP = {"mercator": 2 * math.pi * R}


def step_lengths(f, u, wrap=None):
    """Length on the map of each 1-s step between unit vectors u[i] and u[i+1].
    Steps longer than 20 km on the map (dual-pole rim) are subdivided along the
    great circle until the pieces are short."""
    lat, lon = to_latlon(u)
    x, y = f(lat, lon)
    dx, dy = np.diff(x), np.diff(y)
    if wrap:
        dx = (dx + wrap / 2) % wrap - wrap / 2
    d = np.hypot(dx, dy)
    for i in np.nonzero(d > 20)[0]:
        n = int(min(20000, math.ceil(d[i] / 2)))
        w = np.linspace(0, 1, n + 1)[:, None]
        sub = u[i] * (1 - w) + u[i + 1] * w
        la, lo = to_latlon(sub)
        sx, sy = f(la, lo)
        ddx = np.diff(sx)
        if wrap:
            ddx = (ddx + wrap / 2) % wrap - wrap / 2
        d[i] = np.hypot(ddx, np.diff(sy)).sum()
    return d


def main():
    sn = json.load(open(SATNOGS))
    v_sig = round(sn["v_kms"], 2)     # straight-line reading of the recordings = speed in level flight over a plane
    l1, l2 = [s.strip() for s in open(TLE).read().splitlines()[1:3]]
    sat = Satrec.twoline2rv(l1, l2)
    jd0, fr0 = jday(*START)
    span = int(N_ORBITS * 93 * 60 + 3 * 3600)
    secs = np.arange(span + 1, dtype=float)
    pos, vel = propagate(sat, jd0, fr0, secs)
    lat, lon, alt = geodetic(pos)

    # northbound equator crossings (linear interpolation within the second)
    up = np.nonzero((lat[:-1] < 0) & (lat[1:] >= 0))[0]
    nodes = up + (-lat[up]) / (lat[up + 1] - lat[up])
    assert len(nodes) >= N_ORBITS + 1
    t0 = math.floor(nodes[0])
    nodes = nodes[: N_ORBITS + 1] - t0
    end = int(math.ceil(nodes[-1])) + MARGIN_S
    sl = slice(t0, t0 + end + 1)
    secs, pos, vel, lat, lon, alt = secs[sl] - t0, pos[sl], vel[sl], lat[sl], lon[sl], alt[sl]
    u = unit(lat, lon)
    vmeas = np.linalg.norm(vel, axis=1)                  # globe reading: km/s relative to the ground, at altitude
    orbit_of = np.clip(np.searchsorted(nodes, secs, side="right") - 1, 0, N_ORBITS - 1)
    inorb = secs < nodes[-1]
    spans = [(int(math.ceil(nodes[k])), int(nodes[k + 1])) for k in range(N_ORBITS)]

    # globe: point beneath on the 6371 km sphere (geocentric), and the ground speed the measured speed gives it (x R/r)
    r = np.linalg.norm(pos, axis=1)
    uc = pos / r[:, None]
    vglobe = R * np.arccos(np.clip((uc[:-1] * uc[1:]).sum(1), -1, 1))
    vglobe = np.concatenate([[vglobe[0]], (vglobe[:-1] + vglobe[1:]) / 2, [vglobe[-1]]])
    aglobe = np.zeros_like(vglobe)
    aglobe[5:-5] = (vglobe[10:] - vglobe[:-10]) / 10.0 / G0
    vt = vel - (vel * uc).sum(1)[:, None] * uc
    vg_ref = np.linalg.norm(vt, axis=1) * R / r

    out = {"maps": {}}
    stats = {}
    for name, f in MAPS.items():
        d = step_lengths(f, u, WRAP.get(name))
        S = np.concatenate([[0.0], np.cumsum(d)])                 # map km travelled
        v = np.concatenate([[d[0]], (d[:-1] + d[1:]) / 2, [d[-1]]])
        a = np.zeros_like(v)
        a[5:-5] = (v[10:] - v[:-10]) / 10.0 / G0                    # g, change of speed along the track
        lag, gap = signal(S, secs, nodes, orbit_of, v_sig)
        out["maps"][name] = {"S": S, "v": v, "a": a, "lag": lag, "gap": gap}
        # restarted every minute instead of every orbit
        sync = np.arange(0.0, nodes[-1], 60.0)
        lag60, gap60 = signal(S, secs, sync, np.clip(np.searchsorted(sync, secs, side="right") - 1, 0, len(sync) - 1), v_sig)
        ends = [np.interp(n - 1, secs, lag) for n in nodes[1:]]
        stats[name] = {
            "v_min": float(v[inorb].min()), "v_max": float(v[inorb].max()),
            "a_max": float(np.abs(a[inorb]).max()),
            "ratio_per_orbit": [float(v[a0:b0].max() / v[a0:b0].min()) for a0, b0 in spans],
            "swing_per_orbit": [float(np.abs(np.diff(v[a0:b0])).sum()) for a0, b0 in spans],
            "track_per_orbit": [float(S[b0] - S[a0]) for a0, b0 in spans],
            "lag_end_mean": float(np.mean(ends)), "lag_end_min": float(np.min(ends)), "lag_end_max": float(np.max(ends)),
            "lag_max": float(lag[inorb].max()), "lead_max": float(-lag[inorb].min()),
            "gap_max": float(np.abs(gap[inorb]).max()),
            "lag60_max": float(np.abs(lag60[inorb]).max()), "gap60_max": float(np.abs(gap60[inorb]).max()),
        }
        stats[name].update(one_orbit_after(S, secs, nodes, v_sig))
    # dual-pole rim: closest approach to 0N 180E
    rim = np.degrees(np.arccos(np.clip(-u[inorb][:, 0], -1, 1)))
    stats["dualpole"]["rim_min_deg"] = float(rim.min())

    # cities: station or signal above the horizon (mean altitude), XOR = out of sync
    h_mean = float(alt[inorb].mean())
    lam = math.acos(R / (R + h_mean))
    cities = load_cities()
    cu = unit(np.array([c[1] for c in cities]), np.array([c[2] for c in cities]))
    cl = math.cos(lam)
    idx = np.arange(0, int(nodes[-1]), STEP)
    sees = (u[idx] @ cu.T) >= cl
    worst = np.full(len(idx), 10**9)
    for name in MAPS:
        tau_i = idx - out["maps"][name]["lag"][idx]
        uu = np.stack([np.interp(tau_i, secs, u[:, k]) for k in range(3)], 1)
        bad = sees ^ ((uu @ cu.T) >= cl)
        worst = np.minimum(worst, bad.sum(1))
        stats[name]["cities_now_max"] = int(bad.sum(1).max())

    # opening frame: within the first two orbits, at least 10 min before a restart and with cities out of
    # sync on every map, the moment when the smallest lag over the three maps is largest
    lagmin = np.min([out["maps"][n]["lag"][idx] for n in MAPS], axis=0)
    nxt = nodes[np.minimum(orbit_of[idx] + 1, N_ORBITS)]
    ok = (idx < nodes[2]) & (nxt - idx >= 600) & (worst >= 8)
    t_init = int(idx[np.argmax(np.where(ok, lagmin, -1))])
    stats["globe"] = {"v_min": float(vglobe[inorb].min()), "v_max": float(vglobe[inorb].max()),
                      "a_max": float(np.abs(aglobe[inorb]).max()),
                      "ratio_per_orbit": [float(vglobe[a0:b0].max() / vglobe[a0:b0].min()) for a0, b0 in spans],
                      "swing_per_orbit": [float(np.abs(np.diff(vglobe[a0:b0])).sum()) for a0, b0 in spans],
                      "track_per_orbit": [float(vglobe[a0:b0].sum()) for a0, b0 in spans],
                      "signal_drift_km_max": check_globe(secs, uc, vg_ref, nodes, orbit_of),
                      "ref_mean": float(vg_ref[inorb].mean())}
    stats["measured"] = {"v_sig": v_sig, "v_min": float(vmeas[inorb].min()), "v_max": float(vmeas[inorb].max()),
                         "alt_mean": h_mean, "horizon_deg": math.degrees(lam), "horizon_km": R * lam,
                         "period_min": float(np.diff(nodes).mean() / 60), "signal_per_orbit_km": v_sig * float(np.diff(nodes).mean())}
    stats["recorded"] = recorded_map_speeds(sn, lat, out, inorb)
    stats["one_listener"] = one_listener(u, uc, lat, lon, h_mean)
    stats["satnogs"] = {k: sn[k] for k in ("v_kms", "burst")} | {
        "passes": {k: {x: q[x] for x in ("utc", "stations", "t0_rms_s", "v_kms", "offset_sd_median_hz")} for k, q in sn["passes"].items()}}

    # page data: 10-s samples
    k = np.arange(0, end + 1, STEP)
    epoch = jd0 + fr0 + t0 / 86400.0
    page = {
        "tle": [l1, l2], "epoch_unix": round((epoch - 2440587.5) * 86400.0, 3), "step": STEP,
        "nodes": [round(n, 2) for n in nodes.tolist()], "horizon_deg": round(math.degrees(lam), 3),
        "alt_km": round(h_mean, 1), "t_init": t_init, "v_sig": v_sig, "vg_ref": round(stats["globe"]["ref_mean"], 4),
        "lat": np.round(lat[k], 3).tolist(), "lon": np.round(lon[k], 3).tolist(),
        "vmeas": np.round(vmeas[k] * 1000).astype(int).tolist(),            # m/s
        "vglobe": np.round(vglobe[k] * 1000).astype(int).tolist(),          # m/s
        "aglobe": np.round(aglobe[k] * 1000).astype(int).tolist(),          # milli-g
        "maps": {n: {"v": np.round(m["v"][k] * 1000).astype(int).tolist(),  # m/s
                     "a": np.round(m["a"][k] * 1000).astype(int).tolist(),  # milli-g
                     "S": np.round(m["S"][k], 1).tolist()}                  # map km travelled; the page derives the signal
                 for n, m in out["maps"].items()},
        "cities": [list(c[:4]) for c in cities],
        "stats": stats,
    }
    with open(os.path.join(HERE, "data", "sync.json"), "w") as fh:
        json.dump(page, fh, separators=(",", ":"))
    vals = copy_values(page, stats)
    if not os.path.exists(os.path.join(HERE, "template.html")):
        return
    page["land"] = land()
    tpl = open(os.path.join(HERE, "template.html"), encoding="utf-8").read()
    d3 = open(os.path.join(SRC, "raw", "lib", "d3.min.js"), encoding="utf-8").read()
    assert tpl.count("/*DATA*/null") == 1 and tpl.count("/*D3*/") == 1 and "</script" not in d3
    unused = [k for k in vals if "{{" + k + "}}" not in tpl]
    if unused:
        print("not quoted on the page:", ", ".join(unused))
    for key, val in vals.items():
        tpl = tpl.replace("{{" + key + "}}", val)
    assert "{{" not in tpl, tpl[tpl.index("{{"):][:40]
    html = tpl.replace("/*DATA*/null", json.dumps(page, separators=(",", ":"))).replace("/*D3*/", d3)
    with open(os.path.join(HERE, "site", "index.html"), "w", encoding="utf-8") as fh:
        fh.write(html)


def signal(S, secs, sync, k, v_sig):
    """Source level with the station at each sync time, then moving along the map track at v_sig.
    Returns lag (s since the station was where the signal is) and gap (map km, station ahead > 0)."""
    Ssync = np.interp(sync[k], secs, S)
    sig = Ssync + v_sig * (secs - sync[k])
    return secs - np.interp(sig, S, secs), S - sig


def one_orbit_after(S, secs, nodes, v_sig):
    """Lag and gap one orbit after the signal is lined up, for every start time 1 min apart."""
    P = float(np.diff(nodes).mean())
    starts = np.arange(0.0, nodes[-1] - P, 60.0)
    target = np.interp(starts, secs, S) + v_sig * P
    lag = starts + P - np.interp(target, S, secs)
    gap = np.interp(starts + P, secs, S) - target
    return {"lag1_min": float(lag.min()), "lag1_max": float(lag.max()), "gap1_min": float(gap.min()), "gap1_max": float(gap.max())}


def load_cities():
    """GeoNames cities of MIN_POP or more; places within DEDUP_KM of a larger one are merged into it."""
    rows = []
    for row in open(os.path.join(SRC, "raw", "geonames", "cities15000.txt"), encoding="utf-8"):
        c = row.rstrip("\n").split("\t")
        if int(c[14] or 0) >= MIN_POP:
            rows.append((c[2], round(float(c[4]), 3), round(float(c[5]), 3), c[8], int(c[14])))
    rows.sort(key=lambda c: -c[4])
    kept, ku = [], []
    for c in rows:
        cv = unit(np.array(c[1]), np.array(c[2]))
        if ku and (np.array(ku) @ cv).max() >= math.cos(DEDUP_KM / R):
            continue
        kept.append(c); ku.append(cv)
    return kept


def recorded_map_speeds(sn, lat, out, inorb):
    """Map speed at each SatNOGS receiver's closest approach (3 Oct 2026). On the North-pole AE and Mercator
    maps the speed depends only on latitude and direction (northbound / southbound), so it is read from this
    orbit's track at the same latitude and direction."""
    north = np.gradient(lat) > 0
    res = {}
    for key, p in sn["passes"].items():
        ca = sorted(p["ca"], key=lambda x: x[2])
        nb = ca[-1][0] > ca[0][0] if abs(ca[-1][0] - ca[0][0]) > 1 else None
        row = {"recorded_kms": p["v_kms"], "northbound": nb}
        for name in ("ae", "mercator"):
            v = out["maps"][name]["v"]
            sp = []
            for la, lo, t in ca:
                if nb is None:   # pass near its northern or southern limit: take both branches
                    m = inorb
                else:
                    m = inorb & (north == nb)
                i = np.nonzero(m)[0]
                j = i[np.argmin(np.abs(lat[i] - la))]
                sp.append(float(v[j]))
            row[name] = {"mean": float(np.mean(sp)), "max": float(np.max(sp))}
        res[key] = row
    return res


def one_listener(u, uc, lat, lon, h):
    """The station over the North-pole AE plane, adjusting its own frequency so that one receiver
    (VK4JBE, Brisbane) hears the shift the globe predicts, 2 min before its highest pass in this window.
    Returns the error heard by neighbours at F_UHF."""
    rx = (-27.492, 153.069)
    ball = uc * (R + h)
    dx, dy = ae_xy(lat, lon)
    disc = np.stack([dx, dy, np.full_like(dx, h)], 1)

    def rxpos(la, lo):
        e = unit(np.array(la), np.array(lo))
        q = ae_xy(np.array(la), np.array(lo))
        return e * R, e, np.array([float(q[0]), float(q[1]), 0.0])

    def shifts(la, lo):
        P, e, Q = rxpos(la, lo)
        fb = -F_UHF / C_KMS * np.gradient(np.linalg.norm(ball - P, axis=1))
        fd = -F_UHF / C_KMS * np.gradient(np.linalg.norm(disc - Q, axis=1))
        el = np.degrees(np.arcsin(np.clip(((ball - P) @ e) / np.linalg.norm(ball - P, axis=1), -1, 1)))
        return fb, fd, el

    fb, fd, el = shifts(*rx)
    tca = int(np.argmax(el))
    t = tca - 120
    cor = fb[t] - fd[t]
    res = {"rx": "VK4JBE, Brisbane", "max_el": float(el[tca]), "el": float(el[t]), "correction_hz": float(cor)}
    for dk, brg, nm in ((100, 90, "east100"), (1000, 0, "north1000")):
        p, l, dd, b = math.radians(rx[0]), math.radians(rx[1]), dk / R, math.radians(brg)
        p2 = math.asin(math.sin(p) * math.cos(dd) + math.cos(p) * math.sin(dd) * math.cos(b))
        l2 = l + math.atan2(math.sin(b) * math.sin(dd) * math.cos(p), math.cos(dd) - math.sin(p) * math.sin(p2))
        fb2, fd2, el2 = shifts(math.degrees(p2), math.degrees(l2))
        res[nm] = {"error_hz": float(fd2[t] + cor - fb2[t]), "el": float(el2[t])}
    return res


def mmss(s):
    s = int(round(s))
    return f"{s // 60} min {s % 60:02d} s" if s >= 60 else f"{s} s"


def fmt_v(x):
    return f"{x:,.0f}" if x >= 100 else f"{x:.2f}"


def rng(xs, f):
    lo, hi = f(min(xs)), f(max(xs))
    return lo if lo == hi else f"{lo}–{hi}"


def copy_values(page, st):
    """Numbers quoted in the page text, from this run."""
    import datetime as dt
    utc = lambda t: dt.datetime.fromtimestamp(page["epoch_unix"] + t, dt.timezone.utc).strftime("%Y-%m-%d %H:%M UTC")
    m, ae, dp, me, gl, rec, one, sn = (st[k] for k in ("measured", "ae", "dualpole", "mercator", "globe", "recorded", "one_listener", "satnogs"))
    l1 = page["tle"][0]
    yy, doy = int(l1[18:20]), float(l1[20:32])
    ep = dt.datetime(2000 + yy, 1, 1, tzinfo=dt.timezone.utc) + dt.timedelta(days=doy - 1)
    b = sn["burst"]
    rms = [q["t0_rms_s"] for q in sn["passes"].values()]
    el = lambda km: f"{math.degrees(math.atan2(m['alt_mean'], km)):.1f}"
    v = {
        "V_SIG": f"{m['v_sig']:.2f}", "V_ORBIT_RANGE": f"{m['v_min']:.2f}–{m['v_max']:.2f}",
        "ALT": f"{m['alt_mean']:.0f}", "PERIOD": f"{m['period_min']:.1f}",
        "HORIZON_DEG": f"{m['horizon_deg']:.1f}", "HORIZON_KM": f"{m['horizon_km']:,.0f}",
        "EL_HORIZON": el(m["horizon_km"]), "EL_10000": el(10000),
        "GLOBE_V_RANGE": f"{gl['v_min']:.2f}–{gl['v_max']:.2f}", "GLOBE_DRIFT": f"{gl['signal_drift_km_max']:.2f}",
        "N_CITIES": f"{len(page['cities'])}", "N_ORBITS": f"{len(page['nodes']) - 1}",
        "START_UTC": utc(page["nodes"][0]), "EPOCH_UTC": ep.strftime("%Y-%m-%d %H:%M UTC"),
        "SIG_TRACK": f"{round(m['signal_per_orbit_km'], -2):,.0f}",
        "AE_LAG": f"{ae['lag_end_mean'] / 60:.0f} min", "ME_LAG": f"{me['lag_end_mean'] / 60:.0f} min",
        "DP_LAG": f"{dp['lag_end_min'] / 60:.0f}–{dp['lag_end_max'] / 60:.0f} min",
        "AE_LAG60": mmss(ae["lag60_max"]), "AE_GAP60": f"{round(ae['gap60_max'], -2):,.0f}",
        "SN_RMS": rng(rms, lambda x: f"{x:.1f}"), "SN_N": f"{sum(q['stations'] for q in sn['passes'].values())}",
        "SN_OFFSET": f"{np.median([q['offset_sd_median_hz'] for q in sn['passes'].values()]):.0f}",
        "SN_V": f"{sn['v_kms']:.2f}",
        "BURST_UTC": b["utc"].replace(" UTC", "")[11:] + " UTC on 3 October 2026", "BURST_N": f"{b['receivers']}",
        "BURST_UP": f"{b['up']}", "BURST_DOWN": f"{b['down']}",
        "BURST_MIN": f"−{abs(b['min_hz']) / 1000:.1f}", "BURST_MAX": f"+{b['max_hz'] / 1000:.1f}",
        "REC_AE_SOUTH": rng([rec[k]["ae"]["mean"] for k in ("2026-10-03_SouthAmerica", "2026-10-03_Australia")], lambda x: f"{x:.0f}"),
        "REC_AE_EU": f"{rec['2026-10-03_Europe']['ae']['mean']:.1f}",
        "REC_V_SOUTH": rng([rec[k]["recorded_kms"] for k in ("2026-10-03_SouthAmerica", "2026-10-03_Australia")], lambda x: f"{x:.2f}"),
        "REC_V_EU": f"{rec['2026-10-03_Europe']['recorded_kms']:.2f}",
        "ONE_COR": f"{abs(one['correction_hz']) / 1000:.1f}", "ONE_E100": f"{abs(one['east100']['error_hz']):,.0f}",
        "ONE_N1000": f"{abs(one['north1000']['error_hz']) / 1000:.1f}",
        "DP_RIM": f"{dp['rim_min_deg']:.1f}",
        "DP_VMAX": f"{dp['v_max']:,.0f}",
        "AE_XSIG": f"{ae['v_max'] / m['v_sig']:.1f}",
    }
    for k, s in (("G", gl), ("AE", ae), ("DP", dp), ("ME", me)):
        v[k + "_VR"] = f"{fmt_v(s['v_min'])}–{fmt_v(s['v_max'])}"
        v[k + "_AMAX"] = f"{s['a_max']:,.0f}" if s["a_max"] >= 100 else f"{s['a_max']:.2f}"
        v[k + "_RATIO"] = "×" + rng(s["ratio_per_orbit"], lambda x: f"{x:.2f}" if x < 10 else f"{x:.0f}")
        v[k + "_SWING"] = rng(s["swing_per_orbit"], lambda x: f"{x:,.0f}" if x >= 10 else f"{x:.2f}")
        v[k + "_TRACK"] = rng(s["track_per_orbit"], lambda x: f"{round(x, -2):,.0f}")
        if k != "G":
            v[k + "_LAGEND"] = (mmss(s["lag_end_mean"]) if s["lag_end_max"] - s["lag_end_min"] < 2
                                else f"{mmss(s['lag_end_min'])} – {mmss(s['lag_end_max'])}")
            v[k + "_GAP"] = f"{s['gap_max']:,.0f}"
            v[k + "_LAG1"] = f"{s['lag1_min'] / 60:.0f}–{s['lag1_max'] / 60:.0f} min"
            v[k + "_GAP1"] = f"{round(s['gap1_min'], -2):,.0f}–{round(s['gap1_max'], -2):,.0f}"
            v[k + "_LAG60"] = mmss(s["lag60_max"])
    assert abs(gl["ratio_per_orbit"][0] - 1) < 0.01 and gl["signal_drift_km_max"] < 1
    return v


def check_globe(secs, uc, vg, nodes, orbit_of):
    """Globe: a signal moving at the measured speed scaled to the ground (R/r), started level at each
    crossing, against the station's point beneath on the sphere. Largest separation in km: zero up to
    integration error, since both come from the same orbit."""
    ang = np.arccos(np.clip((uc[:-1] * uc[1:]).sum(1), -1, 1))
    S = np.concatenate([[0.0], np.cumsum(R * ang)])
    M = np.concatenate([[0.0], np.cumsum((vg[:-1] + vg[1:]) / 2)])
    m = secs < nodes[-1]
    sig = np.interp(nodes[orbit_of], secs, S) + M - np.interp(nodes[orbit_of], secs, M)
    return float(np.abs(S - sig)[m].max())


if __name__ == "__main__":
    main()