#!/usr/bin/env python3
"""The numbers behind the 89-year arithmetic, cycle 26 and the cyclone on "S0 News, October 27, 2026".
Fact-Checked Science, sections 4, 5 and 7 ("Our analysis"). Standard library only.

Run it
  python3 -I record_numbers.py [folder]
  The folder (default: the script's own) must hold, or will receive (about 11 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.
    cyclesminmax.html                https://www.sidc.be/SILSO/cyclesminmax
        WDC-SILSO, Royal Observatory of Belgium, Brussels: the dates and smoothed sunspot numbers of every cycle's
        minimum and maximum. CC BY-NC 4.0 (non-commercial use, with credit).
    RGO_NOAA1874_2013.zip            https://solarscience.msfc.nasa.gov/greenwch/RGO_NOAA1874_2013.zip
        Royal Greenwich Observatory (1874-1976) and USAF/NOAA (from 1977) sunspot-group records, compiled by
        D. Hathaway at NASA's Marshall Space Flight Center. Read inside the ZIP; nothing is extracted.

What it prints
  1. Kp: the days whose highest three-hour Kp was 9o (NOAA's G5) or 9- (G4), per solar cycle; the highest Kp of
     1948 and its days; the Kp 9 days of cycle 18; the mean for the complete cycles 17-24. Cycles start in the
     month of the smoothed minimum in SILSO's table.
  2. SILSO: the cycles ranked by smoothed maximum; the mean cycle length and the span of 8 consecutive cycles;
     where the storms of 1770 and 1859, and the year 1948, fall in their cycles.
  3. Greenwich: the five largest sunspot groups of 1874-1976 by peak corrected area (millionths of the Sun's
     hemisphere, MSH), and the largest after 1976, also multiplied by 1.4, Hathaway's correction for the change of
     source.
  4. Typed-in values: Yuan et al.'s (2026) recurrence intervals as chances per decade and per 11-year cycle; the
     cyclone's forecast landfall times in UTC, IST and Arizona time; the days left in 2026.
"""
import datetime as dt
import math
import os
import re
import sys
import urllib.request
import zipfile
from collections import defaultdict

FILES = {
    "Kp_ap_Ap_SN_F107_since_1932.txt": "https://kp.gfz.de/app/files/Kp_ap_Ap_SN_F107_since_1932.txt",
    "cyclesminmax.html": "https://www.sidc.be/SILSO/cyclesminmax",
    "RGO_NOAA1874_2013.zip": "https://solarscience.msfc.nasa.gov/greenwch/RGO_NOAA1874_2013.zip",
}


def fetch(folder, name):
    path = os.path.join(folder, name)
    if not os.path.exists(path):
        print("downloading", FILES[name])
        req = urllib.request.Request(FILES[name], headers={"User-Agent": "factchecked.science record_numbers.py"})
        with urllib.request.urlopen(req, timeout=300) as r, open(path, "wb") as f:
            f.write(r.read())
    return path


def frac_year(y, m, d=15):
    start = dt.date(y, 1, 1)
    return y + (dt.date(y, m, d) - start).days / (dt.date(y + 1, 1, 1) - start).days


def silso_cycles(path):
    raw = open(path, encoding="utf-8", errors="replace").read()
    cells = [re.sub(r"<[^>]+>", "", c).strip() for c in re.findall(r"<td[^>]*>(.*?)</td>", raw, flags=re.S | re.I)]
    cyc, i = {}, 0
    while i + 6 < len(cells):
        ch = cells[i:i + 7]
        if re.fullmatch(r"\d{1,2}", ch[0]) and re.fullmatch(r"\d{4}", ch[1]):
            cyc[int(ch[0])] = {"min": (int(ch[1]), int(ch[2])), "max": (int(ch[4]), int(ch[5])),
                               "sn": float(ch[6].replace(",", "."))}
            i += 7
        else:
            i += 1
    return cyc


def main():
    sys.stdout.reconfigure(encoding="utf-8")
    folder = sys.argv[1] if len(sys.argv) > 1 else os.path.dirname(os.path.abspath(__file__))
    cyc = silso_cycles(fetch(folder, "cyclesminmax.html"))
    starts = sorted((v["min"], k) for k, v in cyc.items())

    def cycle_of(y, m):
        c = None
        for (yy, mm), k in starts:
            if (y, m) >= (yy, mm):
                c = k
        return c

    # 1. Kp
    per = defaultdict(lambda: [0, 0])
    k9_18, y1948, ndays, last = [], [], 0, None
    for line in open(fetch(folder, "Kp_ap_Ap_SN_F107_since_1932.txt"), encoding="utf-8", errors="replace"):
        if line.startswith("#") or not line.strip():
            continue
        p = line.split()
        y, m, d = int(p[0]), int(p[1]), int(p[2])
        kps = [float(x) for x in p[7:15]]
        if any(k < 0 for k in kps):
            continue
        ndays += 1
        last = f"{y}-{m:02d}-{d:02d}"
        k = max(kps)
        c = cycle_of(y, m)
        if k >= 8.99:
            per[c][0] += 1
            if c == 18:
                k9_18.append(f"{y}-{m:02d}-{d:02d}")
        elif k > 8.6:
            per[c][1] += 1
        if y == 1948:
            y1948.append((k, f"{y}-{m:02d}-{d:02d}"))
    print(f"1. Kp, {ndays} days to {last}")
    for c in range(17, max(cyc) + 1):
        print(f"   cycle {c}{' (incomplete)' if c == max(cyc) else ''}: {per[c][0]} at Kp 9o, {per[c][1]} more at 9-")
    full = [c for c in range(17, 25)]
    print(f"   mean for cycles 17-24: {sum(per[c][0] for c in full) / len(full):.2f} Kp 9o days per cycle")
    print(f"   cycle 18's Kp 9o days: {', '.join(k9_18)}")
    top = max(k for k, _ in y1948)
    names = {8.333: "8+", 8.667: "9-", 9.0: "9o", 8.0: "8o", 7.667: "8-"}
    print(f"   1948: highest Kp {names.get(round(top, 3), top)} on {', '.join(dd for k, dd in y1948 if k == top)}; "
          f"days at 9o: {sum(1 for k, _ in y1948 if k >= 8.99)}; at 9-: {sum(1 for k, _ in y1948 if 8.6 < k < 8.99)}")

    # 2. SILSO
    ranked = sorted(cyc, key=lambda k: -cyc[k]["sn"])
    print("\n2. Sunspot cycles by smoothed maximum (SILSO):")
    print("   " + "; ".join(f"{i}. cycle {k} {cyc[k]['sn']:.1f} ({cyc[k]['max'][0]}-{cyc[k]['max'][1]:02d})"
                          for i, k in enumerate(ranked, 1)))
    mins = {k: frac_year(*cyc[k]["min"]) for k in cyc}
    mean = (mins[max(cyc)] - mins[min(cyc)]) / (max(cyc) - min(cyc))
    spans = [mins[k + 8] - mins[k] for k in cyc if k + 8 in cyc]
    print(f"   mean cycle length {mean:.2f} years; 8 cycles = {8 * mean:.1f} years; "
          f"8 consecutive cycles span {min(spans):.1f}-{max(spans):.1f} years")
    for t, label in ((frac_year(1770, 9, 17), "the Qianlong storm, 17 Sep 1770"),
                     (frac_year(1859, 9, 2), "the Carrington storm, 2 Sep 1859"), (1948.5, "mid-1948")):
        k = max(c for c in cyc if mins[c] <= t)
        mx = frac_year(*cyc[k]["max"])
        print(f"   {label}: cycle {k}, {t - mx:+.2f} years from its maximum ({cyc[k]['max'][0]}-{cyc[k]['max'][1]:02d}, "
              f"{cyc[k]['sn']:.1f}, rank {ranked.index(k) + 1} of {len(cyc)})")

    # 3. Greenwich
    best = {}
    with zipfile.ZipFile(fetch(folder, "RGO_NOAA1874_2013.zip")) as z:
        for n in sorted(x for x in z.namelist() if x.lower().endswith(".txt")):
            for raw in z.read(n).decode("latin-1").splitlines():
                if len(raw) < 44 or not raw[:4].strip().isdigit():
                    continue
                try:
                    y, mo, d, g, a = int(raw[0:4]), int(raw[4:6]), int(raw[6:8]), raw[12:20].strip(), int(raw[40:44])
                except ValueError:
                    continue
                key = (y, g) if y <= 1976 else ("NOAA", g, y // 10)
                if a > best.get(key, (0, ""))[0]:
                    best[key] = (a, f"{y}-{mo:02d}-{d:02d}")
    rgo = sorted(((a, w) for k, (a, w) in best.items() if k[0] != "NOAA"), reverse=True)[:5]
    noaa = max((a, w) for k, (a, w) in best.items() if k[0] == "NOAA")
    print("\n3. Largest sunspot groups, RGO 1874-1976 (peak corrected area, MSH): "
          + "; ".join(f"{a:,} on {w}" for a, w in rgo))
    print(f"   largest after 1976 (USAF/NOAA): {noaa[0]:,} MSH on {noaa[1]}, x1.4 = {noaa[0] * 1.4:,.0f}")

    # 4. Typed-in values
    print("\n4. Yuan et al. (2026), recurrence intervals as chances (Poisson):")
    for rec, label in ((14000, "1770-size, central 14,000 years"), (2000, "1770-size, short end about 2,000 years"),
                       (1400, "1859-size, central 1,400 years")):
        print(f"   {label}: {100 * (1 - math.exp(-10 / rec)):.2f}% a decade, {100 * (1 - math.exp(-11 / rec)):.2f}% an 11-year cycle")
    utc, ist = dt.timezone.utc, dt.timezone(dt.timedelta(hours=5, minutes=30))
    az = dt.timezone(dt.timedelta(hours=-7))
    up = dt.datetime(2026, 10, 27, 12, 22, 43, tzinfo=utc)
    print(f"   the upload: {up:%d %b %H:%M} UTC = {up.astimezone(ist):%d %b %H:%M} IST = {up.astimezone(az):%d %b %H:%M} Arizona")
    for label, t in (("IMD, 21:30 IST update: from", dt.datetime(2026, 10, 28, 2, 0, tzinfo=ist)),
                     ("                        to", dt.datetime(2026, 10, 28, 4, 0, tzinfo=ist)),
                     ("RSMC advisory 12: by", dt.datetime(2026, 10, 28, 6, 0, tzinfo=utc))):
        print(f"   {label} {t.astimezone(utc):%d %b %H:%M} UTC, {t.astimezone(ist):%d %b %H:%M} IST, "
              f"{t.astimezone(az):%d %b %H:%M} Arizona ({(t - up).total_seconds() / 3600:.1f} h after the upload)")
    d0 = dt.date(2026, 10, 27)
    print(f"   days after Oct 27 to the end of 2026: {(dt.date(2026, 12, 31) - d0).days}")


if __name__ == "__main__":
    main()