# -*- coding: utf-8 -*-
"""
Qlog Converter V5 — QuSpin Magpie log (.txt) to CSV converter.

Drag one or more Magpie .txt logs onto the exe (or pass them as command-line
arguments). For each input it writes:
    <name>_flat.csv          every parsed record, no processing
    <name>_interpolated.csv  sensor rows with GPS/UTC interpolated onto them
    <name>_dedup.csv         interpolated rows de-duplicated on Timestamp(uS)
    <name>_summary.txt       stats + all '#' status lines from the raw log

V5 changes vs V4 (2026-07-09, Jeff Orton + Claude):
  - Pure standard library (no pandas/numpy) so the packaged exe is ~10 MB
    instead of ~250 MB.
  - FIX: dedup no longer collapses the file to one row when a key column is
    entirely empty (V4 dedup'd on SensorTimestamp/UTCTime; an all-empty column
    made every row a "duplicate"). V5 dedups on Timestamp(uS) only.
  - FIX: lat/lon parsing is now int/1e7 — V4's string slicing returned None
    for any coordinate below 1 degree (e.g. longitudes near Greenwich).
  - FIX: jitter correction is now SEGMENTED — timestamps are uniformly
    re-spaced only within runs of steady cadence; a gap larger than 3x the
    median delta (GPS dropout, sensor restart) starts a new segment instead
    of being smeared across the whole file.
  - FIX: magnetometer jitter correction now also runs on logs with no GPS
    (V4 skipped it unless GPS rows existed).
  - HAE sanity check: heights outside +/-10000 m are treated as not present,
    so old-format logs (fields were ms,sats) degrade to empty Height instead
    of garbage.
  - Headers are written once with the final names (no rewrite pass).
  - Errors and completion pause the console when run as an exe, so
    drag-and-drop users see what happened.
"""

# ------------------------------------------------------------------
# DEV OVERRIDE: set this to test a single .txt file directly.
# Leave as "" to use drag-and-drop / command-line args.
TEST_FILE = r""

# Vector/scalar timing correction. On Survey output the vector component on a
# given row was measured BEFORE the scalar on that same row - see
# VECTOR_SCALAR_TIMING.md. Set to False to write CSVs the old (uncorrected) way.
CORRECT_VECTOR_LAG = True

# Claude 2026-09-02 (Jeff): optional post-processing calibration read from a
# plain text file, normally kept on the Magpie SD card beside the logs.
#
# Why a file and not sensor EEPROM: EEPROM is 512 shared words and every new
# constant costs a firmware flash. A text file holds arbitrary richness, travels
# with the unit, and the converter can gain stages without touching firmware.
#
# Currently used for the SplitVec polarity correction. In split mode the two
# polarity legs of an axis are emitted separately, and a coil-pair DIRECTION
# asymmetry offsets them from each other - on a bench unit, X+ and X- sat 196 nT
# apart. It depends on field direction, so no per-coil scalar (delta or scale)
# can express it, which is why it cannot live in the sensor at all today.
CAL_FILE = "magpie_cal.txt"      # searched next to the log, then next to this script
APPLY_POLARITY_CORRECTION = True


def load_cal(log_path):
    """Read the calibration text file if one is present. Missing file, missing
    keys and unknown keys are all fine - every stage is optional."""
    import os
    cands = [os.path.join(os.path.dirname(os.path.abspath(log_path)), CAL_FILE),
             os.path.join(os.path.dirname(os.path.abspath(__file__)), CAL_FILE)]
    for p in cands:
        if not os.path.isfile(p):
            continue
        cal = {}
        with open(p, encoding="utf-8", errors="replace") as fh:
            for line in fh:
                line = line.split("#", 1)[0].strip()
                if not line or "=" not in line:
                    continue
                k, v = (x.strip() for x in line.split("=", 1))
                cal[k] = v
        cal["_path"] = p
        return cal
    return None


def polarity_coeffs(cal):
    """-> {axis: (a0, ax, ay, az)} or None if the stage is off/absent."""
    if not cal or cal.get("splitvec_polarity", "off").lower() not in ("on", "true", "1"):
        return None
    out = {}
    for ax in "XYZ":
        v = cal.get("polarity_" + ax)
        if not v:
            return None
        parts = v.split()
        if len(parts) != 4:
            return None
        out[ax] = tuple(float(x) for x in parts)
    return out


def apply_polarity_correction(legs, B_avg, coeffs):
    """legs: (plus, minus) arrays for one axis. B_avg: N x 3 field from the
    AVERAGED legs - that combination is immune to the error being corrected,
    which is what keeps this from feeding itself. Returns corrected (plus, minus).

    The alternation is a DC offset between the legs that varies with field
    direction:  alt = a0 + a . B ,  leg+ -= alt , leg- += alt
    """
    a0, ax, ay, az = coeffs
    alt = a0 + B_avg[:, 0] * ax + B_avg[:, 1] * ay + B_avg[:, 2] * az
    plus, minus = legs
    return plus - alt, minus + alt

# Force a mode when auto-detection cannot tell: "survey", "splitvec" or "" (auto)
FORCE_VECTOR_MODE = ""
# ------------------------------------------------------------------

import sys
import os
import re
import math
import textwrap

# ---------------------------------------------
# Column layout (single source of truth)
# ---------------------------------------------
COLUMNS = [
    "Timestamp(uS)", "SensorTimestamp(ms)", "MagData(nT)", "Good S - Data",
    "Data Counter", "ScalarSensitivity", "bX(nT)", "bY(nT)", "bZ(nT)",
    "Good V - Data", "VectorSensitivity", "Latitude", "Longitude", "HAE(m)",
    "UTCDate(YYYY-MM-DD)", "TimeUTC(HH:MM:SS.sss)",
    "AccelerometerX(mg)", "AccelerometerY(mg)", "AccelerometerZ(mg)",
    "GyroX", "GyroY", "GyroZ",
    "IMU_VMagX(uT)", "IMU_VMagY(uT)", "IMU_VMagZ(uT)", "IMU_Temp(C)",
]
IDX = {name: i for i, name in enumerate(COLUMNS)}

# Numeric columns eligible for interpolation onto sensor rows. Only quantities
# that are genuinely continuous belong here.
INTERP_COLS = [
    "SensorTimestamp(ms)", "MagData(nT)",
    "bX(nT)", "bY(nT)", "bZ(nT)",
    "Latitude", "Longitude", "HAE(m)",
    "AccelerometerX(mg)", "AccelerometerY(mg)", "AccelerometerZ(mg)",
    "GyroX", "GyroY", "GyroZ",
    "IMU_VMagX(uT)", "IMU_VMagY(uT)", "IMU_VMagZ(uT)", "IMU_Temp(C)",
]
# Status, string and DISCRETE numeric columns: fill gaps from a neighbour
# rather than interpolating.
#
# Claude 2026-09-21 (Jeff, reported by a customer): Data Counter and the two
# sensitivity indicators used to sit in INTERP_COLS, and linearly interpolating
# them produced values like "Data Counter 386.9881498918329" and
# "ScalarSensitivity 82.00007900072112" on filled rows. Both are discrete: the
# counter ticks once per measurement cycle, and the sensitivity indicators are
# integer quality figures. A value between two of them does not mean anything,
# and on a filled row there was no measurement to have a counter or a quality
# figure for in the first place, so carrying the neighbour's value is the
# honest answer. Data Counter is read by reconstruct_timebase(), which runs
# BEFORE this stage, so the timebase rebuild is unaffected.
FILL_COLS = [
    "Good S - Data", "Good V - Data", "UTCDate(YYYY-MM-DD)",
    "Data Counter", "ScalarSensitivity", "VectorSensitivity",
]

# ---------------------------------------------
# 1) Parsing
# ---------------------------------------------

TOKEN_PATTERN = re.compile(r"""
      (?: (?P<vector>[XYZ]) (?P<vec_val>[+-]?\d+(?:\.\d+)?) (?P<vec_ok>[=?])? )
    | (?:@(?P<counter>\d{1,3}))
    | (?:>(?P<ms_time>\d{1,10}))
    | (?:s(?P<scalar_sens>\d{1,3}))
    | (?:v(?P<vector_sens>\d{1,3}))
""", re.VERBOSE)

IMU_PATTERN = re.compile(r"([abcijkxyzt])([+-]?\d+(?:\.\d+)?)")
IMU_MAP = {
    'a': "AccelerometerX(mg)", 'b': "AccelerometerY(mg)", 'c': "AccelerometerZ(mg)",
    'i': "GyroX", 'j': "GyroY", 'k': "GyroZ",
    'x': "IMU_VMagX(uT)", 'y': "IMU_VMagY(uT)", 'z': "IMU_VMagZ(uT)",
    't': "IMU_Temp(C)",
}
SCALAR_PATTERN = re.compile(r"^([+-]?\d+(?:\.\d+)?)([_*])?")
MAG_LINE_PATTERN = re.compile(r"^(\d+)!(.+)$", re.DOTALL)


def parse_latlon(value_str):
    """Raw lat/lon is an integer of degrees * 1e7. Works for any magnitude."""
    try:
        return int(value_str.strip()) / 1e7
    except (ValueError, AttributeError):
        return None


def parse_mag_line(line):
    m = MAG_LINE_PATTERN.match(line)
    if not m:
        return None
    row = [None] * len(COLUMNS)
    try:
        row[IDX["Timestamp(uS)"]] = int(m.group(1))
    except ValueError:
        return None
    payload = m.group(2).strip()

    ms = SCALAR_PATTERN.match(payload)
    if not ms:
        return None
    try:
        row[IDX["MagData(nT)"]] = float(ms.group(1))
    except ValueError:
        return None
    valid_char = ms.group(2)
    if valid_char == '_':
        row[IDX["Good S - Data"]] = True
    elif valid_char == '*':
        row[IDX["Good S - Data"]] = False

    remainder = payload[ms.end():]

    for t in TOKEN_PATTERN.finditer(remainder):
        if t.group('vector'):
            axis = t.group('vector')
            col = {"X": "bX(nT)", "Y": "bY(nT)", "Z": "bZ(nT)"}[axis]
            try:
                row[IDX[col]] = float(t.group('vec_val'))
            except ValueError:
                pass
            ok = t.group('vec_ok')
            if ok == '=':
                row[IDX["Good V - Data"]] = True
            elif ok == '?':
                row[IDX["Good V - Data"]] = False
        elif t.group('counter'):
            row[IDX["Data Counter"]] = int(t.group('counter'))
        elif t.group('ms_time'):
            row[IDX["SensorTimestamp(ms)"]] = int(t.group('ms_time'))
        elif t.group('scalar_sens'):
            row[IDX["ScalarSensitivity"]] = int(t.group('scalar_sens'))
        elif t.group('vector_sens'):
            row[IDX["VectorSensitivity"]] = int(t.group('vector_sens'))

    for t in IMU_PATTERN.finditer(remainder):
        col = IMU_MAP.get(t.group(1))
        if col:
            try:
                row[IDX[col]] = float(t.group(2))
            except ValueError:
                pass
    return row


def parse_gnss_line(line):
    parts = line.strip().split(',')
    if len(parts) < 12 or parts[1] != "GNSSFIX":
        return None
    row = [None] * len(COLUMNS)
    try:
        row[IDX["Timestamp(uS)"]] = int(parts[0])
    except ValueError:
        return None
    row[IDX["Latitude"]] = parse_latlon(parts[2])
    row[IDX["Longitude"]] = parse_latlon(parts[3])
    row[IDX["UTCDate(YYYY-MM-DD)"]] = f"{parts[4]}-{parts[5]}-{parts[6]}"
    row[IDX["TimeUTC(HH:MM:SS.sss)"]] = f"{parts[7]}:{parts[8]}:{parts[9]}"

    # Height above ellipsoid, split into integer.decimal fields. Old-format
    # logs carry ms,sats in these positions — the sanity bound rejects those
    # so Height degrades to empty instead of garbage.
    hae_int, hae_dec = parts[10].strip(), parts[11].strip()
    hae = None
    if hae_int:
        try:
            hae = float(f"{hae_int}.{hae_dec}") if hae_dec else float(hae_int)
        except ValueError:
            hae = None
    if hae is not None and not (-10000.0 <= hae <= 10000.0):
        hae = None
    row[IDX["HAE(m)"]] = hae
    return row


# Claude 2026-09-15 (Jeff): sensor status lines do NOT reach the SD card bare.
# read_MAG_data() synthesizes a "now" value and prepends it to every line it
# forwards, status lines included, so what is actually on the card looks like
#   78953550#Auto Start Active
# An unanchored startswith('#') therefore never matched a single one of them,
# and because such a line has no '!' and no ',GNSSFIX,' it fell through to
# rec = None and was dropped without trace. That is why "Status Updates" read
# "No status updates found" even on logs that really did contain sensor errors
# (verified on a customer log: 0 lines start with '#', 2 contain one).
# The Magpie's own session header is written straight to the file handle and so
# has no prefix; \d* covers both spellings.
_STATUS_RE = re.compile(r'^\d*#')


# Claude 2026-09-22 (Jeff, from a customer log): a logger timestamp
# that is wildly out of step with the lines either side of it is corruption, not
# data. On that log the SD write moved three characters from the start of one
# line onto the end of the previous one:
#
#   5455461240545,GNSSFIX,184556336,-670658560,2026,08,27,18,34,14,144,8
#   6812336!36379.210_X-23727.138=@167>5419287s078v003
#
# '545' + '6812336' = 5456812336, an ordinary stamp. As written, the GNSS line
# sorts 63 days into the future and the mag line sorts to the very front. One
# such line dragged the timebase fit from a 34 ms residual to 4.8 s, and the
# fitted rate from 999.94 to 999.85 ms/s - about 7 seconds of UTC error
# accumulated across the record.
#
# The test is deliberately LOCAL: the log is written sequentially, so a stamp is
# judged against its neighbours in FILE order, not against a global range. A
# genuine gap (SD stall, sensor restart) is kept, because the lines that follow
# it agree with the new value; only a lone excursion that the file does not
# follow is dropped. Regular small steps are left alone on purpose - this log
# also has 38 backward steps of -5.596 s and -1.449 s which are NOT understood
# yet, and quietly swallowing them would hide a real problem.
TS_JUMP_US = 60_000_000      # 60 s; at 236 Hz nothing legitimate jumps this far
TS_LOOKAHEAD = 8

# Largest timebase fit residual we will build a UTC axis from. A clean log fits
# to tens of milliseconds (Data27 lands at 34 ms once its one corrupt line is
# removed); a second or more means the anchors and the sensor clock disagree.
TIMEBASE_RESID_LIMIT_S = 0.100


def _scan_stream(idxs, ts):
    """Indices (into ts) whose value is a lone excursion among their own kind."""
    if len(idxs) < TS_LOOKAHEAD + 2:
        return []
    seed = sorted(ts[i] for i in idxs[:2 * TS_LOOKAHEAD + 1])
    last = seed[len(seed) // 2]
    out = []
    for p, i in enumerate(idxs):
        t = ts[i]
        if abs(t - last) <= TS_JUMP_US:
            last = t
            continue
        nxt = [ts[j] for j in idxs[p + 1:p + 1 + TS_LOOKAHEAD]]
        if not nxt:
            last = t
            continue
        # Does the file carry on from `last`, or from `t`? Whichever the lines
        # after it are closer to is the real timeline.
        #
        # An earlier version counted how many of the next few landed within
        # TS_JUMP_US of each candidate. That fails when the stream is sparse
        # AND has a real gap: on Data27 the GPS dropped out for 335 s right
        # after the corrupt fix, so nothing landed near EITHER candidate, the
        # vote tied at nothing to nothing, and the corrupt value was kept.
        # Comparing distances has no such blind spot.
        med = sorted(nxt)[len(nxt) // 2]
        if abs(med - t) > abs(med - last):
            out.append((i, t, last))
        else:
            last = t                       # a real discontinuity - keep it
    return out


def drop_corrupt_timestamps(rows):
    """Remove rows whose Timestamp(uS) is a lone excursion. Returns (rows, dropped).

    Claude 2026-09-22: the GNSS and magnetometer lines are scanned as SEPARATE
    streams. They are interleaved in the file but their stamps are not always on
    the same footing - in the older log format (a 2025 log) every
    GNSSFIX line sits about 63 s behind the mag lines around it. Comparing the
    two streams against each other flagged 63 perfectly ordinary GPS fixes on
    that file. Within its own stream each value is consistent, and a genuinely
    corrupt one still stands out.
    """
    ts_col = IDX["Timestamp(uS)"]
    lat_col = IDX["Latitude"]
    mag_col = IDX["MagData(nT)"]
    ts = [r[ts_col] for r in rows]

    # Claude 2026-09-24 (Jeff, from a customer flight log): a timestamp of exactly 0 is
    # never a real reading - it is the Magpie writing a record whose MDP pulse
    # had not arrived yet, which older firmware emitted as "0!...". Left in, the
    # row is placed at the very start of the flight and grids as a spike far
    # from where it was measured.
    #
    # This needs its own rule. The lone-excursion scan below cannot catch them,
    # because they arrive in BURSTS - 519 consecutive records in the first 43 s
    # of that file - and a burst looks to the lookahead like a real
    # discontinuity, so it keeps them. The scan found 3 of 519. Zero is not a
    # judgement call, so test it directly.
    zeros = [(i, 0, None) for i, v in enumerate(ts)
             if v == 0 and (rows[i][mag_col] is not None or rows[i][lat_col] is not None)]

    gnss_idx, mag_idx = [], []
    for i, r in enumerate(rows):
        if ts[i] is None or ts[i] == 0:
            continue
        if r[lat_col] is not None:
            gnss_idx.append(i)
        elif r[mag_col] is not None:
            mag_idx.append(i)

    dropped = zeros + _scan_stream(gnss_idx, ts) + _scan_stream(mag_idx, ts)
    if not dropped:
        return rows, []
    dropped.sort()
    bad = {i for i, _, _ in dropped}
    return [r for i, r in enumerate(rows) if i not in bad], dropped


def parse_file(path):
    """Parse the raw log. Returns (rows, status_lines)."""
    rows, status = [], []
    with open(path, 'r', encoding='utf-8', errors='ignore') as f:
        for line in f:
            line = line.strip()
            if not line:
                continue
            if _STATUS_RE.match(line):
                status.append(line)
                continue
            if ",GNSSFIX," in line:
                rec = parse_gnss_line(line)
            elif "!" in line:
                rec = parse_mag_line(line)
            else:
                rec = None
            if rec is not None:
                rows.append(rec)
    return rows, status


# ---------------------------------------------
# 2) CSV output helpers
# ---------------------------------------------

def fmt_cell(v):
    if v is None:
        return ""
    if isinstance(v, bool):
        return "True" if v else "False"
    if isinstance(v, float):
        if math.isnan(v):
            return ""
        return f"{v:.10g}"
    return str(v)


def write_csv(path, rows):
    with open(path, 'w', encoding='utf-8', newline='') as f:
        f.write(','.join(COLUMNS) + '\n')
        for row in rows:
            f.write(','.join(fmt_cell(v) for v in row) + '\n')


# ---------------------------------------------
# 3) Jitter correction (segmented)
# ---------------------------------------------

def correct_jitter_segmented(rows, indices, label):
    """
    Uniformly re-space Timestamp(uS) within runs of steady cadence.
    `indices` selects rows of one class (mag or GPS) in time order.
    A delta > 3x the median delta ends the current segment, so gaps and
    rate changes are never smeared across the file.
    """
    n = len(indices)
    if n < 3:
        return
    ts_col = IDX["Timestamp(uS)"]
    ts = [rows[i][ts_col] for i in indices]

    deltas = sorted(ts[k + 1] - ts[k] for k in range(n - 1))
    median = deltas[len(deltas) // 2]
    if median <= 0:
        return
    gap_limit = 3 * median

    seg_start = 0
    segments = []
    for k in range(n - 1):
        if ts[k + 1] - ts[k] > gap_limit:
            segments.append((seg_start, k))
            seg_start = k + 1
    segments.append((seg_start, n - 1))

    corrected = 0
    for a, b in segments:
        count = b - a + 1
        if count < 3:
            continue
        t0, t1 = ts[a], ts[b]
        step = (t1 - t0) / (count - 1)
        for j in range(count):
            rows[indices[a + j]][ts_col] = int(round(t0 + j * step))
        corrected += count
    print(f"Jitter correction ({label}): {corrected} rows in "
          f"{len(segments)} segment(s), median dt {median} uS")


# ---------------------------------------------
# 4) Interpolation (linear on Timestamp, ends clamped)
# ---------------------------------------------

def interpolate_column(rows, col_idx, ts_list):
    """Fill None cells of one numeric column by linear interpolation on
    timestamp; values beyond the first/last known point are clamped."""
    known = [(ts_list[i], rows[i][col_idx])
             for i in range(len(rows)) if rows[i][col_idx] is not None]
    if not known:
        return
    if len(known) == 1:
        v = known[0][1]
        for r in rows:
            if r[col_idx] is None:
                r[col_idx] = v
        return
    k = 0
    last = len(known) - 1
    for i, r in enumerate(rows):
        if r[col_idx] is not None:
            continue
        t = ts_list[i]
        while k < last - 1 and known[k + 1][0] <= t:
            k += 1
        t0, v0 = known[k]
        t1, v1 = known[k + 1]
        if t <= known[0][0]:
            r[col_idx] = known[0][1]
        elif t >= known[last][0]:
            r[col_idx] = known[last][1]
        elif t1 == t0:
            r[col_idx] = v0
        else:
            r[col_idx] = v0 + (v1 - v0) * (t - t0) / (t1 - t0)


def fill_column(rows, col_idx):
    """Forward-fill then back-fill a status/string column."""
    last = None
    for r in rows:
        if r[col_idx] is not None:
            last = r[col_idx]
        elif last is not None:
            r[col_idx] = last
    nxt = None
    for r in reversed(rows):
        if r[col_idx] is not None:
            nxt = r[col_idx]
        elif nxt is not None:
            r[col_idx] = nxt


def linear_fit(xs, ys):
    """Least-squares y = a*x + b."""
    n = len(xs)
    mx = sum(xs) / n
    my = sum(ys) / n
    num = den = 0.0
    for x, y in zip(xs, ys):
        dx = x - mx
        num += dx * (y - my)
        den += dx * dx
    if den == 0:
        return None, None
    a = num / den
    return a, my - a * mx


def reconstruct_timebase(rows, mag_idx, gps_idx):
    """
    Build a MEASUREMENT-time axis (continuous seconds) for every row, locked to
    GPS UTC. Returns (axis, description) or (None, reason).

    Claude 2026-08-03: Timestamp(uS) is an ARRIVAL time, not a measurement time,
    and correcting its jitter in place cannot recover measurement time. Two
    separate effects corrupt it:

      1. Burst transmission. The sensor ships samples in groups (typically four),
         which arrive tens of microseconds apart followed by a multi-millisecond
         gap. Measurements are evenly spaced; arrivals are bunched.
      2. Clock drift. The logger's crystal and the sensor's measurement clock are
         different oscillators. Measured on a 500 Hz customer flight log: 13.7 ppm,
         which accumulated 92 ms of skew over 113 minutes. This is systematic,
         grows with recording length, and does not average out.

    correct_jitter_segmented() cannot fix either, because on burst data the
    MEDIAN dt is the within-burst spacing, so its 3x-median gap threshold makes
    every burst its own segment and it merely re-spaces four samples across their
    own arrival window. It is kept only as a last-resort fallback.

    This rebuilds the axis from quantities that are true regardless of what the
    serial link did: the sensor '>' millisecond stamp if present, otherwise the
    '@' counter (exactly one tick per measurement, wrapping once per second),
    fitted by least squares against GPS UTC.
    """
    ctr_col = IDX["Data Counter"]
    ms_col = IDX["SensorTimestamp(ms)"]
    utc_col = IDX["TimeUTC(HH:MM:SS.sss)"]

    if len(gps_idx) < 2:
        return None, "fewer than 2 GPS fixes to anchor absolute time"
    if not mag_idx:
        return None, "no magnetometer rows"

    # --- sensor-side axis: '>' ms stamp preferred, else the '@' counter -------
    n_ms = sum(1 for i in mag_idx if rows[i][ms_col] is not None)
    axis = {}
    if n_ms > 0.99 * len(mag_idx):
        # Claude 2026-09-22 (Jeff, from a customer log): the guard
        # deliberately tolerates up to 1% of mag rows having no '>' stamp, but
        # this loop then called float() on every one of them and died on the
        # first None. Three truncated lines out of 1,158,740 - the UART dropped
        # characters mid-line - killed the entire 63 MB conversion. Rows with
        # no stamp are simply left out of the axis; the neighbour fill at the
        # end of this function already assigns them a time.
        for i in mag_idx:
            v = rows[i][ms_col]
            if v is not None:
                axis[i] = float(v)
        label = "sensor millisecond timestamp ('>' field)"
    else:
        ctrs = [rows[i][ctr_col] for i in mag_idx]
        if any(c is None for c in ctrs):
            return None, "neither the '>' timestamp nor the '@' counter is present"
        period = max(ctrs) + 1
        wraps, prev = 0, None
        for i in mag_idx:
            c = rows[i][ctr_col]
            if prev is not None and c < prev:
                wraps += 1
            axis[i] = float(wraps * period + c)
            prev = c
        label = "'@' data counter (wrap %d, i.e. %d Hz) + GPS" % (period, period)

    # --- GPS UTC, unwrapped across midnight ---------------------------------
    gutc, prev = {}, None
    day = 0.0
    for i in gps_idx:
        s = utc_to_seconds(rows[i][utc_col])
        if s is None:
            continue
        if prev is not None and s < prev - 43200:
            day += 86400.0
        gutc[i] = s + day
        prev = s

    # Anchor each GPS fix to the axis value of the nearest preceding mag row.
    anchor_x, anchor_y, last_axis = [], [], None
    for i in range(len(rows)):
        if i in axis:
            last_axis = axis[i]
        elif i in gutc and last_axis is not None:
            anchor_x.append(last_axis)
            anchor_y.append(gutc[i])
    if len(anchor_x) < 2:
        return None, "could not pair GPS fixes with sensor samples"

    slope, intercept = linear_fit(anchor_x, anchor_y)
    if slope is None or slope <= 0:
        return None, "degenerate fit"

    resid = [abs(y - (slope * x + intercept)) for x, y in zip(anchor_x, anchor_y)]
    rms = (sum(r * r for r in resid) / len(resid)) ** 0.5
    rate = 1.0 / slope
    print("Timebase: %s" % label)
    print("  anchors %d, fitted rate %.6f, fit residual RMS %.1f ms"
          % (len(anchor_x), rate, rms * 1000.0))

    # Claude 2026-09-22 (Jeff): this residual used to be printed and then
    # ignored. On one customer log it read 71,085 ms and the fit was used
    # anyway, which is how a whole survey got UTC that was wrong by seconds and
    # looked fine. A clean log sits at tens of milliseconds; anything near a
    # second means the anchors disagree with the sensor clock and the line
    # through them is not worth trusting. Say so and fall back rather than emit
    # times we do not believe.
    if rms > TIMEBASE_RESID_LIMIT_S:
        worst = sorted(range(len(resid)), key=lambda i: resid[i], reverse=True)
        print("  REJECTED: residual %.0f ms exceeds the %.0f ms limit."
              % (rms * 1000.0, TIMEBASE_RESID_LIMIT_S * 1000.0))
        over = sum(1 for r in resid if r > TIMEBASE_RESID_LIMIT_S)
        print("    %d of %d anchors are off by more than the limit; worst:"
              % (over, len(resid)))
        for i in worst[:3]:
            print("      sensor clock %.0f -> GPS %.3f s, off by %+.1f s"
                  % (anchor_x[i], anchor_y[i],
                     anchor_y[i] - (slope * anchor_x[i] + intercept)))
        print("    A few bad anchors usually means corrupt lines in the log.")
        return None, ("timebase fit residual %.0f ms over the %.0f ms limit"
                      % (rms * 1000.0, TIMEBASE_RESID_LIMIT_S * 1000.0))

    # --- axis for every row --------------------------------------------------
    out = [None] * len(rows)
    for i, a in axis.items():
        out[i] = slope * a + intercept
    for i, y in gutc.items():
        out[i] = y
    # Any row with neither (shouldn't happen) gets filled by neighbours.
    last = None
    for i in range(len(rows)):
        if out[i] is None:
            out[i] = last
        else:
            last = out[i]
    nxt = None
    for i in range(len(rows) - 1, -1, -1):
        if out[i] is None:
            out[i] = nxt
        else:
            nxt = out[i]
    if any(v is None for v in out):
        return None, "could not assign a time to every row"
    return out, label


def utc_to_seconds(utc_str):
    try:
        h, m, s = utc_str.split(':')
        return int(h) * 3600 + int(m) * 60 + float(s)
    except (ValueError, AttributeError):
        return None


def seconds_to_utc(sec):
    # Claude 2026-09-21 (Jeff, reported by a customer): milliseconds, not
    # microseconds. This used to emit 6 decimal places and the extra three were
    # always zeros: when the timebase comes from the sensor's own millisecond
    # stamp or its cycle counter, which is the path we want people on, the
    # reconstructed axis lands on whole milliseconds. Checked over a 61,000 row
    # file, 0 rows had anything in digits 4 to 6. Sub-millisecond digits only
    # appear on the arrival-timestamp fallback, and that path already prints a
    # warning that its UTC carries UART bunching and logger drift, so those
    # digits were false precision exactly where they looked most convincing.
    sec = sec % 86400.0
    h = int(sec // 3600)
    m = int((sec % 3600) // 60)
    s = sec % 60
    return f"{h:02d}:{m:02d}:{s:06.3f}"


# ---------------------------------------------
# 5) Pipeline
# ---------------------------------------------

def process(rows):
    """Trim, jitter-correct, interpolate. Returns sensor-only rows."""
    ts_col = IDX["Timestamp(uS)"]
    mag_col = IDX["MagData(nT)"]
    lat_col = IDX["Latitude"]
    utc_col = IDX["TimeUTC(HH:MM:SS.sss)"]

    rows = sorted(rows, key=lambda r: r[ts_col])

    # Trim to the GPS-covered range so positions are never extrapolated.
    gps_idx = [i for i, r in enumerate(rows) if r[lat_col] is not None]
    if gps_idx:
        rows = rows[gps_idx[0]:gps_idx[-1] + 1]
        print(f"Trimmed to GPS range: {len(rows)} rows remaining")

    mag_idx = [i for i, r in enumerate(rows) if r[mag_col] is not None]
    gps_idx = [i for i, r in enumerate(rows) if r[lat_col] is not None]

    # Tag sensor rows NOW — interpolation fills MagData on every row, so the
    # sensor/GPS distinction must be captured before it runs.
    sensor_flags = [r[mag_col] is not None for r in rows]

    # Jitter correction runs on mag rows regardless of GPS presence (V4 bug).
    # --- Time base -----------------------------------------------------------
    # Preferred: rebuild MEASUREMENT time from the sensor's own clock/counter,
    # locked to GPS. Timestamp(uS) is deliberately left untouched so the raw
    # logger arrival stamp remains in the output for comparison.
    axis, info = reconstruct_timebase(rows, mag_idx, gps_idx)

    if axis is not None:
        ts_list = axis
        for i, r in enumerate(rows):
            r[utc_col] = axis[i]        # numeric seconds; formatted further down
    else:
        print("Timebase: FALLING BACK to arrival timestamps (%s)." % info)
        print("  WARNING: UTC will carry UART arrival bunching and logger clock")
        print("  drift. Enable the sensor timestamp (command 109) when recording.")
        correct_jitter_segmented(rows, mag_idx, "magnetometer")
        if gps_idx:
            correct_jitter_segmented(rows, gps_idx, "GPS")
        # Re-sort after timestamp rewrite, carrying the sensor flags along.
        order = sorted(range(len(rows)), key=lambda i: rows[i][ts_col])
        rows = [rows[i] for i in order]
        sensor_flags = [sensor_flags[i] for i in order]
        ts_list = [r[ts_col] for r in rows]
        utc_secs = [utc_to_seconds(r[utc_col]) if r[utc_col] is not None else None
                    for r in rows]
        for r, s in zip(rows, utc_secs):
            r[utc_col] = s

    for name in INTERP_COLS:
        interpolate_column(rows, IDX[name], ts_list)
    if axis is None:
        interpolate_column(rows, utc_col, ts_list)

    for name in FILL_COLS:
        fill_column(rows, IDX[name])

    for r in rows:
        if r[utc_col] is not None:
            r[utc_col] = seconds_to_utc(r[utc_col])

    # Keep only sensor rows: GPS-only rows served as interpolation anchors.
    sensor_rows = [r for r, is_sensor in zip(rows, sensor_flags) if is_sensor]
    return sensor_rows


def dedup_rows(rows):
    """Drop rows whose Timestamp(uS) repeats. Never dedups on a column that
    can be empty (the V4 bug that collapsed files to one row)."""
    ts_col = IDX["Timestamp(uS)"]
    seen = set()
    out = []
    for r in rows:
        t = r[ts_col]
        if t in seen:
            continue
        seen.add(t)
        out.append(r)
    print(f"De-dup on Timestamp(uS): kept {len(out)} of {len(rows)} rows")
    return out


# ---------------------------------------------
# Vector/scalar timing correction
# ---------------------------------------------
# The scalar and the vector on one output row were NOT measured at the same
# instant. Full derivation in VECTOR_SCALAR_TIMING.md; the short version:
#
# SURVEY (sequential coil order) - a component needs both its polarity slots
# plus their bracketing scalars before it can be computed, so it is printed
# 4 slots after its own measurement midpoint. The printed scalar is whatever
# was measured at print time. With N = data_rate/2 cycles averaged:
#
#     lag_ms  = 4 + 6*(N-1)        dt_ms = 4*N        =>   lag = 1.5*dt - 2
#
#     dt  4 ms (250 Hz)  -> lag  4 ms -> exactly 1 row
#     dt  8 ms (125 Hz)  -> lag 10 ms -> 1.25 rows
#     dt 16 ms (62.5 Hz) -> lag 22 ms -> 1.375 rows
#
# The vector is OLDER than its row's scalar, so the contemporaneous scalar sits
# on the PREVIOUS row: the vector columns move UP.
#
# SPLITVEC (500 Hz mode) - NO shift. split_vec_sample() builds its scalar as the
# mean of the unbiased samples immediately before and after the biased slot, and
# split_vec_emit() prints that stored scalar alongside the vector from the SAME
# FIFO entry. They are paired at source. Applying the Survey shift to a 500 Hz
# file would INTRODUCE a misalignment that was not there.
#
# 2026-09-17 - PAIR-ALIGNED FIRMWARE. Sensor firmware from 2026-09-17 pins the
# Survey scalar averaging to whole polarity pairs (release notes section 11): a
# 250 Hz output is now always the mean of the two unbiased slots that follow one
# axis's + and - pulses, never a straddled pair from two different axes. This only
# applies where the scalar is averaged at all, i.e. 250 Hz and slower; SplitVec
# does not average and is unaffected.
#
# What it means for this shift, at data_rate 2 (250 Hz), sequential order, slots
# numbered in ms within the 12 ms cycle (s x1 s x2 s y1 s y2 s z1 s z2):
#
#   scalar groups close at t = 4, 8, 12   -> (s@2,s@4) (s@6,s@8) (s@10,s@12)
#   ready flags are set at  t = 2, 6, 10  -> Z         X         Y
#     (Filter_functions.ino, Coil_Value == 2 / 6 / 10; one flag pending per row,
#      so the emitted order is the Z X Y rotation seen in real files)
#
#   row t=4 : scalar = X pair (centre 3), vector = Z (midpoint -2)  lag 5 ms
#   row t=8 : scalar = Y pair (centre 7), vector = X (midpoint  2)  lag 5 ms
#   row t=12: scalar = Z pair (centre 11), vector = Y (midpoint 6)  lag 5 ms
#
# So the true lag is 5 ms, not the 4 ms the pre-alignment derivation gave, and the
# integer shift is unchanged at exactly 1 row. The gain is in WHAT that row shift
# now buys: after it, each row's vector axis sits with the scalar average of that
# same axis's own two bracketing unbiased slots. That pairing used to depend on a
# coil-cycle phase that any FPGA parameter write could flip, so it was right only
# by luck. A +1 ms residual remains, the scalar being the newer of the two.
#
# Two consequences handled below:
#   * a pair realign drops one output row (the firmware counts them internally
#     '# pair realigns'), leaving a one-row gap in the timestamps. Shifting the
#     vector across such a gap pairs it with a scalar from the wrong cycle, so
#     rows whose source row is not exactly n rows away in TIME are left empty.
#   * below 250 Hz the vector needs data_rate/2 cycles to accumulate, so three
#     ready flags arrive in one 12 ms burst while rows come every 8 ms or slower.
#     The flags then queue and the X > Y > Z print priority reorders them. No row
#     shift can repair that; see the warning emitted in the summary.

VEC_COLS = ["bX(nT)", "bY(nT)", "bZ(nT)", "Good V - Data", "VectorSensitivity"]


def estimate_dt_ms(rows):
    """Median-ish row spacing in ms, same method write_summary uses."""
    ts_col = IDX["Timestamp(uS)"]
    n = len(rows)
    if n < 2:
        return None
    k = min(10, n - 1)
    try:
        return (rows[k][ts_col] - rows[0][ts_col]) / k / 1000.0
    except (TypeError, ZeroDivisionError):
        return None


def detect_vector_mode(rows, status_lines, dt_ms):
    """Return ("survey"|"splitvec"|"none", why)."""
    if FORCE_VECTOR_MODE:
        return FORCE_VECTOR_MODE, "forced by FORCE_VECTOR_MODE"

    # Does this log even carry vector data?
    bx = IDX["bX(nT)"]
    by = IDX["bY(nT)"]
    bz = IDX["bZ(nT)"]
    if not any(r[bx] is not None or r[by] is not None or r[bz] is not None
               for r in rows):
        return "none", "no vector columns populated"

    # Most reliable: the sensor announces SplitVec in a '#' status line.
    for line in status_lines:
        if "splitvec" in line.lower():
            return "splitvec", "'SplitVec' found in status lines"

    # Fallback on cadence. SplitVec forces data_rate=1 -> ~2 ms rows; Survey with
    # vector never runs there (data_rate 1 + vector is bumped to 2 at startup).
    if dt_ms is not None and dt_ms < 3.0:
        return "splitvec", "dt %.1f ms implies 500 Hz" % dt_ms
    return "survey", "dt %.1f ms" % (dt_ms if dt_ms is not None else -1)


def apply_polarity_stage(rows, mode, cal):
    """SplitVec only: close the DC gap between the two polarity legs.

    In split mode each row carries ONE axis, and consecutive rows of the same
    axis alternate polarity leg. A coil-pair direction asymmetry offsets the two
    legs from each other, so they cannot be used as samples of the same quantity
    - which defeats the point of the mode. Measured on a bench unit the X legs sat
    196 nT apart; corrected, 1.6 nT.

    B is built from the AVERAGE of each axis's legs. That combination is immune
    to the error being corrected, which is what stops this feeding itself.

    Rows are modified in place. Returns a report dict.
    """
    rep = {"applied": False, "note": "", "per_axis": {}}
    if not APPLY_POLARITY_CORRECTION:
        rep["note"] = "disabled (APPLY_POLARITY_CORRECTION = False)"
        return rep
    if mode != "splitvec":
        rep["note"] = "only applies to 500 Hz split vector output"
        return rep
    coeffs = polarity_coeffs(cal)
    if not coeffs:
        rep["note"] = ("no polarity coefficients in %s" % CAL_FILE) if cal                       else ("%s not found" % CAL_FILE)
        return rep

    ts_i = IDX["Timestamp(uS)"]
    col = {"X": IDX["bX(nT)"], "Y": IDX["bY(nT)"], "Z": IDX["bZ(nT)"]}

    # rows carrying each axis, in order
    idx = {a: [] for a in "XYZ"}
    for r, row in enumerate(rows):
        for a in "XYZ":
            v = row[col[a]]
            if v not in (None, ""):
                idx[a].append(r)
    if min(len(idx[a]) for a in "XYZ") < 8:
        rep["note"] = "not enough per-axis samples"
        return rep

    # per-axis averaged value (leg pair mean) on that axis's own timebase,
    # then interpolated to every row so each sample has a full B vector
    avg_t, avg_v = {}, {}
    for a in "XYZ":
        ii = idx[a]
        t = [float(rows[r][ts_i]) for r in ii]
        v = [float(rows[r][col[a]]) for r in ii]
        pt, pv = [], []
        for m in range(0, len(v) - 1, 2):
            pt.append(0.5 * (t[m] + t[m + 1]))
            pv.append(0.5 * (v[m] + v[m + 1]))
        if len(pt) < 4:
            rep["note"] = "not enough leg pairs on %s" % a
            return rep
        avg_t[a], avg_v[a] = pt, pv

    def interp(a, x):
        t, v = avg_t[a], avg_v[a]
        lo, hi = 0, len(t) - 1
        if x <= t[0]:
            return v[0]
        if x >= t[hi]:
            return v[hi]
        while hi - lo > 1:
            mid = (lo + hi) // 2
            if t[mid] <= x:
                lo = mid
            else:
                hi = mid
        f = (x - t[lo]) / (t[hi] - t[lo]) if t[hi] != t[lo] else 0.0
        return v[lo] + f * (v[hi] - v[lo])

    for a in "XYZ":
        a0, ax, ay, az = coeffs[a]
        ii = idx[a]
        alt = []
        for r in ii:
            x = float(rows[r][ts_i])
            alt.append(a0 + interp("X", x) * ax + interp("Y", x) * ay + interp("Z", x) * az)
        vals = [float(rows[r][col[a]]) for r in ii]
        # which subsequence is the "+" leg is not labelled in the output, so try
        # both assignments and keep the one that actually closes the gap
        def gap(sign):
            e = [vals[m] - sign * alt[m] for m in range(0, len(vals), 2)]
            o = [vals[m] + sign * alt[m] for m in range(1, len(vals), 2)]
            n = min(len(e), len(o))
            return abs(sum(e[:n]) / n - sum(o[:n]) / n)
        before = gap(0.0)
        sign = 1.0 if gap(1.0) <= gap(-1.0) else -1.0
        after = gap(sign)
        if after >= before:
            rep["per_axis"][a] = (before, before, False)
            continue
        for m, r in enumerate(ii):
            rows[r][col[a]] = vals[m] - (sign if m % 2 == 0 else -sign) * alt[m]
        rep["per_axis"][a] = (before, after, True)
        rep["applied"] = True
    if rep["applied"]:
        rep["note"] = "from %s" % cal.get("_path", CAL_FILE)
    return rep


def correct_vector_lag(rows, status_lines):
    """Shift vector columns up to align them with their own scalar.

    Returns a report dict; rows are modified in place.
    """
    dt_ms = estimate_dt_ms(rows)
    mode, why = detect_vector_mode(rows, status_lines, dt_ms)
    rep = {"applied": False, "mode": mode, "why": why, "dt_ms": dt_ms,
           "lag_ms": None, "rows_exact": None, "rows_shifted": 0,
           "residual_ms": 0.0, "note": "", "rows_gapped": 0, "pairing": ""}

    if not CORRECT_VECTOR_LAG:
        rep["note"] = "disabled (CORRECT_VECTOR_LAG = False)"
        return rep
    if mode == "none":
        rep["note"] = "no vector data"
        return rep
    if mode == "splitvec":
        rep["lag_ms"] = 0.0
        rep["note"] = ("scalar and vector are paired at source in 500 Hz mode - "
                       "no shift needed or applied")
        return rep
    if not dt_ms or dt_ms <= 0:
        rep["note"] = "could not determine dt - no shift applied"
        return rep

    # Claude 2026-09-17: at 250 Hz the pair-aligned firmware gives exactly 5 ms
    # (derived slot by slot above). The pre-alignment formula 1.5*dt-2 gave 4 ms
    # there; both floor to the same 1-row shift, so this only corrects what the
    # summary reports. Away from 250 Hz the row model itself breaks down, so the
    # old estimate is kept rather than replaced with a guess.
    if abs(dt_ms - 4.0) < 0.3:
        lag_ms = 5.0
    else:
        lag_ms = 1.5 * dt_ms - 2.0
    exact = lag_ms / dt_ms
    # Jeff 2026-08-31: round DOWN, never to nearest. A partial row shift is not a
    # thing we can do by moving rows, and under-correcting leaves the vector
    # slightly late (its original, familiar error direction) rather than pushing
    # it early into a sign flip the user would not expect.
    n = int(math.floor(exact))
    rep["lag_ms"] = lag_ms
    rep["rows_exact"] = exact
    rep["rows_shifted"] = n
    rep["residual_ms"] = lag_ms - n * dt_ms

    if n <= 0:
        rep["note"] = "computed shift rounds to 0 rows - nothing to do"
        return rep

    cols = [IDX[c] for c in VEC_COLS]
    ts_col = IDX["Timestamp(uS)"]
    total = len(rows)
    # Claude 2026-09-17: the shift moves the vector n ROWS up, which is only the
    # right thing when those n rows really span n*dt of time. A pair realign in the
    # sensor drops one output row, and a log can also be stitched or have dropped
    # lines, so check each source row's timestamp before trusting it. Where it does
    # not line up there is no contemporaneous vector for that scalar: leave it empty
    # rather than pairing it with the wrong cycle.
    want_us = n * dt_ms * 1000.0
    tol_us = dt_ms * 1000.0 * 0.5
    src_vals = []
    gapped = 0
    for i in range(total):
        src = i + n
        ok = src < total
        if ok:
            t0, t1 = rows[i][ts_col], rows[src][ts_col]
            if t0 is None or t1 is None or abs((t1 - t0) - want_us) > tol_us:
                ok = False
                gapped += 1
        src_vals.append([rows[src][c] for c in cols] if ok else [None] * len(cols))
    for i in range(total):
        for c, v in zip(cols, src_vals[i]):
            rows[i][c] = v
    rep["applied"] = True
    rep["rows_gapped"] = gapped
    notes = []
    if abs(rep["residual_ms"]) > 0.05:
        notes.append("%.2f ms residual remains (the scalar is the newer of the two); "
                     "correct in the time domain for exactness." % rep["residual_ms"])
    if gapped:
        notes.append("%d row(s) had no source row at the right time - a dropped output, "
                     "most likely a sensor pair realign after an FPGA parameter write - "
                     "so their vector cells were left empty rather than paired with the "
                     "wrong cycle. _flat.csv keeps them empty; _interpolated.csv and "
                     "_dedup.csv fill them by interpolation, as they do any missing "
                     "vector sample." % gapped)
    if dt_ms > 4.3:
        notes.append("BELOW 250 Hz the vector accumulates over data_rate/2 coil cycles, so "
                     "its ready flags arrive in bursts that do not match the row cadence and "
                     "the X>Y>Z print priority reorders them. A row shift cannot fix that; "
                     "treat vector/scalar pairing in this file as approximate.")
    rep["note"] = " ".join(notes)
    return rep


def write_summary(raw_path, rows, status_lines, out_path, vec_report=None,
                  ts_dropped=None):
    ts_col = IDX["Timestamp(uS)"]
    date_col = IDX["UTCDate(YYYY-MM-DD)"]
    utc_col = IDX["TimeUTC(HH:MM:SS.sss)"]

    n = len(rows)
    start_utc = end_utc = "N/A"
    total_sec = 0.0
    dt_ms = "N/A"
    if n:
        first, final = rows[0], rows[-1]
        if first[date_col] is not None and first[utc_col] is not None:
            start_utc = f"{first[date_col]} {first[utc_col]}"
            end_utc = f"{final[date_col]} {final[utc_col]}"
        total_sec = (final[ts_col] - first[ts_col]) / 1e6
        if n > 1:
            k = min(10, n - 1)
            dt_ms = int(round((rows[k][ts_col] - rows[0][ts_col]) / k / 1000.0))

    with open(out_path, 'w', encoding='utf-8') as f:
        f.write("==================================================\n")
        f.write(f"Summary for: {os.path.basename(raw_path)}\n")
        f.write("==================================================\n\n")
        f.write(f"Total Rows: {n}\n")
        f.write(f"Start UTC:  {start_utc}\n")
        f.write(f"End UTC:    {end_utc}\n")
        f.write(f"Total Time: {total_sec:.3f} seconds\n")
        f.write(f"dt:         {dt_ms} ms\n")

        # Claude 2026-09-28 (Jeff, for batch users): the dropped-row detail used to go
        # to the console only, which is no use to anyone converting a season of
        # logs in a batch. Put it in the per-file summary so the output can be
        # QC'd afterwards.
        total_in = n + (len(ts_dropped) if ts_dropped else 0)
        f.write("\n")
        f.write("--------------------------------------------------\n")
        f.write("Corrupt timestamp rows removed:\n")
        f.write("--------------------------------------------------\n")
        if not ts_dropped:
            f.write("None. Every Timestamp(uS) agreed with the rows around it.\n")
        else:
            nz = sum(1 for _, bad, _ in ts_dropped if bad == 0)
            f.write(f"Total removed:   {len(ts_dropped)} of {total_in} rows "
                    f"({100.0 * len(ts_dropped) / max(total_in, 1):.3f}%)\n")
            f.write(f"  exactly zero:  {nz}   (MDP pulse had not arrived; older "
                    f"firmware wrote these as '0!')\n")
            f.write(f"  out of place:  {len(ts_dropped) - nz}   (disagreed with "
                    f"the rows around it by more than 60 s)\n")
            outp = [d for d in ts_dropped if d[2] is not None]
            if outp:
                f.write("\nEvery out-of-place row removed, so each can be checked:\n")
                for _, bad, near in outp:
                    f.write(f"  stamp {bad}  neighbours near {near}  "
                            f"off by {(bad - near) / 1e6:+.1f} s\n")
            if nz:
                f.write(f"\nThe {nz} zero-stamped rows are not listed one by one. "
                        f"They arrive in bursts at\nthe start of a session and are "
                        f"dropped on the value alone, not by comparison.\n")
        if vec_report:
            r = vec_report
            f.write("\n")
            f.write("--------------------------------------------------\n")
            f.write("Vector/scalar timing correction:\n")
            f.write("--------------------------------------------------\n")
            f.write(f"Mode detected:   {r['mode']}  ({r['why']})\n")
            if r['lag_ms'] is not None:
                f.write(f"Vector lag:      {r['lag_ms']:.2f} ms behind its row's scalar\n")
            if r['rows_exact'] is not None:
                f.write(f"Exact shift:     {r['rows_exact']:.3f} rows\n")
            f.write(f"Shift applied:   {r['rows_shifted']} row(s) "
                    f"({'YES' if r['applied'] else 'none'})\n")
            if r['applied'] and abs(r['residual_ms']) > 0.05:
                f.write(f"Residual:        {r['residual_ms']:.2f} ms "
                        f"(the scalar is the newer of the two)\n")
            if r.get('rows_gapped'):
                f.write(f"Rows left empty: {r['rows_gapped']} at the shift stage, no source\n")
                f.write("                 row at the right time (a dropped output). Empty in\n")
                f.write("                 _flat.csv; the interpolation stage then fills them like\n")
                f.write("                 any other missing vector sample.\n")
            # Claude 2026-09-17: state the pairing in words. On pair-aligned sensor
            # firmware this is what the row shift actually buys, and a user reading
            # the CSV has no other way to know it happened.
            if r['applied'] and r['mode'] == 'survey':
                f.write("Pairing:         each row's vector axis now sits with the scalar\n"
                        "                 average of that same axis's two bracketing\n"
                        "                 unbiased slots (pair-aligned firmware, 2026-09-17\n"
                        "                 or later). On older firmware the scalar pairing\n"
                        "                 followed a coil-cycle phase that any FPGA\n"
                        "                 parameter write could flip, so verify against the\n"
                        "                 83.3 Hz scalar line before trusting older files.\n")
            elif r['mode'] == 'splitvec':
                f.write("Pairing:         scalar and vector are paired at source in 500 Hz\n"
                        "                 mode; the pair-alignment change does not apply.\n")
            if r['note']:
                wrapped = textwrap.fill(r['note'], width=70,
                                        initial_indent="", subsequent_indent=" " * 17)
                f.write(f"Note:            {wrapped}\n")
        f.write("\n")
        f.write("--------------------------------------------------\n")
        f.write("Status Updates:\n")
        f.write("--------------------------------------------------\n")
        if status_lines:
            f.write('\n'.join(status_lines) + '\n')
        else:
            f.write("No status updates found.\n")
    print(f"Summary file written to {out_path}")


def process_file(path):
    base, _ = os.path.splitext(path)
    flat_csv = f"{base}_flat.csv"
    interp_csv = f"{base}_interpolated.csv"
    dedup_csv = f"{base}_dedup.csv"
    summary_txt = f"{base}_summary.txt"

    print(f"\n=== {os.path.basename(path)} ===")
    rows, status_lines = parse_file(path)
    print(f"Parsed {len(rows)} data rows")
    if not rows:
        print("No parseable data found — skipping.")
        return

    rows, ts_dropped = drop_corrupt_timestamps(rows)
    if ts_dropped:
        print(f"Corrupt timestamps: dropped {len(ts_dropped)} line(s) whose "
              f"Timestamp(uS) does not belong with the lines around them.")
        nz = sum(1 for _, bad, near in ts_dropped if bad == 0)
        if nz:
            print(f"    {nz} of them had Timestamp(uS) == 0 (MDP pulse had not "
                  f"arrived; older Magpie firmware wrote these as '0!')")
        for _, bad, near in [d for d in ts_dropped if d[2] is not None][:4]:
            print(f"    {bad}  (neighbours are near {near}, "
                  f"off by {(bad - near) / 1e6:+.1f} s)")
        if len(ts_dropped) > 5:
            print(f"    ... and {len(ts_dropped) - 5} more")
        print("    This is SD/UART corruption in the raw log, not a converter "
              "fault. Worth reporting if it is frequent.")

    # Align the vector columns with their own scalar BEFORE anything downstream
    # interpolates or dedups - correcting afterwards would smear the shift.
    vec_report = correct_vector_lag(rows, status_lines)

    # Claude 2026-09-02 (Jeff): close the SplitVec polarity gap. Runs AFTER the
    # lag stage (which is a no-op in split mode) and only when a calibration
    # file supplies coefficients - absent file means this silently does nothing.
    cal = load_cal(path)
    pol_report = apply_polarity_stage(rows, vec_report.get("mode"), cal)
    if pol_report["applied"]:
        bits = ", ".join("%s %.1f->%.1f nT" % (a, b, af)
                         for a, (b, af, ok) in sorted(pol_report["per_axis"].items()) if ok)
        print("Polarity gap corrected: " + bits)
    elif pol_report["note"]:
        print("Polarity correction: " + pol_report["note"])
    if vec_report["applied"]:
        print(f"Vector lag corrected: shifted vector columns up "
              f"{vec_report['rows_shifted']} row(s) "
              f"({vec_report['lag_ms']:.1f} ms, {vec_report['mode']})")
    else:
        print(f"Vector lag: no shift applied ({vec_report['mode']}"
              f"{' - ' + vec_report['note'] if vec_report['note'] else ''})")

    write_csv(flat_csv, rows)
    print(f"Flat file written to {flat_csv}")

    sensor_rows = process(rows)
    write_csv(interp_csv, sensor_rows)
    print(f"Interpolated file written to {interp_csv}")

    write_csv(dedup_csv, dedup_rows(sensor_rows))
    print(f"Dedup file written to {dedup_csv}")

    write_summary(path, sensor_rows, status_lines, summary_txt, vec_report,
                  ts_dropped)


def main():
    if TEST_FILE:
        paths = [TEST_FILE]
    else:
        if len(sys.argv) < 2:
            print("Usage: drag one or more Magpie .txt log files onto this "
                  "program,\nor run: Qlog_ConverterV5.exe <file1.txt> [file2.txt ...]")
            pause_if_frozen()
            sys.exit(1)
        paths = sys.argv[1:]

    for path in paths:
        if not os.path.isfile(path):
            print(f"Skipped missing: {path}")
            continue
        try:
            process_file(path)
        except Exception as e:
            print(f"ERROR processing {path}: {e}")

    print("\nAll done.")
    pause_if_frozen()


def pause_if_frozen():
    # Keep the console open for drag-and-drop users running the packaged exe.
    if getattr(sys, 'frozen', False):
        try:
            input("\nPress Enter to close...")
        except EOFError:
            pass


if __name__ == "__main__":
    try:
        main()
    except Exception as e:
        print(f"FATAL: {e}")
        pause_if_frozen()
        raise