#!/usr/bin/env python3
"""How much has the magnetic field changed where the review's tipping systems sit?

Fact-Checked Science, November 2026, for the paper check /biosphere-tipping-warning/ ("Our analysis").
Standard library only. Conventions follow /field-loss-projections/: Schmidt semi-normalised coefficients,
linear interpolation between IGRF epochs, WGS84 surface, reference radius 6371.2 km.

Rules:
  Sites. 37 coordinates: (a) the four printed in the caption of Figure 8 of Abrams et al. (2026), p. 25; (b) the centres
  of all 33 "Major Systems" polygons in the one map that the authors' script saves, PREDICT_review_map.R (Zenodo
  10.5281/zenodo.17505626), the block it labels "simplified representations"; its plot caption calls them "approximate
  extent". An earlier version of this script used only 13 of the 33 and said the choice wasn't ours. A fresh review
  caught that, and the three largest falls were among the 20 left out.
  Checks on the extremes, reported apart. The saved map's 17 point markers ("smaller or dispersed systems"), and the
  script's two earlier maps, which the saved one replaces: 33 points at "approximate coordinates" and 6 polygons. Those
  coordinates are rough (the Sahel sits at 35 N, 30 E). They reach a larger rise, and a larger fall since 2000, but not
  the largest fall since 1900.
  One name is corrected: the script spells Okavango "Okavanga".
  Epochs. 1900, 1950, 2000 and 2025, with the rate for 2000-2025. IGRF labels its epochs before 1945 non-definitive.
  Scale comparison. The total field F at each site, and the world's range of F in 2025 on a 1-degree grid.
  A check. The 2025 South Atlantic minimum must come out at 22,071 nT at 26.22 S, 60.03 W (the IGRF-14
  paper's figure), within 1 nT, or the run stops.
What this does not show: it compares sizes only. It says nothing on whether any organism responds.

Run it
  python3 field_at_sites.py [folder] [out.json]
  The folder (default: the script's own) must hold, or will receive, igrf14coeffs.txt, IGRF-14's Gauss coefficients,
  from https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt. It is checked against the SHA-256 of the copy
  we used on November 1, 2026. The script prints a summary and writes field-at-sites-numbers.json (or out.json).
"""
import hashlib
import json
import math
import os
import sys
import urllib.request

sys.stdout.reconfigure(encoding="utf-8")
IGRF_URL = "https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt"
IGRF_SHA256 = "8f8d88403028fc4ee92c4f38d97b46e0a87e2cfc496045b43c9e26c1d6b0903c"
NMAX = 13
A_REF = 6371.2
WGS84_A = 6378.137
WGS84_E2 = 0.00669437999014

SITES = [
    # Figure 8 caption, p. 25 of the review
    ("Fig. 8: tropical forest, DR Congo", -1.07, 22.59),
    ("Fig. 8: mangrove forest, Bangladesh", 22.03, 89.43),
    ("Fig. 8: permafrost, Alaska", 65.58, -164.02),
    ("Fig. 8: dryland patterned vegetation, Sudan", 11.12, 28.37),
    # Centres of all 33 "Major Systems" polygons in the authors' PREDICT_review_map.R, in its order
    ("Script: Amazon Basin", -5.0, -60.0),
    ("Script: Congo Basin", 0.0, 20.0),
    ("Script: Southeast Asian Forests", -2.0, 110.0),
    ("Script: North American Boreal", 55.0, -110.0),
    ("Script: Eurasian Taiga", 60.0, 90.0),
    ("Script: Temperate Forests", 40.0, -80.0),
    ("Script: African Savannas", -15.0, 25.0),
    ("Script: Cerrado", -15.0, -50.0),
    ("Script: Northern Australian Savanna", -17.0, 135.0),
    ("Script: Indian Savanna", 20.0, 77.0),
    ("Script: Great Plains", 40.0, -100.0),
    ("Script: Eurasian Steppes", 48.0, 70.0),
    ("Script: Pampas", -35.0, -60.0),
    ("Script: Mediterranean", 38.0, 10.0),
    ("Script: Sahel", 15.0, 10.0),
    ("Script: Central Asian Drylands", 35.0, 70.0),
    ("Script: Australian Drylands", -28.0, 135.0),
    ("Script: Patagonian Steppe", -45.0, -65.0),
    ("Script: Arctic Permafrost", 68.0, -140.0),
    ("Script: Siberian Permafrost", 65.0, 110.0),
    ("Script: Tibetan Permafrost", 35.0, 90.0),
    ("Script: Rocky Mountains", 45.0, -115.0),
    ("Script: European Alps", 46.0, 10.0),
    ("Script: Himalayas/Tibetan Plateau", 30.0, 85.0),
    ("Script: Andes", -20.0, -70.0),
    ("Script: Canadian Peatlands", 55.0, -95.0),
    ("Script: West Siberian Peatlands", 60.0, 75.0),
    ("Script: European Peatlands", 60.0, 25.0),
    ("Script: Tropical Peatlands", 0.0, 105.0),
    ("Script: Prairie Potholes", 48.0, -100.0),
    ("Script: Pantanal", -19.0, -57.0),
    ("Script: African Wetlands", 3.0, 35.0),
    ("Script: Okavango Delta", -19.0, 22.0),
]
# The saved map's 17 point markers, reported apart
POINTS = [
    ("Florida Mangroves", 25.0, -80.0), ("Amazon Mangroves", -5.0, -60.0), ("Caribbean Mangroves", 20.0, -90.0),
    ("West African Mangroves", 5.0, 10.0), ("East African Mangroves", -25.0, 35.0), ("Indian Mangroves", 20.0, 90.0),
    ("Southeast Asian Mangroves", 5.0, 100.0), ("Indonesian Mangroves", -5.0, 130.0), ("Australian Mangroves", -20.0, 150.0),
    ("SE US Forested Wetlands", 32.0, -85.0), ("Amazon Flooded Forests", 5.0, -55.0), ("Congo Swamp Forests", 0.0, 25.0),
    ("Asian Forested Wetlands", 30.0, 100.0), ("Andean Cloud Forests", -10.0, -75.0), ("African Montane", 0.0, -5.0),
    ("Madagascar Montane", -5.0, 35.0), ("Asian Montane", 25.0, 95.0),
]
# The script's two earlier maps, which the saved one replaces, reported apart, with its names and coordinates
EARLIER = [
    ("First map: Amazon Rainforest", -5.0, -60.0), ("First map: North American Boreal", 65.0, -70.0), ("First map: Congo Basin", 0.0, 10.0),
    ("First map: Brazilian Cerrado", -15.0, -60.0), ("First map: Eurasian Taiga", 65.0, 130.0), ("First map: African Savanna", -5.0, 25.0),
    ("First map: Southern African Savanna", -20.0, 15.0), ("First map: Caribbean Mangroves", 10.0, -50.0), ("First map: North American Drylands", 50.0, -100.0),
    ("First map: Central Asian Drylands", 35.0, 110.0), ("First map: Australian Drylands", 15.0, 95.0), ("First map: Mediterranean Drylands", 30.0, -110.0),
    ("First map: Sahel", 35.0, 30.0), ("First map: European Wetlands", 50.0, 20.0), ("First map: Asian Wetlands", 25.0, 95.0),
    ("First map: Australian Wetlands", -20.0, 130.0), ("First map: Arctic Permafrost", 60.0, -105.0), ("First map: European Permafrost", 40.0, 10.0),
    ("First map: Tibetan Permafrost", 30.0, 100.0), ("First map: Prairie Pothole Wetlands", 40.0, -95.0), ("First map: Northern European Wetlands", 55.0, 5.0),
    ("First map: Southeast Asian Mangroves", 20.0, 100.0), ("First map: Australian Mangroves", 10.0, 150.0), ("First map: Central American Mangroves", 15.0, -75.0),
    ("First map: West African Mangroves", 5.0, 20.0), ("First map: Indian Mangroves", 10.0, 105.0), ("First map: Red Sea Mangroves", -5.0, 140.0),
    ("First map: Rocky Mountains", 45.0, -115.0), ("First map: European Alps", 43.0, 15.0), ("First map: Tibetan Plateau", 30.0, 85.0),
    ("First map: Canadian Peatlands", 55.0, -100.0), ("First map: Scandinavian Peatlands", 55.0, 25.0), ("First map: Siberian Peatlands", 60.0, 100.0),
    ("Second map: Amazon Basin", -5.0, -60.0), ("Second map: Congo Basin", 0.0, 20.0), ("Second map: North American Boreal", 60.0, -110.0),
    ("Second map: Eurasian Taiga", 65.0, 90.0), ("Second map: Sahel", 15.0, 10.0), ("Second map: Tibetan Plateau", 35.0, 90.0),
]


def zeros():
    return [[0.0] * (NMAX + 1) for _ in range(NMAX + 1)]


def read_igrf(path):
    epochs, rows = None, {}
    with open(path, encoding="ascii") as f:
        for line in f:
            if line.startswith("g/h"):
                epochs = [float(x) for x in line.split()[3:-1]]
            elif line[:1] in "gh" and not line.startswith("g/h"):
                p = line.split()
                rows[(p[0], int(p[1]), int(p[2]))] = [float(x) for x in p[3:]]
    return epochs, {k: v[:-1] for k, v in rows.items()}, {k: v[-1] for k, v in rows.items()}


def igrf_at(model, t):
    epochs, main, sv = model
    g, h = zeros(), zeros()
    if t >= epochs[-1]:
        i0 = len(epochs) - 1
        for (gh, n, m), v in main.items():
            (g if gh == "g" else h)[n][m] = v[i0] + (t - epochs[-1]) * sv[(gh, n, m)]
    else:
        i0 = max(i for i, e in enumerate(epochs) if e <= t)
        w = (t - epochs[i0]) / (epochs[i0 + 1] - epochs[i0])
        for (gh, n, m), v in main.items():
            (g if gh == "g" else h)[n][m] = v[i0] + w * (v[i0 + 1] - v[i0])
    return g, h


def geodetic_to_geocentric(lat_deg, h_km=0.0):
    phi = math.radians(lat_deg)
    s, c = math.sin(phi), math.cos(phi)
    N = WGS84_A / math.sqrt(1.0 - WGS84_E2 * s * s)
    x = (N + h_km) * c
    z = (N * (1.0 - WGS84_E2) + h_km) * s
    phic = math.atan2(z, x)
    return math.hypot(x, z), phic, phi - phic


def legendre(theta):
    ct, st = math.cos(theta), math.sin(theta)
    P, D = zeros(), zeros()
    P[0][0], P[1][0], D[1][0], P[1][1], D[1][1] = 1.0, ct, -st, st, ct
    for n in range(2, NMAX + 1):
        k = math.sqrt((2.0 * n - 1.0) / (2.0 * n))
        P[n][n] = k * st * P[n - 1][n - 1]
        D[n][n] = k * (ct * P[n - 1][n - 1] + st * D[n - 1][n - 1])
        for m in range(n):
            a2 = math.sqrt((n - 1.0) ** 2 - m * m) if n - 1 >= m else 0.0
            den = math.sqrt(n * n - m * m)
            P2 = P[n - 2][m] if n - 2 >= m else 0.0
            D2 = D[n - 2][m] if n - 2 >= m else 0.0
            P[n][m] = ((2.0 * n - 1.0) * ct * P[n - 1][m] - a2 * P2) / den
            D[n][m] = ((2.0 * n - 1.0) * (ct * D[n - 1][m] - st * P[n - 1][m]) - a2 * D2) / den
    return P, D


class Row:
    def __init__(self, lat_deg, h_km=0.0):
        self.r, phic, self.delta = geodetic_to_geocentric(lat_deg, h_km)
        theta = math.pi / 2.0 - phic
        self.st = math.sin(theta)
        self.P, self.D = legendre(theta)
        self.rr = [(A_REF / self.r) ** (n + 2) for n in range(NMAX + 1)]

    def F(self, g, h, lon_deg):
        lam = math.radians(lon_deg)
        br = bt = bp = 0.0
        for m in range(NMAX + 1):
            c, s = math.cos(m * lam), math.sin(m * lam)
            for n in range(max(1, m), NMAX + 1):
                f, P, D = self.rr[n], self.P[n][m], self.D[n][m]
                gnm, hnm = g[n][m], h[n][m]
                br += (n + 1) * f * (gnm * c + hnm * s) * P
                bt -= f * (gnm * c + hnm * s) * D
                bp -= m * f * (-gnm * s + hnm * c) * P
        bp /= self.st
        xc, zc = -bt, -br
        cd, sd = math.cos(self.delta), math.sin(self.delta)
        x, y, z = xc * cd + zc * sd, bp, zc * cd - xc * sd
        return math.sqrt(x * x + y * y + z * z)


def F_at(model, t, lat, lon):
    g, h = igrf_at(model, t)
    return Row(lat).F(g, h, lon)


def main():
    folder = sys.argv[1] if len(sys.argv) > 1 else os.path.dirname(os.path.abspath(__file__))
    target = sys.argv[2] if len(sys.argv) > 2 else os.path.join(folder, "field-at-sites-numbers.json")
    path = os.path.join(folder, "igrf14coeffs.txt")
    if not os.path.exists(path):
        req = urllib.request.Request(IGRF_URL, headers={"User-Agent": "Mozilla/5.0"})
        with urllib.request.urlopen(req, timeout=120) as r, open(path, "wb") as f:
            f.write(r.read())
    with open(path, "rb") as f:
        digest = hashlib.sha256(f.read()).hexdigest()
    if digest != IGRF_SHA256:
        print("Note: igrf14coeffs.txt differs from the copy we used (SHA-256 %s...); results may differ." % IGRF_SHA256[:16])
    model = read_igrf(path)
    out = {"source": IGRF_URL, "source_sha256": digest, "sites": [], "epochs": [1900.0, 1950.0, 2000.0, 2025.0]}

    # The check.
    g25, h25 = igrf_at(model, 2025.0)
    fmin_check = Row(-26.22).F(g25, h25, -60.03)
    out["check_saa_min_2025_nT"] = round(fmin_check, 1)
    print("Check: F at 26.22S 60.03W, 2025.0 = %.1f nT (target 22,071)" % fmin_check)
    if abs(fmin_check - 22071) > 1.0:
        raise SystemExit("check failed")

    # World range in 2025 on a 1-degree grid.
    lo, hi = (1e9, None, None), (-1e9, None, None)
    for lat10 in range(-89, 90):
        row = Row(float(lat10))
        for lon in range(-180, 180):
            v = row.F(g25, h25, float(lon))
            if v < lo[0]:
                lo = (v, lat10, lon)
            if v > hi[0]:
                hi = (v, lat10, lon)
    out["world_2025_grid"] = {"min_nT": round(lo[0], 0), "min_at": lo[1:], "max_nT": round(hi[0], 0), "max_at": hi[1:]}
    print("World F in 2025, 1-degree grid: min %.0f nT at %s; max %.0f nT at %s; ratio %.2f"
          % (lo[0], lo[1:], hi[0], hi[1:], hi[0] / lo[0]))

    print("%-46s %8s %8s %8s %8s %9s %9s %9s" % ("Site", "1900", "1950", "2000", "2025", "1900-2025", "2000-2025", "%/decade"))
    for name, lat, lon in SITES:
        vals = [F_at(model, t, lat, lon) for t in out["epochs"]]
        ch_all = 100.0 * (vals[3] - vals[0]) / vals[0]
        ch_rec = 100.0 * (vals[3] - vals[2]) / vals[2]
        per_dec = ch_rec / 2.5
        out["sites"].append({"name": name, "lat": lat, "lon": lon, "F_nT": [round(v, 0) for v in vals],
                             "pct_1900_2025": round(ch_all, 1), "pct_2000_2025": round(ch_rec, 1),
                             "pct_per_decade_2000_2025": round(per_dec, 2)})
        print("%-46s %8.0f %8.0f %8.0f %8.0f %8.1f%% %8.1f%% %8.2f%%"
              % (name, vals[0], vals[1], vals[2], vals[3], ch_all, ch_rec, per_dec))

    fs = [s["F_nT"][3] for s in out["sites"]]
    all_pct = [s["pct_1900_2025"] for s in out["sites"]]
    rec_pct = [s["pct_2000_2025"] for s in out["sites"]]
    script = [s for s in out["sites"] if s["name"].startswith("Script: ")]
    ranked = sorted(out["sites"], key=lambda s: s["pct_1900_2025"])
    out["summary"] = {"n_sites": len(out["sites"]), "F2025_min": min(fs), "F2025_max": max(fs), "ratio": round(max(fs) / min(fs), 2),
                      "pct_1900_2025_range": [min(all_pct), max(all_pct)], "pct_2000_2025_range": [min(rec_pct), max(rec_pct)],
                      "n_rose_1900_2025": sum(p > 0 for p in all_pct),
                      "n_script_sites": len(script), "n_script_rose_1900_2025": sum(s["pct_1900_2025"] > 0 for s in script),
                      "largest_falls_1900_2025": [[s["name"], s["pct_1900_2025"]] for s in ranked[:5]],
                      "largest_rises_1900_2025": [[s["name"], s["pct_1900_2025"]] for s in ranked[::-1][:3]]}
    print("Across these %d sites in 2025: F from %.0f to %.0f nT, ratio %.2f" % (len(fs), min(fs), max(fs), max(fs) / min(fs)))
    print("1900-2025: %.1f%% to %.1f%%; rose at %d of %d (%d of the script's %d). 2000-2025: %.1f%% to %.1f%%"
          % (min(all_pct), max(all_pct), sum(p > 0 for p in all_pct), len(all_pct), out["summary"]["n_script_rose_1900_2025"], len(script), min(rec_pct), max(rec_pct)))
    print("Largest falls:", out["summary"]["largest_falls_1900_2025"])
    print("Largest rises:", out["summary"]["largest_rises_1900_2025"])
    for key, rows, label in (("point_markers", POINTS, "The saved map's 17 point markers"), ("earlier_maps", EARLIER, "The script's two earlier maps")):
        pts = []
        for name, lat, lon in rows:
            vals = [F_at(model, t, lat, lon) for t in (1900.0, 2000.0, 2025.0)]
            pts.append({"name": name, "lat": lat, "lon": lon, "F_nT_1900_2000_2025": [round(v, 0) for v in vals],
                        "pct_1900_2025": round(100.0 * (vals[2] - vals[0]) / vals[0], 1), "pct_2000_2025": round(100.0 * (vals[2] - vals[1]) / vals[1], 1)})
        out[key] = pts
        p_all = [q["pct_1900_2025"] for q in pts]
        p_rec = [q["pct_2000_2025"] for q in pts]
        p_f = [q["F_nT_1900_2000_2025"][2] for q in pts]
        print("%s, apart: 1900-2025 %.1f%% to %.1f%%; 2000-2025 %.1f%% to %.1f%%; F 2025 %.0f to %.0f nT"
              % (label, min(p_all), max(p_all), min(p_rec), max(p_rec), min(p_f), max(p_f)))
    with open(target, "w", encoding="utf-8", newline="\n") as f:
        json.dump(out, f, ensure_ascii=False, indent=1)


if __name__ == "__main__":
    main()