"""SatNOGS facts quoted on the page, from the "Listening to the ISS" tab of
The Shape of Distance (src/handoff/The_Shape_of_Distance_PRESENTATION.html,
page 3; recordings of 3 October 2026 at 437.55 / 437.80 MHz).

  python3 satnogs.py   -> data/satnogs_20261003.json
"""
import base64, datetime as dt, json, os, re
import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
PRES = os.path.join(HERE, "..", "src", "handoff", "The_Shape_of_Distance_PRESENTATION.html")
BURST_T = -10.0      # s from the Europe pass epoch: the instant compared across receivers


def page3():
    s = open(PRES, encoding="utf-8").read()
    blobs = re.findall(r"[A-Za-z0-9+/]{20000,}={0,2}", s)
    html = base64.b64decode(blobs[3]).decode("utf-8")
    assert "<title>Listening to the ISS</title>" in html
    js = re.findall(r"<script[^>]*>(.*?)</script>", html, flags=re.S)[0]
    return json.loads(js[js.index("{"): js.rindex("}") + 1])


def main():
    D = page3()
    out = {"source": "The Shape of Distance, Listening to the ISS (https://distance.aethercosmology.com/#doppler)", "passes": {}}
    for key, p in D["passes"].items():
        st = p["stations"]
        dt0 = np.array([s["t0"] - s["pub"]["t_ca"] for s in st])
        V = np.array([s["V"] for s in st]); E = np.array([s["V_err"] for s in st]); w = 1 / E**2
        down = [s["t0"] for s in st]
        out["passes"][key] = {
            "utc": dt.datetime.fromtimestamp(p["t_epoch"], dt.timezone.utc).strftime("%Y-%m-%d %H:%M UTC"),
            "stations": len(st),
            "t0_rms_s": float(np.sqrt((dt0**2).mean())),
            "v_kms": float((V * w).sum() / w.sum() / 1000),
            # point beneath the station at each receiver's closest approach (published orbit)
            "ca": [[s["pub"]["lat"], s["pub"]["lon"], s["t0"]] for s in st],
            "offset_sd_median_hz": float(np.median([s["offset_sd"] for s in st])),
        }
    allV = np.array([s["V"] for p in D["passes"].values() for s in p["stations"]])
    allE = np.array([s["V_err"] for p in D["passes"].values() for s in p["stations"]])
    out["v_kms"] = float((allV / allE**2).sum() / (1 / allE**2).sum() / 1000)
    # one instant of the Europe pass: carrier-corrected shift at every receiver with data within 3 s
    p = D["passes"]["2026-10-03_Europe"]
    shifts = []
    for s in p["stations"]:
        pts = np.array(s["pts"]); t = pts[:, 0]; i = np.searchsorted(t, BURST_T)
        if 0 < i < len(t) and t[i] - t[i - 1] < 6:
            shifts.append((s["name"], float(np.interp(BURST_T, t, pts[:, 1]) - (s["f"] - s["f_nom"])), s["lat"], s["lon"]))
    v = np.array([x[1] for x in shifts])
    lo, hi = min(shifts, key=lambda x: x[1]), max(shifts, key=lambda x: x[1])
    out["burst"] = {
        "utc": dt.datetime.fromtimestamp(p["t_epoch"] + BURST_T, dt.timezone.utc).strftime("%Y-%m-%d %H:%M:%S UTC"),
        "receivers": len(shifts), "up": int((v > 0).sum()), "down": int((v < 0).sum()),
        "min_hz": float(v.min()), "max_hz": float(v.max()), "min_name": lo[0], "max_name": hi[0],
        "lat_range": [min(x[2] for x in shifts), max(x[2] for x in shifts)],
        "lon_range": [min(x[3] for x in shifts), max(x[3] for x in shifts)],
    }
    with open(os.path.join(HERE, "data", "satnogs_20261003.json"), "w") as fh:
        json.dump(out, fh, indent=1)
    print(json.dumps({k: v for k, v in out.items() if k != "passes"}, indent=1))
    for k, q in out["passes"].items():
        print(k, q["utc"], q["stations"], "t0 rms %.2f s" % q["t0_rms_s"], "v %.3f" % q["v_kms"], "offset sd %.0f Hz" % q["offset_sd_median_hz"])


if __name__ == "__main__":
    main()