#!/usr/bin/env python3
"""The counts behind section 6 of https://factchecked.science/qianlong-1770-storm/ (Fact-Checked Science, 27 October
2026): days with Kp 9 around the episode's 1948 mark, the ranking of solar cycles by their peaks, and the largest
sunspot groups in the Greenwich record. Standard library only.

It fetches three public files by URL (about 11 MB in all, over verified HTTPS) and reads them in memory as data,
nothing in them is run; or give it your own copies of the three files, in this order:
  python3 -I counts.py [Kp_ap_Ap_SN_F107_since_1932.txt cyclesminmax.html RGO_NOAA1874_2013.zip] [out.json]

Data:
  Kp: GFZ Helmholtz Centre for Geosciences, Geomagnetic Observatory Niemegk. Matzka, J., et al. (2021), Space
      Weather 19, e2020SW002641, https://doi.org/10.1029/2020SW002641; data https://doi.org/10.5880/Kp.0001, CC BY 4.0.
      File https://kp.gfz.de/app/files/Kp_ap_Ap_SN_F107_since_1932.txt. Its sunspot columns (CC BY-NC) are not used.
  Cycles: WDC-SILSO, Royal Observatory of Belgium, Brussels, table of cycle minima and maxima (13-month smoothed
      sunspot number, version 2.0), https://www.sidc.be/SILSO/cyclesminmax, CC BY-NC 4.0, used with credit on a
      non-commercial site.
  Groups: NASA Marshall Space Flight Center's compilation of the Royal Greenwich Observatory (1874-1976) and USAF/NOAA
      (1977-2013) sunspot-group records, https://solarscience.msfc.nasa.gov/greenwch.shtml, file
      https://solarscience.msfc.nasa.gov/greenwch/RGO_NOAA1874_2013.zip, format
      https://solarscience.msfc.nasa.gov/greenwch/format.txt. A US government work.

Rules, as in our verification (fixed before the counts were first run):
  Kp is given in thirds: 9- = 8.667, 9o = 9.000. NOAA's G5 ("Extreme") is Kp 9o; "a Kp 9 day" is a UT day with at
      least one 3-hour value of 9o. Days whose highest value is 9- are counted separately (NOAA puts 9- in G4).
      Only definitive values are used (the file's last column is 2); days with missing values are skipped.
  Solar cycles run from the month of one minimum in SILSO's table to the month before the next. The cycles are
      ranked by the maximum in the same table.
  A sunspot group's size is its largest daily corrected whole-spot area (columns 41-44, millionths of the visible
      hemisphere), the group being keyed by its year file and Greenwich number. The ranking uses the RGO years
      1874-1976; NASA notes that the USAF/NOAA areas after 1976 run low and multiplies them by 1.4, so those are
      listed apart, for comparison.
It prints the numbers the page uses, and writes them as JSON if a file name is given last.
"""
import hashlib
import html
import io
import json
import re
import sys
import urllib.request
import zipfile
from collections import defaultdict

URLS = ["https://kp.gfz.de/app/files/Kp_ap_Ap_SN_F107_since_1932.txt",
        "https://www.sidc.be/SILSO/cyclesminmax",
        "https://solarscience.msfc.nasa.gov/greenwch/RGO_NOAA1874_2013.zip"]
KP9, KP9M = 8.99, 8.66


def get(src):
    if src.startswith("https://"):
        req = urllib.request.Request(src, headers={"User-Agent": "factchecked.science counts.py"})
        with urllib.request.urlopen(req, timeout=300) as r:
            data = r.read()
    else:
        with open(src, "rb") as f:
            data = f.read()
    print(f"read {src}: {len(data):,} bytes, SHA-256 {hashlib.sha256(data).hexdigest()}")
    return data


def cycles_table(raw):
    t = raw.decode("utf-8", errors="replace")
    cells = [html.unescape(re.sub(r"<[^>]+>", "", c)).strip() for c in re.findall(r"(?is)<td[^>]*>(.*?)</td>", t)]
    rows = []
    for i in range(0, len(cells) - 6, 7):
        c = cells[i:i + 7]
        if c[0].isdigit():
            rows.append({"cycle": int(c[0]), "min": (int(c[1]), int(c[2])), "max": (int(c[4]), int(c[5])),
                         "max_SN": float(c[6].replace(",", "."))})
    if [r["cycle"] for r in rows] != list(range(1, len(rows) + 1)) or len(rows) < 25:
        sys.exit("SILSO's table did not parse as expected; stopping")
    for k, r in enumerate(sorted(rows, key=lambda r: -r["max_SN"])):
        r["rank"] = k + 1
    return rows


def cycle_of(y, m, rows):
    c = None
    for r in rows:
        if (y, m) >= r["min"]:
            c = r["cycle"]
    return c


def kp_counts(raw, rows):
    by_year9, by_year9m, year_max = defaultdict(int), defaultdict(int), {}
    by_cyc9, by_cyc9m = defaultdict(int), defaultdict(int)
    days9 = []
    last = None
    for line in raw.decode("ascii", errors="replace").splitlines():
        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 int(p[-1]) < 2 or any(k < 0 for k in kps):
            continue
        last = f"{y:04d}-{m:02d}-{d:02d}"
        mx = max(kps)
        year_max[y] = max(year_max.get(y, 0.0), mx)
        c = cycle_of(y, m, rows)
        if mx >= KP9:
            by_year9[y] += 1
            by_cyc9[c] += 1
            days9.append(last)
        elif mx >= KP9M:
            by_year9m[y] += 1
            by_cyc9m[c] += 1
    return by_year9, by_year9m, year_max, by_cyc9, by_cyc9m, days9, last


def thirds(v):
    """Kp in the usual notation: 8.333 -> 8+, 8.667 -> 9-, 9.0 -> 9o."""
    n = round(v * 3)
    base, rem = divmod(n, 3)
    return f"{base}o" if rem == 0 else f"{base}+" if rem == 1 else f"{base + 1}-"


def groups(raw):
    best = {}
    with zipfile.ZipFile(io.BytesIO(raw)) as z:
        names = sorted(n for n in z.namelist() if n.lower().endswith(".txt"))
        for n in names:
            for line in z.read(n).decode("ascii", errors="replace").splitlines():
                if len(line) < 44:
                    continue
                try:
                    y, m, d = int(line[0:4]), int(line[4:6]), int(line[6:8])
                    grp, area = line[12:20].strip(), int(line[40:44])
                except ValueError:
                    continue
                if grp and area > 0 and ((n, grp) not in best or area > best[(n, grp)][0]):
                    best[(n, grp)] = (area, f"{y:04d}-{m:02d}-{d:02d}", y)
    rgo = sorted(((a, dt) for (a, dt, y) in best.values() if y <= 1976), reverse=True)
    noaa = sorted(((a, dt) for (a, dt, y) in best.values() if y >= 1977), reverse=True)
    return len(names), rgo, noaa


def main():
    sys.stdout.reconfigure(encoding="utf-8")
    args = sys.argv[1:]
    outp = args.pop() if len(args) in (1, 4) else None
    srcs = args if len(args) == 3 else URLS
    kp_raw, cyc_raw, zip_raw = (get(s) for s in srcs)
    out = {}

    rows = cycles_table(cyc_raw)
    rank = {r["cycle"]: r["rank"] for r in rows}
    print(f"\nSolar cycles ranked by their smoothed maxima ({len(rows)} cycles, SILSO):")
    for c, what in ((2, "the 1770 storm"), (10, "the 1859 storm"), (14, "the 1903 and 1909 storms"), (18, "1948")):
        r = rows[c - 1]
        print(f"  cycle {c:2d} ({what}): maximum {r['max_SN']} in {r['max'][0]}-{r['max'][1]:02d}, {rank[c]} of {len(rows)}")
    out["cycle_rank"] = {str(r["cycle"]): [r["rank"], r["max_SN"], list(r["max"])] for r in rows}

    y9, y9m, ymax, c9, c9m, days9, last = kp_counts(kp_raw, rows)
    complete = [c for c in range(17, 25)]                      # complete cycles since the record began in 1932
    mean9 = sum(c9[c] for c in complete) / len(complete)
    both = {c: c9[c] + c9m[c] for c in complete}
    order = sorted(complete, key=lambda c: -both[c])
    print(f"\nKp (GFZ, definitive values to {last}):")
    print(f"  1948: {y9[1948]} days with Kp 9, {y9m[1948]} more at 9-; its highest value was {thirds(ymax[1948])}")
    print(f"  cycle 18: {c9[18]} days with Kp 9 ({', '.join(d for d in days9 if cycle_of(int(d[:4]), int(d[5:7]), rows) == 18)});"
          f" {both[18]} counting the days at 9-")
    print(f"  cycles 17-24, days with Kp 9: " + ", ".join(f"{c}: {c9[c]}" for c in complete) + f"; mean {mean9:.2f}")
    print(f"  cycles 17-24, days at 9- or 9: " + ", ".join(f"{c}: {both[c]}" for c in order)
          + f"; mean {sum(both.values()) / len(both):.2f}")
    out["kp"] = {"definitive_to": last, "1948": {"kp9_days": y9[1948], "kp9minus_only_days": y9m[1948],
                 "highest": ymax[1948]}, "kp9_days_by_cycle": {str(c): c9[c] for c in complete},
                 "kp9_or_9minus_days_by_cycle": {str(c): both[c] for c in complete}, "mean_kp9_days_17_24": mean9,
                 "kp9_days": days9}

    nfiles, rgo, noaa = groups(zip_raw)
    c18 = [r for r in rows if r["cycle"] in (18, 19)]
    lo, hi = c18[0]["min"], c18[1]["min"]
    print(f"\nLargest sunspot groups, RGO 1874-1976 ({nfiles} year files), millionths of the visible hemisphere:")
    for i, (a, dt) in enumerate(rgo[:8]):
        inside = lo <= (int(dt[:4]), int(dt[5:7])) < hi
        print(f"  {i + 1}. {a:5d} on {dt}" + ("  (cycle 18)" if inside else ""))
    print("  largest after 1976 (USAF/NOAA), as reported and x1.4: "
          + ", ".join(f"{a} ({a * 1.4:.0f}) on {dt}" for a, dt in noaa[:3]))
    out["groups"] = {"rgo_top8": rgo[:8], "usaf_noaa_top3": noaa[:3]}
    if outp:
        with open(outp, "w", encoding="utf-8", newline="\n") as f:
            json.dump(out, f, indent=1)
        print("wrote", outp)


if __name__ == "__main__":
    main()