#!/usr/bin/env python3
"""Where is Earth's magnetic field weakening fastest, and how far is it from the laboratory?

A pure-Python (standard library only) spherical-harmonic evaluator for IGRF-14 and WMM2025,
written for the Fact-Checked Science check of S0 News, 20 October 2026 ("Our analysis").

    python3 geomag_sh.py                      inputs in ./geomag_inputs (downloaded from NOAA if missing), results in ./geomag_out
    python3 geomag_sh.py <inputs> <out>       your own folders

What it does
  1. Validates the evaluator against NOAA's twelve published WMM2025 test values (X, Y, Z, F and
     their yearly rates), so the IGRF-14 numbers can be trusted.
  2. Evaluates the IGRF-14 total intensity F at 2025.0 and its forecast rate dF/dt (2025-2030 secular
     variation) on a 2-degree global grid and at named cities, to find where the field is weakening fastest.
  3. Cross-checks the 2020-2025 mean rate from the DGRF-2020 and IGRF-2025 main-field models, and runs
     WMM2025 as a second model.
  4. Computes the dipole moment from IGRF-14 (1900-2025, forecast 2030), its present rate, and how long
     a straight-line continuation would take to reach laboratory "hypomagnetic" levels (5 uT, 1 uT) and
     the Laschamp low (5% of today).
  5. Computes the north dip pole's drift speeds from NOAA's yearly positions (NP.xy).

Inputs (read-only; downloaded to <inputs> if absent; the IGRF file's SHA-256 is checked, the others noted)
  igrf14/igrf14coeffs.txt          IGRF-14 Gauss coefficients   https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt
  wmm2025/WMM2025COF.zip           WMM2025 package              https://www.ncei.noaa.gov/sites/default/files/2025-01/WMM2025COF.zip
      -> wmm2025/WMM2025.COF and wmm2025/WMM2025_TEST_VALUES.txt (extracted from the zip)
  ncei-poles/NP.xy                 north dip pole positions     https://www.ngdc.noaa.gov/geomag/data/poles/NP.xy

Outputs (<out>)
  validation_wmm2025.txt, dFdt_summary.txt, dFdt_grid_igrf14_2025.csv, dipole_igrf14.txt, pole_speeds.txt,
  and s0-news-2026-10-20-numbers.json with every number used on the page.

Conventions
  Schmidt quasi-normalised associated Legendre functions; geocentric spherical coordinates; reference
  radius a = 6371.2 km; WGS84 ellipsoid for geodetic positions; X north, Y east, Z down;
  F = sqrt(X^2 + Y^2 + Z^2); dF/dt = (X Xdot + Y Ydot + Z Zdot) / F.
"""
import csv
import hashlib
import json
import math
import os
import sys
import urllib.request
import zipfile

A_REF = 6371.2            # km, IGRF/WMM reference radius
WGS84_A = 6378.137        # km
WGS84_F = 1.0 / 298.257223563
WGS84_E2 = WGS84_F * (2.0 - WGS84_F)
MU0_OVER_4PI = 1.0e-7     # T m / A

# relative path -> (URL, SHA-256 of the copy we used on 2026-10-21, whether a mismatch stops the run)
INPUTS = {
    "igrf14/igrf14coeffs.txt": ("https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt",
                                "8f8d88403028fc4ee92c4f38d97b46e0a87e2cfc496045b43c9e26c1d6b0903c", True),
    "wmm2025/WMM2025COF.zip": ("https://www.ncei.noaa.gov/sites/default/files/2025-01/WMM2025COF.zip",
                               "999080a12c95122c2031aeac8be77229d97584530e87e0f1a6df7d946ecf90c5", False),
    "ncei-poles/NP.xy": ("https://www.ngdc.noaa.gov/geomag/data/poles/NP.xy",
                         "7c329dcebbc2dc40dce7a8b109a27bcd5b80a439b296d304d16efc12a0f913e7", False),
}


# ----------------------------------------------------------------------------- inputs
def sha256(path):
    h = hashlib.sha256()
    with open(path, "rb") as f:
        for chunk in iter(lambda: f.read(1 << 20), b""):
            h.update(chunk)
    return h.hexdigest()


def fetch_inputs(dl):
    """Download any missing input into <dl>, extract the WMM package, check the hashes."""
    for rel, (url, want, strict) in INPUTS.items():
        path = os.path.join(dl, *rel.split("/"))
        if not os.path.exists(path):
            os.makedirs(os.path.dirname(path), exist_ok=True)
            print("downloading", url)
            req = urllib.request.Request(url, headers={"User-Agent": "Mozilla/5.0 (geomag_sh.py; factchecked.science)"})
            with urllib.request.urlopen(req, timeout=60) as r, open(path, "wb") as f:
                f.write(r.read())
        if sha256(path) != want:
            msg = "%s differs from the copy we used (SHA-256 %s)" % (path, want)
            if strict:
                raise SystemExit(msg + "; delete it to download again, or check NOAA's IGRF page")
            print("note:", msg, "(NOAA updates this file; the numbers may differ slightly)")
    zpath = os.path.join(dl, "wmm2025", "WMM2025COF.zip")
    for member, name in (("WMM2025COF/WMM.COF", "WMM2025.COF"), ("WMM2025COF/WMM2025_TEST_VALUES.txt", "WMM2025_TEST_VALUES.txt")):
        target = os.path.join(dl, "wmm2025", name)
        if not os.path.exists(target):
            with zipfile.ZipFile(zpath) as z:
                data = z.read(member)
            with open(target, "wb") as f:
                f.write(data)


# ----------------------------------------------------------------------------- coefficients
class Model:
    """A set of Gauss coefficients g[n][m], h[n][m] (nT) with secular variation (nT/yr)."""

    def __init__(self, nmax, nmax_sv, epoch):
        self.nmax = nmax
        self.nmax_sv = nmax_sv
        self.epoch = epoch
        z = lambda: [[0.0] * (nmax + 1) for _ in range(nmax + 1)]
        self.g, self.h, self.gd, self.hd = z(), z(), z(), z()

    def at(self, t):
        """Coefficients linearly extrapolated to decimal year t (as WMM/IGRF prescribe)."""
        m = Model(self.nmax, self.nmax_sv, t)
        dt = t - self.epoch
        for n in range(1, self.nmax + 1):
            for k in range(0, n + 1):
                m.g[n][k] = self.g[n][k] + dt * self.gd[n][k]
                m.h[n][k] = self.h[n][k] + dt * self.hd[n][k]
                m.gd[n][k] = self.gd[n][k]
                m.hd[n][k] = self.hd[n][k]
        return m


def read_igrf(path):
    """Return (epochs, {epoch: Model}) from the IGRF-14 coefficient file.

    The last column is the 2025-2030 secular variation (nT/yr), degrees 1-8; main-field
    columns run 1900.0 ... 2025.0 (degree 10 to 1995, degree 13 from 2000)."""
    rows = []
    epochs = None
    with open(path, encoding="utf-8") as f:
        for line in f:
            if line.startswith("#"):
                continue
            parts = line.split()
            if not parts or parts[0] == "c/s":
                continue
            if parts[0] == "g/h":
                epochs = [float(x) for x in parts[3:-1]]
                continue
            rows.append(parts)
    models = {}
    for ep in epochs:
        models[ep] = Model(13 if ep >= 2000 else 10, 8, ep)
    for parts in rows:
        gh, n, m = parts[0], int(parts[1]), int(parts[2])
        vals = [float(x) for x in parts[3:]]
        sv = vals[-1]
        for i, ep in enumerate(epochs):
            mod = models[ep]
            if n > mod.nmax:
                continue
            if gh == "g":
                mod.g[n][m] = vals[i]
            else:
                mod.h[n][m] = vals[i]
            if ep == epochs[-1] and n <= 8:
                if gh == "g":
                    mod.gd[n][m] = sv
                else:
                    mod.hd[n][m] = sv
    return epochs, models


def read_wmm(path):
    """Return a Model from a WMM .COF file (n m g h gdot hdot; header line first)."""
    mod = None
    with open(path, encoding="utf-8") as f:
        for line in f:
            parts = line.split()
            if not parts:
                continue
            if mod is None:
                mod = Model(12, 12, float(parts[0]))
                continue
            if parts[0].startswith("9999"):
                break
            n, m = int(parts[0]), int(parts[1])
            mod.g[n][m], mod.h[n][m] = float(parts[2]), float(parts[3])
            mod.gd[n][m], mod.hd[n][m] = float(parts[4]), float(parts[5])
    return mod


# ----------------------------------------------------------------------------- evaluation
def schmidt_legendre(nmax, theta):
    """Schmidt quasi-normalised P[n][m](cos theta) and dP/dtheta, theta = colatitude (rad)."""
    ct, st = math.cos(theta), math.sin(theta)
    P = [[0.0] * (nmax + 1) for _ in range(nmax + 1)]
    dP = [[0.0] * (nmax + 1) for _ in range(nmax + 1)]
    P[0][0] = 1.0
    for n in range(1, nmax + 1):
        for m in range(0, n + 1):
            if n == m:
                if n == 1:
                    P[1][1], dP[1][1] = st, ct
                else:
                    k = math.sqrt((2.0 * n - 1.0) / (2.0 * n))
                    P[n][n] = k * st * P[n - 1][n - 1]
                    dP[n][n] = k * (st * dP[n - 1][n - 1] + ct * P[n - 1][n - 1])
            else:
                d = math.sqrt(float(n * n - m * m))
                k1 = (2.0 * n - 1.0) / d
                k2 = math.sqrt(float((n - 1) * (n - 1) - m * m)) / d
                p2 = P[n - 2][m] if n - 2 >= m else 0.0
                dp2 = dP[n - 2][m] if n - 2 >= m else 0.0
                P[n][m] = k1 * ct * P[n - 1][m] - k2 * p2
                dP[n][m] = k1 * (ct * dP[n - 1][m] - st * P[n - 1][m]) - k2 * dp2
    return P, dP


def geodetic_to_geocentric(lat_deg, lon_deg, h_km):
    """WGS84 geodetic (lat, lon, height) -> (r km, colatitude rad, longitude rad, geocentric lat rad)."""
    phi = math.radians(lat_deg)
    sphi, cphi = math.sin(phi), math.cos(phi)
    rc = WGS84_A / math.sqrt(1.0 - WGS84_E2 * sphi * sphi)
    p = (rc + h_km) * cphi
    z = (rc * (1.0 - WGS84_E2) + h_km) * sphi
    r = math.hypot(p, z)
    phic = math.atan2(z, p)
    return r, math.pi / 2.0 - phic, math.radians(lon_deg), phic


def field_components(g, h, nmax, r, theta, lon):
    """Geocentric X' (north), Y (east), Z' (down) in nT from coefficient arrays."""
    P, dP = schmidt_legendre(nmax, theta)
    st = math.sin(theta)
    X = Y = Z = 0.0
    cm = [math.cos(m * lon) for m in range(nmax + 1)]
    sm = [math.sin(m * lon) for m in range(nmax + 1)]
    for n in range(1, nmax + 1):
        ar = (A_REF / r) ** (n + 2)
        for m in range(0, n + 1):
            gc = g[n][m] * cm[m] + h[n][m] * sm[m]
            X += ar * gc * dP[n][m]
            Y += ar * m * (g[n][m] * sm[m] - h[n][m] * cm[m]) * P[n][m]
            Z -= ar * (n + 1) * gc * P[n][m]
    Y = Y / st if abs(st) > 1e-10 else 0.0
    return X, Y, Z


def evaluate(model, lat, lon, h_km=0.0, t=None):
    """Field and its rate at a geodetic point: X, Y, Z, H, F, I, D and their rates (Xd ... Fd)."""
    mod = model if t is None else model.at(t)
    r, theta, lam, phic = geodetic_to_geocentric(lat, lon, h_km)
    Xp, Y, Zp = field_components(mod.g, mod.h, mod.nmax, r, theta, lam)
    Xdp, Yd, Zdp = field_components(mod.gd, mod.hd, mod.nmax_sv, r, theta, lam)
    psi = phic - math.radians(lat)      # geocentric minus geodetic latitude
    cps, sps = math.cos(psi), math.sin(psi)
    X = Xp * cps - Zp * sps
    Z = Xp * sps + Zp * cps
    Xd = Xdp * cps - Zdp * sps
    Zd = Xdp * sps + Zdp * cps
    H = math.hypot(X, Y)
    F = math.sqrt(X * X + Y * Y + Z * Z)
    return dict(X=X, Y=Y, Z=Z, H=H, F=F, I=math.degrees(math.atan2(Z, H)), D=math.degrees(math.atan2(Y, X)),
                Xd=Xd, Yd=Yd, Zd=Zd, Hd=(X * Xd + Y * Yd) / H if H > 0 else 0.0,
                Fd=(X * Xd + Y * Yd + Z * Zd) / F if F > 0 else 0.0)


# ----------------------------------------------------------------------------- tasks
def validate_wmm(wmm, test_path, out, J):
    """Compare against NOAA's WMM2025 test values; report the largest differences."""
    worst = {}
    lines = ["Validation of the evaluator against NOAA's WMM2025 test values",
             "(difference = ours minus NOAA; NOAA rounds to 0.1 nT and 0.1 nT/yr)", "",
             "date   h(km)   lat     lon  |  dX     dY     dZ     dF   | dXdot dYdot dZdot dFdot"]
    n_tests = 0
    with open(test_path, encoding="utf-8") as f:
        for line in f:
            if line.startswith("#") or not line.strip():
                continue
            p = line.split()
            t, hk, lat, lon = (float(x) for x in p[0:4])
            ref = dict(X=float(p[4]), Y=float(p[5]), Z=float(p[6]), F=float(p[8]),
                       Xd=float(p[12]), Yd=float(p[13]), Zd=float(p[14]), Fd=float(p[16]))
            v = evaluate(wmm, lat, lon, hk, t)
            diffs = {k: v[k] - ref[k] for k in ref}
            for k, d in diffs.items():
                worst[k] = max(worst.get(k, 0.0), abs(d))
            lines.append("%6.1f %6.1f %6.1f %7.1f | %5.2f %5.2f %5.2f %5.2f | %5.2f %5.2f %5.2f %5.2f" % (
                t, hk, lat, lon, diffs["X"], diffs["Y"], diffs["Z"], diffs["F"], diffs["Xd"], diffs["Yd"], diffs["Zd"], diffs["Fd"]))
            n_tests += 1
    ok = all(worst[k] <= 0.06 for k in worst)
    lines += ["", "Largest absolute difference over %d test points: " % n_tests +
              ", ".join("%s %.3f" % (k, worst[k]) for k in ("X", "Y", "Z", "F", "Xd", "Yd", "Zd", "Fd")),
              "Result: %s (all within rounding of NOAA's values)" % ("PASS" if ok else "FAIL")]
    with open(os.path.join(out, "validation_wmm2025.txt"), "w", encoding="utf-8") as f:
        f.write("\n".join(lines) + "\n")
    print("\n".join(lines))
    J["validation"] = {"test_points": n_tests, "largest_abs_difference_nT_or_nT_per_yr": {k: round(v, 3) for k, v in worst.items()}, "pass": ok}
    return ok


CITIES = [
    # name, geodetic latitude, longitude (east positive)
    ("San Antonio, USA", 29.4241, -98.4936),
    ("Shanghai, China", 31.2304, 121.4737),
    ("Sao Paulo, Brazil", -23.5505, -46.6333),
    ("Johannesburg, South Africa", -26.2041, 28.0473),
    ("Ottawa, Canada", 45.4215, -75.6972),
    ("Washington DC, USA", 38.9072, -77.0369),
    ("Los Angeles, USA", 34.0522, -118.2437),
    ("Miami, USA", 25.7617, -80.1918),
    ("Anchorage, USA", 61.2181, -149.9003),
    ("Honolulu, USA", 21.3099, -157.8581),
    ("Mexico City, Mexico", 19.4326, -99.1332),
    ("Bogota, Colombia", 4.7110, -74.0721),
    ("Rio de Janeiro, Brazil", -22.9068, -43.1729),
    ("Buenos Aires, Argentina", -34.6037, -58.3816),
    ("Ascension Island", -7.95, -14.36),
    ("Cape Town, South Africa", -33.9249, 18.4241),
    ("London, UK", 51.5074, -0.1278),
    ("Moscow, Russia", 55.7558, 37.6173),
    ("Delhi, India", 28.6139, 77.2090),
    ("Tokyo, Japan", 35.6762, 139.6503),
    ("Sydney, Australia", -33.8688, 151.2093),
]

BOXES = {
    "North America (lat 25..60, lon -130..-60)": (25, 60, -130, -60),
    "Contiguous USA approx. (lat 25..49, lon -125..-67)": (25, 49, -125, -67),
    "Canada approx. (lat 49..70, lon -140..-55)": (49, 70, -140, -55),
    "Central America and Caribbean (lat 7..25, lon -105..-60)": (7, 25, -105, -60),
    "South America (lat -55..12, lon -82..-34)": (-55, 12, -82, -34),
    "South Atlantic (lat -50..0, lon -50..20)": (-50, 0, -50, 20),
    "Southern Africa (lat -35..-15, lon 10..40)": (-35, -15, 10, 40),
    "East Asia (lat 20..50, lon 100..145)": (20, 50, 100, 145),
    "Europe (lat 35..70, lon -10..40)": (35, 70, -10, 40),
}


def grid(model, step=2.0):
    rows = []
    lat = -89.0
    while lat <= 89.0 + 1e-9:
        lon = -180.0
        while lon < 180.0 - 1e-9:
            v = evaluate(model, lat, lon)
            rows.append((lat, lon, v["F"], v["Fd"]))
            lon += step
        lat += step
    return rows


def grid_and_cities(igrf2025, dgrf2020, wmm, out, J):
    """dF/dt from IGRF-14 at 2025.0 on a 2-degree grid and at cities; WMM2025 and 2020-2025 cross-checks."""
    rows = grid(igrf2025)
    fmin = min(rows, key=lambda r: r[2])
    dmin = min(rows, key=lambda r: r[3])
    dmax = max(rows, key=lambda r: r[3])
    area_total = sum(math.cos(math.radians(r[0])) for r in rows)
    share = {thr: 100.0 * sum(math.cos(math.radians(r[0])) for r in rows if r[3] < thr) / area_total for thr in (-50, -75, -100)}
    with open(os.path.join(out, "dFdt_grid_igrf14_2025.csv"), "w", newline="", encoding="utf-8") as f:
        wr = csv.writer(f)
        wr.writerow(["lat_deg", "lon_deg", "F_nT_2025.0", "dFdt_nT_per_yr_2025-2030"])
        for r in rows:
            wr.writerow(["%.1f" % r[0], "%.1f" % r[1], "%.1f" % r[2], "%.2f" % r[3]])
    lines = ["IGRF-14 total intensity F at 2025.0 and its forecast rate dF/dt (2025-2030 secular variation),",
             "evaluated on the WGS84 surface (height 0) on a 2-degree grid (lat -89..89, lon -180..178).", "",
             "Weakest field on the grid: F = %.0f nT at lat %.0f, lon %.0f (South Atlantic Anomaly)" % (fmin[2], fmin[0], fmin[1]),
             "Fastest decline on the grid: dF/dt = %.1f nT/yr at lat %.0f, lon %.0f" % (dmin[3], dmin[0], dmin[1]),
             "Fastest rise on the grid:    dF/dt = %+.1f nT/yr at lat %.0f, lon %.0f" % (dmax[3], dmax[0], dmax[1])]
    for thr in sorted(share):
        lines.append("Share of Earth's surface with dF/dt below %d nT/yr: %.1f%%" % (thr, share[thr]))
    lines += ["", "Regional dF/dt (nT/yr): area-weighted mean, and the fastest decline in the box"]
    regions = {}
    for name, (la0, la1, lo0, lo1) in BOXES.items():
        sel = [r for r in rows if la0 <= r[0] <= la1 and lo0 <= r[1] <= lo1]
        sw = sum(math.cos(math.radians(r[0])) for r in sel)
        mean = sum(math.cos(math.radians(r[0])) * r[3] for r in sel) / sw
        mn = min(sel, key=lambda r: r[3])
        regions[name] = {"mean_dFdt": round(mean, 1), "min_dFdt": round(mn[3], 1), "min_at": [mn[0], mn[1]]}
        lines.append("  %-62s mean %7.1f   min %7.1f at (%.0f, %.0f)" % (name, mean, mn[3], mn[0], mn[1]))
    lines += ["", "Cities (geodetic position, height 0): F in 2025.0, forecast dF/dt 2025-2030, per cent per year,",
              "the 2020-2025 mean rate from the DGRF-2020 and IGRF-2025 main-field models (degree 13), and WMM2025's rate",
              "%-28s %8s %9s %9s %10s %9s %6s" % ("city", "F 2025", "dF/dt", "%/yr", "dF/dt20-25", "WMM2025", "rank")]
    cities = []
    for name, la, lo in CITIES:
        v = evaluate(igrf2025, la, lo)
        v20 = evaluate(dgrf2020, la, lo)
        w = evaluate(wmm, la, lo)
        cities.append({"city": name, "lat": la, "lon": lo, "F_2025_nT": round(v["F"]), "dFdt_2025_2030_nT_per_yr": round(v["Fd"], 1),
                       "percent_per_yr": round(100.0 * v["Fd"] / v["F"], 3), "dFdt_2020_2025_nT_per_yr": round((v["F"] - v20["F"]) / 5.0, 1),
                       "wmm2025_dFdt_nT_per_yr": round(w["Fd"], 1), "wmm2025_F_nT": round(w["F"])})
    order = sorted(cities, key=lambda c: c["dFdt_2025_2030_nT_per_yr"])
    for i, c in enumerate(order):
        c["rank_fastest_decline"] = i + 1
    for c in cities:
        lines.append("%-28s %8d %9.1f %9.3f %10.1f %9.1f %6d" % (c["city"], c["F_2025_nT"], c["dFdt_2025_2030_nT_per_yr"], c["percent_per_yr"],
                                                                 c["dFdt_2020_2025_nT_per_yr"], c["wmm2025_dFdt_nT_per_yr"], c["rank_fastest_decline"]))
    lines += ["(rank 1 = fastest decline among the listed cities)", "",
              "Share of Earth's surface declining faster than each of S0's named places (2025-2030 forecast):"]
    for c in cities[:5]:
        below = 100.0 * sum(math.cos(math.radians(r[0])) for r in rows if r[3] < c["dFdt_2025_2030_nT_per_yr"]) / area_total
        c["share_of_surface_declining_faster_percent"] = round(below, 1)
        lines.append("  %-28s %.1f%% of the surface has a more negative dF/dt than %.1f nT/yr" % (c["city"], below, c["dFdt_2025_2030_nT_per_yr"]))
    # relative rates: per cent of the local field per year
    rmin = min(rows, key=lambda r: r[3] / r[2])
    la0, la1, lo0, lo1 = BOXES["North America (lat 25..60, lon -130..-60)"]
    na_rmin = min((r for r in rows if la0 <= r[0] <= la1 and lo0 <= r[1] <= lo1), key=lambda r: r[3] / r[2])
    lines += ["", "In proportion to the local field: fastest relative decline %.2f%%/yr at lat %.0f, lon %.0f (F = %.0f nT, dF/dt = %.1f nT/yr)" % (
        100.0 * rmin[3] / rmin[2], rmin[0], rmin[1], rmin[2], rmin[3]),
        "  North America's fastest relative decline %.2f%%/yr at lat %.0f, lon %.0f (F = %.0f nT, dF/dt = %.1f nT/yr)" % (
        100.0 * na_rmin[3] / na_rmin[2], na_rmin[0], na_rmin[1], na_rmin[2], na_rmin[3])]
    wrows = grid(wmm)
    wmin = min(wrows, key=lambda r: r[3])
    v = evaluate(igrf2025, 85.76, 139.30)
    lines += ["", "WMM2025 as a second model: fastest decline %.1f nT/yr at lat %.0f, lon %.0f" % (wmin[3], wmin[0], wmin[1]),
              "Cross-check: at WMM2025's north dip pole (85.76 N, 139.30 E) IGRF-14 gives H = %.0f nT, F = %.0f nT (H should be near 0)" % (v["H"], v["F"])]
    text = "\n".join(lines)
    with open(os.path.join(out, "dFdt_summary.txt"), "w", encoding="utf-8") as f:
        f.write(text + "\n")
    print(text)
    J["field_2025"] = {
        "grid_step_deg": 2.0, "weakest_F_nT": round(fmin[2]), "weakest_at": [fmin[0], fmin[1]],
        "fastest_decline_nT_per_yr": round(dmin[3], 1), "fastest_decline_at": [dmin[0], dmin[1]],
        "fastest_rise_nT_per_yr": round(dmax[3], 1), "fastest_rise_at": [dmax[0], dmax[1]],
        "share_of_surface_below_threshold_percent": {str(k): round(v, 1) for k, v in share.items()},
        "wmm2025_fastest_decline_nT_per_yr": round(wmin[3], 1), "wmm2025_fastest_decline_at": [wmin[0], wmin[1]],
        "regions": regions, "cities": cities,
        "fastest_relative_decline": {"percent_per_yr": round(100.0 * rmin[3] / rmin[2], 2), "at": [rmin[0], rmin[1]], "F_nT": round(rmin[2]), "dFdt_nT_per_yr": round(rmin[3], 1)},
        "north_america_fastest_relative_decline": {"percent_per_yr": round(100.0 * na_rmin[3] / na_rmin[2], 2), "at": [na_rmin[0], na_rmin[1]], "F_nT": round(na_rmin[2]), "dFdt_nT_per_yr": round(na_rmin[3], 1)},
        "H_at_wmm2025_north_dip_pole_nT": round(v["H"]),
    }


def city_history(models, J):
    """F at four cities through the century, from the IGRF-14 main-field models."""
    pts = {n: (la, lo) for n, la, lo in CITIES}
    names = ["San Antonio, USA", "Shanghai, China", "Ottawa, Canada", "Sao Paulo, Brazil"]
    hist = {}
    for ep in (1900.0, 1925.0, 1950.0, 1975.0, 2000.0, 2025.0):
        hist[str(int(ep))] = {n: round(evaluate(models[ep], *pts[n])["F"]) for n in names}
    J["F_history_nT"] = hist
    sa = hist["1900"]["San Antonio, USA"], hist["2025"]["San Antonio, USA"]
    J["san_antonio_1900_2025"] = {"F_1900_nT": sa[0], "F_2025_nT": sa[1], "change_percent": round(100.0 * (sa[1] - sa[0]) / sa[0], 1),
                                  "mean_rate_nT_per_yr": round((sa[1] - sa[0]) / 125.0, 1)}
    v = evaluate(models[2025.0], *pts["San Antonio, USA"])
    J["san_antonio_local_straight_line_years"] = {"to_5_uT": round((v["F"] - 5000.0) / -v["Fd"]), "to_1_uT": round((v["F"] - 1000.0) / -v["Fd"])}
    J["moscow_field_over_350_nT"] = round(evaluate(models[2025.0], 55.7558, 37.6173)["F"] / 350.0)


def dipole(epochs, models, out, J):
    """Dipole moment from the degree-1 coefficients, its rate and straight-line projections."""
    a_m = A_REF * 1000.0
    lines = ["Dipole moment from IGRF-14 degree-1 coefficients: m = (4 pi / mu0) a^3 sqrt(g10^2 + g11^2 + h11^2)",
             "a = 6371.2 km; 1 nT = 1e-9 T. Axial dipole g10 given for comparison.", "",
             "%8s %10s %10s %10s %14s %12s" % ("epoch", "g10 (nT)", "g11 (nT)", "h11 (nT)", "|g1| (nT)", "m (1e22 Am2)")]
    series = []
    for ep in epochs:
        mod = models[ep]
        g10, g11, h11 = mod.g[1][0], mod.g[1][1], mod.h[1][1]
        g1 = math.sqrt(g10 * g10 + g11 * g11 + h11 * h11)
        m = a_m ** 3 * g1 * 1e-9 / MU0_OVER_4PI
        series.append({"epoch": ep, "g10_nT": g10, "g1_nT": round(g1, 1), "moment_1e22_Am2": round(m / 1e22, 3)})
        lines.append("%8.1f %10.1f %10.1f %10.1f %14.1f %12.3f" % (ep, g10, g11, h11, g1, m / 1e22))
    last = models[epochs[-1]]
    g10, g11, h11 = last.g[1][0], last.g[1][1], last.h[1][1]
    gd10, gd11, hd11 = last.gd[1][0], last.gd[1][1], last.hd[1][1]
    g1 = math.sqrt(g10 * g10 + g11 * g11 + h11 * h11)
    g1d = (g10 * gd10 + g11 * gd11 + h11 * hd11) / g1
    g1_2030 = math.sqrt((g10 + 5 * gd10) ** 2 + (g11 + 5 * gd11) ** 2 + (h11 + 5 * hd11) ** 2)
    m_2030 = a_m ** 3 * g1_2030 * 1e-9 / MU0_OVER_4PI
    lines.append("%8.1f %10.1f %10.1f %10.1f %14.1f %12.3f  (forecast)" % (2030.0, g10 + 5 * gd10, g11 + 5 * gd11, h11 + 5 * hd11, g1_2030, m_2030 / 1e22))
    g1_1900 = series[0]["g1_nT"]
    lines += ["", "2025-2030 forecast rates: d|g1|/dt = %.1f nT/yr (%.3f%% per year, %.1f%% per century); dg10/dt = %+.1f nT/yr" % (g1d, 100.0 * g1d / g1, 10000.0 * g1d / g1, gd10),
              "Change 1900-2025: |g1| %.0f -> %.0f nT, %.1f%%; mean %.1f nT/yr" % (g1_1900, g1, 100.0 * (g1 - g1_1900) / g1_1900, (g1 - g1_1900) / 125.0),
              "Five-year paces of |g1| (nT/yr):"]
    paces = []
    for i in range(1, len(series)):
        pace = (series[i]["g1_nT"] - series[i - 1]["g1_nT"]) / 5.0
        paces.append({"from": series[i - 1]["epoch"], "to": series[i]["epoch"], "nT_per_yr": round(pace, 1)})
        lines.append("  %s-%s: %6.1f" % (int(series[i - 1]["epoch"]), int(series[i]["epoch"]), pace))
    # The axial dipole g10 alone, the measure of Finlay et al. 2016 (9% since 1840, 16 nT/yr, zero within 1,900 years)
    g10_1900 = series[0]["g10_nT"]
    axial_paces = [{"from": series[i - 1]["epoch"], "to": series[i]["epoch"], "nT_per_yr": round((abs(series[i]["g10_nT"]) - abs(series[i - 1]["g10_nT"])) / 5.0, 1)} for i in range(1, len(series))]
    axial_fast = max((p for p in axial_paces if p["from"] >= 1945), key=lambda p: abs(p["nT_per_yr"]))
    fastest = max((p for p in paces if p["from"] >= 1945), key=lambda p: abs(p["nT_per_yr"]))
    lines += ["", "Axial dipole g10 alone (Finlay et al. 2016's measure): %.0f nT in 1900, %.0f in 2025 (%.1f%%); forecast %+.1f nT/yr; fastest five-year pace since 1945 %.1f nT/yr (%d-%d)" % (
        abs(g10_1900), abs(g10), 100.0 * (abs(g10) - abs(g10_1900)) / abs(g10_1900), gd10, abs(axial_fast["nT_per_yr"]), axial_fast["from"], axial_fast["to"]),
        "Fastest five-year pace of the full dipole |g1| since 1945: %.1f nT/yr (%d-%d)" % (abs(fastest["nT_per_yr"]), fastest["from"], fastest["to"])]
    lines += ["", "Straight-line continuation of the dipole from 2025 (|g1| = %.0f nT):" % g1]
    targets = (("half of today (where past excursions began)", "half", 0.5 * g1), ("Laschamp's low, 5% of today", "laschamp_5pct", 0.05 * g1),
               ("5 uT equatorial dipole field (the 'hypomagnetic' threshold)", "5_uT", 5000.0), ("1 uT (the mole-rat study's residual field)", "1_uT", 1000.0), ("zero", "zero", 0.0))
    proj = {}
    for label, key, rate in (("2025-2030 forecast pace of the full dipole, %.1f nT/yr" % abs(g1d), "forecast_2025_2030", abs(g1d)),
                             ("1840-2010 average pace of the axial dipole, 16 nT/yr (Finlay et al. 2016)", "average_1840_2010_axial", 16.0),
                             ("fastest five-year pace of the full dipole since 1945, %.1f nT/yr (%d-%d)" % (abs(fastest["nT_per_yr"]), fastest["from"], fastest["to"]), "fastest_since_1945", abs(fastest["nT_per_yr"]))):
        lines.append("  at the %s:" % label)
        proj[key] = {"rate_nT_per_yr": round(rate, 1)}
        for name, tkey, target in targets:
            yrs = (g1 - target) / rate
            proj[key][tkey] = {"years": round(yrs), "year_AD": round(2025 + yrs)}
            lines.append("    %-62s %6.0f years, about AD %4.0f" % (name, yrs, 2025 + yrs))
    text = "\n".join(lines)
    with open(os.path.join(out, "dipole_igrf14.txt"), "w", encoding="utf-8") as f:
        f.write(text + "\n")
    print(text)
    J["dipole"] = {"series": series, "forecast_2030": {"g1_nT": round(g1_2030, 1), "moment_1e22_Am2": round(m_2030 / 1e22, 3)},
                   "g1_2025_nT": round(g1, 1), "rate_2025_2030_nT_per_yr": round(g1d, 1), "rate_percent_per_yr": round(100.0 * g1d / g1, 3),
                   "rate_percent_per_century": round(10000.0 * g1d / g1, 1),
                   "change_1900_2025_percent": round(100.0 * (g1 - g1_1900) / g1_1900, 1), "five_year_paces": paces,
                   "fastest_five_year_pace_since_1945": {"nT_per_yr": round(abs(fastest["nT_per_yr"]), 1), "from": fastest["from"], "to": fastest["to"]},
                   "axial_dipole_g10": {"g10_1900_nT": abs(g10_1900), "g10_2025_nT": abs(g10), "change_1900_2025_percent": round(100.0 * (abs(g10) - abs(g10_1900)) / abs(g10_1900), 1),
                                        "forecast_rate_2025_2030_nT_per_yr": gd10, "five_year_paces": axial_paces,
                                        "fastest_five_year_pace_since_1945": {"nT_per_yr": round(abs(axial_fast["nT_per_yr"]), 1), "from": axial_fast["from"], "to": axial_fast["to"]},
                                        "note": "Finlay et al. 2016's 9% since 1840, 16 nT/yr average and 'zero within 1,900 years' refer to this axial dipole"},
                   "straight_line_projections": proj}


def pole_speeds(np_path, out, J):
    """Year-to-year great-circle speeds of the north dip pole from NOAA's NP.xy (lon lat year)."""
    pts = []
    with open(np_path, encoding="utf-8") as f:
        for line in f:
            p = line.split()
            if len(p) < 3:
                continue
            try:
                lon, lat, yr = float(p[0]), float(p[1]), float(p[2])
            except ValueError:
                continue
            pts.append((yr, lat, lon))
    pts.sort()
    R = 6371.0

    def gc(p, q):
        la1, lo1, la2, lo2 = map(math.radians, (p[1], p[2], q[1], q[2]))
        d = math.sin((la2 - la1) / 2) ** 2 + math.cos(la1) * math.cos(la2) * math.sin((lo2 - lo1) / 2) ** 2
        return 2 * R * math.asin(math.sqrt(d))

    by_year = {p[0]: p for p in pts}
    lines = ["North dip pole from NOAA's NP.xy (model positions; lon east, lat north):", "", "%6s %8s %9s %10s" % ("year", "lat", "lon", "km/yr to next year")]
    yearly = []
    for i in range(len(pts) - 1):
        p, q = pts[i], pts[i + 1]
        if p[0] >= 1985 and q[0] - p[0] == 1:
            yearly.append({"year": int(p[0]), "lat": p[1], "lon": p[2], "km_per_yr": round(gc(p, q), 1)})
            lines.append("%6.0f %8.3f %9.3f %10.1f" % (p[0], p[1], p[2], gc(p, q)))
    lines.append("")
    means = {}
    for (y0, y1) in ((1990, 1995), (1995, 2000), (2000, 2005), (2005, 2010), (2010, 2015), (2015, 2020), (2020, 2025)):
        if y0 in by_year and y1 in by_year:
            d = sum(gc(by_year[y], by_year[y + 1]) for y in range(y0, y1) if y in by_year and y + 1 in by_year)
            means["%d-%d" % (y0, y1)] = round(d / (y1 - y0), 1)
            lines.append("Mean speed %d-%d: %.1f km/yr" % (y0, y1, d / (y1 - y0)))
    dist = None
    if 1831 in by_year and 2025 in by_year:
        dist = round(gc(by_year[1831], by_year[2025]))
        lines.append("Straight-line distance 1831 -> 2025: %.0f km" % dist)
    text = "\n".join(lines)
    with open(os.path.join(out, "pole_speeds.txt"), "w", encoding="utf-8") as f:
        f.write(text + "\n")
    print(text)
    J["north_dip_pole"] = {"yearly_speeds": yearly, "five_year_mean_speeds_km_per_yr": means, "straight_line_1831_to_2025_km": dist}


def main():
    here = os.path.dirname(os.path.abspath(__file__))
    dl = sys.argv[1] if len(sys.argv) > 1 else os.path.join(here, "geomag_inputs")
    out = sys.argv[2] if len(sys.argv) > 2 else os.path.join(here, "geomag_out")
    os.makedirs(out, exist_ok=True)
    fetch_inputs(dl)
    J = {"what": "Numbers behind factchecked.science/s0-news-2026-10-20/, from geomag_sh.py (IGRF-14, WMM2025, NOAA pole positions)",
         "inputs_sha256": {rel: sha256(os.path.join(dl, *rel.split("/"))) for rel in INPUTS}}
    wmm = read_wmm(os.path.join(dl, "wmm2025", "WMM2025.COF"))
    ok = validate_wmm(wmm, os.path.join(dl, "wmm2025", "WMM2025_TEST_VALUES.txt"), out, J)
    print()
    epochs, models = read_igrf(os.path.join(dl, "igrf14", "igrf14coeffs.txt"))
    igrf2025, dgrf2020 = models[epochs[-1]], models[2020.0]
    grid_and_cities(igrf2025, dgrf2020, wmm, out, J)
    city_history(models, J)
    print()
    dipole(epochs, models, out, J)
    print()
    pole_speeds(os.path.join(dl, "ncei-poles", "NP.xy"), out, J)
    with open(os.path.join(out, "s0-news-2026-10-20-numbers.json"), "w", encoding="utf-8") as f:
        json.dump(J, f, indent=1)
    print("\nwrote", os.path.join(out, "s0-news-2026-10-20-numbers.json"))
    if not ok:
        print("\nWARNING: validation against WMM2025 test values FAILED; do not use the numbers above.")
        sys.exit(1)


if __name__ == "__main__":
    main()