#!/usr/bin/env python3
"""Cycle 25's largest geomagnetic storms and proton events, behind "S0 News, October 28, 2026".
Fact-Checked Science, section 5 ("Our analysis"). Standard library only.

Run it
  python3 -I storm_numbers.py [folder] [--until YYYY-MM-DD]
  --until (default 2026-10-27, the last day of the Kp record the page used) ignores everything after that day, so
  later storms and proton events do not change the page's counts and ranks; pass a later date to bring them in.
  The folder (default: the script's own) must hold, or will receive (about 6 MB in all):
    Kp_ap_Ap_SN_F107_since_1932.txt  https://kp.gfz.de/app/files/Kp_ap_Ap_SN_F107_since_1932.txt
        Kp from Geomagnetic Observatory Niemegk, GFZ Helmholtz Centre for Geosciences, CC BY 4.0: Matzka, Stolle,
        Yamazaki, Bronkalla and Morschhauser (2021), Space Weather, https://doi.org/10.1029/2020SW002641; data
        https://doi.org/10.5880/Kp.0001. The file grows daily; its newest values are nowcasts. Its sunspot numbers
        (CC BY-NC 4.0) are not used.
    dst_YYYYMM.html                  https://wdc.kugi.kyoto-u.ac.jp/dst_provisional/YYYYMM/index.html
        Hourly Dst from the World Data Center for Geomagnetism, Kyoto, for every month of cycle 25 with a day of
        Kp 7+ or more: the provisional values the page used; where a month has none, its final values (dst_final),
        and failing those its real-time values (dst_realtime). Final values can differ from provisional ones by a
        few nanotesla.
    spe_list.html                    https://www.ngdc.noaa.gov/stp/space-weather/interplanetary-data/solar-proton-events/SEP%20page%20code.html
        NOAA SWPC's "Solar Proton Events Affecting the Earth Environment, 1976 - Present", as hosted by NCEI.

What it prints
  1. Kp: every cycle-25 storm (consecutive days merged) whose highest three-hour Kp reached 8- or more, NOAA's G4
     (Kp 8, including 9-) or G5 (Kp 9o), with its highest daily Ap, that day's rank among all days since 1932, and
     the lowest hourly Dst on its days and the day after. Cycle 25 starts in December 2019 (WDC-SILSO's minimum).
  2. For comparison, the number of storms reaching Kp 9- or more in each cycle since 1933.
  3. The lowest hourly Dst of each month checked, deepest first.
  4. The proton events since December 2019, ranked by their peak flux above 10 MeV, with NOAA's S level
     (S1 >= 10, S2 >= 100, S3 >= 1,000, S4 >= 10,000, S5 >= 100,000 particle flux units).
  Kp thirds: 7.667 = 8-, 8.000 = 8o, 8.333 = 8+, 8.667 = 9-, 9.000 = 9o.
"""
import datetime as dt
import html
import os
import re
import sys
import urllib.error
import urllib.request

KP_URL = "https://kp.gfz.de/app/files/Kp_ap_Ap_SN_F107_since_1932.txt"
DST_URLS = ["https://wdc.kugi.kyoto-u.ac.jp/dst_provisional/{ym}/index.html",
            "https://wdc.kugi.kyoto-u.ac.jp/dst_final/{ym}/index.html",
            "https://wdc.kugi.kyoto-u.ac.jp/dst_realtime/{ym}/index.html"]
SPE_URL = "https://www.ngdc.noaa.gov/stp/space-weather/interplanetary-data/solar-proton-events/SEP%20page%20code.html"
UNTIL = dt.date(2026, 10, 27)
CYCLE_STARTS = [(17, dt.date(1933, 9, 1)), (18, dt.date(1944, 2, 1)), (19, dt.date(1954, 4, 1)), (20, dt.date(1964, 10, 1)),
                (21, dt.date(1976, 3, 1)), (22, dt.date(1986, 9, 1)), (23, dt.date(1996, 8, 1)), (24, dt.date(2008, 12, 1)),
                (25, dt.date(2019, 12, 1))]
KPNAME = {7.667: "8-", 8.0: "8o", 8.333: "8+", 8.667: "9-", 9.0: "9o"}


def kpname(k):
    return next((n for v, n in KPNAME.items() if abs(k - v) < 0.01), f"{k:.3f}")


def fetch(folder, name, urls):
    path = os.path.join(folder, name)
    if not os.path.exists(path):
        for url in urls:
            try:
                print("downloading", url)
                req = urllib.request.Request(url, headers={"User-Agent": "factchecked.science storm_numbers.py"})
                with urllib.request.urlopen(req, timeout=300) as r:
                    data = r.read()
                with open(path, "wb") as f:
                    f.write(data)
                break
            except urllib.error.HTTPError as e:
                print("  not there:", e.code)
    return path


def read_kp(path, until):
    rows = []
    for line in open(path, encoding="utf-8", errors="replace"):
        if line.startswith("#") or not line.strip():
            continue
        f = line.split()
        day = dt.date(int(f[0]), int(f[1]), int(f[2]))
        kps = [float(x) for x in f[7:15]]
        if day <= until and max(kps) >= 0:
            rows.append((day, max(kps), int(f[23])))
    return rows


def storms(rows, thr):
    out = []
    for r in (r for r in rows if r[1] >= thr):
        if out and (r[0] - out[-1][-1][0]).days <= 1:
            out[-1].append(r)
        else:
            out.append([r])
    return out


def read_dst(path):
    """{date: [24 hourly values]} from a Kyoto monthly page, and the version the page names."""
    txt = re.sub(r"<[^>]+>", "", open(path, encoding="utf-8", errors="replace").read())
    m = re.search(r"\n\s+([A-Z]+)\s+(\d{4})\s*\n", txt)
    month = dt.datetime.strptime(m.group(1).title() + " " + m.group(2), "%B %Y").date()
    ver = re.search(r"Hourly Equatorial Dst Values\s*\(([^)]*)\)", txt)
    days = {}
    for line in txt.splitlines():
        mm = re.match(r"^\s*(\d{1,2})\s+(-?\d.*)$", line)
        if mm:
            vals = [int(v) for v in re.findall(r"-?\d+", mm.group(2).replace("-", " -"))]
            if len(vals) >= 24:
                days[month.replace(day=int(mm.group(1)))] = [v for v in vals[:24] if v != 9999]
    return days, (ver.group(1) if ver else "?")


def main():
    sys.stdout.reconfigure(encoding="utf-8")
    args = sys.argv[1:]
    until = UNTIL
    if "--until" in args:
        i = args.index("--until")
        until = dt.date.fromisoformat(args[i + 1])
        del args[i:i + 2]
    folder = args[0] if args else os.path.dirname(os.path.abspath(__file__))
    rows = read_kp(fetch(folder, "Kp_ap_Ap_SN_F107_since_1932.txt", [KP_URL]), until)
    rank = {r[0]: i + 1 for i, r in enumerate(sorted(rows, key=lambda r: -r[2]))}
    c25 = [r for r in rows if r[0] >= CYCLE_STARTS[-1][1]]

    months = sorted({(r[0].year, r[0].month) for r in c25 if r[1] >= 7.33})
    dst, versions = {}, {}
    for y, m in months:
        ym = f"{y:04d}{m:02d}"
        p = fetch(folder, f"dst_{ym}.html", [u.format(ym=ym) for u in DST_URLS])
        if os.path.exists(p):
            d, ver = read_dst(p)
            dst.update(d)
            versions[ym] = ver

    print(f"GFZ Kp record: {rows[0][0]} to {rows[-1][0]} (cut at {until}); cycle 25 from {CYCLE_STARTS[-1][1]}\n")
    print("1. Cycle-25 storms reaching Kp 8- or more (G4 or G5), consecutive days merged")
    print("   dates                     highest Kp  highest daily Ap (rank since 1932)  lowest Dst, nT")
    big = storms(c25, 7.66)
    for s in big:
        top = max(s, key=lambda r: (r[1], r[2]))
        apday = max(s, key=lambda r: r[2])
        span = [s[0][0] + dt.timedelta(days=i) for i in range((s[-1][0] - s[0][0]).days + 2)]
        vals = [v for d in span for v in dst.get(d, [])]
        low = min(vals) if vals else None
        flag = "  <- Kp 9- or more" if top[1] >= 8.66 else ""
        print(f"   {s[0][0]} to {s[-1][0]}   {kpname(top[1]):>3}         {apday[2]:3d} on {apday[0]} ({rank[apday[0]]:3d})"
              f"         {low if low is not None else 'n/a':>5}{flag}")
    print(f"   {len(big)} storms; {sum(1 for s in big if max(r[1] for r in s) >= 8.66)} reached Kp 9- or more; "
          f"{sum(1 for s in big if max(r[1] for r in s) >= 8.99)} reached Kp 9o (G5)")

    print("\n2. Storms reaching Kp 9- or more, per cycle (cycle 25 to date)")
    for i, (cyc, start) in enumerate(CYCLE_STARTS):
        end = CYCLE_STARTS[i + 1][1] if i + 1 < len(CYCLE_STARTS) else rows[-1][0] + dt.timedelta(days=1)
        print(f"   cycle {cyc}: {len(storms([r for r in rows if start <= r[0] < end], 8.66))}")

    print("\n3. Lowest hourly Dst of each month checked, deepest first (Kyoto; version in brackets)")
    lows = []
    for ym, ver in versions.items():
        y, m = int(ym[:4]), int(ym[4:])
        vals = [(v, d) for d, vs in dst.items() if (d.year, d.month) == (y, m) and d <= until for v in vs]
        if vals:
            v, d = min(vals)
            lows.append((v, d, ver))
    for v, d, ver in sorted(lows):
        print(f"   {v:5d} nT on {d} ({ver.lower()})")

    print("\n4. Proton events since December 2019, by peak flux above 10 MeV (NOAA's list)")
    raw = open(fetch(folder, "spe_list.html", [SPE_URL]), encoding="utf-8", errors="replace").read()
    events = []
    for tr in re.findall(r"(?is)<tr[^>]*>(.*?)</tr>", raw):
        cells = [html.unescape(re.sub(r"\s+", " ", re.sub(r"<[^>]+>", " ", c))).strip()
                 for c in re.findall(r"(?is)<t[dh][^>]*>(.*?)</t[dh]>", tr)]
        if len(cells) >= 6 and re.match(r"^\d{4} \d\d/\d\d", cells[0]):
            began = dt.date(int(cells[0][:4]), int(cells[0][5:7]), int(cells[0][8:10]))
            m = re.match(r"[\d,\.]+", cells[2])
            if m and dt.date(2019, 12, 1) <= began <= until:
                events.append((float(m.group(0).replace(",", "")), cells))
    for pfu, c in sorted(events, key=lambda e: -e[0])[:8]:
        level = next(s for t, s in ((1e5, "S5"), (1e4, "S4"), (1e3, "S3"), (1e2, "S2"), (10, "S1"), (0, "-")) if pfu >= t)
        print(f"   {c[2]:>8} pfu  {level}  began {c[0]}, peak {c[1]}, region {c[3]}, flare {c[5]}")
    print(f"   {len(events)} events in all; {sum(1 for p, _ in events if p >= 1000)} at S3 or above")


if __name__ == "__main__":
    main()