#!/usr/bin/env python3
"""
process_gps.py - Post-processing + diagnostics pipeline for u-blox ZED-X20P captures.

Pipeline (Mode 2 of the GPS project):
    1. inventory   - parse the .ubx (pyubx2): RXM-RAWX/SFRBX, NAV-PVT/SAT/SIG, MON-RF
    2. rinex       - convert to RINEX 3.04 obs+nav (RTKLIB demo5 convbin)
    3. single      - standalone single-point solution (rnx2rtkp -p 0), ~1 m, offline-safe
    4. ppk         - PPK vs auto-fetched NGS CORS base (cm-level when fixed)
    5. ppp         - optional (--ppp) PPP vs open IGS orbit/clock products (BKG/ESA)
    6. track       - track.csv (decimal-degree lat/lon @ --interval) + track.kml
    7. diagnostics - diagnostics.json + report.md (constellation, iono, error, timing, RF)

Outputs land in processed/<capture-stem>/.  External stages (rinex/single/ppk/ppp)
are idempotent: outputs are reused if present, so the script can be re-run later to
pick up CORS/IGS data that had not been published yet.  Run setup.sh first.

Motion policy (buoy deployment): processing is ALWAYS kinematic by default - static
is never inferred from the data, because a slowly drifting platform is
indistinguishable from measurement noise over short windows. Pass --static only
when the antenna is independently known to have been fixed.

Usage:
    venv/bin/python process_gps.py <capture.ubx> [--interval N] [--static]
                                   [--no-ppk] [--ppp] [--base SSSS] [--max-base-km KM]
                                   [--force-stage {rinex,single,ppk,ppp,track}]
"""

import argparse
import gzip
import json
import math
import re
import shutil
import subprocess
import sys
import urllib.error
import urllib.request
from collections import Counter
from datetime import datetime, timedelta, timezone
from pathlib import Path

import numpy as np
import pandas as pd
from pyubx2 import UBXReader, UBX_PROTOCOL

SCRIPT_DIR = Path(__file__).resolve().parent
TOOLS = SCRIPT_DIR / "tools"
BIN = TOOLS / "bin"
ANTEX_DIR = TOOLS / "antex"
CACHE = TOOLS / "cache"
# combine_rinex.sh lives in notes/ since the repo reorganization; accept either spot
COMBINE_SH = next((p for p in (SCRIPT_DIR / "notes" / "combine_rinex.sh",
                               SCRIPT_DIR / "combine_rinex.sh") if p.is_file()),
                  SCRIPT_DIR / "notes" / "combine_rinex.sh")

GPS_EPOCH = datetime(1980, 1, 6, tzinfo=timezone.utc)
CORS_MIRRORS = [
    "https://geodesy.noaa.gov/corsdata",
    "https://noaa-cors-pds.s3.amazonaws.com",
]
IGS_PRODUCT_MIRRORS = [  # (base url, filename prefix regex)
    ("https://igs.bkg.bund.de/root_ftp/IGS/products", r"IGS0OPS(FIN|RAP)"),
    ("http://navigation-office.esa.int/products/gnss-products", r"ESA0OPS(FIN|RAP)"),
]

GNSS_NAME = {0: "GPS", 1: "SBAS", 2: "Galileo", 3: "BeiDou", 5: "QZSS", 6: "GLONASS", 7: "NavIC"}
QUALITY_NAME = {1: "fix", 2: "float", 3: "sbas", 4: "dgps", 5: "single", 6: "ppp"}
IONO_MODEL = {0: "none/default", 1: "Klobuchar (GPS)", 2: "SBAS", 3: "Klobuchar (BDS)",
              4: "NTCM (GAL)", 8: "dual-frequency"}
CORR_SOURCE = {0: "none", 1: "SBAS", 2: "BDS", 3: "RTCM2", 4: "RTCM3 OSR", 5: "RTCM3 SSR",
               6: "QZSS SLAS", 7: "SPARTN", 9: "CLAS", 10: "LPP OSR", 11: "LPP SSR", 12: "GAL HAS"}
ANT_STATUS = {0: "init", 1: "unknown", 2: "OK", 3: "short", 4: "open"}

# Carrier frequency [Hz] by (gnssId, sigId) - u-blox F9/X20 signal identifiers.
SIG_FREQ = {
    (0, 0): ("L1C/A", 1575.42e6), (0, 3): ("L2CL", 1227.60e6), (0, 4): ("L2CM", 1227.60e6),
    (0, 6): ("L5I", 1176.45e6), (0, 7): ("L5Q", 1176.45e6),
    (1, 0): ("L1C/A", 1575.42e6),
    (2, 0): ("E1C", 1575.42e6), (2, 1): ("E1B", 1575.42e6),
    (2, 3): ("E5aI", 1176.45e6), (2, 4): ("E5aQ", 1176.45e6),
    (2, 5): ("E5bI", 1207.14e6), (2, 6): ("E5bQ", 1207.14e6),
    (2, 8): ("E6B", 1278.75e6), (2, 9): ("E6C", 1278.75e6),
    (3, 0): ("B1I D1", 1561.098e6), (3, 1): ("B1I D2", 1561.098e6),
    (3, 2): ("B2I D1", 1207.14e6), (3, 3): ("B2I D2", 1207.14e6),
    (3, 4): ("B3I", 1268.52e6),
    (3, 5): ("B1Cp", 1575.42e6), (3, 6): ("B1Cd", 1575.42e6),
    (3, 7): ("B2ap", 1176.45e6), (3, 8): ("B2ad", 1176.45e6),
    (5, 0): ("L1C/A", 1575.42e6), (5, 1): ("L1S", 1575.42e6),
    (5, 4): ("L2CM", 1227.60e6), (5, 5): ("L2CL", 1227.60e6),
    (5, 8): ("L5I", 1176.45e6), (5, 9): ("L5Q", 1176.45e6),
    (7, 0): ("L5A", 1176.45e6),
}


def sig_info(gnss, sig, freq_slot):
    """(signal name, carrier freq Hz) - GLONASS FDMA needs the frequency slot."""
    if gnss == 6:
        k = freq_slot - 7
        if sig == 0:
            return ("L1OF", 1602.0e6 + k * 562.5e3)
        if sig == 2:
            return ("L2OF", 1246.0e6 + k * 437.5e3)
        return (f"GLO sig{sig}", 0.0)
    return SIG_FREQ.get((gnss, sig), (f"sig{sig}", 0.0))


def band_of(freq_hz):
    if freq_hz > 1.5e9:
        return "L1"
    if freq_hz > 1.26e9:
        return "L6/E6"
    if freq_hz > 1.2e9:
        return "L2"
    if freq_hz > 0:
        return "L5"
    return "?"


def log(msg):
    print(f"[process_gps] {msg}", flush=True)


def fetch(url, dest=None, timeout=120):
    """Download url; return bytes (and write dest if given). None on failure."""
    try:
        req = urllib.request.Request(url, headers={"User-Agent": "process_gps.py"})
        with urllib.request.urlopen(req, timeout=timeout) as r:
            data = r.read()
        if dest is not None:
            Path(dest).write_bytes(data)
        return data
    except (urllib.error.URLError, urllib.error.HTTPError, OSError):
        return None


def run(cmd, **kw):
    """subprocess.run with logging; returns CompletedProcess."""
    log("  $ " + " ".join(str(c) for c in cmd))
    return subprocess.run([str(c) for c in cmd], capture_output=True, text=True, **kw)


def gps_to_utc(week, tow, leap_s):
    return GPS_EPOCH + timedelta(weeks=week, seconds=tow - leap_s)


def haversine_km(lat1, lon1, lat2, lon2):
    p = math.pi / 180
    a = (math.sin((lat2 - lat1) * p / 2) ** 2
         + math.cos(lat1 * p) * math.cos(lat2 * p) * math.sin((lon2 - lon1) * p / 2) ** 2)
    return 12742 * math.asin(math.sqrt(a))


def llh_to_enu(lat, lon, h, lat0, lon0, h0):
    """Small-extent geodetic -> local ENU (m) about (lat0, lon0, h0). Arrays OK."""
    a, f = 6378137.0, 1 / 298.257223563
    e2 = f * (2 - f)
    phi = math.radians(lat0)
    n_rad = a / math.sqrt(1 - e2 * math.sin(phi) ** 2)            # prime vertical
    m_rad = a * (1 - e2) / (1 - e2 * math.sin(phi) ** 2) ** 1.5   # meridian
    d_n = np.radians(np.asarray(lat) - lat0) * m_rad
    d_e = np.radians(np.asarray(lon) - lon0) * n_rad * math.cos(phi)
    d_u = np.asarray(h) - h0
    return d_n, d_e, d_u


# ----------------------------------------------------------------------------
# Stage 1: UBX inventory
# ----------------------------------------------------------------------------

def parse_ubx(ubx_path):
    """Single streaming pass over the capture; returns dict of DataFrames + summary."""
    log(f"parsing {ubx_path.name} with pyubx2 ...")
    counts = Counter()
    rawx_hdr, rawx_meas = [], []
    pvt, sat, sig_rows, rf = [], [], [], []
    sfrbx_gnss = Counter()

    with open(ubx_path, "rb") as f:
        ubr = UBXReader(f, protfilter=UBX_PROTOCOL, quitonerror=0)
        for _, m in ubr:
            if m is None:
                continue
            ident = m.identity
            counts[ident] += 1
            if ident == "RXM-RAWX":
                ei = len(rawx_hdr)
                rawx_hdr.append((m.week, m.rcvTow, m.leapS, m.leapSec, m.numMeas))
                for i in range(1, m.numMeas + 1):
                    s = f"_{i:02d}"
                    gnss = getattr(m, "gnssId" + s)
                    sid = getattr(m, "sigId" + s)
                    fslot = getattr(m, "freqId" + s)
                    name, freq = sig_info(gnss, sid, fslot)
                    rawx_meas.append((ei, gnss, getattr(m, "svId" + s), sid, name, freq,
                                      getattr(m, "cno" + s), getattr(m, "locktime" + s),
                                      getattr(m, "prStd" + s), getattr(m, "cpStd" + s),
                                      getattr(m, "prValid" + s), getattr(m, "cpValid" + s),
                                      getattr(m, "halfCyc" + s), getattr(m, "prMes" + s)))
            elif ident == "RXM-SFRBX":
                sfrbx_gnss[m.gnssId] += 1
            elif ident == "NAV-PVT":
                pvt.append((datetime(m.year, m.month, m.day, m.hour, m.min, m.second,
                                     tzinfo=timezone.utc),
                            m.fixType, m.carrSoln, m.diffSoln, m.numSV,
                            m.lat, m.lon, m.height / 1000.0, m.hMSL / 1000.0,
                            m.hAcc / 1000.0, m.vAcc / 1000.0, m.tAcc, m.pDOP,
                            m.gSpeed / 1000.0))
            elif ident == "NAV-SAT":
                for i in range(1, m.numSvs + 1):
                    s = f"_{i:02d}"
                    sat.append((getattr(m, "gnssId" + s), getattr(m, "svId" + s),
                                getattr(m, "cno" + s), getattr(m, "elev" + s),
                                getattr(m, "azim" + s), getattr(m, "prRes" + s),
                                getattr(m, "qualityInd" + s), getattr(m, "svUsed" + s),
                                getattr(m, "health" + s)))
            elif ident == "NAV-SIG":
                for i in range(1, m.numSigs + 1):
                    s = f"_{i:02d}"
                    sig_rows.append((getattr(m, "gnssId" + s), getattr(m, "svId" + s),
                                     getattr(m, "sigId" + s), getattr(m, "cno" + s),
                                     getattr(m, "qualityInd" + s), getattr(m, "corrSource" + s),
                                     getattr(m, "ionoModel" + s), getattr(m, "prUsed" + s),
                                     getattr(m, "prRes" + s)))
            elif ident == "MON-RF":
                for i in range(1, m.nBlocks + 1):
                    s = f"_{i:02d}"
                    rf.append((getattr(m, "blockId" + s), getattr(m, "jammingState" + s),
                               getattr(m, "antStatus" + s), getattr(m, "antPower" + s),
                               getattr(m, "noisePerMS" + s), getattr(m, "agcCnt" + s),
                               getattr(m, "jamInd" + s)))

    data = {
        "counts": counts,
        "sfrbx_gnss": sfrbx_gnss,
        "rawx_hdr": pd.DataFrame(rawx_hdr, columns=["week", "rcvTow", "leapS", "leapSecValid", "numMeas"]),
        "rawx": pd.DataFrame(rawx_meas, columns=["epoch", "gnss", "sv", "sig", "signame",
                                                 "freq", "cno", "locktime", "prStd", "cpStd",
                                                 "prValid", "cpValid", "halfCyc", "prMes"]),
        "pvt": pd.DataFrame(pvt, columns=["utc", "fixType", "carrSoln", "diffSoln", "numSV",
                                          "lat", "lon", "height", "hMSL", "hAcc", "vAcc",
                                          "tAcc_ns", "pDOP", "gSpeed"]),
        "sat": pd.DataFrame(sat, columns=["gnss", "sv", "cno", "elev", "azim", "prRes",
                                          "quality", "used", "health"]),
        "sig": pd.DataFrame(sig_rows, columns=["gnss", "sv", "sig", "cno", "quality",
                                               "corrSource", "ionoModel", "prUsed", "prRes"]),
        "rf": pd.DataFrame(rf, columns=["blockId", "jammingState", "antStatus", "antPower",
                                        "noisePerMS", "agcCnt", "jamInd"]),
    }

    # session window in UTC from RAWX (receiver time is GPS-aligned)
    hdr = data["rawx_hdr"]
    if len(hdr):
        leap = int(hdr["leapS"].iloc[0])
        data["t_start"] = gps_to_utc(int(hdr["week"].iloc[0]), float(hdr["rcvTow"].iloc[0]), leap)
        data["t_end"] = gps_to_utc(int(hdr["week"].iloc[-1]), float(hdr["rcvTow"].iloc[-1]), leap)
    # uc2x sidecar metadata (JSON header at start of the u-center 2 index file)
    data["meta"] = {}
    uc2x = ubx_path.with_suffix(".uc2x")
    if uc2x.is_file():
        blob = uc2x.read_bytes()[:4096]
        m = re.search(rb'\{"time".*?\}', blob)
        if m:
            try:
                data["meta"] = json.loads(m.group(0))
            except json.JSONDecodeError:
                pass
    return data


# ----------------------------------------------------------------------------
# Stage 2: RINEX conversion
# ----------------------------------------------------------------------------

def stage_rinex(ubx_path, outdir, n_rawx_epochs, warnings):
    obs, nav = outdir / "rover.obs", outdir / "rover.nav"
    if obs.is_file() and nav.is_file():
        log("rinex: reusing existing rover.obs / rover.nav")
    else:
        r = run([BIN / "convbin", "-r", "ubx", "-v", "3.04", "-od", "-os", "-oi", "-ot", "-ol",
                 "-f", "4", "-o", obs, "-n", nav, ubx_path])
        if not obs.is_file() or not nav.is_file():
            raise RuntimeError(f"convbin failed:\n{r.stderr[-2000:]}")
    if (TOOLS / "STOCK_RTKLIB").exists():
        warnings.append("STOCK RTKLIB in use: BeiDou and L5/E5 observations are MISSING "
                        "from the RINEX output (~30% of this receiver's measurements). "
                        "Re-run setup.sh where the demo5 build can succeed.")

    header, sys_lines, iono = [], [], {}
    with open(obs, errors="replace") as f:
        for line in f:
            header.append(line)
            if "SYS / # / OBS TYPES" in line and line[0] != " ":
                sys_lines.append(line[0])
            if "END OF HEADER" in line:
                break
    n_epochs = sum(1 for line in open(obs, errors="replace") if line.startswith(">"))
    with open(nav, errors="replace") as f:
        for line in f:
            if "IONOSPHERIC CORR" in line:
                iono[line[:4].strip()] = [float(x.replace("D", "E"))
                                          for x in line[5:60].split()]
            if "END OF HEADER" in line:
                break
    nav_counts = Counter(line[0] for line in open(nav, errors="replace")
                         if re.match(r"^[GRECJSI][0-9]{2} ", line))

    if n_epochs != n_rawx_epochs:
        warnings.append(f"RINEX epoch count ({n_epochs}) != RAWX epoch count ({n_rawx_epochs})")
    if "C" not in sys_lines:
        warnings.append("BeiDou missing from RINEX obs (stock convbin limitation?)")
    log(f"rinex: {n_epochs} epochs, systems {sys_lines}, "
        f"nav ephemerides {dict(nav_counts)}, iono keys {list(iono)}")
    return {"obs": obs, "nav": nav, "n_epochs": n_epochs, "systems": sys_lines,
            "iono": iono, "nav_counts": dict(nav_counts)}


# ----------------------------------------------------------------------------
# RTKLIB config + solution helpers
# ----------------------------------------------------------------------------

def build_merged_atx(names, dest):
    """Extract ANTEX blocks for `names` from ngs20/igs20 into one small file."""
    found = []
    with open(dest, "w") as out:
        out.write("     1.4            M                                       "
                  "ANTEX VERSION / SYST\n")
        out.write("A                                                           "
                  "PCV TYPE / REFANT\n")
        out.write("                                                            "
                  "END OF HEADER\n")
        remaining = {n for n in names if n}
        for atx in ("ngs20.atx", "igs20.atx"):
            p = ANTEX_DIR / atx
            if not p.is_file() or not remaining:
                continue
            block, keep = [], False
            for line in open(p, encoding="utf-8", errors="replace"):
                if "START OF ANTENNA" in line:
                    block, keep = [line], False
                    continue
                if not block:
                    continue
                block.append(line)
                if "TYPE / SERIAL" in line:
                    name = line[:20].rstrip()
                    if name in remaining:
                        keep = True
                if "END OF ANTENNA" in line:
                    if keep:
                        out.write("                                                            "
                                  "START OF ANTENNA\n")
                        out.writelines(block)
                        found.append(line and name)
                        remaining.discard(name)
                    block = []
    return found


def rover_antenna():
    aj = TOOLS / "antenna.json"
    if aj.is_file():
        return json.loads(aj.read_text()).get("antex_name")
    return None


def write_conf(path, ant1=None, ant2=None, atx=None):
    lines = [
        "pos1-elmask =15",
        "pos1-navsys =45",          # GPS+GLO+GAL+BDS
        "pos2-armode =continuous",
        "pos2-gloarmode =off",      # mixed-manufacturer base/rover: GLONASS AR off
        "out-outstat =residual",
    ]
    if ant1:
        lines.append(f"ant1-anttype ={ant1}")
    if ant2:
        lines.append(f"ant2-anttype ={ant2}")
    if atx:
        lines.append(f"file-rcvantfile ={atx}")
        lines.append(f"file-satantfile ={ANTEX_DIR / 'igs20.atx'}")
    path.write_text("\n".join(lines) + "\n")
    return path


def read_pos(pos_path):
    """Parse an RTKLIB .pos file (-t -u => UTC, decimal degrees) into a DataFrame."""
    rows = []
    if not Path(pos_path).is_file():
        return pd.DataFrame()
    for line in open(pos_path, errors="replace"):
        if line.startswith("%") or not line.strip():
            continue
        p = line.split()
        if len(p) < 9:
            continue
        try:
            t = datetime.strptime(p[0] + " " + p[1], "%Y/%m/%d %H:%M:%S.%f"
                                  ).replace(tzinfo=timezone.utc)
            rows.append((t, float(p[2]), float(p[3]), float(p[4]), int(p[5]), int(p[6]),
                         float(p[7]), float(p[8]), float(p[9]),
                         float(p[13]) if len(p) > 13 else 0.0))
        except (ValueError, IndexError):
            continue
    return pd.DataFrame(rows, columns=["utc", "lat", "lon", "height", "Q", "ns",
                                       "sdn", "sde", "sdu", "ratio"])


# ----------------------------------------------------------------------------
# Stage 3: single-point solution
# ----------------------------------------------------------------------------

def stage_single(outdir, obs, nav, conf):
    pos = outdir / "single.pos"
    if not pos.is_file():
        r = run([BIN / "rnx2rtkp", "-k", conf, "-p", "0", "-f", "3", "-t", "-u",
                 "-o", pos, obs, nav])
        if not pos.is_file():
            raise RuntimeError(f"rnx2rtkp single failed:\n{r.stderr[-1500:]}")
    df = read_pos(pos)
    log(f"single: {len(df)} epochs")
    return df


# ----------------------------------------------------------------------------
# Stage 4: PPK vs NGS CORS
# ----------------------------------------------------------------------------

def load_cors_catalog():
    cat = CACHE / "itrf2020_geo.comp.txt"
    if not cat.is_file():
        CACHE.mkdir(parents=True, exist_ok=True)
        if fetch("https://geodesy.noaa.gov/corsdata/coord/coord_20/itrf2020_geo.comp.txt",
                 cat) is None:
            return []
    stations = []
    for line in open(cat, errors="replace"):
        p = line.split()
        if len(p) >= 17 and p[1] == "2020.00":
            try:
                lat = int(p[2]) + int(p[3]) / 60 + float(p[4]) / 3600
                lat = lat if p[5] == "N" else -lat
                lon = int(p[6]) + int(p[7]) / 60 + float(p[8]) / 3600
                lon = lon if p[9] == "E" else -lon
                stations.append({"id": p[0], "lat": lat, "lon": lon,
                                 "height": float(p[10]), "state": p[15], "status": p[16]})
            except ValueError:
                continue
    return stations


def cors_station_xyz(ssss, epoch_year):
    """ITRF2020 ECEF of a CORS station, velocity-propagated to epoch_year."""
    txt = fetch(f"https://geodesy.noaa.gov/corsdata/coord/coord_20/{ssss.lower()}_20.coord.txt")
    if txt is None:
        return None
    text = txt.decode(errors="replace")
    itrf = text[text.find("ITRF2020 POSITION"):]
    try:
        x = float(re.search(r"X\s*=\s*(-?[\d.]+)\s*m", itrf).group(1))
        y = float(re.search(r"Y\s*=\s*(-?[\d.]+)\s*m", itrf).group(1))
        z = float(re.search(r"Z\s*=\s*(-?[\d.]+)\s*m", itrf).group(1))
        vx = float(re.search(r"VX\s*=\s*(-?[\d.]+)\s*m/yr", itrf).group(1))
        vy = float(re.search(r"VY\s*=\s*(-?[\d.]+)\s*m/yr", itrf).group(1))
        vz = float(re.search(r"VZ\s*=\s*(-?[\d.]+)\s*m/yr", itrf).group(1))
    except AttributeError:
        return None
    dt = epoch_year - 2020.0
    return (x + vx * dt, y + vy * dt, z + vz * dt)


def download_base_hours(station, t0, t1, dest_dir, warnings):
    """Fetch CORS RINEX obs covering [t0, t1]; returns list of decompressed files."""
    dest_dir.mkdir(parents=True, exist_ok=True)
    ssss = station.lower()
    files = []
    hours = []
    t = t0.replace(minute=0, second=0, microsecond=0)
    while t <= t1:
        hours.append(t)
        t += timedelta(hours=1)
    age_days = (datetime.now(timezone.utc) - t1).days

    for t in hours:
        yy, ddd = t.strftime("%y"), t.strftime("%j")
        letter = chr(ord("a") + t.hour)
        candidates = [f"{ssss}{ddd}{letter}.{yy}o.gz", f"{ssss}{ddd}{letter}.{yy}d.gz"]
        if age_days >= 2:  # hourly files expire after ~2 days; daily file fallback
            candidates += [f"{ssss}{ddd}0.{yy}o.gz", f"{ssss}{ddd}0.{yy}d.gz"]
        got = None
        for cand in candidates:
            local_gz = dest_dir / cand
            local = dest_dir / cand[:-3]
            if local.is_file():
                got = local
                break
            for mirror in CORS_MIRRORS:
                url = f"{mirror}/rinex/{t.year}/{ddd}/{ssss}/{cand}"
                if fetch(url, local_gz) is not None:
                    local.write_bytes(gzip.decompress(local_gz.read_bytes()))
                    if cand.endswith(f"d.gz"):  # Hatanaka -> RINEX
                        crx2rnx = BIN / "CRX2RNX"
                        rnx = local.with_suffix(local.suffix[:-1] + "o")
                        r = run([crx2rnx, "-f", local])
                        local = rnx if rnx.is_file() else None
                    got = local
                    break
            if got:
                break
        if got and got not in files:
            files.append(got)
        elif not got:
            warnings.append(f"CORS {station}: no file for {t:%Y-%m-%d %H}:00 UTC yet")
    return files


def stage_ppk(outdir, obs, nav, data, args, warnings):
    result = {"status": "skipped"}
    t0 = data["t_start"] - timedelta(minutes=5)
    t1 = data["t_end"] + timedelta(minutes=5)
    mean_lat, mean_lon = data["pvt"]["lat"].mean(), data["pvt"]["lon"].mean()

    stations = load_cors_catalog()
    if not stations:
        warnings.append("CORS catalog unavailable (offline?) - PPK skipped")
        return result
    for s in stations:
        s["km"] = haversine_km(mean_lat, mean_lon, s["lat"], s["lon"])
    usable = sorted((s for s in stations if s["status"] == "Operational"
                     and s["km"] <= args.max_base_km), key=lambda s: s["km"])
    if args.base:
        usable = [s for s in stations if s["id"].upper() == args.base.upper()] + usable
    if not usable:
        warnings.append(f"no operational CORS station within {args.max_base_km} km")
        return result

    base_obs = None
    base_station = None
    for s in usable[:3]:
        log(f"ppk: trying CORS {s['id']} ({s['km']:.1f} km, {s['state']})")
        base_dir = outdir / "base" / s["id"].lower()
        files = download_base_hours(s["id"], t0, t1, base_dir, warnings)
        if not files:
            continue
        if len(files) == 1:
            base_obs = files[0]
        else:
            base_obs = base_dir / "base_combined.obs"
            if not base_obs.is_file():
                # combine_rinex.sh word-splits its args: run in base_dir with bare
                # filenames so paths never contain spaces
                r = run(["bash", COMBINE_SH, base_obs.name, *[f.name for f in files]],
                        cwd=base_dir)
                if not base_obs.is_file() or r.returncode != 0:
                    warnings.append(f"combine_rinex.sh failed for {s['id']}: {r.stdout[-500:]}")
                    base_obs = None
                    continue
        base_station = s
        break
    if base_obs is None:
        result["status"] = "pending"
        result["message"] = ("CORS base data not yet published (typical lag ~1 h). "
                             "Re-run this script later; completed stages are cached.")
        log("ppk: " + result["message"])
        return result

    # base antenna type from its RINEX header; base ECEF from NGS coord file
    ant2 = None
    for line in open(base_obs, errors="replace"):
        if "ANT # / TYPE" in line:
            ant2 = line[20:40].rstrip()
        if "END OF HEADER" in line:
            break
    epoch_year = data["t_start"].year + data["t_start"].timetuple().tm_yday / 365.25
    xyz = cors_station_xyz(base_station["id"], epoch_year)
    if xyz is None:
        warnings.append(f"NGS coord file for {base_station['id']} unavailable - "
                        "using base RINEX header position (NAD83, ~1.5 m frame offset)")

    ant1 = rover_antenna()
    merged = outdir / "merged.atx"
    build_merged_atx([ant1, ant2], merged)
    conf = write_conf(outdir / "ppk.conf", ant1=ant1, ant2=ant2, atx=merged)

    ts = ["-ts", t0.strftime("%Y/%m/%d"), t0.strftime("%H:%M:%S"),
          "-te", t1.strftime("%Y/%m/%d"), t1.strftime("%H:%M:%S")]
    rpos = ["-r", *map(str, xyz)] if xyz else []

    # Motion policy: NEVER infer static (deployment platform is a buoy - slow drift
    # within metres looks like noise). Static only when explicitly flagged.
    dn, de, _ = llh_to_enu(data["pvt"]["lat"], data["pvt"]["lon"], data["pvt"]["height"],
                           mean_lat, mean_lon, 0)
    spread95 = float(np.percentile(np.hypot(dn, de), 95))
    static = args.static
    if static and spread95 > 5.0:
        warnings.append(f"--static flagged but onboard track spread is {spread95:.1f} m "
                        f"(95th pct) - the static solution may be averaging real motion")
    elif not static and spread95 < 5.0:
        log(f"ppk: track spread {spread95:.2f} m would pass as static, but processing "
            f"kinematic (default policy; use --static only for a verified fixed antenna)")
    log(f"ppk: base={base_station['id']} ({base_station['km']:.1f} km) "
        f"ant2={ant2!r} static={static} (95% spread {spread95:.2f} m)")

    out = {}
    runs = [("kinematic", "2")] + ([("static", "3")] if static else [])
    for name, mode in runs:
        pos = outdir / f"ppk_{name}.pos"
        if not pos.is_file():
            r = run([BIN / "rnx2rtkp", "-k", conf, "-p", mode, "-f", "3", "-t", "-u",
                     *ts, *rpos, "-o", pos, obs, base_obs, nav])
            if not pos.is_file() or not read_pos(pos).shape[0]:
                warnings.append(f"rnx2rtkp ppk {name} produced no solution: {r.stderr[-800:]}")
                continue
        out[name] = read_pos(pos)
        qc = Counter(out[name]["Q"])
        log(f"ppk {name}: {len(out[name])} epochs, Q={dict(qc)}")

    result.update({"status": "ok" if out else "failed", "station": base_station,
                   "static": static, "solutions": out, "base_obs": str(base_obs),
                   "ant2": ant2, "base_xyz": xyz, "spread95": spread95})
    return result


# ----------------------------------------------------------------------------
# Stage 5: PPP vs open IGS products (optional)
# ----------------------------------------------------------------------------

def stage_ppp(outdir, obs, nav, data, args, warnings):
    day = data["t_start"]
    week = int((day - GPS_EPOCH).days // 7)
    yyyyddd = day.strftime("%Y%j")
    prod_dir = outdir / "igs"
    prod_dir.mkdir(exist_ok=True)

    sp3 = clk = None
    for base_url, prefix in IGS_PRODUCT_MIRRORS:
        listing = fetch(f"{base_url}/{week}/", timeout=60)
        if listing is None:
            continue
        names = set(re.findall(r'href="([^"]+\.(?:SP3|CLK)\.gz)"', listing.decode(errors="replace")))
        for kind in ("FIN", "RAP"):
            cand_sp3 = [n for n in names if re.match(prefix, n) and kind in n
                        and yyyyddd + "0000" in n and n.endswith("ORB.SP3.gz")]
            cand_clk = [n for n in names if re.match(prefix, n) and kind in n
                        and yyyyddd + "0000" in n and n.endswith("CLK.CLK.gz")]
            if cand_sp3 and cand_clk:
                fs, fc = prod_dir / cand_sp3[0], prod_dir / cand_clk[0]
                if ((fs.with_suffix("").is_file() or fetch(f"{base_url}/{week}/{cand_sp3[0]}", fs))
                        and (fc.with_suffix("").is_file() or fetch(f"{base_url}/{week}/{cand_clk[0]}", fc))):
                    for fgz in (fs, fc):
                        if fgz.is_file() and not fgz.with_suffix("").is_file():
                            fgz.with_suffix("").write_bytes(gzip.decompress(fgz.read_bytes()))
                    sp3, clk = fs.with_suffix(""), fc.with_suffix("")
                    log(f"ppp: using {sp3.name} + {clk.name}")
                    break
        if sp3:
            break
    if sp3 is None:
        msg = ("IGS precise products for this date not yet published "
               "(rapid ~1 day lag). Re-run with --ppp later.")
        log("ppp: " + msg)
        return {"status": "pending", "message": msg}

    ant1 = rover_antenna()
    merged = outdir / "merged_ppp.atx"
    build_merged_atx([ant1], merged)
    conf = write_conf(outdir / "ppp.conf", ant1=ant1, atx=merged)
    pos = outdir / "ppp.pos"
    ppp_mode = "7" if args.static else "6"      # ppp-static only when flagged static
    if pos.is_file() and not read_pos(pos).shape[0]:
        pos.unlink()                             # stale empty output - retry
    if not pos.is_file():
        run([BIN / "rnx2rtkp", "-k", conf, "-p", ppp_mode, "-f", "3", "-t", "-u",
             "-o", pos, obs, nav, sp3, clk])
        if not pos.is_file() or not read_pos(pos).shape[0]:
            pos.unlink(missing_ok=True)          # empty output: let future runs retry
            mode_name = "static" if args.static else "kinematic"
            warnings.append(f"PPP ({mode_name}) produced no solution - RTKLIB's "
                            "PPP-kinematic engine is a known weak point. PPK kinematic "
                            "remains the post-processed reference solution. "
                            "(PPP-static is available only with --static.)")
            return {"status": "failed"}
    df = read_pos(pos)
    log(f"ppp: {len(df)} epochs")
    warnings.append("PPP note: 36-minute session is short for full PPP convergence; "
                    "treat PPP as a cross-check, not truth.")
    return {"status": "ok", "df": df, "sp3": sp3.name, "clk": clk.name}


# ----------------------------------------------------------------------------
# Stage 6: track outputs
# ----------------------------------------------------------------------------

def stage_track(outdir, solutions, pvt, interval):
    # best available source, in order of accuracy
    for source in ("ppk_kinematic", "ppp", "single"):
        if source in solutions and len(solutions[source]):
            best = solutions[source]
            break
    else:
        return None, None

    df = best.copy()
    df["utc_s"] = df["utc"].dt.round("1s")
    ob = pvt.copy()
    ob["utc_s"] = ob["utc"].dt.round("1s")
    ob = ob[["utc_s", "lat", "lon", "height", "hAcc", "fixType"]].rename(
        columns={"lat": "onboard_lat", "lon": "onboard_lon", "height": "onboard_height",
                 "hAcc": "onboard_hAcc_m", "fixType": "onboard_fixType"})
    df = df.merge(ob, on="utc_s", how="left").drop(columns=["utc_s"])
    df["source"] = source
    if interval > 1:
        epoch_s = (df["utc"] - pd.Timestamp(0, tz="utc")).dt.total_seconds().round()
        df = df[epoch_s % interval == 0]
    csv_path = outdir / "track.csv"
    with open(csv_path, "w") as f:
        f.write(f"# GPS track - solution source: {source}\n")
        f.write("# lat/lon: decimal degrees; heights: ellipsoidal metres\n")
        if source.startswith("ppk"):
            f.write("# datum note: PPK positions are in the CORS base frame "
                    "(ITRF2020 current epoch when -r was applied)\n")
        cols = ["utc", "lat", "lon", "height", "Q", "ns", "sdn", "sde", "sdu",
                "source", "onboard_lat", "onboard_lon", "onboard_height",
                "onboard_hAcc_m", "onboard_fixType"]
        df.to_csv(f, index=False, columns=cols,
                  float_format="%.9f", date_format="%Y-%m-%dT%H:%M:%S.%fZ")
    # KML for the chosen solution
    kml_src = {"ppk_kinematic": "ppk_kinematic.pos", "ppp": "ppp.pos",
               "single": "single.pos"}[source]
    run([BIN / "pos2kml", outdir / kml_src])
    kml = (outdir / kml_src).with_suffix(".kml")
    if kml.is_file():
        shutil.move(kml, outdir / "track.kml")
    log(f"track: {len(df)} rows from {source} -> track.csv"
        + (" + track.kml" if (outdir / 'track.kml').is_file() else ""))
    return df, source


# ----------------------------------------------------------------------------
# Figures (dataviz reference palette; static light-mode PNGs)
# ----------------------------------------------------------------------------
# Categorical slots (validated adjacent-pairs, light mode): fixed constellation order.
SERIES = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100", "#e87ba4"]  # blue orange aqua yellow magenta
CONST_ORDER = ["GPS", "Galileo", "BeiDou", "GLONASS", "SBAS"]
CONST_COLOR = dict(zip(CONST_ORDER, SERIES))
SEQ_RAMP = ["#cde2fb", "#9ec5f4", "#6da7ec", "#3987e5", "#256abf", "#184f95", "#0d366b"]
SURFACE, INK, INK2, MUTED = "#fcfcfb", "#0b0b0b", "#52514e", "#898781"
GRID, BASELINE = "#e1e0d9", "#c3c2b7"


def _style_ax(ax, xlabel=None, ylabel=None, title=None):
    ax.set_facecolor(SURFACE)
    for side in ("top", "right"):
        ax.spines[side].set_visible(False)
    for side in ("left", "bottom"):
        ax.spines[side].set_color(BASELINE)
    ax.tick_params(colors=MUTED, labelsize=8)
    ax.grid(True, color=GRID, linewidth=0.7, alpha=0.9)
    ax.set_axisbelow(True)
    if xlabel:
        ax.set_xlabel(xlabel, color=INK2, fontsize=9)
    if ylabel:
        ax.set_ylabel(ylabel, color=INK2, fontsize=9)
    if title:
        ax.set_title(title, color=INK, fontsize=10, loc="left", fontweight="bold")


def stage_figures(outdir, data, solutions, track_source, diag):
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from matplotlib.lines import Line2D
    from matplotlib.colors import LinearSegmentedColormap

    figdir = outdir / "figures"
    figdir.mkdir(exist_ok=True)
    figs = []
    pvt, rawx, sat = data["pvt"], data["rawx"], data["sat"]
    ref = diag["position_error"]["reference"]
    refpos = (ref["lat"], ref["lon"], ref["ellip_height_m"])
    best = solutions.get(track_source) if track_source else None
    seq_cmap = LinearSegmentedColormap.from_list("seq", SEQ_RAMP)

    def newfig(w=8.0, h=4.5):
        f = plt.figure(figsize=(w, h), facecolor=SURFACE, dpi=150)
        return f

    def save(f, name, title):
        f.savefig(figdir / name, facecolor=SURFACE, bbox_inches="tight")
        plt.close(f)
        figs.append((name, title))

    # ---- 1. horizontal position scatter (track map) with CEP circles --------
    f = newfig(6.4, 6.0)
    ax = f.add_subplot(111)
    series = []
    if best is not None and len(best):
        dn, de, _ = llh_to_enu(best["lat"], best["lon"], best["height"], *refpos)
        series.append((track_source, de, dn, SERIES[0]))
    dn_o, de_o, _ = llh_to_enu(pvt["lat"], pvt["lon"], pvt["height"], *refpos)
    series.append(("onboard (NAV-PVT)", de_o, dn_o, SERIES[1]))
    for name, de, dn, color in series:
        ax.scatter(de, dn, s=6, c=color, alpha=0.5, linewidths=0, label=name)
    for name, de, dn, color in series:            # mean markers with surface ring
        ax.scatter([np.mean(de)], [np.mean(dn)], s=90, c=color, edgecolors=SURFACE,
                   linewidths=2, zorder=5)
    is_static = diag["session"].get("motion_mode", "").startswith("static")
    radius_label = "CEP" if is_static else "R"     # kinematic: radii include motion
    if best is not None and len(best):
        dn, de = series[0][2], series[0][1]
        r = np.hypot(dn - dn.mean(), de - de.mean())
        for pct, ls in ((50, "-"), (95, "--")):
            rad = np.percentile(r, pct)
            circ = plt.Circle((de.mean(), dn.mean()), rad, fill=False,
                              color=MUTED, linewidth=1.0, linestyle=ls)
            ax.add_patch(circ)
            ax.annotate(f"{radius_label}{pct} {rad:.2f} m", (de.mean(), dn.mean() + rad),
                        color=INK2, fontsize=8, ha="center", va="bottom")
    ax.set_aspect("equal")
    ax.legend(loc="upper right", fontsize=8, frameon=False, labelcolor=INK2)
    scatter_word = "scatter" if is_static else "dispersion (incl. platform motion)"
    _style_ax(ax, "East (m)", "North (m)",
              f"Horizontal position {scatter_word}")
    save(f, "fig1_track_scatter.png",
         f"Horizontal position {scatter_word} with {radius_label}50/{radius_label}95 radii")

    # ---- 2. ENU error time series ------------------------------------------
    f = newfig(8.5, 6.0)
    axes = f.subplots(3, 1, sharex=True)
    comps = ["North", "East", "Up"]
    on_enu = llh_to_enu(pvt["lat"], pvt["lon"], pvt["height"], *refpos)
    best_enu = (llh_to_enu(best["lat"], best["lon"], best["height"], *refpos)
                if best is not None and len(best) else None)
    for i, ax in enumerate(axes):
        if best_enu is not None:
            ax.plot(best["utc"], best_enu[i], color=SERIES[0], linewidth=1.4,
                    label=track_source)
        ax.plot(pvt["utc"], on_enu[i], color=SERIES[1], linewidth=1.4,
                label="onboard")
        ax.axhline(0, color=BASELINE, linewidth=0.8)
        _style_ax(ax, None, f"{comps[i]} (m)")
    axes[0].legend(loc="lower right", bbox_to_anchor=(1, 1.02), fontsize=8,
                   frameon=False, ncols=2, labelcolor=INK2)
    axes[0].set_title("Position error vs reference, by component", color=INK,
                      fontsize=10, loc="left", fontweight="bold")
    axes[-1].set_xlabel("UTC", color=INK2, fontsize=9)
    f.autofmt_xdate()
    save(f, "fig2_enu_timeseries.png", "Per-component position error time series")

    # ---- 3. satellites: counts over time + C/N0 heatmap --------------------
    f = newfig(8.5, 6.4)
    gs = f.add_gridspec(2, 1, height_ratios=[1.15, 1], hspace=0.42)
    ax = f.add_subplot(gs[0])
    hdr = data["rawx_hdr"]
    times = [gps_to_utc(int(w), t, int(hdr["leapS"].iloc[0]))
             for w, t in zip(hdr["week"], hdr["rcvTow"])]
    end_vals = {}
    for name in CONST_ORDER:
        gid = {v: k for k, v in GNSS_NAME.items()}[name]
        per_epoch = rawx[rawx["gnss"] == gid].groupby("epoch")["sv"].nunique()
        counts = per_epoch.reindex(range(len(hdr)), fill_value=0)
        ax.plot(times, counts, color=CONST_COLOR[name], linewidth=1.4, label=name)
        end_vals[name] = float(counts.iloc[-1])
    # selective direct labels: only when the line end is clearly separated
    for name, v in end_vals.items():
        if all(abs(v - o) > 1.0 for n, o in end_vals.items() if n != name):
            ax.annotate(name, (times[-1], v), xytext=(4, 0),
                        textcoords="offset points", color=CONST_COLOR[name],
                        fontsize=8, va="center")
    # legend sits in the title-pad band, directly above the axes; title above it
    ax.legend(loc="lower left", fontsize=8, frameon=False, ncols=5,
              bbox_to_anchor=(0, 1.01), labelcolor=INK2)
    _style_ax(ax, None, "satellites tracked")
    ax.set_title("Satellites tracked per constellation", color=INK, fontsize=10,
                 loc="left", fontweight="bold", pad=26)

    axc = f.add_subplot(gs[1])
    bands = ["L1", "L2", "L5", "L6/E6"]
    mat = np.full((len(CONST_ORDER), len(bands)), np.nan)
    for i, cname in enumerate(CONST_ORDER):
        gid = {v: k for k, v in GNSS_NAME.items()}[cname]
        g = rawx[rawx["gnss"] == gid]
        for j, b in enumerate(bands):
            sel = g[g["freq"].map(band_of) == b]["cno"]
            if len(sel):
                mat[i, j] = sel.mean()
    im = axc.imshow(mat, cmap=seq_cmap, aspect="auto", vmin=25, vmax=45)
    axc.set_xticks(range(len(bands)), bands)
    axc.set_yticks(range(len(CONST_ORDER)), CONST_ORDER)
    for i in range(mat.shape[0]):
        for j in range(mat.shape[1]):
            if not np.isnan(mat[i, j]):
                axc.text(j, i, f"{mat[i, j]:.1f}", ha="center", va="center",
                         color=INK if mat[i, j] < 38 else SURFACE, fontsize=9)
    axc.set_title("Mean C/N0 by constellation and band (dB-Hz)", color=INK,
                  fontsize=10, loc="left", fontweight="bold")
    axc.tick_params(colors=MUTED, labelsize=9)
    for sp in axc.spines.values():
        sp.set_visible(False)
    save(f, "fig3_satellites.png", "Satellite counts and C/N0 by band")

    # ---- 4. skyplot ---------------------------------------------------------
    f = newfig(6.4, 6.4)
    ax = f.add_subplot(111, projection="polar")
    ax.set_facecolor(SURFACE)
    ax.set_theta_zero_location("N")
    ax.set_theta_direction(-1)
    sky = (sat[(sat["elev"] > 0) & (sat["cno"] > 0)]
           .groupby(["gnss", "sv"]).agg(elev=("elev", "mean"), azim=("azim", "mean"),
                                        cno=("cno", "mean")).reset_index())
    sc = ax.scatter(np.radians(sky["azim"]), 90 - sky["elev"], c=sky["cno"],
                    cmap=seq_cmap, s=42, vmin=25, vmax=45,
                    edgecolors=SURFACE, linewidths=1)
    ax.set_rlim(0, 90)
    ax.set_rgrids([30, 60], ["60°", "30°"], color=MUTED, fontsize=8)
    ax.set_thetagrids([0, 90, 180, 270], ["N", "E", "S", "W"], color=INK2)
    ax.grid(color=GRID, linewidth=0.7)
    cb = f.colorbar(sc, ax=ax, shrink=0.7, pad=0.08)
    cb.set_label("mean C/N0 (dB-Hz)", color=INK2, fontsize=9)
    cb.ax.tick_params(colors=MUTED, labelsize=8)
    cb.outline.set_visible(False)
    ax.set_title("Sky plot (session-mean position per satellite)", color=INK,
                 fontsize=10, fontweight="bold")
    save(f, "fig4_skyplot.png", "Sky plot colored by C/N0")

    # ---- 5. accuracy metrics over time -------------------------------------
    f = newfig(8.5, 5.4)
    ax1, ax2 = f.subplots(2, 1, sharex=True)
    ax1.plot(pvt["utc"], pvt["hAcc"], color=SERIES[1], linewidth=1.4,
             label="onboard hAcc")
    ax1.plot(pvt["utc"], pvt["vAcc"], color=SERIES[1], linewidth=1.0,
             linestyle="--", label="onboard vAcc")
    if best is not None and "sdn" in best and best["sdn"].abs().sum() > 0:
        sdh = np.hypot(best["sdn"], best["sde"])
        ax1.plot(best["utc"], sdh, color=SERIES[0], linewidth=1.4,
                 label=f"{track_source} sd horiz")
        ax1.plot(best["utc"], best["sdu"], color=SERIES[0], linewidth=1.0,
                 linestyle="--", label=f"{track_source} sd up")
    ax1.legend(loc="upper right", fontsize=8, frameon=False, ncols=2, labelcolor=INK2)
    _style_ax(ax1, None, "1-σ accuracy (m)", "Reported accuracy estimates")
    ax2.plot(pvt["utc"], pvt["pDOP"], color=SERIES[2], linewidth=1.4)
    ax2.annotate("pDOP", (pvt["utc"].iloc[-1], pvt["pDOP"].iloc[-1]),
                 xytext=(4, 0), textcoords="offset points", color=SERIES[2],
                 fontsize=8, va="center")
    _style_ax(ax2, "UTC", "pDOP", None)
    f.autofmt_xdate()
    save(f, "fig5_accuracy.png", "Accuracy estimates and DOP")

    # ---- 6. ionosphere + RF -------------------------------------------------
    f = newfig(8.5, 6.2)
    gs = f.add_gridspec(2, 2, height_ratios=[1.2, 1], hspace=0.5, wspace=0.3)
    axi = f.add_subplot(gs[0, :])
    gf = diag["ionosphere"]["geometry_free_per_satellite"]
    if gf:
        labels = [f"{r['gnss'][:3]}{r['sv']:02d}" for r in gf]
        vals = [r["iono_L1_m_median"] for r in gf]
        cols = [CONST_COLOR.get(r["gnss"], MUTED) for r in gf]
        x = np.arange(len(gf))
        axi.bar(x, vals, color=cols, width=0.7)
        axi.set_xticks(x, labels, rotation=90, fontsize=6.5)
        handles = [Line2D([0], [0], marker="s", linestyle="", color=CONST_COLOR[c],
                          label=c) for c in CONST_ORDER if any(r["gnss"] == c for r in gf)]
        axi.legend(handles=handles, loc="upper right", fontsize=8, frameon=False,
                   ncols=len(handles), labelcolor=INK2)
    _style_ax(axi, None, "slant delay @L1 (m)",
              "Dual-frequency geometry-free ionospheric delay (session median per satellite)")

    axr = f.add_subplot(gs[1, 0])
    rfagg = data["rf"].groupby("blockId")["jamInd"].agg(["mean", "max"])
    x = np.arange(len(rfagg))
    axr.bar(x, rfagg["mean"], color=SERIES[0], width=0.55, label="mean")
    axr.scatter(x, rfagg["max"], color=INK2, s=24, zorder=5, label="max")
    axr.set_xticks(x, [f"block {b}" for b in rfagg.index])
    axr.set_ylim(0, 255)
    axr.axhline(255 * 0.6, color=MUTED, linewidth=0.8, linestyle="--")
    axr.annotate("interference concern >153", (0, 158), color=MUTED, fontsize=7)
    axr.legend(loc="upper right", fontsize=8, frameon=False, labelcolor=INK2)
    _style_ax(axr, None, "jamInd (0-255)", "RF interference indicator")

    axj = f.add_subplot(gs[1, 1])
    jitter_ms = (np.diff(hdr["rcvTow"].to_numpy()) - 1.0) * 1000
    axj.hist(jitter_ms, bins=40, color=SERIES[0])
    _style_ax(axj, "epoch interval error (ms)", "epochs", "Measurement epoch jitter")
    save(f, "fig6_iono_rf.png", "Ionospheric delay, RF health, timing jitter")

    log(f"figures: {len(figs)} PNGs -> {figdir}")
    return figs


# ----------------------------------------------------------------------------
# Stage 7: diagnostics + report
# ----------------------------------------------------------------------------

def solution_stats(df, ref):
    """Error metrics for one solution DataFrame about reference (lat, lon, h)."""
    if df is None or not len(df):
        return None
    dn, de, du = llh_to_enu(df["lat"], df["lon"], df["height"], *ref)
    horiz = np.hypot(dn - dn.mean(), de - de.mean())      # scatter about own mean
    stats = {
        "epochs": int(len(df)),
        "mean_lat": float(df["lat"].mean()), "mean_lon": float(df["lon"].mean()),
        "mean_height": float(df["height"].mean()),
        "offset_from_ref_NEU_m": [float(dn.mean()), float(de.mean()), float(du.mean())],
        "rms_NEU_m": [float(np.sqrt(np.mean((dn - dn.mean()) ** 2))),
                      float(np.sqrt(np.mean((de - de.mean()) ** 2))),
                      float(np.sqrt(np.mean((du - du.mean()) ** 2)))],
        "cep50_m": float(np.percentile(horiz, 50)),
        "cep95_m": float(np.percentile(horiz, 95)),
        "2drms_m": float(2 * np.sqrt(np.var(dn) + np.var(de))),
    }
    if "sdn" in df and df["sdn"].abs().sum() > 0:
        stats["formal_sd_NEU_m_mean"] = [float(df["sdn"].mean()), float(df["sde"].mean()),
                                         float(df["sdu"].mean())]
    if "Q" in df:
        qc = Counter(df["Q"])
        stats["quality_histogram"] = {QUALITY_NAME.get(k, str(k)): int(v)
                                      for k, v in sorted(qc.items())}
        if qc.get(1) or qc.get(2):          # fix rate only meaningful for RTK/PPK
            stats["fix_rate"] = float(qc.get(1, 0) / len(df))
    return stats


def geometry_free_iono(rawx):
    """Per-satellite slant iono estimate (m at L1) from dual-frequency pseudoranges."""
    valid = rawx[(rawx["prValid"] == 1) & (rawx["freq"] > 0)]
    out = []
    for (gnss, sv), g in valid.groupby(["gnss", "sv"]):
        freqs = g.groupby("freq")["prMes"].count()
        if len(freqs) < 2:
            continue
        f1 = max(freqs.index)                    # highest frequency (L1-ish)
        f2 = min(freqs.index)                    # lowest (L5/L2)
        if f1 - f2 < 50e6:
            continue
        a = g[g["freq"] == f1].set_index("epoch")["prMes"]
        b = g[g["freq"] == f2].set_index("epoch")["prMes"]
        both = pd.concat([a.groupby(level=0).first(), b.groupby(level=0).first()],
                         axis=1, keys=["p1", "p2"]).dropna()
        if len(both) < 30:
            continue
        gamma = (f1 / f2) ** 2
        i_f1 = (both["p2"] - both["p1"]) / (gamma - 1)       # iono delay at f1, metres
        i_l1 = i_f1 * (f1 / 1575.42e6) ** 2                  # scale to L1
        tecu = i_f1 * f1 ** 2 / 40.3 / 1e16                  # TECU (1e16 el/m^2)
        out.append({"gnss": GNSS_NAME.get(gnss, gnss), "sv": int(sv),
                    "pair": f"{band_of(f1)}/{band_of(f2)}", "n": int(len(both)),
                    "iono_L1_m_median": round(float(i_l1.median()), 3),
                    "iono_L1_m_iqr": round(float(i_l1.quantile(.75) - i_l1.quantile(.25)), 3),
                    "TECU_median": round(float(tecu.median()), 2)})
    return sorted(out, key=lambda r: (r["gnss"], r["sv"]))


def stage_diagnostics(outdir, data, rinex_info, solutions, ppk_result, track_source,
                      warnings, args):
    rawx, pvt, sat, sig, rf = data["rawx"], data["pvt"], data["sat"], data["sig"], data["rf"]

    # --- reference position: PPK static end-state (only if --static was flagged)
    #     > PPK kinematic mean > onboard mean. In kinematic mode the reference is a
    #     dispersion anchor only - "scatter" includes any real platform motion.
    ref_source = "onboard session mean"
    ref = (float(pvt["lat"].mean()), float(pvt["lon"].mean()), float(pvt["height"].mean()))
    if "ppk_static" in solutions and len(solutions["ppk_static"]):
        last = solutions["ppk_static"].iloc[-1]
        ref = (float(last["lat"]), float(last["lon"]), float(last["height"]))
        ref_source = f"PPK static final epoch (Q={QUALITY_NAME.get(int(last['Q']))})"
    elif "ppk_kinematic" in solutions and len(solutions["ppk_kinematic"]):
        d = solutions["ppk_kinematic"]
        ref = (float(d["lat"].mean()), float(d["lon"].mean()), float(d["height"].mean()))
        ref_source = "PPK kinematic session mean (dispersion anchor; includes platform motion)"

    # --- constellation metrics
    per_gnss = {}
    for gnss, g in rawx.groupby("gnss"):
        name = GNSS_NAME.get(gnss, str(gnss))
        bands = {}
        for band, gb in g.groupby(g["freq"].map(band_of)):
            bands[band] = {"signals": int(gb.groupby(["sv", "sig"]).ngroups),
                           "meas": int(len(gb)),
                           "cno_mean": round(float(gb["cno"].mean()), 1),
                           "cno_min": int(gb["cno"].min()), "cno_max": int(gb["cno"].max())}
        # cycle-slip proxy: locktime decreases per (sv, sig) stream
        slips = 0
        for _, gs in g.groupby(["sv", "sig"]):
            lt = gs.sort_values("epoch")["locktime"].to_numpy()
            slips += int((np.diff(lt) < 0).sum())
        used = sat[(sat["gnss"] == gnss) & (sat["used"] == 1)]
        prres = sat[(sat["gnss"] == gnss) & (sat["used"] == 1) & (sat["prRes"] != 0)]["prRes"]
        per_gnss[name] = {
            "satellites_tracked": int(g["sv"].nunique()),
            "satellites_used_avg": round(len(used) / max(1, data["counts"]["NAV-SAT"]), 1),
            "signals_tracked": int(g.groupby(["sv", "sig"]).ngroups),
            "measurements": int(len(g)),
            "bands": bands,
            "cycle_slip_events": slips,
            "sfrbx_frames": int(data["sfrbx_gnss"].get(gnss, 0)),
            "prRes_rms_m": round(float(np.sqrt((prres ** 2).mean())), 2) if len(prres) else None,
        }

    # --- ionosphere
    iono_gf = geometry_free_iono(rawx)
    iono_section = {
        "broadcast_klobuchar": rinex_info.get("iono", {}),
        "ionoModel_census": {IONO_MODEL.get(k, str(k)): int(v)
                             for k, v in sig["ionoModel"].value_counts().items()},
        "corrSource_census": {CORR_SOURCE.get(k, str(k)): int(v)
                              for k, v in sig["corrSource"].value_counts().items()},
        "geometry_free_per_satellite": iono_gf,
        "note": "geometry-free estimates are slant (not vertical) delays derived from "
                "dual-frequency code; median over session, metres at L1 frequency",
    }

    # --- positional error metrics
    pos_section = {"reference": {"source": ref_source, "lat": ref[0], "lon": ref[1],
                                 "ellip_height_m": ref[2]}}
    onboard_df = pvt.rename(columns={})[["utc", "lat", "lon", "height"]].copy()
    for name, df in [("onboard", onboard_df)] + list(solutions.items()):
        s = solution_stats(df, ref)
        if s:
            pos_section[name] = s
    pos_section["onboard_reported_accuracy"] = {
        "hAcc_m": {"mean": round(float(pvt["hAcc"].mean()), 3),
                   "min": round(float(pvt["hAcc"].min()), 3),
                   "max": round(float(pvt["hAcc"].max()), 3)},
        "vAcc_m": {"mean": round(float(pvt["vAcc"].mean()), 3)},
        "pDOP": {"mean": round(float(pvt["pDOP"].mean()), 2),
                 "max": round(float(pvt["pDOP"].max()), 2)},
    }
    if ppk_result.get("station"):
        st = ppk_result["station"]
        pos_section["ppk_base"] = {"id": st["id"], "distance_km": round(st["km"], 1),
                                   "antenna": ppk_result.get("ant2"),
                                   "ecef_itrf2020": ppk_result.get("base_xyz")}

    # --- timing metrics
    jitter_ms = (np.diff(data["rawx_hdr"]["rcvTow"].to_numpy()) - 1.0) * 1000
    timing = {
        "epochs": int(len(data["rawx_hdr"])),
        "nominal_interval_s": 1.0,
        "epoch_jitter_ms": {"mean": round(float(jitter_ms.mean()), 4),
                            "std": round(float(jitter_ms.std()), 4),
                            "max_abs": round(float(np.abs(jitter_ms).max()), 4)},
        "tAcc_ns": {"mean": round(float(pvt["tAcc_ns"].mean()), 1),
                    "min": int(pvt["tAcc_ns"].min()), "max": int(pvt["tAcc_ns"].max())},
        "leap_seconds": {"rawx": int(data["rawx_hdr"]["leapS"].iloc[0]),
                         "valid_flag": bool(data["rawx_hdr"]["leapSecValid"].iloc[0])},
    }

    # --- validation of computed solutions against the receiver's own embedded
    #     position solution (NAV-PVT), matched epoch-by-epoch
    validation = {}
    on = pvt[["utc", "lat", "lon", "height"]].copy()
    on["utc_s"] = on["utc"].dt.round("1s")
    for name, df in solutions.items():
        if df is None or not len(df):
            continue
        d = df.copy()
        d["utc_s"] = d["utc"].dt.round("1s")
        j = d.merge(on, on="utc_s", suffixes=("", "_ob"))
        if not len(j):
            continue
        dn, de, du = llh_to_enu(j["lat"], j["lon"], j["height"],
                                float(j["lat_ob"].mean()), float(j["lon_ob"].mean()),
                                float(j["height_ob"].mean()))
        dn_o, de_o, du_o = llh_to_enu(j["lat_ob"], j["lon_ob"], j["height_ob"],
                                      float(j["lat_ob"].mean()), float(j["lon_ob"].mean()),
                                      float(j["height_ob"].mean()))
        ddn, dde, ddu = dn - dn_o, de - de_o, du - du_o
        horiz = float(np.hypot(ddn.mean(), dde.mean()))
        # agreement thresholds vs the onboard ~1 m (1-sigma) SBAS solution, scaled
        # by solution quality: ambiguity-fixed PPK must agree tightly; float PPK
        # carries metre-level bias uncertainty; single-point is the loosest
        if name.startswith(("ppk", "ppp")):
            fixed = (df["Q"] == 1).mean() > 0.5 if "Q" in df else False
            thresh_h = 1.5 if fixed else 4.0
        else:
            thresh_h = 5.0
        validation[name] = {
            "matched_epochs": int(len(j)),
            "mean_offset_NEU_m": [round(float(ddn.mean()), 3), round(float(dde.mean()), 3),
                                  round(float(ddu.mean()), 3)],
            "rms_diff_NEU_m": [round(float(np.sqrt((ddn ** 2).mean())), 3),
                               round(float(np.sqrt((dde ** 2).mean())), 3),
                               round(float(np.sqrt((ddu ** 2).mean())), 3)],
            "horizontal_offset_m": round(horiz, 3),
            "threshold_m": thresh_h,
            "verdict": "PASS" if horiz < thresh_h else "FAIL",
        }
        log(f"validate {name} vs onboard: horiz offset {horiz:.3f} m "
            f"(threshold {thresh_h} m) -> {validation[name]['verdict']}")

    # --- RF health
    rf_section = {}
    for block, g in rf.groupby("blockId"):
        rf_section[f"block_{block}"] = {
            "jamInd_mean": round(float(g["jamInd"].mean()), 1),
            "jamInd_max": int(g["jamInd"].max()),
            "jammingState_nonzero": int((g["jammingState"] != 0).sum()),
            "noisePerMS_mean": round(float(g["noisePerMS"].mean()), 1),
            "agcCnt_mean": round(float(g["agcCnt"].mean()), 0),
            "antStatus": ANT_STATUS.get(int(g["antStatus"].mode()[0]), "?"),
        }

    diagnostics = {
        "session": {
            "file": str(args.capture),
            "motion_mode": "static (explicitly flagged)" if args.static
                           else "kinematic (default)",
            "receiver": data["meta"].get("extensions", []),
            "module_fw": data["meta"].get("swVersion"),
            "start_utc": data["t_start"].isoformat(),
            "end_utc": data["t_end"].isoformat(),
            "duration_min": round((data["t_end"] - data["t_start"]).total_seconds() / 60, 1),
            "gps_week": int(data["rawx_hdr"]["week"].iloc[0]),
            "message_counts": {k: int(v) for k, v in data["counts"].most_common()},
            "track_solution_source": track_source,
        },
        "constellations": per_gnss,
        "ionosphere": iono_section,
        "position_error": pos_section,
        "validation_vs_onboard": validation,
        "timing": timing,
        "rf_health": rf_section,
        "warnings": warnings,
    }
    (outdir / "diagnostics.json").write_text(json.dumps(diagnostics, indent=2, default=str))
    log(f"diagnostics: diagnostics.json written ({len(warnings)} warnings)")
    return diagnostics


def write_report(outdir, d, figs=()):
    L = []
    s = d["session"]
    L.append(f"# GPS Session Report - {Path(s['file']).name}\n")
    L.append(f"- **Receiver**: {s.get('module_fw')} "
             f"({', '.join(x for x in (s.get('receiver') or []) if 'MOD=' in x)})")
    L.append(f"- **Window**: {s['start_utc']} → {s['end_utc']} ({s['duration_min']} min, "
             f"GPS week {s['gps_week']})")
    L.append(f"- **Track solution**: `{s['track_solution_source']}`\n")

    ref = d["position_error"]["reference"]
    L.append("## Position summary\n")
    L.append(f"Motion mode: **{s.get('motion_mode', '?')}**\n")
    L.append(f"Reference ({ref['source']}): **{ref['lat']:.9f}°, {ref['lon']:.9f}°**, "
             f"h = {ref['ellip_height_m']:.3f} m\n")
    if not s.get("motion_mode", "").startswith("static"):
        L.append("> Kinematic mode: RMS/CEP dispersion figures include any real platform "
                 "motion, not just measurement error.\n")
    L.append("| Solution | Epochs | Offset from ref N/E/U (m) | RMS N/E/U (m) | CEP50 | CEP95 | Fix rate |")
    L.append("|---|---|---|---|---|---|---|")
    for name in ("ppk_kinematic", "ppk_static", "ppp", "single", "onboard"):
        st = d["position_error"].get(name)
        if not st:
            continue
        off = "/".join(f"{v:+.3f}" for v in st["offset_from_ref_NEU_m"])
        rms = "/".join(f"{v:.3f}" for v in st["rms_NEU_m"])
        fr = f"{st['fix_rate']*100:.1f}%" if "fix_rate" in st else "-"
        L.append(f"| {name} | {st['epochs']} | {off} | {rms} | {st['cep50_m']:.3f} | "
                 f"{st['cep95_m']:.3f} | {fr} |")
    oa = d["position_error"]["onboard_reported_accuracy"]
    L.append(f"\nOnboard self-reported: hAcc {oa['hAcc_m']['mean']} m mean "
             f"({oa['hAcc_m']['min']}-{oa['hAcc_m']['max']}), pDOP {oa['pDOP']['mean']}\n")
    if "ppk_base" in d["position_error"]:
        b = d["position_error"]["ppk_base"]
        L.append(f"PPK base: **{b['id']}** at {b['distance_km']} km, antenna `{b['antenna']}`\n")

    if d.get("validation_vs_onboard"):
        L.append("## Validation against onboard (NAV-PVT) positions\n")
        L.append("Computed solutions cross-checked epoch-by-epoch against the receiver's "
                 "own embedded position solution:\n")
        L.append("| Solution | Matched epochs | Mean offset N/E/U (m) | RMS diff N/E/U (m) | Horiz offset | Threshold | Verdict |")
        L.append("|---|---|---|---|---|---|---|")
        for name, v in d["validation_vs_onboard"].items():
            off = "/".join(f"{x:+.3f}" for x in v["mean_offset_NEU_m"])
            rms = "/".join(f"{x:.3f}" for x in v["rms_diff_NEU_m"])
            L.append(f"| {name} | {v['matched_epochs']} | {off} | {rms} | "
                     f"{v['horizontal_offset_m']} m | < {v['threshold_m']} m | "
                     f"**{v['verdict']}** |")
        L.append("")

    if figs:
        L.append("## Figures\n")
        for name, title in figs:
            L.append(f"### {title}\n")
            L.append(f"![{title}](figures/{name})\n")

    L.append("## Constellations\n")
    L.append("| GNSS | Sats | Signals | Meas | C/N0 by band (mean dBHz) | Slips | SFRBX | prRes RMS |")
    L.append("|---|---|---|---|---|---|---|---|")
    for name, c in d["constellations"].items():
        bands = ", ".join(f"{b}:{v['cno_mean']}" for b, v in sorted(c["bands"].items()))
        L.append(f"| {name} | {c['satellites_tracked']} | {c['signals_tracked']} | "
                 f"{c['measurements']} | {bands} | {c['cycle_slip_events']} | "
                 f"{c['sfrbx_frames']} | {c['prRes_rms_m'] or '-'} |")

    L.append("\n## Ionosphere\n")
    for k, v in d["ionosphere"]["broadcast_klobuchar"].items():
        L.append(f"- Broadcast `{k}`: {v}")
    L.append(f"- Signal iono-model census: {d['ionosphere']['ionoModel_census']}")
    gf = d["ionosphere"]["geometry_free_per_satellite"]
    if gf:
        med = np.median([r["iono_L1_m_median"] for r in gf])
        L.append(f"- Dual-frequency geometry-free slant delay, session median across "
                 f"{len(gf)} satellites: **{med:.2f} m at L1** "
                 f"(~{np.median([r['TECU_median'] for r in gf]):.1f} TECU); "
                 f"per-satellite detail in diagnostics.json")

    t = d["timing"]
    L.append("\n## Timing\n")
    L.append(f"- Epoch interval jitter: mean {t['epoch_jitter_ms']['mean']} ms, "
             f"σ {t['epoch_jitter_ms']['std']} ms, max |{t['epoch_jitter_ms']['max_abs']}| ms")
    L.append(f"- Receiver time accuracy estimate (tAcc): mean {t['tAcc_ns']['mean']} ns "
             f"({t['tAcc_ns']['min']}-{t['tAcc_ns']['max']} ns)")
    L.append(f"- Leap seconds: {t['leap_seconds']['rawx']} (valid={t['leap_seconds']['valid_flag']})")

    L.append("\n## RF health\n")
    L.append("| Block | jamInd mean/max | jamState≠0 | noise/ms | AGC | Antenna |")
    L.append("|---|---|---|---|---|---|")
    for b, v in d["rf_health"].items():
        L.append(f"| {b} | {v['jamInd_mean']}/{v['jamInd_max']} | {v['jammingState_nonzero']} "
                 f"| {v['noisePerMS_mean']} | {v['agcCnt_mean']:.0f} | {v['antStatus']} |")

    if d["warnings"]:
        L.append("\n## Warnings\n")
        for w in d["warnings"]:
            L.append(f"- ⚠ {w}")
    (outdir / "report.md").write_text("\n".join(L) + "\n")


# ----------------------------------------------------------------------------

def write_html_report(outdir, d, figs=()):
    """Self-contained report.html: summary cards, tables, embedded figures."""
    import base64
    import html as H

    def esc(x):
        return H.escape(str(x))

    def table(headers, rows):
        th = "".join(f"<th>{esc(h)}</th>" for h in headers)
        trs = "".join("<tr>" + "".join(f"<td>{c}</td>" for c in r) + "</tr>"
                      for r in rows)
        return f"<table><thead><tr>{th}</tr></thead><tbody>{trs}</tbody></table>"

    s = d["session"]
    pe = d["position_error"]
    ref = pe["reference"]
    best_name = s.get("track_solution_source") or "onboard"
    best = pe.get(best_name, pe.get("onboard", {}))
    val = d.get("validation_vs_onboard", {})
    val_ok = all(v["verdict"] == "PASS" for v in val.values()) if val else None

    cards = [
        ("Best solution", esc(best_name)),
        ("CEP50", f"{best.get('cep50_m', float('nan')):.2f} m" if best else "-"),
        ("Epochs", f"{s.get('message_counts', {}).get('NAV-PVT', best.get('epochs', 0))}"),
        ("Duration", f"{s['duration_min']} min"),
        ("Validation", "-" if val_ok is None else ("all PASS" if val_ok else "FAIL")),
    ]
    if "fix_rate" in best:
        cards.append(("PPK fix rate", f"{best['fix_rate']*100:.0f}%"))

    mod = ", ".join(x for x in (s.get("receiver") or []) if "MOD=" in x) or "u-blox"
    parts = [f"""<style>
:root {{ --surface:#fcfcfb; --page:#f9f9f7; --ink:#0b0b0b; --ink2:#52514e;
        --muted:#898781; --grid:#e1e0d9; --accent:#2a78d6; --border:rgba(11,11,11,.10); }}
@media (prefers-color-scheme: dark) {{ :root:not([data-theme="light"]) {{
  --surface:#1a1a19; --page:#0d0d0d; --ink:#fff; --ink2:#c3c2b7;
  --grid:#2c2c2a; --border:rgba(255,255,255,.10); }} }}
:root[data-theme="dark"] {{ --surface:#1a1a19; --page:#0d0d0d; --ink:#fff;
  --ink2:#c3c2b7; --grid:#2c2c2a; --border:rgba(255,255,255,.10); }}
body {{ background:var(--page); color:var(--ink); margin:0 auto; max-width:1000px;
       padding:2rem 1.2rem; font:15px/1.5 system-ui,-apple-system,"Segoe UI",sans-serif; }}
h1 {{ font-size:1.5rem; margin:.2rem 0 .3rem; }}
h2 {{ font-size:1.1rem; margin:2rem 0 .6rem; }}
.sub {{ color:var(--ink2); font-size:.9rem; margin-bottom:1.2rem; }}
.cards {{ display:flex; flex-wrap:wrap; gap:.7rem; margin:1rem 0; }}
.card {{ background:var(--surface); border:1px solid var(--border); border-radius:10px;
        padding:.7rem 1rem; min-width:8.5rem; }}
.card .k {{ color:var(--muted); font-size:.75rem; text-transform:uppercase;
           letter-spacing:.04em; }}
.card .v {{ font-size:1.35rem; font-weight:650; }}
table {{ border-collapse:collapse; width:100%; margin:.6rem 0 1rem;
        font-variant-numeric:tabular-nums; font-size:.87rem; }}
th {{ text-align:left; color:var(--ink2); font-weight:600; }}
th,td {{ padding:.35rem .6rem; border-bottom:1px solid var(--grid); }}
.pass {{ color:#006300; font-weight:650; }} .fail {{ color:#d03b3b; font-weight:650; }}
figure {{ background:#fcfcfb; border:1px solid var(--border); border-radius:10px;
         padding:.8rem; margin:1rem 0; overflow-x:auto; }}
figure img {{ max-width:100%; height:auto; display:block; margin:0 auto; }}
figcaption {{ color:var(--ink2); font-size:.85rem; padding-top:.5rem; }}
.warn {{ background:var(--surface); border-left:3px solid #fab219; padding:.5rem .8rem;
        margin:.4rem 0; border-radius:0 8px 8px 0; font-size:.9rem; }}
code {{ color:var(--ink2); }}
</style>
<h1>GPS Session Report</h1>
<div class="sub">{esc(Path(s['file']).name)} &middot; {esc(mod)} &middot;
{esc(s.get('module_fw'))}<br>{esc(s['start_utc'])} &rarr; {esc(s['end_utc'])}
(GPS week {s['gps_week']}) &middot; track source: <code>{esc(best_name)}</code></div>"""]

    parts.append('<div class="cards">' + "".join(
        f'<div class="card"><div class="k">{esc(k)}</div><div class="v">{esc(v)}</div></div>'
        for k, v in cards) + "</div>")

    motion = s.get("motion_mode", "?")
    motion_note = ("" if motion.startswith("static") else
                   " &middot; dispersion figures include any real platform motion")
    parts.append(f"<h2>Position summary</h2><div class='sub'>Motion mode: "
                 f"<b>{esc(motion)}</b>{motion_note}<br>Reference "
                 f"({esc(ref['source'])}): <b>{ref['lat']:.9f}&deg;, {ref['lon']:.9f}&deg;"
                 f"</b>, h = {ref['ellip_height_m']:.3f} m</div>")
    rows = []
    for name in ("ppk_kinematic", "ppk_static", "ppp", "single", "onboard"):
        st = pe.get(name)
        if not st:
            continue
        rows.append([esc(name), st["epochs"],
                     "/".join(f"{v:+.3f}" for v in st["offset_from_ref_NEU_m"]),
                     "/".join(f"{v:.3f}" for v in st["rms_NEU_m"]),
                     f"{st['cep50_m']:.3f}", f"{st['cep95_m']:.3f}",
                     f"{st['fix_rate']*100:.1f}%" if "fix_rate" in st else "-"])
    parts.append(table(["Solution", "Epochs", "Offset from ref N/E/U (m)",
                        "RMS N/E/U (m)", "CEP50 (m)", "CEP95 (m)", "Fix rate"], rows))

    if val:
        parts.append("<h2>Validation against onboard (NAV-PVT)</h2>")
        rows = [[esc(n), v["matched_epochs"],
                 "/".join(f"{x:+.3f}" for x in v["mean_offset_NEU_m"]),
                 f"{v['horizontal_offset_m']} m", f"&lt; {v['threshold_m']} m",
                 f"<span class='{v['verdict'].lower()}'>{v['verdict']}</span>"]
                for n, v in val.items()]
        parts.append(table(["Solution", "Matched epochs", "Mean offset N/E/U (m)",
                            "Horiz offset", "Threshold", "Verdict"], rows))

    parts.append("<h2>Constellations</h2>")
    rows = []
    for name, c in d["constellations"].items():
        bands = ", ".join(f"{b}: {v['cno_mean']}" for b, v in sorted(c["bands"].items()))
        rows.append([esc(name), c["satellites_tracked"], c["signals_tracked"],
                     c["measurements"], esc(bands), c["cycle_slip_events"],
                     c["sfrbx_frames"], c["prRes_rms_m"] or "-"])
    parts.append(table(["GNSS", "Sats", "Signals", "Meas", "C/N0 by band (dB-Hz)",
                        "Slips", "SFRBX", "prRes RMS (m)"], rows))

    io = d["ionosphere"]
    parts.append("<h2>Ionosphere</h2><ul>")
    for k, v in io["broadcast_klobuchar"].items():
        parts.append(f"<li>Broadcast <code>{esc(k)}</code>: {esc(v)}</li>")
    parts.append(f"<li>Signal iono-model census: {esc(io['ionoModel_census'])}</li>")
    gf = io["geometry_free_per_satellite"]
    if gf:
        med_m = float(np.median([r["iono_L1_m_median"] for r in gf]))
        med_t = float(np.median([r["TECU_median"] for r in gf]))
        parts.append(f"<li>Dual-frequency geometry-free slant delay, median across "
                     f"{len(gf)} satellites: <b>{med_m:.2f} m at L1</b> "
                     f"(~{med_t:.0f} TECU)</li>")
    parts.append("</ul>")

    t = d["timing"]
    parts.append("<h2>Timing</h2><ul>")
    parts.append(f"<li>Epoch jitter: mean {t['epoch_jitter_ms']['mean']} ms, "
                 f"&sigma; {t['epoch_jitter_ms']['std']} ms, "
                 f"max |{t['epoch_jitter_ms']['max_abs']}| ms</li>")
    parts.append(f"<li>Receiver time accuracy (tAcc): mean {t['tAcc_ns']['mean']} ns "
                 f"({t['tAcc_ns']['min']}&ndash;{t['tAcc_ns']['max']} ns)</li>")
    parts.append(f"<li>Leap seconds: {t['leap_seconds']['rawx']} "
                 f"(valid={t['leap_seconds']['valid_flag']})</li></ul>")

    parts.append("<h2>RF health</h2>")
    rows = [[esc(b), f"{v['jamInd_mean']}/{v['jamInd_max']}", v["jammingState_nonzero"],
             v["noisePerMS_mean"], f"{v['agcCnt_mean']:.0f}", esc(v["antStatus"])]
            for b, v in d["rf_health"].items()]
    parts.append(table(["Block", "jamInd mean/max", "jamState&ne;0", "noise/ms",
                        "AGC", "Antenna"], rows))

    if d["warnings"]:
        parts.append("<h2>Warnings</h2>")
        for w in d["warnings"]:
            parts.append(f"<div class='warn'>{esc(w)}</div>")

    if figs:
        parts.append("<h2>Figures</h2>")
        for name, title in figs:
            p = outdir / "figures" / name
            if p.is_file():
                b64 = base64.b64encode(p.read_bytes()).decode()
                parts.append(f"<figure><img src='data:image/png;base64,{b64}' "
                             f"alt='{esc(title)}'><figcaption>{esc(title)}"
                             f"</figcaption></figure>")

    (outdir / "report.html").write_text(
        f"<!doctype html><html><head><meta charset='utf-8'>"
        f"<meta name='viewport' content='width=device-width,initial-scale=1'>"
        f"<title>GPS Report - {H.escape(Path(s['file']).stem)}</title></head>"
        f"<body>{''.join(parts)}</body></html>")


def main():
    ap = argparse.ArgumentParser(description=__doc__,
                                 formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("capture", type=Path, help="input .ubx capture")
    ap.add_argument("--interval", type=int, default=1, help="track.csv output interval, s")
    ap.add_argument("--static", action="store_true",
                    help="antenna was VERIFIABLY stationary: enables static solutions and "
                         "static error framing. Never set for buoy/drifting deployments, "
                         "slow motion within a few metres is indistinguishable from noise "
                         "and would be averaged away.")
    ap.add_argument("--no-ppk", action="store_true")
    ap.add_argument("--ppp", action="store_true", help="also run PPP vs open IGS products")
    ap.add_argument("--base", help="force a CORS station id (e.g. MSEV)")
    ap.add_argument("--max-base-km", type=float, default=150.0)
    ap.add_argument("--force-stage", action="append", default=[],
                    choices=["rinex", "single", "ppk", "ppp", "track"])
    args = ap.parse_args()

    if not args.capture.is_file():
        sys.exit(f"no such file: {args.capture}")
    for tool in ("convbin", "rnx2rtkp", "pos2kml"):
        if not (BIN / tool).is_file():
            sys.exit(f"missing tools/bin/{tool} - run setup.sh first")

    outdir = SCRIPT_DIR / "processed" / args.capture.stem
    outdir.mkdir(parents=True, exist_ok=True)
    # --force-stage: remove that stage's outputs so it re-runs
    force_files = {"rinex": ["rover.obs", "rover.nav"], "single": ["single.pos"],
                   "ppk": ["ppk_kinematic.pos", "ppk_static.pos"],
                   "ppp": ["ppp.pos"], "track": ["track.csv", "track.kml"]}
    for st in args.force_stage:
        for fn in force_files[st]:
            (outdir / fn).unlink(missing_ok=True)

    warnings = []
    data = parse_ubx(args.capture)
    n_rawx = len(data["rawx_hdr"])
    log(f"inventory: {sum(data['counts'].values())} UBX messages, {n_rawx} RAWX epochs, "
        f"{len(data['pvt'])} PVT epochs")

    if n_rawx == 0:
        warnings.append("NO RAW OBSERVATIONS (RXM-RAWX) in this capture - it is a "
                        "debug/NMEA-style log. RINEX conversion and corrected positioning "
                        "are impossible; diagnostics below cover onboard data only.")
        solutions = {}
        if len(data["pvt"]):
            diag = stage_diagnostics(outdir, data, {"iono": {}}, solutions, {}, None,
                                     warnings, args)
            write_report(outdir, diag)
            write_html_report(outdir, diag)
        else:
            log("no NAV-PVT either; nothing to report")
        sys.exit(0)

    rinex_info = stage_rinex(args.capture, outdir, n_rawx, warnings)
    solutions = {}
    solutions["single"] = stage_single(outdir, rinex_info["obs"], rinex_info["nav"],
                                       write_conf(outdir / "single.conf",
                                                  ant1=rover_antenna(),
                                                  atx=None))
    ppk_result = {}
    if not args.no_ppk:
        ppk_result = stage_ppk(outdir, rinex_info["obs"], rinex_info["nav"], data, args, warnings)
        for name, df in (ppk_result.get("solutions") or {}).items():
            solutions[f"ppk_{name}"] = df
    if args.ppp:
        ppp_result = stage_ppp(outdir, rinex_info["obs"], rinex_info["nav"], data, args,
                               warnings)
        if ppp_result.get("status") == "ok":
            solutions["ppp"] = ppp_result["df"]

    track_df, track_source = stage_track(outdir, solutions, data["pvt"], args.interval)
    diag = stage_diagnostics(outdir, data, rinex_info, solutions, ppk_result, track_source,
                             warnings, args)
    figs = stage_figures(outdir, data, solutions, track_source, diag)
    write_report(outdir, diag, figs)
    write_html_report(outdir, diag, figs)
    log("report: report.md + report.html written")

    log(f"done -> {outdir}")
    if ppk_result.get("status") == "pending":
        log("NOTE: " + ppk_result["message"])


if __name__ == "__main__":
    main()