commit 65a483660822e85d68d3f9a9417d07f858c84a0e
Author: = <=>
Date:   Sat Jul 25 14:36:41 2026 -0400

    Initial commit of fully working GPS processing pipeline. See README for more information.

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..e9fb985
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,18 @@
+# Python environment
+venv/
+__pycache__/
+*.pyc
+
+# Bootstrapped tooling: RTKLIB build, ANTEX models, CORS catalog cache
+# (fully reproducible with setup.sh)
+tools/
+
+# Pipeline outputs (regenerate with process_gps.py)
+processed/
+
+# Live streaming session outputs (stream_gps.py)
+stream_*.csv
+
+# Raw GPS capture data
+*.ubx
+*.uc2x
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..9818468
--- /dev/null
+++ b/README.md
@@ -0,0 +1,193 @@
+# u-blox GPS Module Data Procesing (u-blox ZED-X20P)
+
+Tooling for evaluating u-blox ZED-X20P all-band GNSS modules, targeting **buoy
+deployments**. Converts raw u-center capture files (`.ubx`) into corrected
+positions and produces a full diagnostics report, and supports live streaming
+with real-time corrections.
+
+Two operating modes:
+
+| Mode | Script | Accuracy | Correction source |
+|---|---|---|---|
+| Post-processing | `process_gps.py` | dm-cm (PPK) | NGS CORS base data, auto-fetched (free, open) |
+| Real-time streaming | `stream_gps.py` | 1-3 cm (RTK fix) | GCGC RTN NTRIP (free registration) or Galileo HAS (FW ≥ HPG 2.10) |
+
+Both consume/produce standard formats (RINEX 3, `.pos`, CSV, KML) via RTKLIB
+(demo5 fork) and widely used open-source Python libraries (pyubx2, pygnssutils,
+numpy, pandas, matplotlib).
+
+## Motion policy (important)
+
+**Processing never assumes the antenna was stationary.** Buoys drift slowly and
+drift within a few metres is statistically indistinguishable from measurement
+noise. Everything defaults to kinematic; pass `--static` to `process_gps.py` 
+only when the antenna is independently known to have been fixed (e.g. bench/roof test).
+
+## Quick start
+
+```bash
+bash setup.sh                                    # one-time environment bootstrap
+venv/bin/python process_gps.py <capture.ubx>     # process a capture
+# results in processed/<capture-stem>/  (report.html is the human summary)
+```
+
+## Requirements
+
+Linux with `python3` (3.10+), `gcc`/`make`, `git`, `curl`, and network access
+for `setup.sh` downloads. No sudo required and everything installs into the
+project directory (`venv/`, `tools/`).
+
+---
+
+## setup.sh
+
+Environment bootstrap (safe to re-run; each step skips if already satisfied):
+
+1. **Python venv** (`venv/`): pyubx2, pygnssutils, pyserial, numpy, pandas,
+   pyproj, matplotlib.
+2. **RTKLIB demo5 v2.5.1** built from source into `tools/bin/`
+   (`convbin`, `rnx2rtkp`, `pos2kml`). The demo5 fork is required: stock
+   RTKLIB 2.4.3 silently drops all BeiDou and GPS L5 / Galileo E5 observations
+   from this receiver (~30% of the measurements). If the build fails, setup
+   falls back to the stock apt package and creates a `tools/STOCK_RTKLIB`
+   marker — `process_gps.py` then adds a prominent data-loss warning to every
+   report until a demo5 build succeeds.
+3. **Antenna calibration models** into `tools/antex/`: `igs20.atx` (IGS, covers
+   CORS base-station antennas) and `ngs20.atx` (NGS absolute calibrations,
+   covers the rover antenna). The rover antenna is resolved automatically and
+   written to `tools/antenna.json` (see *Antenna configuration* below).
+4. **CRX2RNX** (Hatanaka RINEX decompressor) — fallback for CORS stations that
+   only publish compressed `.d` observation files.
+5. **NGS CORS station catalog** cache (used for nearest-base auto-selection).
+6. Self-check summary of every tool and the resolved antenna.
+
+## process_gps.py
+
+```
+venv/bin/python process_gps.py <capture.ubx> [options]
+```
+
+Stages (idempotent — outputs are cached, so re-running later only fills in
+what was previously unavailable, e.g. CORS/IGS data that had not been
+published yet):
+
+1. **inventory** — parses the `.ubx` (RXM-RAWX/SFRBX raw observations,
+   NAV-PVT/SAT/SIG, MON-RF) and the `.uc2x` sidecar's session metadata.
+   A capture without raw observations (e.g. u-center debug-mode logs) gets a
+   diagnostics-only report and a clear notice.
+2. **rinex** — RINEX 3.04 obs+nav via demo5 `convbin`; verifies epoch count and
+   BDS/L5 presence; extracts broadcast ionosphere parameters.
+3. **single** — standalone single-point solution (~1 m, works offline).
+4. **ppk** — differential post-processing against the nearest operational NGS
+   CORS station (auto-selected by distance from the capture's own mean
+   position; MSEV at 34 km for the Hattiesburg test site). Hourly base files
+   are fetched automatically (~1 h publication lag; if not yet published the
+   stage reports "re-run later" and everything else proceeds). Base antenna
+   type is read from the base RINEX header; base coordinates are ITRF2020,
+   velocity-propagated to the capture epoch.
+5. **ppp** (`--ppp`) — precise point positioning against IGS rapid/final
+   orbit+clock products (anonymous download from BKG/ESA, ~1 day lag).
+   Kinematic PPP by default per the motion policy; note RTKLIB's PPP-kinematic
+   engine is weak and may yield no solution — PPK is the reference method.
+6. **track** — `track.csv` (UTC + decimal-degree lat/lon + per-epoch quality
+   and 1-σ sd columns + onboard comparison columns) and `track.kml`.
+7. **diagnostics** — `diagnostics.json`, `report.md`, `report.html`
+   (self-contained, figures embedded) covering: constellation metrics (per-band
+   C/N0, cycle slips, residuals), ionosphere (broadcast Klobuchar + measured
+   dual-frequency slant delay/TECU), positional error metrics (RMS/R50/R95
+   dispersion, formal σ, inter-solution offsets), **validation of every
+   computed solution against the receiver's own onboard NAV-PVT positions**
+   (epoch-matched, PASS/FAIL), timing (epoch jitter, tAcc, leap seconds), and
+   RF health (jamming indicator, AGC, antenna supervisor state).
+
+### Options
+
+| Flag | Meaning |
+|---|---|
+| `--interval N` | `track.csv` output interval in seconds (default 1 = native rate) |
+| `--static` | antenna verifiably stationary: enables static PPK/PPP solutions and static error framing. **Never for buoys.** |
+| `--no-ppk` | skip base-station download/processing |
+| `--ppp` | additionally run PPP against IGS precise products |
+| `--base SSSS` | force a specific CORS station id (e.g. `MSEV`) |
+| `--max-base-km KM` | base-station search radius (default 150) |
+| `--force-stage S` | delete and re-run one stage (`rinex`,`single`,`ppk`,`ppp`,`track`) |
+
+## Antenna configuration
+
+Receiver-antenna phase-center corrections (PCO/PCV) matter at the cm level and
+are applied during PPK/PPP via ANTEX calibration files.
+
+- **Automatic**: `setup.sh` downloads `ngs20.atx`/`igs20.atx` and searches for
+  the project's antenna (ArduSimple calibrated quad-band, NGS calibration
+  AS-ANT3BCAL01). The resolved 20-character ANTEX type string is written to
+  `tools/antenna.json`:
+
+  ```json
+  { "antex_name": "AS-ANT3BCAL     NONE", "source": "ngs20.atx" }
+  ```
+
+- **Different antenna**: edit `antex_name` in `tools/antenna.json` to the exact
+  ANTEX type string of your antenna (padding/spacing matters, copy it verbatim
+  from the `TYPE / SERIAL NO` line of an antenna block in
+  `tools/antex/ngs20.atx` or `igs20.atx`). Browse available calibrations at
+  https://geodesy.noaa.gov/ANTCAL/ . If your antenna has an NGS calibration but
+  a different model string, also update the `ANT3B` search pattern in `setup.sh`
+  so future setups re-resolve it automatically.
+- **No calibration available**: set `antex_name` to `null` (or leave it
+  unresolved) — processing continues without rover PCO/PCV and the report
+  carries a warning. Expect a few cm of systematic error, mostly in height.
+- **Base antenna** needs no configuration: its type is read from the CORS RINEX
+  header and looked up in the same ANTEX files.
+- **Antenna height**: assumed 0 (positions refer to the antenna reference
+  point). For a surveyed monument, subtract your ARP offset downstream.
+
+## Other important settings and inputs
+
+**Capture requirements.** Input `.ubx` files must contain raw observations:
+UBX-RXM-RAWX + UBX-RXM-SFRBX enabled (plus NAV-PVT; NAV-SAT/NAV-SIG/MON-RF
+enrich diagnostics). In u-center, do **not** record with the "debug messages"
+option — it floods the log with undocumented TRK/TUN/DBG/SEC messages and, in
+past sessions here, raw output was configured in flash but not active in RAM.
+The optional `.uc2x` sidecar (u-center 2 index) is only used for its session
+metadata header.
+
+**Real-time credentials (`stream_gps.py`).** Free registration at
+http://rtn.usm.edu/RegisterAccount.aspx (Mississippi GCGC RTN), then:
+
+```bash
+export GCGC_USER=<username> GCGC_PASS=<password>
+venv/bin/python stream_gps.py --port /dev/ttyUSB0     # or COM7 on Windows
+```
+
+Key streaming defaults (buoy-oriented): receiver dynamic-platform model `sea`
+(`--dynmodel` to override), **live** NMEA GGA uplink so the network-RTK virtual
+reference follows platform drift (`--gga-fixed` for bench tests), raw
+RXM-RAWX/SFRBX always logged so every session is post-processable, and
+`--replay <file.ubx>` offline test mode. Galileo HAS (~20 cm, no internet)
+activates automatically as fallback when firmware ≥ HPG 2.10 is detected.
+
+**Datums.** PPK output is in the CORS base frame (ITRF2020, current epoch);
+GCGC real-time corrections are NAD83(2011) epoch 2010.0 — a constant ~1.5 m
+offset from ITRF/WGS84 in CONUS. Each CSV records its frame in the header;
+comparisons in the report account for the reference used. Onboard receiver
+positions are WGS84-aligned.
+
+**Data-source timing.** CORS hourly files: ~1 h lag, expire after ~2 days
+(daily files remain, decimated to 30 s after 30 days). IGS rapid products:
+~1 day lag. The pipeline caches everything it fetches under
+`processed/<stem>/` and tells you when a re-run will find more data.
+
+## Repository layout
+
+```
+setup.sh                 environment bootstrap (run first)
+process_gps.py           post-processing pipeline (Mode 2)
+stream_gps.py            real-time streaming client (Mode 1)
+run.sh                   convenience wrapper
+notes/                   manual workflow notes + combine_rinex.sh helper
+tools/                   built binaries, ANTEX models, caches   (generated)
+venv/                    Python environment                     (generated)
+processed/<capture>/     per-capture outputs: RINEX, .pos solutions,
+                         track.csv, track.kml, figures/, diagnostics.json,
+                         report.md, report.html                 (generated)
+```
diff --git a/notes/RTK Processing Notes.txt b/notes/RTK Processing Notes.txt
new file mode 100755
index 0000000..6e0dca1
--- /dev/null
+++ b/notes/RTK Processing Notes.txt	
@@ -0,0 +1,42 @@
+### THESE ARE NOTES FROM PREVIOUS WORK, KEPT FOR REFERENCE. NONE OF THESE STEPS ARE REQUIRED TO USE THE GPS DATA PROCESSING SCRIPTS IN THE PARENT DIRECTORY.
+
+## Download the CRX2RNX tool precompiled binary
+
+# Download the pre-compiled Linux binary
+cd /tmp
+wget https://terras.gsi.go.jp/ja/crx2rnx/RNXCMP_4.2.0_Linux_gcc_64bit.tar.gz
+
+# Extract
+tar -xzf RNXCMP_4.2.0_Linux_gcc_64bit.tar.gz
+
+# Install to system path
+sudo cp RNXCMP_4.2.0_Linux_gcc_64bit/CRX2RNX /usr/local/bin/
+
+# Make executable
+sudo chmod +x /usr/local/bin/CRX2RNX
+
+# Create lowercase symlink (rnx2rtkp looks for lowercase)
+sudo ln -s /usr/local/bin/CRX2RNX /usr/local/bin/crx2rnx
+
+# Test
+crx2rnx -h
+
+
+## Download CORS corrections from here: https://www.ngs.noaa.gov/CORS/
+
+## If base station CORS files are compressed (.25d format), decompress them
+ls msin*.25d | xargs -n 1 -P 4 crx2rnx
+
+## Combine base station observations into one file
+./combine_rinex.sh msin_combined.obs msin*.25o
+
+## Perform corrections on a .obs file using base station files as reference
+cd CORS\ STATION\ MSIN\ NOV\ 2025/
+
+rnx2rtkp \ -t -u \
+  -o ../../A/DATA/result.pos \
+  ../../A/DATA/20251120_033544.006_UTC_GPS_A_combined.obs \
+  msin_combined.obs \
+  ../../A/DATA/20251120_033544.006_UTC_GPS_A_combined.nav
+
+## The result.pos file now contains PPK corrected GPS data in LLH format
\ No newline at end of file
diff --git a/notes/combine_rinex.sh b/notes/combine_rinex.sh
new file mode 100755
index 0000000..52e9837
--- /dev/null
+++ b/notes/combine_rinex.sh
@@ -0,0 +1,141 @@
+#!/bin/bash
+# combine_rinex.sh - Combine RINEX observation files (supports both 2.x and 3.x)
+
+if [ $# -lt 2 ]; then
+    echo "Usage: $0 output.obs input_files..."
+    echo "Example: $0 combined.obs msin*.25o"
+    echo "Example: $0 combined.obs msin3240.25o msin3250.25o msin3260.25o"
+    exit 1
+fi
+
+output="$1"
+shift
+
+# Sort files chronologically (RINEX naming is sortable)
+files=$(printf '%s\n' "$@" | sort)
+
+if [ -z "$files" ]; then
+    echo "Error: No input files specified"
+    exit 1
+fi
+
+# Get first file
+first_file=$(echo "$files" | head -1)
+
+# Verify all files exist
+echo "Verifying input files..."
+for file in $files; do
+    if [ ! -f "$file" ]; then
+        echo "Error: File not found: $file"
+        exit 1
+    fi
+done
+
+file_count=$(echo "$files" | wc -w)
+echo "Found $file_count files to combine"
+echo ""
+
+# Get header from first file
+header_lines=$(grep -n "END OF HEADER" "$first_file" | cut -d: -f1)
+
+if [ -z "$header_lines" ]; then
+    echo "Error: No RINEX header found in $first_file"
+    exit 1
+fi
+
+# Detect RINEX version from first line
+rinex_version=$(head -1 "$first_file" | awk '{print substr($1,1,1)}')
+
+if [ "$rinex_version" = "3" ]; then
+    echo "Detected: RINEX 3.x format"
+    epoch_pattern="^>"
+elif [ "$rinex_version" = "2" ]; then
+    echo "Detected: RINEX 2.x format"
+    epoch_pattern="^ [0-9][0-9] [ 0-9][0-9] [ 0-9][0-9]"
+else
+    echo "Warning: Could not detect RINEX version, assuming 2.x"
+    epoch_pattern="^ [0-9][0-9] [ 0-9][0-9] [ 0-9][0-9]"
+fi
+
+echo ""
+
+# Copy header from first file
+head -n $header_lines "$first_file" > "$output"
+echo "✓ Added header from $first_file"
+
+# Append data from all files (skip their headers)
+total_lines=0
+for file in $files; do
+    data_start=$(($(grep -n "END OF HEADER" "$file" | cut -d: -f1) + 1))
+    data_lines=$(tail -n +$data_start "$file" | wc -l)
+    tail -n +$data_start "$file" >> "$output"
+    total_lines=$((total_lines + data_lines))
+    echo "✓ Added data from $file ($data_lines lines)"
+done
+
+echo ""
+
+# Count epochs based on detected version
+if [ "$rinex_version" = "3" ]; then
+    epochs=$(grep -c "^>" "$output" 2>/dev/null || echo "0")
+else
+    epochs=$(grep -c "^ [0-9][0-9] [ 0-9][0-9] [ 0-9][0-9]" "$output" 2>/dev/null || echo "0")
+fi
+
+# Get file size
+file_size=$(ls -lh "$output" | awk '{print $5}')
+
+echo "========================================="
+echo "✓ Combined file created: $output"
+echo "  RINEX Version: $rinex_version.x"
+echo "  Size: $file_size"
+echo "  Files combined: $file_count"
+echo "  Total data lines: $total_lines"
+echo "  Total epochs: $epochs"
+echo "========================================="
+
+# Verify reasonable results
+if [ $epochs -eq 0 ]; then
+    echo ""
+    echo "⚠ WARNING: No epochs detected!"
+    echo "  This may indicate a problem with the combined file."
+    echo "  Checking first few lines of output..."
+    echo ""
+    head -40 "$output"
+    exit 1
+fi
+
+# Show time range if epochs found
+if [ $epochs -gt 0 ]; then
+    echo ""
+    echo "Time range:"
+
+    if [ "$rinex_version" = "3" ]; then
+        # RINEX 3.x format: > YYYY MM DD HH MM SS.SSSSSSS
+        first_epoch=$(grep "^>" "$output" | head -1)
+        last_epoch=$(grep "^>" "$output" | tail -1)
+        first_time=$(echo "$first_epoch" | awk '{printf "%s-%02d-%02d %02d:%02d:%06.3f", $2, $3, $4, $5, $6, $7}')
+        last_time=$(echo "$last_epoch" | awk '{printf "%s-%02d-%02d %02d:%02d:%06.3f", $2, $3, $4, $5, $6, $7}')
+    else
+        # RINEX 2.x format:  YY MM DD HH MM SS.SSSSSSS
+        first_epoch=$(grep "^ [0-9][0-9] [ 0-9][0-9] [ 0-9][0-9]" "$output" | head -1)
+        last_epoch=$(grep "^ [0-9][0-9] [ 0-9][0-9] [ 0-9][0-9]" "$output" | tail -1)
+        first_time=$(echo "$first_epoch" | awk '{printf "20%02d-%02d-%02d %02d:%02d:%011.7f", $1, $2, $3, $4, $5, $6}')
+        last_time=$(echo "$last_epoch" | awk '{printf "20%02d-%02d-%02d %02d:%02d:%011.7f", $1, $2, $3, $4, $5, $6}')
+    fi
+
+    echo "  First: $first_time"
+    echo "  Last:  $last_time"
+fi
+
+# Verify single header (no corruption)
+header_count=$(grep -c "END OF HEADER" "$output")
+if [ $header_count -ne 1 ]; then
+    echo ""
+    echo "⚠ WARNING: Found $header_count headers (expected 1)"
+    echo "  The file may be corrupted!"
+    exit 1
+fi
+
+echo ""
+echo "✓ File validation passed"
diff --git a/process_gps.py b/process_gps.py
new file mode 100644
index 0000000..dcf854b
--- /dev/null
+++ b/process_gps.py
@@ -0,0 +1,1584 @@
+#!/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()
diff --git a/run.sh b/run.sh
new file mode 100755
index 0000000..91950c2
--- /dev/null
+++ b/run.sh
@@ -0,0 +1,3 @@
+source ./venv/bin/activate
+python3 process_gps.py -h
+python3 process_gps.py 2026-7-23_124722_serial-COM7.ubx
diff --git a/setup.sh b/setup.sh
new file mode 100755
index 0000000..c6db532
--- /dev/null
+++ b/setup.sh
@@ -0,0 +1,193 @@
+#!/usr/bin/env bash
+# setup.sh - One-time environment bootstrap for the GPS processing pipeline.
+#
+# Installs (all local to this directory, no sudo required):
+#   venv/          Python env: pyubx2, pygnssutils, pyserial, numpy, pandas, pyproj
+#   tools/bin/     RTKLIB demo5 v2.5.1 (convbin, rnx2rtkp, pos2kml)  [stock apt fallback]
+#   tools/antex/   igs20.atx + ngs20.atx antenna models; tools/antenna.json (rover antenna)
+#   tools/bin/CRX2RNX   Hatanaka RINEX decompressor (fallback for .d CORS files)
+#   tools/cache/   NGS CORS station catalog
+#
+# Idempotent: re-running skips anything already installed.
+
+set -uo pipefail
+
+SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+cd "$SCRIPT_DIR"
+
+TOOLS="$SCRIPT_DIR/tools"
+VENV="$SCRIPT_DIR/venv"
+DEMO5_TAG="v2.5.1"
+DEMO5_REPO="https://github.com/rtklibexplorer/RTKLIB"
+
+PASS=()
+WARN=()
+
+say()  { printf '\n=== %s ===\n' "$*"; }
+ok()   { printf '  [ok] %s\n' "$*"; PASS+=("$*"); }
+warn() { printf '  [WARN] %s\n' "$*"; WARN+=("$*"); }
+skip() { printf '  [skip] %s (already present)\n' "$*"; }
+
+fetch() { # fetch <url> <dest>
+    curl -fsSL --max-time 300 -o "$2" "$1"
+}
+
+# ---------------------------------------------------------------- 1. Python venv
+say "1/6 Python virtual environment"
+if [ -x "$VENV/bin/python" ]; then
+    skip "venv"
+else
+    python3 -m venv "$VENV" || { echo "FATAL: venv creation failed"; exit 1; }
+fi
+if "$VENV/bin/python" -c "import pyubx2, pygnssutils, serial, numpy, pandas, pyproj, matplotlib" 2>/dev/null; then
+    skip "python packages"
+else
+    "$VENV/bin/pip" install --quiet --upgrade pip
+    "$VENV/bin/pip" install --quiet pyubx2 pygnssutils pyserial numpy pandas pyproj matplotlib \
+        || { echo "FATAL: pip install failed"; exit 1; }
+fi
+"$VENV/bin/python" -c "import pyubx2, pygnssutils, serial, numpy, pandas, pyproj, matplotlib" \
+    && ok "venv with pyubx2/pygnssutils/pyserial/numpy/pandas/pyproj/matplotlib" \
+    || { echo "FATAL: package import check failed"; exit 1; }
+
+mkdir -p "$TOOLS/bin" "$TOOLS/antex" "$TOOLS/cache" "$TOOLS/src"
+
+# ---------------------------------------------------------------- 2. RTKLIB demo5
+say "2/6 RTKLIB demo5 $DEMO5_TAG (convbin, rnx2rtkp, pos2kml)"
+build_demo5_app() { # build_demo5_app <appname>
+    local app="$1"
+    local gccdir="$TOOLS/src/RTKLIB/app/consapp/$app/gcc"
+    [ -d "$gccdir" ] || return 1
+    # -lgfortran in some makefiles is vestigial; try without it first.
+    if ! make -C "$gccdir" -j"$(nproc)" LDLIBS="-lm" >"$TOOLS/src/build_$app.log" 2>&1; then
+        make -C "$gccdir" clean >/dev/null 2>&1
+        make -C "$gccdir" -j"$(nproc)" >>"$TOOLS/src/build_$app.log" 2>&1 || return 1
+    fi
+    cp "$gccdir/$app" "$TOOLS/bin/$app"
+}
+
+if [ -x "$TOOLS/bin/convbin" ] && [ -x "$TOOLS/bin/rnx2rtkp" ] && [ -x "$TOOLS/bin/pos2kml" ]; then
+    skip "RTKLIB binaries"
+    [ -f "$TOOLS/STOCK_RTKLIB" ] && warn "using STOCK RTKLIB (BDS/L5 signals will be missing from RINEX)"
+else
+    if [ ! -d "$TOOLS/src/RTKLIB" ]; then
+        git clone --depth 1 --branch "$DEMO5_TAG" "$DEMO5_REPO" "$TOOLS/src/RTKLIB" \
+            || warn "git clone of demo5 RTKLIB failed"
+    fi
+    built=true
+    for app in convbin rnx2rtkp pos2kml; do
+        if build_demo5_app "$app"; then
+            printf '  built %s\n' "$app"
+        else
+            warn "demo5 build failed for $app (see tools/src/build_$app.log)"
+            built=false
+        fi
+    done
+    if $built; then
+        rm -f "$TOOLS/STOCK_RTKLIB"
+        ok "RTKLIB demo5 $DEMO5_TAG built from source"
+    else
+        # Fallback: stock RTKLIB 2.4.3 from the Ubuntu archive (works, but silently
+        # drops BeiDou + L5/E5 observations from this receiver's raw data).
+        echo "  falling back to stock RTKLIB 2.4.3 from apt archive..."
+        tmpd="$(mktemp -d)"
+        if (cd "$tmpd" && apt-get download rtklib >/dev/null 2>&1 && dpkg -x rtklib_*.deb x); then
+            cp "$tmpd"/x/usr/bin/{convbin,rnx2rtkp,pos2kml} "$TOOLS/bin/"
+            touch "$TOOLS/STOCK_RTKLIB"
+            warn "STOCK RTKLIB installed - BDS/L5 will be MISSING from RINEX output"
+        else
+            warn "stock RTKLIB fallback also failed - pipeline cannot run"
+        fi
+        rm -rf "$tmpd"
+    fi
+fi
+for app in convbin rnx2rtkp pos2kml; do
+    [ -x "$TOOLS/bin/$app" ] || warn "$app missing from tools/bin"
+done
+
+# ---------------------------------------------------------------- 3. Antenna models
+say "3/6 Antenna calibration models (ANTEX)"
+if [ -s "$TOOLS/antex/igs20.atx" ]; then skip "igs20.atx"; else
+    fetch "https://files.igs.org/pub/station/general/igs20.atx" "$TOOLS/antex/igs20.atx" \
+        && ok "igs20.atx (IGS antenna models incl. CORS base antennas)" \
+        || warn "igs20.atx download failed"
+fi
+if [ -s "$TOOLS/antex/ngs20.atx" ]; then skip "ngs20.atx"; else
+    fetch "https://geodesy.noaa.gov/ANTCAL/LoadFile?file=ngs20.atx" "$TOOLS/antex/ngs20.atx" \
+        && ok "ngs20.atx (NGS absolute calibrations incl. ArduSimple)" \
+        || warn "ngs20.atx download failed"
+fi
+# Resolve the rover antenna (NGS calibration AS-ANT3BCAL01, ArduSimple) to its
+# exact 20-char ANTEX type string, consumed by process_gps.py as ant1-anttype.
+"$VENV/bin/python" - "$TOOLS" <<'EOF'
+import json, re, sys, pathlib
+tools = pathlib.Path(sys.argv[1])
+result = {"antex_name": None, "source": None, "candidates": []}
+for atx in ("ngs20.atx", "igs20.atx"):
+    p = tools / "antex" / atx
+    if not p.is_file():
+        continue
+    for line in p.open(encoding="utf-8", errors="replace"):
+        if "TYPE / SERIAL" in line and "ANT3B" in line[:20].upper():
+            name = line[:20].rstrip()
+            result["candidates"].append({"name": name, "file": atx})
+            if result["antex_name"] is None:
+                result["antex_name"], result["source"] = name, atx
+(tools / "antenna.json").write_text(json.dumps(result, indent=2))
+if result["antex_name"]:
+    print(f'  [ok] rover antenna resolved: "{result["antex_name"]}" ({result["source"]})')
+else:
+    print("  [WARN] no ANT3B antenna found in ANTEX files - PPK will run without rover PCO/PCV")
+EOF
+
+# ---------------------------------------------------------------- 4. CRX2RNX
+say "4/6 CRX2RNX (Hatanaka RINEX decompressor, fallback for .d CORS files)"
+if [ -x "$TOOLS/bin/CRX2RNX" ]; then
+    skip "CRX2RNX"
+else
+    tmpd="$(mktemp -d)"
+    if fetch "https://terras.gsi.go.jp/ja/crx2rnx/RNXCMP_4.2.0_Linux_gcc_64bit.tar.gz" "$tmpd/rnxcmp.tar.gz" \
+        && tar -xzf "$tmpd/rnxcmp.tar.gz" -C "$tmpd" \
+        && cp "$tmpd"/RNXCMP_*/bin/CRX2RNX "$TOOLS/bin/CRX2RNX" 2>/dev/null \
+            || cp "$tmpd"/RNXCMP_*/CRX2RNX "$TOOLS/bin/CRX2RNX" 2>/dev/null; then
+        chmod +x "$TOOLS/bin/CRX2RNX"
+        ok "CRX2RNX 4.2.0"
+    else
+        warn "CRX2RNX install failed (only needed for CORS stations without plain .o files)"
+    fi
+    rm -rf "$tmpd"
+fi
+
+# ---------------------------------------------------------------- 5. CORS catalog
+say "5/6 NGS CORS station catalog"
+if [ -s "$TOOLS/cache/itrf2020_geo.comp.txt" ]; then
+    skip "station catalog"
+else
+    fetch "https://geodesy.noaa.gov/corsdata/coord/coord_20/itrf2020_geo.comp.txt" \
+          "$TOOLS/cache/itrf2020_geo.comp.txt" \
+        && ok "CORS station catalog (ITRF2020 coordinates)" \
+        || warn "CORS catalog download failed (process_gps.py will retry at run time)"
+fi
+
+# ---------------------------------------------------------------- 6. Self-check
+say "6/6 Self-check summary"
+for app in convbin rnx2rtkp pos2kml; do
+    if [ -x "$TOOLS/bin/$app" ]; then
+        printf '  %-10s %s\n' "$app" "$("$TOOLS/bin/$app" -h 2>&1 | grep -m1 -iE 'synops|usage|version' || echo present)"
+    fi
+done
+[ -x "$TOOLS/bin/CRX2RNX" ] && printf '  %-10s %s\n' "CRX2RNX" "present"
+[ -s "$TOOLS/cache/itrf2020_geo.comp.txt" ] && printf '  %-10s %d stations\n' "catalog" "$(wc -l < "$TOOLS/cache/itrf2020_geo.comp.txt")"
+[ -s "$TOOLS/antenna.json" ] && printf '  %-10s %s\n' "antenna" "$(cat "$TOOLS/antenna.json" | "$VENV/bin/python" -c 'import json,sys; print(json.load(sys.stdin)["antex_name"])')"
+
+echo
+if [ ${#WARN[@]} -gt 0 ]; then
+    echo "Completed with ${#WARN[@]} warning(s):"
+    printf '  - %s\n' "${WARN[@]}"
+else
+    echo "Setup complete - no warnings."
+fi
+echo
+echo "Mode 1 (real-time RTK) additionally requires free registration at"
+echo "  http://rtn.usm.edu/RegisterAccount.aspx  (Mississippi GCGC RTN)"
+echo "then:  export GCGC_USER=<user> GCGC_PASS=<pass>"
diff --git a/stream_gps.py b/stream_gps.py
new file mode 100644
index 0000000..25a6bd5
--- /dev/null
+++ b/stream_gps.py
@@ -0,0 +1,446 @@
+#!/usr/bin/env python3
+"""
+stream_gps.py - Mode 1: continuous real-time positioning for the u-blox ZED-X20P.
+
+Best-available correction tier is selected automatically:
+  1. Network RTK via NTRIP  (default: Mississippi GCGC RTN, free registration at
+     http://rtn.usm.edu/RegisterAccount.aspx) -> onboard RTK engine, 1-3 cm.
+     Credentials: env GCGC_USER / GCGC_PASS (or --ntrip-user/--ntrip-pass).
+  2. Galileo HAS (E6-B, no internet needed)  -> ~20-25 cm.  Requires firmware
+     HPG >= 2.10 (auto-probed via MON-VER; this project's module reports 2.10).
+  3. Neither available -> the onboard multiband+SBAS solution (~1 m) cannot be
+     improved in real time; per project policy the enhanced mode is SKIPPED
+     (run with --log-only to still record a .ubx for post-processing).
+
+Every session simultaneously logs raw RXM-RAWX/SFRBX to a .ubx file, so a
+streaming session can always be post-processed later with process_gps.py.
+
+Buoy deployment defaults: the receiver's dynamic platform model is set to "sea"
+(--dynmodel to override) and the NTRIP GGA uplink sends the LIVE position so the
+VRS reference follows platform drift (--gga-fixed to disable for a fixed antenna).
+
+Outputs: live console status, stream_<ts>.csv (decimal-degree lat/lon),
+         stream_<ts>.ubx (raw capture).
+
+Usage:
+    venv/bin/python stream_gps.py --port /dev/ttyUSB0            (Linux)
+    venv/bin/python stream_gps.py --port COM7                    (Windows)
+    venv/bin/python stream_gps.py --replay capture.ubx           (offline test)
+
+Datum note: NTRIP corrections define the output frame. GCGC RTN broadcasts
+NAD83(2011) epoch 2010.00 - constant ~1.5 m from WGS84/ITRF in CONUS; the CSV
+header records which tier (and therefore frame) was active.
+"""
+
+import argparse
+import csv
+import os
+import queue
+import signal
+import socket
+import sys
+import threading
+import time
+from collections import Counter
+from datetime import datetime, timezone
+from pathlib import Path
+
+from pyubx2 import UBXMessage, UBXReader, UBX_PROTOCOL, SET_LAYER_RAM
+
+FIX_NAME = {0: "none", 1: "DR", 2: "2D", 3: "3D", 4: "GNSS+DR", 5: "time-only"}
+CARR_NAME = {0: "", 1: "RTK-FLOAT", 2: "RTK-FIXED"}
+
+# UBX-CFG keys (u-blox X20 HPG interface descriptions; NAVCOR from HPG 2.10 JSON)
+K = {
+    "UART1_BAUD":        0x40520001,
+    "UART1_IN_RTCM3":    0x10730004,
+    "USB_IN_RTCM3":      0x10770004,
+    "PVT_UART1":         0x20910007, "PVT_USB":       0x20910009,
+    "HPPOSLLH_UART1":    0x20910034, "HPPOSLLH_USB":  0x20910036,
+    "RAWX_UART1":        0x209102a5, "RAWX_USB":      0x209102a7,
+    "SFRBX_UART1":       0x20910232, "SFRBX_USB":     0x20910234,
+    "NAVCOR_ENABLE_HOST":    0x100d0001,   # HPG >= 2.10
+    "NAVCOR_ENABLE_GAL_HAS": 0x100d0002,   # HPG >= 2.10
+    "DYNMODEL":          0x20110021,       # CFG-NAVSPG-DYNMODEL
+}
+
+# CFG-NAVSPG-DYNMODEL values; default "sea" for the buoy deployment
+DYNMODEL = {"portable": 0, "stationary": 2, "pedestrian": 3, "automotive": 4,
+            "sea": 5, "airborne1g": 6, "airborne2g": 7, "airborne4g": 8, "wrist": 9}
+
+
+def log(msg):
+    print(f"[stream_gps] {msg}", flush=True)
+
+
+# ----------------------------------------------------------------------------
+# Receiver configuration
+# ----------------------------------------------------------------------------
+
+def valset(ser, ubr, cfg, label):
+    """Send CFG-VALSET (RAM) and wait briefly for ACK; returns True on ACK-ACK."""
+    msg = UBXMessage.config_set(layers=SET_LAYER_RAM, transaction=0, cfgData=cfg)
+    ser.write(msg.serialize())
+    t_end = time.time() + 2.0
+    while time.time() < t_end:
+        try:
+            _, parsed = ubr.read()
+        except Exception:
+            continue
+        if parsed is None:
+            continue
+        if parsed.identity == "ACK-ACK":
+            return True
+        if parsed.identity == "ACK-NAK":
+            log(f"  NAK for {label}")
+            return False
+    log(f"  no ACK for {label} (continuing)")
+    return False
+
+
+def poll_fw_version(ser, ubr):
+    """Poll MON-VER; returns firmware string like 'HPG 2.10' or None."""
+    ser.write(UBXMessage("MON", "MON-VER", 2).serialize())  # 2 = POLL msgmode
+    t_end = time.time() + 3.0
+    while time.time() < t_end:
+        try:
+            _, parsed = ubr.read()
+        except Exception:
+            continue
+        if parsed is not None and parsed.identity == "MON-VER":
+            sw = getattr(parsed, "swVersion", b"")
+            sw = sw.decode(errors="replace") if isinstance(sw, bytes) else str(sw)
+            return sw.strip("\x00 ")
+    return None
+
+
+def configure_receiver(ser, ubr, want_has, dynmodel):
+    """Enable RTCM input, HP output, raw logging; optionally HAS. RAM layer only."""
+    base_cfg = [
+        (K["UART1_IN_RTCM3"], 1), (K["USB_IN_RTCM3"], 1),
+        (K["PVT_UART1"], 1), (K["PVT_USB"], 1),
+        (K["HPPOSLLH_UART1"], 1), (K["HPPOSLLH_USB"], 1),
+        (K["RAWX_UART1"], 1), (K["RAWX_USB"], 1),
+        (K["SFRBX_UART1"], 1), (K["SFRBX_USB"], 1),
+        (K["DYNMODEL"], DYNMODEL[dynmodel]),
+    ]
+    ok = valset(ser, ubr, base_cfg, f"base output/input config (dynmodel={dynmodel})")
+    if want_has:
+        ok_has = valset(ser, ubr, [(K["NAVCOR_ENABLE_HOST"], 0),
+                                   (K["NAVCOR_ENABLE_GAL_HAS"], 1)], "GAL HAS enable")
+        return ok, ok_has
+    return ok, False
+
+
+# ----------------------------------------------------------------------------
+# NTRIP
+# ----------------------------------------------------------------------------
+
+def caster_reachable(server, port, timeout=8):
+    try:
+        with socket.create_connection((server, port), timeout=timeout):
+            return True
+    except OSError:
+        return False
+
+
+class SerialRTCMSink:
+    """Write-only wrapper handed to GNSSNTRIPClient as its output stream."""
+    def __init__(self, ser):
+        self._ser = ser
+        self.bytes_out = 0
+
+    def write(self, data):
+        # GNSSNTRIPClient may hand us (raw, parsed) tuples or raw bytes
+        if isinstance(data, tuple):
+            data = data[0]
+        if isinstance(data, (bytes, bytearray)):
+            self._ser.write(bytes(data))
+            self.bytes_out += len(data)
+
+    # queue-style interface (pygnssutils uses .put on Queue outputs)
+    def put(self, item):
+        self.write(item)
+
+
+def start_ntrip(args, sink, reflat, reflon, app=None):
+    """app (with .get_coordinates()) => live GGA that follows a drifting platform;
+    without app, falls back to a fixed reference position."""
+    from pygnssutils import GNSSNTRIPClient
+    client = GNSSNTRIPClient(app=app)
+    kwargs = dict(
+        server=args.ntrip_caster, port=args.ntrip_port, https=0,
+        mountpoint=args.ntrip_mount, datatype="RTCM",
+        ntripuser=args.ntrip_user, ntrippassword=args.ntrip_pass,
+        ggainterval=args.gga_interval,
+        ggamode=0 if app is not None else 1,    # 0 = live from receiver, 1 = fixed
+        reflat=reflat, reflon=reflon, refalt=0.0, refsep=0.0,
+        output=sink,
+    )
+    ok = client.run(**kwargs)
+    return client if ok else None
+
+
+# ----------------------------------------------------------------------------
+# Main streaming loop
+# ----------------------------------------------------------------------------
+
+class Session:
+    def __init__(self, args):
+        self.args = args
+        self.q = queue.Queue()
+        self.stop = threading.Event()
+        self.fix_hist = Counter()
+        self.hacc = []
+        self.tier = "none"
+        # latest fix, served to GNSSNTRIPClient.get_coordinates() for live GGA
+        self.live = {"lat": 0.0, "lon": 0.0, "alt": 0.0, "sep": 0.0}
+        ts = datetime.now(timezone.utc).strftime("%Y%m%d_%H%M%S")
+        self.csv_path = Path(f"stream_{ts}.csv")
+        self.ubx_path = Path(f"stream_{ts}.ubx")
+
+    def get_coordinates(self):
+        """pygnssutils app interface: live position for the NTRIP GGA uplink,
+        so the VRS reference follows a drifting platform (buoy)."""
+        return dict(self.live)
+
+    def reader_thread(self, stream, raw_log):
+        """Serial/file -> raw .ubx tee -> parsed NAV messages -> queue."""
+        ubr = UBXReader(stream, protfilter=UBX_PROTOCOL, quitonerror=0)
+        while not self.stop.is_set():
+            try:
+                raw, parsed = ubr.read()
+            except Exception:
+                if self.args.replay:
+                    break
+                continue
+            if raw is None:
+                if self.args.replay:
+                    break
+                continue
+            if raw_log:
+                raw_log.write(raw)
+            if parsed is None:
+                continue
+            if parsed.identity in ("NAV-PVT", "NAV-HPPOSLLH"):
+                self.q.put(parsed)
+        self.q.put(None)
+
+    def run(self):
+        args = self.args
+        # ---------------- input: serial or replay file ----------------------
+        if args.replay:
+            stream = open(args.replay, "rb")
+            ser = None
+            log(f"REPLAY mode from {args.replay} (no config writes, no NTRIP)")
+        else:
+            import serial
+            from serial.tools import list_ports
+            if not args.port:
+                ports = [p.device for p in list_ports.comports()]
+                sys.exit(f"--port required. Detected ports: {ports or 'none'}")
+            ser = serial.Serial(args.port, args.baud, timeout=1)
+            stream = ser
+            log(f"connected {args.port} @ {args.baud}")
+
+            cfg_ubr = UBXReader(ser, protfilter=UBX_PROTOCOL, quitonerror=0)
+            fw = poll_fw_version(ser, cfg_ubr)
+            log(f"firmware: {fw or 'unknown'}")
+            fw_has = bool(fw and any(
+                fw.split("HPG")[-1].strip().startswith(v)
+                for v in ("2.1", "2.2", "3.")))
+
+            # ------------- correction tier selection -------------------------
+            ntrip_ok = False
+            if not args.log_only and args.ntrip_user and args.ntrip_pass:
+                ntrip_ok = caster_reachable(args.ntrip_caster, args.ntrip_port)
+                if not ntrip_ok:
+                    log(f"NTRIP caster {args.ntrip_caster}:{args.ntrip_port} unreachable")
+            elif not args.log_only:
+                log("no NTRIP credentials (set GCGC_USER/GCGC_PASS or --ntrip-user/-pass)")
+
+            want_has = not ntrip_ok and not args.log_only and fw_has
+            if ntrip_ok:
+                self.tier = f"NTRIP RTK ({args.ntrip_caster}/{args.ntrip_mount})"
+            elif want_has:
+                self.tier = "Galileo HAS (E6-B, ~20 cm)"
+            elif args.log_only:
+                self.tier = "log-only (onboard solution)"
+            else:
+                log("no correction source available that improves on the onboard "
+                    "solution (NTRIP unreachable/no credentials; HAS needs FW >= "
+                    "HPG 2.10). Skipping enhanced streaming per project policy.")
+                log("use --log-only to record raw data for post-processing anyway.")
+                ser.close()
+                return 2
+
+            ok_base, ok_has = configure_receiver(ser, cfg_ubr, want_has, args.dynmodel)
+            log(f"receiver configured (base={'ok' if ok_base else 'partial'}"
+                + (f", HAS={'ok' if ok_has else 'failed'}" if want_has else "") + ")")
+
+        log(f"correction tier: {self.tier}")
+
+        # ---------------- outputs -------------------------------------------
+        raw_log = None if args.replay else open(self.ubx_path, "wb")
+        csv_f = open(self.csv_path, "w", newline="")
+        csv_f.write(f"# stream_gps.py session - correction tier: {self.tier}\n")
+        csv_f.write("# datum: NAD83(2011) if GCGC NTRIP tier, else WGS84/ITRF\n")
+        writer = csv.writer(csv_f)
+        writer.writerow(["utc", "lat_deg", "lon_deg", "ellip_height_m",
+                         "hAcc_m", "vAcc_m", "fixType", "carrSoln", "numSV"])
+
+        # ---------------- threads -------------------------------------------
+        t_read = threading.Thread(target=self.reader_thread, args=(stream, raw_log),
+                                  daemon=True)
+        t_read.start()
+
+        ntrip_client, sink = None, None
+        if ser is not None and self.tier.startswith("NTRIP"):
+            # need a first fix for the GGA reference position
+            reflat, reflon = args.lat, args.lon
+            if reflat is None:
+                log("waiting for first fix to seed NTRIP GGA position ...")
+                first = self._wait_first_fix(timeout=60)
+                if first is None:
+                    log("no fix within 60 s; using 0/0 GGA (VRS may reject)")
+                    reflat = reflon = 0.0
+                else:
+                    reflat, reflon = first
+            sink = SerialRTCMSink(ser)
+            app = None if args.gga_fixed else self
+            log("GGA uplink: " + ("fixed reference" if app is None
+                                  else "live (follows platform drift)"))
+            ntrip_client = start_ntrip(args, sink, reflat, reflon, app=app)
+            if ntrip_client is None:
+                log("NTRIP connection failed (check credentials/mountpoint) - "
+                    "continuing with onboard solution only")
+                self.tier += " [FAILED - onboard only]"
+
+        signal.signal(signal.SIGINT, lambda *_: self.stop.set())
+
+        # ---------------- consume + display ---------------------------------
+        last_pvt = {}
+        try:
+            while not self.stop.is_set():
+                try:
+                    m = self.q.get(timeout=1.0)
+                except queue.Empty:
+                    continue
+                if m is None:
+                    break
+                if m.identity == "NAV-PVT":
+                    last_pvt = dict(
+                        utc=datetime(m.year, m.month, m.day, m.hour, m.min, m.second,
+                                     tzinfo=timezone.utc),
+                        fixType=m.fixType, carrSoln=m.carrSoln, numSV=m.numSV,
+                        lat=m.lat, lon=m.lon, height=m.height / 1000.0,
+                        hAcc=m.hAcc / 1000.0, vAcc=m.vAcc / 1000.0)
+                elif m.identity == "NAV-HPPOSLLH" and last_pvt:
+                    # high-precision refinement: lat/lon + Hp components (deg)
+                    last_pvt["lat"] = m.lat + m.latHp
+                    last_pvt["lon"] = m.lon + m.lonHp
+                    last_pvt["height"] = (m.height + m.heightHp) / 1000.0
+                    last_pvt["hAcc"] = m.hAcc / 1000.0
+                    last_pvt["vAcc"] = m.vAcc / 1000.0
+                    continue
+                if not last_pvt:
+                    continue
+                p = last_pvt
+                self.live.update(lat=p["lat"], lon=p["lon"], alt=p["height"])
+                self.fix_hist[(p["fixType"], p["carrSoln"])] += 1
+                self.hacc.append(p["hAcc"])
+                writer.writerow([p["utc"].isoformat(), f"{p['lat']:.9f}",
+                                 f"{p['lon']:.9f}", f"{p['height']:.3f}",
+                                 f"{p['hAcc']:.3f}", f"{p['vAcc']:.3f}",
+                                 p["fixType"], p["carrSoln"], p["numSV"]])
+                n_epochs = sum(self.fix_hist.values())
+                if not sys.stdout.isatty() and n_epochs % 60 != 1:
+                    continue                      # avoid log spam when piped
+                fixlab = FIX_NAME.get(p["fixType"], "?")
+                if p["carrSoln"]:
+                    fixlab += "/" + CARR_NAME[p["carrSoln"]]
+                rtcm = f"  rtcm {sink.bytes_out//1024} kB" if sink else ""
+                end = "   " if sys.stdout.isatty() else "\n"
+                status = (f"\r{p['utc']:%H:%M:%S}  {fixlab:<12s} "
+                          f"{p['lat']:+.9f} {p['lon']:+.9f}  h={p['height']:8.3f} m  "
+                          f"hAcc={p['hAcc']:6.3f} m  sv={p['numSV']:2d}{rtcm}{end}")
+                sys.stdout.write(status)
+                sys.stdout.flush()
+        finally:
+            self.stop.set()
+            print()
+            if ntrip_client is not None:
+                try:
+                    ntrip_client.stop()
+                except Exception:
+                    pass
+            csv_f.close()
+            if raw_log:
+                raw_log.close()
+            if ser:
+                ser.close()
+            self._summary()
+        return 0
+
+    def _wait_first_fix(self, timeout):
+        t_end = time.time() + timeout
+        while time.time() < t_end:
+            try:
+                m = self.q.get(timeout=1.0)
+            except queue.Empty:
+                continue
+            if m is not None and m.identity == "NAV-PVT" and m.fixType >= 2:
+                self.live.update(lat=m.lat, lon=m.lon, alt=m.height / 1000.0)
+                return (m.lat, m.lon)
+        return None
+
+    def _summary(self):
+        n = sum(self.fix_hist.values())
+        if not n:
+            log("no position epochs received")
+            return
+        log(f"session summary: {n} epochs, tier = {self.tier}")
+        for (fix, carr), c in sorted(self.fix_hist.items(), key=lambda kv: -kv[1]):
+            lab = FIX_NAME.get(fix, "?") + ("/" + CARR_NAME[carr] if carr else "")
+            log(f"  {lab:<14s} {c:6d}  ({100*c/n:.1f}%)")
+        h = sorted(self.hacc)
+        log(f"  hAcc: mean {sum(h)/len(h):.3f} m, median {h[len(h)//2]:.3f} m, "
+            f"95% {h[int(len(h)*0.95)]:.3f} m")
+        if not self.args.replay:
+            log(f"  outputs: {self.csv_path}, {self.ubx_path}")
+            log(f"  post-process the raw log: venv/bin/python process_gps.py {self.ubx_path}")
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__,
+                                 formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--port", help="serial port (COM7, /dev/ttyUSB0, ...)")
+    ap.add_argument("--baud", type=int, default=38400)
+    ap.add_argument("--replay", type=Path, help="read a .ubx file instead of serial "
+                                                "(offline parser/CSV test mode)")
+    ap.add_argument("--log-only", action="store_true",
+                    help="skip corrections; just log raw data + onboard positions")
+    ap.add_argument("--lat", type=float, help="initial latitude for NTRIP GGA (else first fix)")
+    ap.add_argument("--lon", type=float, help="initial longitude for NTRIP GGA")
+    ap.add_argument("--gga-fixed", action="store_true",
+                    help="send a FIXED GGA reference instead of live position "
+                         "(only for a truly stationary antenna; default is live, "
+                         "so the VRS reference follows a drifting buoy)")
+    ap.add_argument("--dynmodel", default="sea", choices=sorted(DYNMODEL),
+                    help="receiver dynamic platform model, set in RAM at startup "
+                         "(default: sea, for buoy deployment)")
+    ap.add_argument("--ntrip-caster", default="rtn.usm.edu")
+    ap.add_argument("--ntrip-port", type=int, default=2101)
+    ap.add_argument("--ntrip-mount", default="VRS_GNSS_RTCM34_NAD83")
+    ap.add_argument("--ntrip-user", default=os.environ.get("GCGC_USER"))
+    ap.add_argument("--ntrip-pass", default=os.environ.get("GCGC_PASS"))
+    ap.add_argument("--gga-interval", type=int, default=15)
+    args = ap.parse_args()
+
+    if not args.replay and not args.port:
+        ap.error("one of --port or --replay is required")
+    sys.exit(Session(args).run())
+
+
+if __name__ == "__main__":
+    main()