#!/usr/bin/env python3
"""Do extreme geomagnetic storms recur every ~89 years? A test of the S0 News arithmetic of 27 October 2026
("1770 to 1859 is 89 years. It's the Glyceberg cycle ... about 88.8 years ... 1948 ... 2037 2038"), in the manner
of Fact-Checked Science's /excursion-cycles/ test, for https://factchecked.science/qianlong-1770-storm/.
Standard library only; fixed seeds, so every run gives the same numbers; it reads no file and needs no network.

THE RULES, fixed before the first run
-------------------------------------
Claim under test: "killer events" recur on an 88.8-year (or 89-year) cycle counted from the 1770 storm, so the
marks fall at 1770.7, 1859.5, 1948.3 and 2037.1 (88.8 years) or 1770.7, 1859.7, 1948.7, 2037.7 (89 years).

Catalogue A, "extreme storms": every storm from 1770 to 2025 that a published reconstruction or the Dst record puts
at a minimum Dst of about -500 nT or below. This is the threshold of Cliver and Dietrich (2013, JSWSC 3 A31,
note 2: "Events with Dst < -500 nT occur approximately every 50 years (1859, 1909(?), 1921, and 1989)"),
updated with the reconstructions published since. Dates are the storm's main day, as decimal years
(the middle of the day).
    1770-09-17  -1319 nT from the aurora (Yuan et al. 2026); "comparable to, or slightly larger than" 1859
                (Kataoka and Iwahashi 2017, Space Weather 15, 1314, abstract)
    1859-09-02  -949 +/- 31 nT (Hayakawa et al. 2022, as cited in Hayakawa et al. 2023, ApJ 959:23); ~-900 nT
                (Cliver and Dietrich 2013)
    1872-02-04  <= -834 nT at Colaba (Hayakawa et al. 2023, abstract; about -1402 nT at Tiflis, one station, §6)
    1903-10-31  ~-531 nT (Hayakawa et al. 2020, ApJL 897 L10, abstract)
    1909-09-25  -595 nT (Love, Hayakawa and Cliver 2019, Space Weather 17, 37, abstract)
    1921-05-15  -907 +/- 132 nT (Love, Hayakawa and Cliver 2019, Space Weather 17, 1281)
    1946-03-28  <= -512 nT (Hayakawa et al. 2020, MNRAS 497, 5507, abstract)
    1989-03-14  -589 nT (WDC Kyoto Dst, final)
  Not in A (above -500 nT): 1938 January (Dcx ~ -336 nT), 1941 March (<= -464 nT), 1957-1959 (-427 to -429 nT),
  2003 November (-422 nT), 2024 May (-406 nT). Before 1859 the record is incomplete: storms between 1770 and 1859
  may be missing, which can only add events that the comb does not predict.
Catalogue B, "Carrington class": storms put at about -800 nT or below, or with an auroral oval as low as 1859's:
  1770, 1859, 1872, 1921 (sources as above; Hayakawa et al. 2023 rank 1859, 1872 and 1921 as the three largest
  storms observed).

Test 1, out-of-sample hits. The comb was fitted to 1770 and 1859, so only the other storms can test it. For the
  storms of A (and B) after 1859, count those within +/- w years of a mark (w = 2, 3, 5) and take their mean
  distance to the nearest mark. Chance level: the same number of storms placed uniformly at random over
  1860.0-2025.99, 200,000 times (seed 89). p = share of random catalogues with at least as many hits (or a mean
  distance at least as small).
Test 2, any period near the Gleissberg band. For each period P from 50 to 150 years (step 0.1), the phase
  bunching of the storm dates, R(P) = |mean of exp(2 pi i t / P)| (1 = all at one phase). The statistic is the
  largest R in the band; chance level from 5,000 random catalogues (seed 1770): the same number of dates uniform
  over the catalogue's window, judged the same way (so the search over P is paid for). Run for A (1770-2025, eight
  dates) and for the six storms of A after the two anchors (window 1860-2025, less than two 89-year cycles, so
  its power is very low).
Test 3, pairs as combs. For every pair of storms in A at least 20 years apart, project the next mark (later date
  plus the gap) and count the pairs whose mark falls within +/- 3 years of another storm in A, for marks up to
  2025. This shows how often an arbitrary pair "predicts" a later storm.
Also reported: the gaps between successive storms.
Run:  python3 -I comb_test.py [out.json]   (about 30 seconds; prints the results, and writes them as JSON if a file
      name is given)
"""
import cmath
import json
import math
import random
import sys

A = [("1770-09-17", 1770.71), ("1859-09-02", 1859.67), ("1872-02-04", 1872.09), ("1903-10-31", 1903.83),
     ("1909-09-25", 1909.73), ("1921-05-15", 1921.37), ("1946-03-28", 1946.24), ("1989-03-14", 1989.20)]
B = [x for x in A if x[0][:4] in ("1770", "1859", "1872", "1921")]
END = 2025.99


def marks(period, start=1770.71, end=2040.0):
    m, k = [], 0
    while start + k * period <= end:
        m.append(start + k * period)
        k += 1
    return m


def dist(t, ms):
    return min(abs(t - x) for x in ms)


def test1(cat, period, ws=(2, 3, 5), n=200_000, seed=89):
    ms = marks(period)
    after = [t for _, t in cat if t > 1860.0]
    obs_d = sum(dist(t, ms) for t in after) / len(after)
    obs_h = {w: sum(dist(t, ms) <= w for t in after) for w in ws}
    rng = random.Random(seed)
    ge_h = {w: 0 for w in ws}
    le_d = 0
    exp_h = {w: 0.0 for w in ws}
    for _ in range(n):
        r = [rng.uniform(1860.0, END) for _ in after]
        d = [dist(t, ms) for t in r]
        if sum(d) / len(d) <= obs_d:
            le_d += 1
        for w in ws:
            h = sum(x <= w for x in d)
            exp_h[w] += h / n
            if h >= obs_h[w]:
                ge_h[w] += 1
    return {"period": period, "marks": [round(x, 1) for x in ms], "storms_tested": len(after),
            "distances": [round(dist(t, ms), 1) for t in after], "mean_distance": obs_d,
            "p_mean_distance": le_d / n,
            "hits": {w: {"observed": obs_h[w], "expected": exp_h[w], "p": ge_h[w] / n} for w in ws}}


def bunching(ts, p):
    return abs(sum(cmath.exp(2j * math.pi * t / p) for t in ts)) / len(ts)


def test2(cat, start, n=5_000, seed=1770):
    ts = [t for _, t in cat if t >= start]
    grid = [50 + 0.1 * i for i in range(1001)]
    rs = [bunching(ts, p) for p in grid]
    best = max(range(len(grid)), key=lambda i: rs[i])
    r888 = bunching(ts, 88.8)
    rng = random.Random(seed)
    lo = start if start > 1770 else 1770.0
    ge_best = ge_888 = 0
    for _ in range(n):
        r = [rng.uniform(lo, END) for _ in ts]
        if max(bunching(r, p) for p in grid) >= rs[best]:
            ge_best += 1
        if bunching(r, 88.8) >= r888:
            ge_888 += 1
    return {"dates": len(ts), "window": [lo, END], "random_catalogues": n, "best_period": grid[best], "best_R": rs[best],
            "R_at_88.8": r888, "p_R_at_88.8_alone": ge_888 / n, "p_best_anywhere_50_150": ge_best / n}


def test3(cat, tol=3.0, mingap=20.0):
    rows = []
    for i in range(len(cat)):
        for j in range(i + 1, len(cat)):
            (ni, ti), (nj, tj) = cat[i], cat[j]
            g = tj - ti
            nxt = tj + g
            if g < mingap or nxt > END:
                continue
            others = [(n, t) for n, t in cat if n not in (ni, nj)]
            near = min(others, key=lambda x: abs(x[1] - nxt))
            rows.append({"pair": [ni[:4], nj[:4]], "gap": round(g, 1), "next_mark": round(nxt, 1),
                         "nearest_storm": near[0][:4], "off_by": round(abs(near[1] - nxt), 1),
                         "hit": abs(near[1] - nxt) <= tol})
    return {"pairs": len(rows), "hits": sum(r["hit"] for r in rows), "rows": rows}


def main():
    sys.stdout.reconfigure(encoding="utf-8")
    out = {"catalogue_A": A, "catalogue_B": B}
    gaps = [round(A[i + 1][1] - A[i][1], 1) for i in range(len(A) - 1)]
    gapsB = [round(B[i + 1][1] - B[i][1], 1) for i in range(len(B) - 1)]
    out["gaps_A"], out["gaps_B"] = gaps, gapsB
    print("gaps between storms, A:", gaps, "| after 1859: mean", round(sum(gaps[1:]) / len(gaps[1:]), 1))
    print("gaps between storms, B:", gapsB)
    for period in (88.8, 89.0):
        for name, cat in (("A", A), ("B", B)):
            r = test1(cat, period)
            out[f"test1_{name}_{period}"] = r
            print(f"Test 1, {name}, P = {period}: marks {r['marks']}; storms after 1859 {r['storms_tested']},"
                  f" distances {r['distances']}, mean {r['mean_distance']:.1f} yr (p = {r['p_mean_distance']:.2f})")
            for w, h in r["hits"].items():
                print(f"    within +/-{w} yr: {h['observed']} (chance {h['expected']:.2f}; p = {h['p']:.2f})")
    for start in (1770.0, 1860.0):
        r = test2(A, start)
        out[f"test2_A_from_{int(start)}"] = r
        print(f"Test 2, A from {int(start)} ({r['dates']} dates): best period {r['best_period']:.1f} yr, R = {r['best_R']:.2f}"
              f" (global p = {r['p_best_anywhere_50_150']:.2f}); R at 88.8 yr = {r['R_at_88.8']:.2f} (p = {r['p_R_at_88.8_alone']:.2f})")
    r = test3(A)
    out["test3_A"] = r
    print(f"Test 3: {r['hits']} of {r['pairs']} pairs 20+ years apart project a mark within 3 years of another storm:")
    for row in r["rows"]:
        print(f"    {row['pair'][0]}-{row['pair'][1]} gap {row['gap']:5.1f} -> {row['next_mark']:.1f}: nearest {row['nearest_storm']}"
              f" ({row['off_by']} yr){'  HIT' if row['hit'] else ''}")
    if len(sys.argv) > 1:
        with open(sys.argv[1], "w", encoding="utf-8", newline="\n") as f:
            json.dump(out, f, indent=1)
        print("wrote", sys.argv[1])


if __name__ == "__main__":
    main()