"""Our analysis for the page on Baird et al. (2026), GLE 77 and flight radiation, from public records.

Usage:  python3 -I gle77_records.py [FILE ...] [OUT.json]

With no FILE arguments it fetches its eight inputs (about 1.2 MB) and prints the SHA-256 of each; any input you
already have can be given as a path instead, matched by file name. It treats every input as data and never runs it;
the 2025 event lists are read from their tarball in memory, never extracted. Standard library only.

Inputs
  alerts_202511.html, alerts_202601.html  NOAA SWPC's alerts, warnings and summaries for November 2025 and
                                           January 2026, ftp://ftp.swpc.noaa.gov/pub/alerts/
  2025_events.tar.gz                       NOAA SWPC's daily event lists for 2025, ftp://ftp.swpc.noaa.gov/pub/warehouse/2025/
  kp_2019-12-01_2026-10-28.json            GFZ's Kp index (CC BY 4.0; Matzka et al. 2021, doi:10.1029/2020SW002641;
                                           data doi:10.5880/Kp.0001), from https://kp.gfz.de/app/json/
  dst_provisional_YYYYMM.html (4 months)   WDC for Geomagnetism, Kyoto, provisional Dst (doi:10.17593/14515-74000)

It prints
  1. SWPC's proton, X-ray and storm-watch products of 10-13 November 2025, and the minutes to the first 100 MeV
     proton product from GLE 77's onset: 10:10 UTC at Newark (Mishev et al. 2026, Advances in Space Research 77,
     p. 9770) and 10:15 UTC (Baird et al. 2026, Space Weather 24, e2026SW005279, p. 4).
  2. The X-class flares in SWPC's 2025 event lists, largest first.
  3. January 2026: SWPC's largest proton summary and whether a 100 MeV product was issued (18-22 January).
  4. The storms of solar cycle 25 (from December 2019) whose highest three-hour Kp reached 8- (NOAA's G4) or more,
     consecutive days merged, and Kyoto's lowest provisional Dst in the months of the four that reached Kp 9-.
  5. Hours at ICAO's thresholds (30 and 80 microsieverts an hour, effective dose rate) as doses.
"""
import datetime
import hashlib
import html
import io
import json
import os
import re
import sys
import tarfile
import urllib.request

KP_URL = 'https://kp.gfz.de/app/json/?start=2019-12-01T00:00:00Z&end=2026-10-28T00:00:00Z&index=Kp'
INPUTS = {
    'alerts_202511.html': 'ftp://ftp.swpc.noaa.gov/pub/alerts/alerts_202511.html',
    'alerts_202601.html': 'ftp://ftp.swpc.noaa.gov/pub/alerts/alerts_202601.html',
    '2025_events.tar.gz': 'ftp://ftp.swpc.noaa.gov/pub/warehouse/2025/2025_events.tar.gz',
    'kp_2019-12-01_2026-10-28.json': KP_URL,
    'dst_provisional_202405.html': 'https://wdc.kugi.kyoto-u.ac.jp/dst_provisional/202405/index.html',
    'dst_provisional_202410.html': 'https://wdc.kugi.kyoto-u.ac.jp/dst_provisional/202410/index.html',
    'dst_provisional_202511.html': 'https://wdc.kugi.kyoto-u.ac.jp/dst_provisional/202511/index.html',
    'dst_provisional_202601.html': 'https://wdc.kugi.kyoto-u.ac.jp/dst_provisional/202601/index.html',
}
ONSETS = {'10:10 (Mishev et al. 2026, p. 9770)': datetime.datetime(2025, 11, 11, 10, 10),
          '10:15 (Baird et al. 2026, p. 4)': datetime.datetime(2025, 11, 11, 10, 15)}
MOD, SEV = 30.0, 80.0                                  # ICAO Doc 10100 thresholds, microsieverts an hour


def load(args):
    given = {os.path.basename(a): a for a in args}
    data, sums = {}, {}
    for name, url in INPUTS.items():
        if name in given:
            raw = open(given[name], 'rb').read()
        else:
            req = urllib.request.Request(url, headers={'User-Agent': 'Mozilla/5.0 (fact-check script)'})
            with urllib.request.urlopen(req, timeout=120) as r:
                raw = r.read()
        sums[name] = hashlib.sha256(raw).hexdigest()
        data[name] = raw if name.endswith('.gz') else raw.decode('utf-8', 'replace')
    return data, sums


def swpc_messages(text):
    out = []
    for block in re.split(r'(?i)<hr>', text):
        t = html.unescape(re.sub(r'<[^>]+>', '', re.sub(r'(?i)<br\s*/?>', '\n', block))).strip()
        code = re.search(r'Space Weather Message Code:\s*(\S+)', t)
        issued = re.search(r'Issue Time:\s*(\d{4} \w{3} \d{2} \d{4}) UTC', t)
        if code and issued:
            when = datetime.datetime.strptime(issued.group(1), '%Y %b %d %H%M')
            lines = [ln.strip() for ln in t.splitlines() if ln.strip()]
            out.append((when, code.group(1), lines))
    return sorted(out)


def field(lines, key):
    for ln in lines:
        if ln.startswith(key):
            return ln[len(key):].strip()
    return ''


def xflares(raw):
    rows = []
    with tarfile.open(fileobj=io.BytesIO(raw), mode='r:gz') as tf:
        for m in tf.getmembers():
            if not m.isfile() or not m.name.endswith('events.txt'):
                continue
            day = re.search(r'(\d{8})events\.txt$', m.name).group(1)
            if not day.startswith('2025'):
                continue
            for line in tf.extractfile(m).read().decode('ascii', 'replace').splitlines():
                p = line.split()
                if 'XRA' in p and len(p) > p.index('XRA') + 2:
                    i = p.index('XRA')
                    cls = p[i + 2]
                    if cls.startswith('X'):
                        rows.append((float(cls[1:]), day, cls, p[i + 4] if len(p) > i + 4 else ''))
    return sorted(rows, key=lambda r: -r[0])


def label(k):
    base, rem = divmod(round(k * 3), 3)
    return {0: f'{base}o', 1: f'{base}+', 2: f'{base + 1}-'}[rem]


def dst_min(text):
    t = html.unescape(re.sub(r'<[^>]+>', '', text))
    best, in_table = None, False
    for line in t.splitlines():
        if line.strip() == 'DAY':
            in_table = True
            continue
        m = re.match(r'^\s*(\d{1,2})\s+(.*)$', line) if in_table else None
        if m:
            vals = [int(v) for v in re.findall(r'-?\d+', m.group(2))]   # 4-character fields, e.g. '-98-116'
            if len(vals) == 24:
                for h, v in enumerate(vals):
                    if best is None or v < best[0]:
                        best = (v, int(m.group(1)), h + 1)
    return best


def main():
    sys.stdout.reconfigure(encoding='utf-8')
    args = [a for a in sys.argv[1:] if not a.endswith('.json') or os.path.basename(a) in INPUTS]
    out_json = next((a for a in sys.argv[1:] if a.endswith('.json') and os.path.basename(a) not in INPUTS), None)
    data, sums = load(args)
    res = {'inputs_sha256': sums}
    for name, s in sums.items():
        print(f'{name}: sha256 {s}')

    print('\n1. SWPC products, 10-13 November 2025 (UTC issue times)')
    keep = ('ALTXMF', 'SUMX01', 'WARPC0', 'ALTPC0', 'SUMPC0', 'WARPX1', 'ALTPX1', 'ALTPX2', 'ALTPX3', 'SUMPX2',
            'SUMPX3', 'WATA99')
    nov = [m for m in swpc_messages(data['alerts_202511.html'])
           if datetime.date(2025, 11, 10) <= m[0].date() <= datetime.date(2025, 11, 13) and m[1] in keep]
    res['nov_2025'] = []
    for when, code, lines in nov:
        head = next((ln for ln in lines if ln.split(':')[0] in ('ALERT', 'WARNING', 'EXTENDED WARNING', 'SUMMARY',
                                                                   'WATCH', 'CONTINUED ALERT')), lines[3] if len(lines) > 3 else '')
        extra = '; '.join(x for x in (field(lines, 'X-ray Class:'), field(lines, 'NOAA Scale:'),
                                      field(lines, 'Maximum 10MeV Flux:'), field(lines, 'Maximum 100MeV Flux:')) if x)
        print(f'  {when:%d %b %H:%M}  {code:7s} {head}' + (f'  [{extra}]' if extra else ''))
        res['nov_2025'].append({'issued': f'{when:%Y-%m-%d %H:%M}', 'code': code, 'text': head, 'details': extra})
    res['minutes_onset_to_first_100MeV_product'] = {}
    for name, onset in ONSETS.items():
        first = min((w for w, c, _ in nov if c in ('WARPC0', 'ALTPC0') and w >= onset), default=None)
        if first:
            mins = int((first - onset).total_seconds() // 60)
            print(f'  Minutes from the onset at {name} to SWPC\'s first 100 MeV proton product: {mins}')
            res['minutes_onset_to_first_100MeV_product'][name] = mins

    print('\n2. X-class flares in SWPC\'s 2025 event lists, largest first')
    fl = xflares(data['2025_events.tar.gz'])
    for size, day, cls, reg in fl[:5]:
        print(f'  {cls:6s} {day[:4]}-{day[4:6]}-{day[6:]}  region {reg}')
    print(f'  {len(fl)} X-class flares in 2025')
    res['xflares_2025'] = [{'class': c, 'day': d, 'region': r} for _, d, c, r in fl]

    print('\n3. January 2026')
    jan = [m for m in swpc_messages(data['alerts_202601.html'])
           if datetime.date(2026, 1, 18) <= m[0].date() <= datetime.date(2026, 1, 22)]
    s4 = [m for m in jan if m[1] == 'SUMPX4']
    for when, code, lines in s4:
        print(f'  {when:%d %b %H:%M} {code}: maximum {field(lines, "Maximum 10MeV Flux:")} at '
              f'{field(lines, "Maximum Time:")}, {field(lines, "NOAA Scale:")}')
    p100 = [m for m in jan if m[1] in ('WARPC0', 'ALTPC0', 'SUMPC0')]
    print(f'  100 MeV proton warnings, alerts or summaries issued 18-22 January: {len(p100)}')
    res['jan_2026'] = {'s4_summaries': [field(m[2], 'Maximum 10MeV Flux:') for m in s4], 'n_100MeV_products': len(p100)}

    print('\n4. Solar cycle 25 storms at G4 or more, by GFZ Kp')
    kp = json.loads(data['kp_2019-12-01_2026-10-28.json'])
    days = {}
    for t, k in zip(kp['datetime'], kp['Kp']):
        days.setdefault(t[:10], []).append(k)
    big = sorted((d, max(v)) for d, v in days.items() if max(v) >= 7.66)
    storms = []
    for d, mx in big:
        dt = datetime.date.fromisoformat(d)
        if storms and (dt - storms[-1]['last']).days == 1:
            storms[-1]['last'], storms[-1]['max'] = dt, max(storms[-1]['max'], mx)
        else:
            storms.append({'first': dt, 'last': dt, 'max': mx})
    for s in storms:
        print(f'  {s["first"]} to {s["last"]}: highest Kp {label(s["max"])}')
    top = [s for s in storms if s['max'] >= 8.6]
    print(f'  {len(storms)} storms reached Kp 8- or more; {len(top)} reached Kp 9- or more '
          f'({sum(1 for s in storms if s["max"] >= 8.99)} reached Kp 9o, NOAA\'s G5)')
    res['cycle25'] = [{'first': str(s['first']), 'last': str(s['last']), 'max_kp': label(s['max'])} for s in storms]
    for name in ('dst_provisional_202405.html', 'dst_provisional_202410.html', 'dst_provisional_202511.html',
                 'dst_provisional_202601.html'):
        v, d, h = dst_min(data[name])
        print(f'  {name[-11:-5]}: lowest provisional Dst {v} nT on day {d}, hour ending {h:02d} UT')
        res.setdefault('dst_min', {})[name[-11:-5]] = v

    print('\n5. Hours at ICAO\'s thresholds, as doses')
    for hours in (7, 8):
        print(f'  {hours} h at {MOD:.0f} uSv/h = {MOD * hours / 1000:.2f} mSv; at {SEV:.0f} uSv/h = {SEV * hours / 1000:.2f} mSv')
    res['dose'] = {'7h_mod_mSv': MOD * 7 / 1000, '7h_sev_mSv': SEV * 7 / 1000}
    if out_json:
        with open(out_json, 'w', encoding='utf-8', newline='\n') as f:
            json.dump(res, f, indent=1, ensure_ascii=False)
        print(f'\nSaved {out_json}')


if __name__ == '__main__':
    main()