"""Check Baird et al. (2026)'s headline numbers against its own data ("Our analysis").

Data: Baird, F., et al. (2026). Observations of GLE 77 from Radiosondes, Aircraft, and UK Neutron Monitors,
Zenodo, version 2, https://doi.org/10.5281/zenodo.23082794, CC BY 4.0 (the record says version 1's aircraft and
radiosonde files used a faulty calibration and that version 2 matches the paper). Read as data only.

Usage: python3 -I gle77_data_check.py DATA_DIR [OUT.json]
DATA_DIR holds the record's CSV files under their Zenodo names (the radiosonde files end in _GLE77.csv or
_20251017.csv; FileDescriptions.csv calls the GLE files _GLE.csv). Standard library only.

What it does
- Radiosondes: the ascent only (up to the highest GPS altitude). Altitudes are the files' "Flight Level" column,
  pressure altitude in hundreds of feet, the measure ICAO's advisories use. One-minute samples of the ambient dose
  equivalent rate H*(10) are noisy, so we also take a running median of three samples ("smoothed"). For each flight:
  the smoothed rate at flight levels 250-600, the lowest flight level above which the smoothed rate stays above
  ICAO's moderate (30) and severe (80) thresholds, in microsieverts an hour, up to FL600, and the first run of three
  raw samples above each. The three quiet-day files are flights of 17 October 2025, as their names and
  FileDescriptions.csv say, although their timestamps read 2025-11-11.
- Flights: the total H*(10) from the cumulative column; the one-minute samples above 30 uSv/h; a trailing ten-minute
  mean (ours); the cruise mean. The "Altitude (m)" column holds feet (cruise values of 37,000).
- Neutron monitors: the percentage increase over the mean of 08:00-10:00 UTC, the baseline Mishev et al. (2026) use,
  with a ten-minute running mean, and its peak and time.
ICAO's thresholds are effective dose rates; the files give H*(10). We compare them as the paper's abstract does.
"""
import csv
import json
import os
import statistics as st
import sys
from datetime import datetime, timedelta

MOD, SEV = 30.0, 80.0
LEVELS = [250, 300, 350, 400, 450, 500, 550, 600]
SONDES = [('Lerwick_Radiosonde_Flight_1_GLE77.csv', 'Lerwick, launched 11:12 UTC'),
          ('Lerwick_Radiosonde_Flight_2_GLE77.csv', 'Lerwick, launched 15:45 UTC'),
          ('Lerwick_Radiosonde_GCR_20251017.csv', 'Lerwick, quiet day, 17 October 2025'),
          ('Camborne_Radiosonde_GLE77.csv', 'Camborne, launched 15:45 UTC'),
          ('Camborne_Radiosonde_GCR_20251017.csv', 'Camborne, quiet day, 17 October 2025'),
          ('DeBilt_Radiosonde_GLE77.csv', 'De Bilt, launched 14:15 UTC'),
          ('DeBilt_Radiosonde_GCR_20251017.csv', 'De Bilt, quiet day, 17 October 2025')]
FLIGHTS = [('LHR_JFK.csv', 'Heathrow–JFK'), ('LHR_EWR.csv', 'Heathrow–Newark')]
MONITORS = ['Lerwick.csv', 'Lancaster.csv', 'Camborne.csv', 'Guildford.csv']


def rows(path):
    with open(path, newline='', encoding='utf-8-sig', errors='replace') as f:
        r = list(csv.reader(f))
    return r[0], r[1:]


def col(header, *keys, exclude=()):
    for i, h in enumerate(header):
        hl = h.lower()
        if all(k in hl for k in keys) and not any(x in hl for x in exclude):
            return i
    raise KeyError(keys)


def t(s):
    return datetime.strptime(s.strip()[:19], '%Y-%m-%d %H:%M:%S')


def median3(v):
    return [st.median(v[max(0, i - 1):i + 2]) for i in range(len(v))]


def interp(xs, ys, x):
    pts = sorted(zip(xs, ys))
    for (x0, y0), (x1, y1) in zip(pts, pts[1:]):
        if x0 <= x <= x1 and x1 > x0:
            return y0 + (y1 - y0) * (x - x0) / (x1 - x0)
    return None


def sonde(path):
    h, rs = rows(path)
    it, ia, ifl = col(h, 'datetime'), col(h, 'gps altitude'), col(h, 'flight level')
    irc, ih = col(h, 'rigidity'), col(h, 'ambient dose')
    data = [(t(r[it]), float(r[ia]), float(r[ifl]), float(r[irc]), float(r[ih])) for r in rs if len(r) > ih]
    top = max(range(len(data)), key=lambda i: data[i][1])
    asc = data[:top + 1]
    fl = [d[2] for d in asc]
    raw = [d[4] for d in asc]
    sm = median3(raw)
    out = {'samples': len(data), 'ascent_samples': len(asc), 'start': f'{asc[0][0]:%H:%M}',
           'top_fl': round(asc[-1][2]), 'top_time': f'{asc[-1][0]:%H:%M}',
           'cutoff_GV': round(st.median(d[3] for d in data), 2), 'at_level': {}, 'time_at_level': {},
           'profile': [[round(f, 1), round(s, 2), round(r, 2)] for f, s, r in zip(fl, sm, raw) if f <= 650]}
    for L in LEVELS:
        v = interp(fl, sm, L)
        out['at_level'][L] = None if v is None else round(v, 1)
        tl = interp(fl, [(d[0] - asc[0][0]).total_seconds() for d in asc], L)
        out['time_at_level'][L] = None if tl is None else f'{asc[0][0] + timedelta(seconds=tl):%H:%M}'
    band = [(f, s, r) for f, s, r in zip(fl, sm, raw) if 250 <= f <= 600]
    for name, thr in (('mod', MOD), ('sev', SEV)):
        lvl = None   # lowest flight level above which the smoothed rate stays above thr, up to FL600
        for i in range(len(band) - 1, -1, -1):
            if band[i][1] > thr:
                lvl = band[i][0]
            else:
                break
        out[f'{name}_stays_above_from_FL'] = None if lvl is None else round(lvl)
        k, first = 0, None
        for f, s, r in band:
            k = k + 1 if r > thr else 0
            if k == 3:
                first = f
                break
        out[f'{name}_first_run3_FL'] = None if first is None else round(first)
    out['max_raw_FL250_600'] = max(((r, round(f)) for f, s, r in band), default=(None, None))
    out['max_smoothed_FL250_600'] = max(((round(s, 1), round(f)) for f, s, r in band), default=(None, None))
    return out


def flight(path):
    h, rs = rows(path)
    it, ia = col(h, 'datetime'), col(h, 'altitude')
    ic, ir = col(h, 'ambient dose equivalent (usv)'), col(h, 'ambient dose equivalent rate')
    irc = col(h, 'rigidity')
    data = [(t(r[it]), float(r[ia]), float(r[ic]), float(r[ir]), float(r[irc])) for r in rs if len(r) > ir]
    rate = [d[3] for d in data]
    roll = [st.mean(rate[max(0, i - 9):i + 1]) for i in range(len(rate))]
    cruise = [d for d in data if d[1] >= 30000]
    above = [(f'{d[0]:%H:%M}', d[3], round(d[1])) for d in data if d[3] > MOD]
    imax = max(range(len(roll)), key=lambda i: roll[i])
    total = data[-1][2] - data[0][2]
    hours = (data[-1][0] - data[0][0]).total_seconds() / 3600
    return {'samples': len(data), 'start': f'{data[0][0]:%H:%M}', 'end': f'{data[-1][0]:%Y-%m-%d %H:%M}',
            'hours': round(hours, 2), 'total_uSv': round(total, 2), 'sum_of_rates_uSv': round(sum(rate) / 60, 2),
            'max_altitude_ft': max(d[1] for d in data), 'cruise_mean_uSv_h': round(st.mean(d[3] for d in cruise), 2),
            'cruise_minutes': len(cruise), 'samples_above_30': above,
            'roll10_max': [round(roll[imax], 2), f'{data[imax][0]:%H:%M}'],
            'roll10_minutes_above_30': sum(1 for v in roll if v > MOD),
            'series': [[f'{d[0]:%Y-%m-%d %H:%M}', round(d[3], 2), round(r10, 2), round(d[1])] for d, r10 in zip(data, roll)]}


def monitor(path):
    h, rs = rows(path)
    data = [(t(r[0]), float(r[1])) for r in rs if len(r) > 1 and r[1].strip()]
    base = st.mean(c for d, c in data if 8 <= d.hour < 10)
    pct = [(d, (c / base - 1) * 100) for d, c in data]
    sm = [(pct[i][0], st.mean(p for _, p in pct[max(0, i - 4):i + 6])) for i in range(len(pct))]
    after = [x for x in sm if x[0].hour >= 10]
    peak = max(after, key=lambda x: x[1])
    later = [x for x in after if x[0] >= peak[0] + timedelta(minutes=40) and x[0].hour < 14]
    second = max(later, key=lambda x: x[1]) if later else None
    return {'samples': len(data), 'baseline_08_10': round(base, 4), 'peak_pct': round(peak[1], 1),
            'peak_time': f'{peak[0]:%H:%M}', 'later_max_pct': None if not second else round(second[1], 1),
            'later_max_time': None if not second else f'{second[0]:%H:%M}',
            'series': [[f'{d:%H:%M}', round(p, 2)] for d, p in sm if 9 <= d.hour < 17]}


def main():
    sys.stdout.reconfigure(encoding='utf-8')
    d = sys.argv[1]
    out = {'sondes': {}, 'flights': {}, 'monitors': {}}
    print('== Radiosondes: smoothed H*(10) rate (uSv/h) at flight levels, ascent only')
    for f, label in SONDES:
        r = sonde(os.path.join(d, f))
        r['label'] = label
        out['sondes'][f] = r
        lv = '  '.join(f'FL{L}:{r["at_level"][L]}' for L in LEVELS)
        print(f'\n{f} ({label}): {r["samples"]} samples, ascent {r["ascent_samples"]} from {r["start"]} to FL{r["top_fl"]} '
              f'at {r["top_time"]}; cutoff {r["cutoff_GV"]} GV')
        print(f'  {lv}')
        print(f'  smoothed stays above 30 from FL{r["mod_stays_above_from_FL"]}, above 80 from FL{r["sev_stays_above_from_FL"]}; '
              f'first run of 3 raw samples above 30 at FL{r["mod_first_run3_FL"]}, above 80 at FL{r["sev_first_run3_FL"]}')
        print(f'  within FL250-FL600: max raw {r["max_raw_FL250_600"]}, max smoothed {r["max_smoothed_FL250_600"]}; '
              f'times at FL400 {r["time_at_level"][400]}, FL550 {r["time_at_level"][550]}, FL600 {r["time_at_level"][600]}')
    for g, q in (('Lerwick_Radiosonde_Flight_1_GLE77.csv', 'Lerwick_Radiosonde_GCR_20251017.csv'),
                 ('Lerwick_Radiosonde_Flight_2_GLE77.csv', 'Lerwick_Radiosonde_GCR_20251017.csv'),
                 ('Camborne_Radiosonde_GLE77.csv', 'Camborne_Radiosonde_GCR_20251017.csv'),
                 ('DeBilt_Radiosonde_GLE77.csv', 'DeBilt_Radiosonde_GCR_20251017.csv')):
        a, b = out['sondes'][g]['at_level'], out['sondes'][q]['at_level']
        ratios = {L: round(a[L] / b[L], 1) for L in LEVELS if a[L] and b[L]}
        out['sondes'][g]['ratio_to_quiet'] = ratios
        print(f'  {g} / quiet day: ' + '  '.join(f'FL{L} x{v}' for L, v in ratios.items()))
    print('\n== Flights')
    for f, label in FLIGHTS:
        r = flight(os.path.join(d, f))
        out['flights'][f] = r
        print(f'\n{f} ({label}): {r["start"]} to {r["end"]} UTC, {r["hours"]} h; max altitude {r["max_altitude_ft"]:.0f} ft')
        print(f'  total {r["total_uSv"]} uSv (sum of rates {r["sum_of_rates_uSv"]}); cruise mean {r["cruise_mean_uSv_h"]} uSv/h '
              f'over {r["cruise_minutes"]} min')
        print(f'  one-minute samples above 30: {r["samples_above_30"]}')
        print(f'  ten-minute mean: max {r["roll10_max"]}; minutes above 30: {r["roll10_minutes_above_30"]}')
    jfk = out['flights']['LHR_JFK.csv']['total_uSv']
    out['jfk_if_426pct'] = {'typical_uSv': round(jfk / 1.426, 1), 'excess_uSv': round(jfk - jfk / 1.426, 1)}
    print(f'  JFK: if 42.6% above typical, typical {jfk / 1.426:.1f} uSv and the excess {jfk - jfk / 1.426:.1f} uSv')
    print('\n== Neutron monitors: % over the 08:00-10:00 UTC mean, ten-minute running mean')
    for f in MONITORS:
        r = monitor(os.path.join(d, f))
        out['monitors'][f] = r
        print(f'  {f}: peak {r["peak_pct"]}% at {r["peak_time"]}; later maximum {r["later_max_pct"]}% at {r["later_max_time"]}')
    if len(sys.argv) > 2:
        with open(sys.argv[2], 'w', encoding='utf-8', newline='\n') as fh:
            json.dump(out, fh, ensure_ascii=False, indent=1)
        print(f'\nSaved {sys.argv[2]}')


if __name__ == '__main__':
    main()