#!/usr/bin/env python3
"""The electrojet over India and its drifting magnetic equator: the numbers we derive from the field models and the
observatories' records on Fact-Checked Science's page https://factchecked.science/equatorial-electrojet-trend/ .
Standard library only.

Run it
  python3 electrojet_numbers.py [folder [output.json]]
  The folder (default: the script's own) is searched, with its subfolders, for the input files below. Any that is
  missing is downloaded into the folder from its source, and every file's SHA-256 is compared with the copy we used
  on October 30, 2026; a different hash means the source has changed the file, so numbers may differ. GFZ's file
  grows by a day every day, and WDC Edinburgh stamps each file with the day it was made, so those two hashes always
  differ; the values we use do not change. It prints a summary and writes electrojet-numbers.json (into the folder,
  or to the path given). Once the inputs are in place a run takes about a quarter of a minute; downloading them first
  takes longer, depending on the connection.

Inputs (about 75 MB in all; read as data, never executed)
  igrf14coeffs.txt   IGRF-14 Gauss coefficients 1900-2025 and secular variation 2025-30 (IAGA V-MOD), NOAA
  WMM2025COF.zip     World Magnetic Model 2025 with NOAA's 12 official test values (checks our field synthesis)
  Kp_ap_Ap_SN_F107_since_1932.txt   GFZ Potsdam: daily Ap and F10.7 (CC BY 4.0; Matzka et al. 2021)
  hour_<obs>_<years>.csv   hourly values from the WDC for Geomagnetism, Kyoto (HAPI, "best available"), of the
                     Indian Institute of Geomagnetism's observatories Trivandrum (TRD), Alibag (ABG), Kodaikanal
                     (KOD), Annamalainagar (ANN), Ettaiyapuram (ETT) and Tirunelveli (TIR, for a cross-check)
  bgs-hour-tir-1999-2014-hdzf-csv.zip   Tirunelveli's hourly values from the WDC for Geomagnetism, Edinburgh
  No observatory data are distributed with the page: the Indian Institute of Geomagnetism does not allow its data
  to be passed on without its permission. The World Data Centres provide them for scientific use; if you use them,
  acknowledge the Indian Institute of Geomagnetism and the WDC you took them from.

What it computes (WGS84 surface; 1 degree of latitude = 110.57 km, the meridian arc at 8.5 N)
  1. A check of our spherical-harmonic field synthesis against NOAA's WMM2025 test values.
  2. The dip equator (inclination 0, Z = 0) along meridians 72-80 E every year 1960-2030 in IGRF-14, by
     bisection: its drift over spans and five-year steps; how far each observatory lies from it; when it passed
     south of Kanyakumari, India's southern tip; how much of the drift is the centred dipole's equator, and the
     dipole pole's latitude, geocentric and geodetic; the field on the dip equator; WMM2025's forecast; and the dip equator over Tatuoca,
     Brazil, for contrast.
  3. The observatories' own annual means of H and Z (TRD, TIR, ABG) against IGRF-14: the slopes and offsets.
  4. Our re-run of a long-term trend in the electrojet (EEJ), after Zhang et al. (2026, Geophysical Research
     Letters, 10.1029/2026GL126510), with the textbook index; the paper's own method differs in its details (its
     nightly level, its storm filter, its fit and its errors; the page's section 4 lists them):
       local time LT = UT + 5 h; dH = H minus the mean of that day's 00, 01 and 02 LT values;
       EEJ = dH(TRD or TIR) - dH(ABG); monthly medians for each LT hour 07-18 (at least 15 days);
       per LT hour, least squares on the month's mean F10.7 (observed), annual and semiannual harmonics and a
       linear trend (nT per decade), with the trend's standard error inflated by sqrt((1+r1)/(1-r1)) for lag-1
       autocorrelation r1 of the residuals, and a two-sided 95% test (|t| > 1.96).
     A: TRD-ABG, 1967-01 to 1999-09 (TRD's record ends there); variants: F and F^2, two-step, days with Ap > 20
        dropped. B: TRD-ABG to 1998-12, then TIR-ABG 1999-01 to 2014-12 (the public record we found ends in 2014;
        the paper runs to 2015), TIR scaled to TRD per LT hour by their ratio in the Jan-Sep 1999 overlap;
        variants: unscaled, ending in 2011, Ap > 20 dropped, F and F^2. C: TIR-ABG alone, 1999-2014.
     And the solar cycle's swing at 11 LT in A beside its trend: the range of the 13-month running mean of F10.7
     times the fitted nT per unit of F10.7.
  5. How fast the EEJ's ground signal falls off away from the dip equator: the ratio of KOD's, ANN's and ETT's EEJ
     to TRD's at 10-12 LT, in months when both ran, against their distances from the dip equator; a Gaussian
     exp(-ln2 (x/s)^2) fitted for the half-width s; and what each station's drift alone did to its signal.
  6. The trends of 4 again, with each month divided by that profile factor for the station's distance that month.
"""
import csv
import datetime as dt
import hashlib
import json
import math
import os
import sys
import urllib.request
import zipfile

HAPI = ("https://wdc.kugi.kyoto-u.ac.jp/hapi/data?id=hour_{0}&time.min={1}&time.max={2}"
        "&parameters=geomagneticFieldH,geomagneticFieldX,geomagneticFieldY,geomagneticFieldZ,versionCode&format=csv")
INPUTS = {  # file name: (where it comes from, SHA-256 of the copy we used on 2026-10-30)
    "igrf14coeffs.txt": ("https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt",
                         "8f8d88403028fc4ee92c4f38d97b46e0a87e2cfc496045b43c9e26c1d6b0903c"),
    "WMM2025COF.zip": ("https://www.ncei.noaa.gov/sites/default/files/2025-01/WMM2025COF.zip",
                       "999080a12c95122c2031aeac8be77229d97584530e87e0f1a6df7d946ecf90c5"),
    # grows by a day every day, so the hash differs; its values for 1967-2014 do not
    "Kp_ap_Ap_SN_F107_since_1932.txt": ("https://kp.gfz.de/app/files/Kp_ap_Ap_SN_F107_since_1932.txt",
                                        "13ec3dcfce8aed6b7356ce341aeaca9409a85a016d200936c282cff038cc033d"),
    "hour_trd_1957-1999.csv": (HAPI.format("trd", "1957-10-01T00:29:30Z", "1999-10-01T00:29:30Z"),
                               "53b2c41d2992335acde0075cb92fd9a72b13c79b958143f53fa4fe3ae58ccb02"),
    "hour_abg_1957-1990.csv": (HAPI.format("abg", "1957-01-01T00:00:00Z", "1990-01-01T00:00:00Z"),
                               "cc7b96b4ede039119c5ecd935ba7d386f2a6ee4c746d89d1a9c17037e4fbac1b"),
    "hour_abg_1990-2023.csv": (HAPI.format("abg", "1990-01-01T00:00:00Z", "2024-01-01T00:29:30Z"),
                               "768636ffc3f7fb823c98534d5d190e8f3e81194bc2fa93081b529882e54771ed"),
    "hour_kod_1970-1980.csv": (HAPI.format("kod", "1970-01-01T00:00:00Z", "1980-07-01T00:29:30Z"),
                               "e40f28f35759a35f58ed64093e58acc168e15b7e9e6c1c1df2fc1d2b15bb99ce"),
    "hour_ann_1970-1980.csv": (HAPI.format("ann", "1970-01-01T00:00:00Z", "1981-01-01T00:00:00Z"),
                               "b03919274c63f7fefd71bbb92fac1298f6da179a7e364634941c26f8bf0267c0"),
    "hour_ett_1996-1998.csv": (HAPI.format("ett", "1996-01-01T00:29:30Z", "1998-05-01T00:29:30Z"),
                               "06ab79e1025c24a53cdab501db52581f6fd243ae0fd4792f70741cda0866ad19"),
    "hour_tir_1999-2011.csv": (HAPI.format("tir", "1999-01-01T00:29:30Z", "2012-01-01T00:29:30Z"),
                               "2ba5e137a86e94d66a05f1e1fd97a1c7c6873c0d2ce81a320d6d3dbf0675b9ba"),
    # each CSV in the zip is stamped with the day it was generated, so the hash differs on every download
    "bgs-hour-tir-1999-2014-hdzf-csv.zip": (
        "https://wdcapi.bgs.ac.uk/data/zipped-data/hour/tir?start_dt=1999-01-01&end_dt=2015-01-01&file_format=csv&orientation=HDZF",
        "42083ba3e4eeb2e39cd28827917e38a42a19a419617555312961eec1f2b4a9d1"),
}
CHANGES_ANYWAY = {"Kp_ap_Ap_SN_F107_since_1932.txt", "bgs-hour-tir-1999-2014-hdzf-csv.zip"}

A_REF = 6371.2
WGS84_A = 6378.137
WGS84_F = 1.0 / 298.257223563
WGS84_E2 = WGS84_F * (2.0 - WGS84_F)
NMAX = 13
FILL = 99000.0
LT_OFFSET = 5
HOURS = list(range(7, 19))
STATIONS = {  # WDC Kyoto's coordinates for the observatories
    "TRD": ("Trivandrum", 8.48, 76.95), "TIR": ("Tirunelveli", 8.70, 77.80), "ETT": ("Ettaiyapuram", 9.167, 78.017),
    "KOD": ("Kodaikanal", 10.233, 77.467), "ANN": ("Annamalainagar", 11.367, 79.683), "ABG": ("Alibag", 18.638, 72.872),
}
KANYAKUMARI = (8.078, 77.541)
TATUOCA = (-1.205, 311.487)


# ------------------------------------------------------------------ inputs
def find_or_fetch(folder, name):
    for dp, _, fns in os.walk(folder):
        if name in fns:
            path = os.path.join(dp, name)
            break
    else:
        url = INPUTS[name][0]
        path = os.path.join(folder, name)
        print("downloading", name, "from", url.split("?")[0])
        req = urllib.request.Request(url, headers={"User-Agent": "factchecked.science electrojet_numbers.py"})
        with urllib.request.urlopen(req, timeout=300) as r, open(path, "wb") as f:
            f.write(r.read())
    with open(path, "rb") as f:
        digest = hashlib.sha256(f.read()).hexdigest()
    ours = INPUTS[name][1]
    same = digest == ours
    if not same and name not in CHANGES_ANYWAY:
        print("Note: %s differs from the copy we used (SHA-256 %s...); results may differ." % (name, ours[:16]))
    return path, {"url": INPUTS[name][0], "sha256": digest, "same_as_ours": same}


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 zeros():
    return [[0.0] * (NMAX + 1) for _ in range(NMAX + 1)]


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 read_wmm(zpath):
    z = zipfile.ZipFile(zpath)
    names = {os.path.basename(n): n for n in z.namelist()}
    cof = z.read(names["WMM.COF"]).decode("ascii")
    tests = z.read(names["WMM2025_TEST_VALUES.txt"]).decode("ascii")
    g, h, gd, hd = zeros(), zeros(), zeros(), zeros()
    lines = cof.splitlines()
    epoch = float(lines[0].split()[0])
    for line in lines[1:]:
        p = line.split()
        if len(p) < 6 or p[0].startswith("9999"):
            continue
        n, m = int(p[0]), int(p[1])
        g[n][m], h[n][m], gd[n][m], hd[n][m] = map(float, p[2:6])
    return epoch, g, h, gd, hd, tests


# ------------------------------------------------------------------ field synthesis (as in the site's field_decline.py)
def geodetic_to_geocentric(lat_deg, h_km):
    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


def field(g, h, lat, lon, hkm=0.0):
    """X (north), Y (east), Z (down), nT, at geodetic latitude and longitude (degrees) and height hkm."""
    r, phic, delta = geodetic_to_geocentric(lat, hkm)
    theta = math.pi / 2.0 - phic
    st = math.sin(theta)
    P, D = legendre(theta)
    lam = math.radians(lon)
    br = bt = bp = 0.0
    for n in range(1, NMAX + 1):
        f = (A_REF / r) ** (n + 2)
        for m in range(n + 1):
            c, s = math.cos(m * lam), math.sin(m * lam)
            gnm, hnm = g[n][m], h[n][m]
            br += (n + 1) * f * (gnm * c + hnm * s) * P[n][m]
            bt -= f * (gnm * c + hnm * s) * D[n][m]
            bp -= m * f * (-gnm * s + hnm * c) * P[n][m]
    bp /= st
    xc, zc = -bt, -br
    cd, sd = math.cos(delta), math.sin(delta)
    return xc * cd + zc * sd, bp, zc * cd - xc * sd


def validate(zpath):
    epoch, g, h, gd, hd, tests = read_wmm(zpath)
    worst = {k: 0.0 for k in ("X", "Y", "Z", "F")}
    n = 0
    for line in tests.splitlines():
        if line.startswith("#") or not line.strip():
            continue
        p = [float(x) for x in line.split()]
        t, hkm, lat, lon = p[0], p[1], p[2], p[3]
        gg = [[g[i][j] + (t - epoch) * gd[i][j] for j in range(NMAX + 1)] for i in range(NMAX + 1)]
        hh = [[h[i][j] + (t - epoch) * hd[i][j] for j in range(NMAX + 1)] for i in range(NMAX + 1)]
        X, Y, Z = field(gg, hh, lat, lon, hkm)
        F = math.sqrt(X * X + Y * Y + Z * Z)
        for key, mine, ref in (("X", X, p[4]), ("Y", Y, p[5]), ("Z", Z, p[6]), ("F", F, p[8])):
            worst[key] = max(worst[key], abs(mine - ref))
        n += 1
    return {"test_points": n, "largest_abs_difference_nT": {k: round(v, 3) for k, v in worst.items()},
            "pass_within_0.1_nT": all(v <= 0.1 for v in worst.values())}


def dip_equator_lat(g, h, lon, lo=-15.0, hi=25.0):
    """Geodetic latitude where Z = 0 on meridian lon, by bisection (one sign change in the bracket)."""
    zlo = field(g, h, lo, lon)[2]
    if zlo * field(g, h, hi, lon)[2] > 0:
        raise ValueError("no sign change of Z between %s and %s at %s E" % (lo, hi, lon))
    for _ in range(60):
        mid = 0.5 * (lo + hi)
        zm = field(g, h, mid, lon)[2]
        if zm * zlo > 0:
            lo, zlo = mid, zm
        else:
            hi = mid
    return 0.5 * (lo + hi)


def km_per_degree(lat):
    s = math.sin(math.radians(lat))
    return WGS84_A * (1 - WGS84_E2) / (1 - WGS84_E2 * s * s) ** 1.5 * math.pi / 180.0


KM = km_per_degree(8.5)


# ------------------------------------------------------------------ 2. the dip equator
def dip_equator_numbers(model, wmm_path):
    out = {}
    years = [float(y) for y in range(1960, 2031)]
    lons = [72.0, 76.0, 76.95, 77.0, 77.541, 77.8, 78.0, 80.0]
    deq = {}
    for t in years:
        g, h = igrf_at(model, t)
        deq[t] = {lon: dip_equator_lat(g, h, lon) for lon in lons}
    out["latitude_by_year"] = {"%d" % t: {"%gE" % lon: round(deq[t][lon], 4) for lon in lons} for t in years}

    def rate(lon, a, b):
        d = deq[b][lon] - deq[a][lon]
        return {"from": int(a), "to": int(b), "lat_start": round(deq[a][lon], 3), "lat_end": round(deq[b][lon], 3),
                "moved_km": round(d * KM, 1), "km_per_year": round(d * KM / (b - a), 2)}
    spans = [(1967.0, 1999.0), (1999.0, 2015.0), (1967.0, 2015.0), (2015.0, 2025.0), (2020.0, 2025.0), (1967.0, 2025.0),
             (2025.0, 2030.0)]
    out["drift_77E"] = [rate(77.0, a, b) for a, b in spans]
    out["drift_76.95E_TRD"] = [rate(76.95, a, b) for a, b in spans[:3]]
    out["drift_77.8E_TIR"] = [rate(77.8, a, b) for a, b in spans[:3]]
    out["five_year_paces_77E"] = [rate(77.0, float(y), float(y + 5)) for y in range(1960, 2030, 5)]
    # Kanyakumari: when the dip equator passed south of India's southern tip on its meridian (linear within the year)
    lat_k, lon_k = KANYAKUMARI
    t1 = next(t for t in years if deq[t][lon_k] < lat_k)
    a, b = deq[t1 - 1.0][lon_k], deq[t1][lon_k]
    out["passed_south_of_kanyakumari"] = round(t1 - 1.0 + (a - lat_k) / (a - b), 2)
    # stations
    st = {}
    for code in ("TRD", "TIR", "ABG", "KOD", "ANN", "ETT"):
        name, lat, lon = STATIONS[code]
        rows = {}
        for t in (1967.0, 1980.0, 1985.0, 1986.0, 1987.0, 1990.0, 1999.0, 2005.0, 2010.0, 2015.0, 2025.0):
            g, h = igrf_at(model, t)
            X, Y, Z = field(g, h, lat, lon)
            inc = math.degrees(math.atan2(Z, math.hypot(X, Y)))
            rows["%d" % t] = {"km_north_of_dip_equator": round((lat - dip_equator_lat(g, h, lon)) * KM, 1),
                              "inclination_deg": round(inc, 3),
                              "dip_latitude_deg": round(math.degrees(math.atan(math.tan(math.radians(inc)) / 2.0)), 3)}
        st[code] = {"name": name, "lat": lat, "lon": lon, "by_year": rows}
    out["stations"] = st
    # the centred dipole, the field on the dip equator, at 77 E
    fe = {}
    for t in (1967.0, 1999.0, 2015.0, 2020.0, 2025.0):
        g, h = igrf_at(model, t)
        g10, g11, h11 = g[1][0], g[1][1], h[1][1]
        b0 = math.sqrt(g10 ** 2 + g11 ** 2 + h11 ** 2)
        plat = 90.0 - math.degrees(math.acos(-g10 / b0))
        plon = math.degrees(math.atan2(-h11, -g11))
        pl, pn, lo = math.radians(plat), math.radians(plon), math.radians(77.0)
        dipole_eq = math.degrees(math.atan(-math.cos(pl) * math.cos(lo - pn) / math.sin(pl)))  # geocentric
        lat = deq[t][77.0]
        X, Y, Z = field(g, h, lat, 77.0)
        # the pole's latitude is geocentric (90 degrees minus the tilt); the geodetic value, as on the site's guide,
        # converts it on the WGS84 ellipsoid
        plat_gd = math.degrees(math.atan(math.tan(math.radians(plat)) / (1.0 - WGS84_E2)))
        fe["%d" % t] = {"dip_equator_lat": round(lat, 4), "dipole_equator_lat_geocentric": round(dipole_eq, 4),
                        "dipole_north_pole_lat": round(plat, 3), "dipole_north_pole_lat_geodetic": round(plat_gd, 3),
                        "dipole_tilt_deg": round(90.0 - plat, 2),
                        "F_on_dip_equator_nT": round(math.sqrt(X * X + Y * Y + Z * Z), 1)}
    out["dipole_and_field_77E"] = fe

    def share(a, b):
        dd = fe[b]["dipole_equator_lat_geocentric"] - fe[a]["dipole_equator_lat_geocentric"]
        de = fe[b]["dip_equator_lat"] - fe[a]["dip_equator_lat"]
        return {"from": int(a), "to": int(b), "dipole_equator_moved_km": round(dd * KM, 1),
                "dip_equator_moved_km": round(de * KM, 1), "dipole_share": round(dd / de, 2),
                "pole_lat_deg_per_year": round((fe[b]["dipole_north_pole_lat"] - fe[a]["dipole_north_pole_lat"]) / (int(b) - int(a)), 4),
                "pole_lat_geodetic_deg_per_year": round((fe[b]["dipole_north_pole_lat_geodetic"] - fe[a]["dipole_north_pole_lat_geodetic"]) / (int(b) - int(a)), 4)}
    out["dipole_share_of_drift"] = [share("1967", "1999"), share("1999", "2015"), share("1967", "2015"), share("2015", "2025")]
    f67 = fe["1967"]["F_on_dip_equator_nT"]
    out["F_on_dip_equator_change_pct"] = {k: round(100 * (fe[k]["F_on_dip_equator_nT"] - f67) / f67, 2) for k in ("1999", "2015", "2025")}
    # WMM2025's forecast at 77 E
    epoch, gw, hw, gdw, hdw, _ = read_wmm(wmm_path)
    wm = {}
    for t in (2025.0, 2030.0):
        gg = [[gw[i][j] + (t - epoch) * gdw[i][j] for j in range(NMAX + 1)] for i in range(NMAX + 1)]
        hh = [[hw[i][j] + (t - epoch) * hdw[i][j] for j in range(NMAX + 1)] for i in range(NMAX + 1)]
        wm[t] = dip_equator_lat(gg, hh, 77.0)
    out["wmm2025_77E_km_per_year_2025_2030"] = round((wm[2030.0] - wm[2025.0]) * KM / 5.0, 2)
    # Tatuoca, Brazil
    tt = {}
    for t in (1957.0, 1967.0, 1999.0, 2012.0, 2013.0, 2014.0, 2015.0, 2025.0):
        g, h = igrf_at(model, t)
        tt["%d" % t] = round(dip_equator_lat(g, h, TATUOCA[1], lo=-30.0, hi=20.0), 3)
    out["tatuoca_dip_equator_lat"] = tt
    out["tatuoca_km_per_year_1967_2015"] = round((tt["2015"] - tt["1967"]) * KM / 48.0, 2)
    out["tatuoca_passed_observatory_between"] = [[int(t), int(t) + 1] for t in (2012.0, 2013.0, 2014.0)
                                                 if tt["%d" % t] <= TATUOCA[0] < tt["%d" % (t + 1)]]
    return out


# ------------------------------------------------------------------ observatory data
def read_kyoto(paths):
    """WDC Kyoto HAPI CSV (Time, H, X, Y, Z, version) -> {UT datetime: (H, Z)}; H from X, Y where H is not given."""
    rec = {}
    for p in paths:
        with open(p, encoding="ascii", newline="") as f:
            for row in csv.reader(f):
                if not row or not row[0][:4].isdigit():
                    continue
                t = dt.datetime(int(row[0][:4]), int(row[0][5:7]), int(row[0][8:10]), int(row[0][11:13]))
                H, X, Y, Z = (float(v) for v in row[1:5])
                if H > FILL and X < FILL and Y < FILL:
                    H = math.hypot(X, Y)
                rec[t] = (None if H > FILL else H, None if Z > FILL else Z)
    return rec


def read_bgs(path):
    """WDC Edinburgh zipped hourly CSV (H, D, Z, F, datetime at the start of the UT hour) -> {UT datetime: (H, Z)}."""
    rec = {}
    z = zipfile.ZipFile(path)
    for name in sorted(z.namelist()):
        if not name.endswith(".csv"):
            continue
        for line in z.read(name).decode("ascii", "replace").splitlines():
            if not line or line.startswith("#") or line.startswith("H,"):
                continue
            p = line.split(",")
            H, Z, s = float(p[0]), float(p[2]), p[4]
            t = dt.datetime(int(s[:4]), int(s[5:7]), int(s[8:10]), int(s[11:13]))
            rec[t] = (None if H > FILL else H, None if Z > FILL else Z)
    return rec


def lt_days(rec):
    """{LT date: [H or None] * 24}, LT = UT + 5 h."""
    days = {}
    for t, (H, _) in rec.items():
        if H is None:
            continue
        lt = t + dt.timedelta(hours=LT_OFFSET)
        days.setdefault(lt.date(), [None] * 24)[lt.hour] = H
    return days


def dH(days):
    out = {}
    for d, hs in days.items():
        if any(v is None for v in hs[0:3]):
            continue
        b = sum(hs[0:3]) / 3.0
        out[d] = [None if v is None else v - b for v in hs]
    return out


def read_gfz(path):
    f107, ap = {}, {}
    with open(path, encoding="ascii") as f:
        for line in f:
            if line.startswith("#") or not line.strip():
                continue
            p = line.split()
            d = dt.date(int(p[0]), int(p[1]), int(p[2]))
            ap[d] = int(p[23])
            if float(p[25]) > 0:
                f107[d] = float(p[25])
    return f107, ap


# ------------------------------------------------------------------ 3. observatories against IGRF-14
def obs_vs_igrf(model, recs):
    res = {}
    for code, rec in recs.items():
        acc = {}
        for t, (H, Z) in rec.items():
            if H is None or Z is None:
                continue
            a = acc.setdefault(t.year, [0, 0.0, 0.0])
            a[0] += 1
            a[1] += H
            a[2] += Z
        _, lat, lon = STATIONS[code]
        rows = {}
        for y, (n, sh, sz) in sorted(acc.items()):
            hours = 8784 if (y % 4 == 0 and (y % 100 != 0 or y % 400 == 0)) else 8760
            if n < 0.8 * hours:
                continue
            g, h = igrf_at(model, y + 0.5)
            X, Y, Z = field(g, h, lat, lon)
            rows[y] = (sh / n, sz / n, math.hypot(X, Y), Z)
        ys = sorted(rows)

        def slope(i):
            vs = [rows[y][i] for y in ys]
            mx, mv = sum(ys) / len(ys), sum(vs) / len(vs)
            return sum((x - mx) * (v - mv) for x, v in zip(ys, vs)) / sum((x - mx) ** 2 for x in ys)
        dz = [rows[y][1] - rows[y][3] for y in ys]
        res[code] = {"years": [ys[0], ys[-1]], "n_years": len(ys),
                     "Z_obs_nT_per_year": round(slope(1), 1), "Z_igrf_nT_per_year": round(slope(3), 1),
                     "H_obs_nT_per_year": round(slope(0), 2), "H_igrf_nT_per_year": round(slope(2), 2),
                     "Z_obs_minus_igrf_mean_nT": round(sum(dz) / len(dz), 0),
                     "Z_obs_minus_igrf_range_nT": [round(min(dz), 0), round(max(dz), 0)],
                     "inclination_obs_first_last_deg": [round(math.degrees(math.atan2(rows[y][1], rows[y][0])), 2)
                                                        for y in (ys[0], ys[-1])]}
    return res


# ------------------------------------------------------------------ 4. the trend re-run
def median(v):
    s = sorted(v)
    n = len(s)
    return s[n // 2] if n % 2 else 0.5 * (s[n // 2 - 1] + s[n // 2])


def monthly_medians(eq, off, start, end, ap=None, ap_max=None, scale=None):
    acc = {}
    for d, a in eq.items():
        ym = (d.year, d.month)
        if ym < start or ym > end or d not in off:
            continue
        if ap_max is not None and ap.get(d, 999) > ap_max:
            continue
        b = off[d]
        for h in HOURS:
            if a[h] is not None and b[h] is not None:
                acc.setdefault(ym, {}).setdefault(h, []).append(a[h] - b[h])
    out = {}
    for ym, byh in acc.items():
        out[ym] = {h: median(v) * (scale[h] if scale else 1.0) for h, v in byh.items() if len(v) >= 15}
    return out


def solve(A, b):
    n = len(b)
    M = [row[:] + [bb] for row, bb in zip(A, b)]
    for i in range(n):
        p = max(range(i, n), key=lambda k: abs(M[k][i]))
        M[i], M[p] = M[p], M[i]
        for k in range(n):
            if k != i:
                f = M[k][i] / M[i][i]
                M[k] = [x - f * y for x, y in zip(M[k], M[i])]
    return [M[i][n] / M[i][i] for i in range(n)]


def inverse(A):
    n = len(A)
    M = [row[:] + [1.0 if i == j else 0.0 for j in range(n)] for i, row in enumerate(A)]
    for i in range(n):
        p = max(range(i, n), key=lambda k: abs(M[k][i]))
        M[i], M[p] = M[p], M[i]
        piv = M[i][i]
        M[i] = [x / piv for x in M[i]]
        for k in range(n):
            if k != i:
                f = M[k][i]
                M[k] = [x - f * y for x, y in zip(M[k], M[i])]
    return [row[n:] for row in M]


def ols(X, y):
    k = len(X[0])
    XtX = [[sum(r[i] * r[j] for r in X) for j in range(k)] for i in range(k)]
    Xty = [sum(r[i] * yy for r, yy in zip(X, y)) for i in range(k)]
    beta = solve(XtX, Xty)
    res = [yy - sum(bi * xi for bi, xi in zip(beta, r)) for r, yy in zip(X, y)]
    s2 = sum(e * e for e in res) / (len(y) - k)
    return beta, [[s2 * c for c in row] for row in inverse(XtX)], res


def lag1(res, keys):
    """Lag-1 autocorrelation of the residuals over consecutive months only."""
    idx = {k: i for i, k in enumerate(keys)}
    m = sum(res) / len(res)
    num, pairs = 0.0, 0
    for i, (y, mo) in enumerate(keys):
        nxt = (y, mo + 1) if mo < 12 else (y + 1, 1)
        if nxt in idx:
            num += (res[i] - m) * (res[idx[nxt]] - m)
            pairs += 1
    den = sum((e - m) ** 2 for e in res)
    return (num / pairs) / (den / len(res))


def design(keys, fmon, t0, quad=False, trend=True):
    X = []
    for (y, m) in keys:
        w = 2 * math.pi * (m - 0.5) / 12.0
        F = fmon[(y, m)]
        row = [1.0, F] + ([F * F / 100.0] if quad else []) + [math.cos(w), math.sin(w), math.cos(2 * w), math.sin(2 * w)]
        if trend:
            row.append((y + (m - 0.5) / 12.0 - t0) / 10.0)
        X.append(row)
    return X


def fit_hour(med, fmon, h, quad=False, two_step=False):
    keys = sorted(k for k in med if h in med[k] and k in fmon)
    y = [med[k][h] for k in keys]
    t0 = 0.5 * (keys[0][0] + (keys[0][1] - 0.5) / 12 + keys[-1][0] + (keys[-1][1] - 0.5) / 12)
    if two_step:
        _, _, r0 = ols(design(keys, fmon, t0, quad, trend=False), y)
        beta, cov, res = ols([[1.0, (yy + (mm - 0.5) / 12.0 - t0) / 10.0] for (yy, mm) in keys], r0)
        b, se, solar = beta[1], math.sqrt(cov[1][1]), None
    else:
        beta, cov, res = ols(design(keys, fmon, t0, quad), y)
        b, se, solar = beta[-1], math.sqrt(cov[-1][-1]), round(beta[1], 3)
    r1 = lag1(res, keys)
    se_adj = se * math.sqrt((1 + r1) / (1 - r1))
    mean_y = sum(y) / len(y)
    return {"LT": h, "n_months": len(keys), "first": "%d-%02d" % keys[0], "last": "%d-%02d" % keys[-1],
            "mean_EEJ_nT": round(mean_y, 1), "trend_nT_per_decade": round(b, 2), "se_adj": round(se_adj, 2),
            "t": round(b / se_adj, 2), "passes_95pct": abs(b / se_adj) > 1.96, "lag1_r": round(r1, 3),
            "pct_of_mean_per_decade": round(100 * b / mean_y, 1), "solar_nT_per_sfu": solar}


def run(med, fmon, label, **kw):
    return {"label": label, "rows": [fit_hour(med, fmon, h, **kw) for h in HOURS]}


def overlap_ratio(trd_m, tir_m):
    scale = {}
    for h in HOURS:
        ks = [k for k in trd_m if h in trd_m[k] and k in tir_m and h in tir_m[k]]
        scale[h] = sum(trd_m[k][h] for k in ks) / sum(tir_m[k][h] for k in ks)
    return scale


# ------------------------------------------------------------------ main
def main():
    folder = os.path.abspath(sys.argv[1]) if len(sys.argv) > 1 else os.path.dirname(os.path.abspath(__file__))
    out_path = os.path.abspath(sys.argv[2]) if len(sys.argv) > 2 else os.path.join(folder, "electrojet-numbers.json")
    os.makedirs(folder, exist_ok=True)
    paths, inputs = {}, {}
    for name in INPUTS:
        paths[name], inputs[name] = find_or_fetch(folder, name)
    out = {"page": "https://factchecked.science/equatorial-electrojet-trend/", "inputs": inputs}

    print("1. field synthesis against WMM2025's test values")
    out["check_wmm2025"] = validate(paths["WMM2025COF.zip"])
    model = read_igrf(paths["igrf14coeffs.txt"])
    print("2. the dip equator in IGRF-14")
    out["dip_equator"] = dip_equator_numbers(model, paths["WMM2025COF.zip"])

    print("3. the observatories against IGRF-14")
    rec_trd = read_kyoto([paths["hour_trd_1957-1999.csv"]])
    rec_abg = read_kyoto([paths["hour_abg_1957-1990.csv"], paths["hour_abg_1990-2023.csv"]])
    rec_tir_k = read_kyoto([paths["hour_tir_1999-2011.csv"]])
    rec_tir = read_bgs(paths["bgs-hour-tir-1999-2014-hdzf-csv.zip"])
    out["observatories_vs_igrf14"] = obs_vs_igrf(model, {"TRD": rec_trd, "TIR": rec_tir_k, "ABG": rec_abg})
    common = [t for t in rec_tir_k if t in rec_tir and rec_tir_k[t][0] is not None and rec_tir[t][0] is not None]
    out["tir_kyoto_vs_edinburgh"] = {"hours_in_both": len(common),
                                     "largest_H_difference_nT": max(abs(rec_tir_k[t][0] - rec_tir[t][0]) for t in common)}

    print("4. the electrojet's trend by local time")
    f107, ap = read_gfz(paths["Kp_ap_Ap_SN_F107_since_1932.txt"])
    acc = {}
    for d, v in f107.items():
        acc.setdefault((d.year, d.month), []).append(v)
    fmon = {k: sum(v) / len(v) for k, v in acc.items()}
    trd, abg, tir = dH(lt_days(rec_trd)), dH(lt_days(rec_abg)), dH(lt_days(rec_tir))
    blocks = {}
    mA = monthly_medians(trd, abg, (1967, 1), (1999, 9))
    blocks["A"] = run(mA, fmon, "TRD-ABG, 1967-01 to 1999-09")
    blocks["A_F2"] = run(mA, fmon, "A with F and F^2", quad=True)
    blocks["A_two_step"] = run(mA, fmon, "A in two steps (season and F removed, then the residual's trend)", two_step=True)
    blocks["A_quiet"] = run(monthly_medians(trd, abg, (1967, 1), (1999, 9), ap=ap, ap_max=20), fmon, "A without days of Ap > 20")
    scale = overlap_ratio(monthly_medians(trd, abg, (1999, 1), (1999, 9)), monthly_medians(tir, abg, (1999, 1), (1999, 9)))
    out["tir_to_trd_ratio_1999_overlap"] = {str(h): round(v, 3) for h, v in scale.items()}

    def spliced(ap_max=None, sc=scale, end=(2014, 12)):
        m = dict(monthly_medians(trd, abg, (1967, 1), (1998, 12), ap=ap, ap_max=ap_max))
        m.update(monthly_medians(tir, abg, (1999, 1), end, ap=ap, ap_max=ap_max, scale=sc))
        return m
    mB = spliced()
    blocks["B"] = run(mB, fmon, "TRD-ABG 1967-01 to 1998-12, then TIR-ABG (scaled to TRD) 1999-01 to 2014-12")
    blocks["B_unscaled"] = run(spliced(sc=None), fmon, "B with TIR not scaled")
    blocks["B_to_2011"] = run(spliced(end=(2011, 12)), fmon, "B ending 2011-12")
    blocks["B_quiet"] = run(spliced(ap_max=20), fmon, "B without days of Ap > 20")
    blocks["B_F2"] = run(mB, fmon, "B with F and F^2", quad=True)
    mC = monthly_medians(tir, abg, (1999, 1), (2014, 12))
    blocks["C"] = run(mC, fmon, "TIR-ABG alone, 1999-01 to 2014-12")
    out["trends"] = blocks
    # the solar cycle's swing at 11 LT beside the trend, 1967-1999: the range of the 13-month running mean of the
    # monthly F10.7 over the months of period A, times the fitted nT per sfu, against the trend over the same span
    ks = sorted(k for k in mA if k in fmon)
    run13 = [sum(fmon[x] for x in ks[i - 6:i + 7]) / 13.0 for i in range(6, len(ks) - 6)]
    a11 = next(r for r in blocks["A"]["rows"] if r["LT"] == 11)
    span = (ks[-1][0] + (ks[-1][1] - 0.5) / 12.0) - (ks[0][0] + (ks[0][1] - 0.5) / 12.0)
    out["solar_swing_11LT_A"] = {"F107_13_month_mean_range_sfu": round(max(run13) - min(run13), 1),
                                 "solar_nT_per_sfu": a11["solar_nT_per_sfu"],
                                 "swing_nT": round((max(run13) - min(run13)) * a11["solar_nT_per_sfu"], 1),
                                 "years": round(span, 2),
                                 "trend_over_the_years_nT": round(a11["trend_nT_per_decade"] * span / 10.0, 1),
                                 "mean_EEJ_nT": a11["mean_EEJ_nT"]}

    print("5. the fall-off of the electrojet's signal with distance from the dip equator")
    dist_cache = {}

    def dist_km(code, t):
        key = (code, round(t, 4))
        if key not in dist_cache:
            _, lat, lon = STATIONS[code]
            g, h = igrf_at(model, t)
            dist_cache[key] = (lat - dip_equator_lat(g, h, lon)) * KM
        return dist_cache[key]
    pts = []
    for code, fn, a, b, tmid in (("KOD", "hour_kod_1970-1980.csv", (1980, 1), (1980, 6), 1980.25),
                                 ("ANN", "hour_ann_1970-1980.csv", (1970, 1), (1980, 12), 1975.5),
                                 ("ETT", "hour_ett_1996-1998.csv", (1996, 1), (1998, 4), 1997.2)):
        st = dH(lt_days(read_kyoto([paths[fn]])))
        ms, mt = monthly_medians(st, abg, a, b), monthly_medians(trd, abg, a, b)
        num = den = 0.0
        shared = set()
        for k in ms:
            for h in (10, 11, 12):
                if k in mt and h in ms[k] and h in mt[k]:
                    num += ms[k][h]
                    den += mt[k][h]
                    shared.add(k)
        pts.append({"station": code, "months": "%d-%02d to %d-%02d" % (a + b), "months_shared": len(shared),
                    "years_shared": sorted({y for y, _ in shared}), "ratio_to_TRD_10_12LT": round(num / den, 3),
                    "km_from_dip_equator": round(dist_km(code, tmid), 1), "TRD_km_from_dip_equator": round(dist_km("TRD", tmid), 1)})
    u = [(p["km_from_dip_equator"] ** 2 - p["TRD_km_from_dip_equator"] ** 2, math.log(p["ratio_to_TRD_10_12LT"])) for p in pts]
    k = -sum(a * b for a, b in u) / sum(a * a for a, _ in u)
    s_fit = math.sqrt(math.log(2) / k)
    per_point = {p["station"]: round(math.sqrt(-math.log(2) * a / b), 0) for p, (a, b) in zip(pts, u)}

    def change(x0, x1, s):
        return round(100 * (math.exp(-math.log(2) * ((x1 / s) ** 2 - (x0 / s) ** 2)) - 1), 1)
    trd0, trd1 = dist_km("TRD", 1967.0), dist_km("TRD", 1999.7)
    tir0, tir1 = dist_km("TIR", 1999.0), dist_km("TIR", 2015.0)
    out["profile"] = {"points": pts, "half_width_fitted_km": round(s_fit, 0), "half_width_each_point_km": per_point,
                      "TRD_km_1967_and_1999.7": [round(trd0, 1), round(trd1, 1)],
                      "TRD_signal_change_1967_1999_pct": change(trd0, trd1, s_fit),
                      "TIR_km_1999_and_2015": [round(tir0, 1), round(tir1, 1)],
                      "TIR_signal_change_1999_2015_pct": change(tir0, tir1, s_fit),
                      "TIR_signal_change_1999_2015_pct_for_half_widths": {
                          "%d" % w: change(tir0, tir1, w) for w in (min(per_point.values()), max(per_point.values()), 250, 400)}}

    print("6. the trends with the drift taken out")

    def corrected(med, code):
        out_m = {}
        for (y, m), byh in med.items():
            x = dist_km(code, y + (m - 0.5) / 12.0)
            fac = math.exp(-math.log(2) * (x / s_fit) ** 2)
            out_m[(y, m)] = {h: v / fac for h, v in byh.items()}
        return out_m
    gA = corrected(mA, "TRD")
    gscale = overlap_ratio(corrected(monthly_medians(trd, abg, (1999, 1), (1999, 9)), "TRD"),
                           corrected(monthly_medians(tir, abg, (1999, 1), (1999, 9)), "TIR"))
    gB = corrected(monthly_medians(trd, abg, (1967, 1), (1998, 12)), "TRD")
    for key, v in corrected(mC, "TIR").items():
        gB[key] = {h: x * gscale[h] for h, x in v.items()}
    out["tir_to_trd_ratio_1999_overlap_drift_removed"] = {str(h): round(v, 3) for h, v in gscale.items()}
    out["trends_drift_removed"] = {"A": run(gA, fmon, "A with the drift removed"),
                                   "B": run(gB, fmon, "B with the drift removed"),
                                   "C": run(corrected(mC, "TIR"), fmon, "C with the drift removed")}

    with open(out_path, "w", encoding="utf-8", newline="\n") as f:
        json.dump(out, f, indent=1, ensure_ascii=False)
        f.write("\n")

    d = out["dip_equator"]
    print("\nWMM2025 test values:", out["check_wmm2025"])
    for r in d["drift_77E"]:
        print("dip equator at 77 E, %(from)d-%(to)d: %(lat_start).2f to %(lat_end).2f N, %(moved_km)+.1f km, %(km_per_year)+.2f km a year" % r)
    print("passed south of Kanyakumari in", d["passed_south_of_kanyakumari"], "| dipole share:", [r["dipole_share"] for r in d["dipole_share_of_drift"]])
    print("observatories:", json.dumps(out["observatories_vs_igrf14"]))
    for key in ("A", "B", "C"):
        print(key, " ".join("%d:%+.2f%s" % (r["LT"], r["trend_nT_per_decade"], "*" if r["passes_95pct"] else "")
                            for r in blocks[key]["rows"]))
    print("profile: half-width %d km; TIR %.0f%% over 1999-2015" % (out["profile"]["half_width_fitted_km"],
                                                                   out["profile"]["TIR_signal_change_1999_2015_pct"]))
    print("wrote", out_path)


if __name__ == "__main__":
    sys.stdout.reconfigure(encoding="utf-8")
    main()