#!/usr/bin/env python3 """The South Atlantic Anomaly's stronger ground response: every number of ours on Fact-Checked Science's page https://factchecked.science/sama-sudden-commencement/ , recomputed from public data. Standard library only. Run it python3 -I sama_numbers.py [data-folder [output-folder]] Inputs are read from the data folder (default: ./sama_numbers_data) and downloaded there first if missing. They are read as data, parsed as numbers and never executed. A full run needs about 4,800 files and 490 MB, which took about an hour to fetch here, and then takes a few minutes. Observatories revise their files as data become definitive, and EMBRACE's archive grows, so a newer download may differ slightly from the numbers on the page. Outputs, written to the output folder (default: the folder holding the script), derived numbers only: sama-numbers.json every result below sama-event-metrics.csv A, V, Pc5, PA and the step changes in H and Z, per station and event No magnetometer series is written out and none is hosted with the page: EMBRACE's conditions bar passing on its data to third parties, and INTERMAGNET asks for written permission for bulk distribution. The tests. Tests 1-9 follow test-rules.md, written before any event other than 10 May 2024 was computed and hosted with this script. run_added() holds the checks added after their results were read, and the page marks them so: sections A-G as first built, H-J after the fresh review of 1 November 2026. That review also found a bug in this script, now fixed: a fill value in Z alone used to drop a minute whose H was good. Rule 6's last bullet, SJG against the other Northern band stations, was missing until the same review. Reproduction Tables 4, 5 and 6 of Piassi et al. (2026), 10.22541/essoar.15011353/v1, on raw H, with T = 17:04 UT. The preprint does not state its window, baseline, wavelet or filter. These settings fit its tables, and were frozen before any other event was computed: A = largest H from T to T+12 min minus the mean of the 60 minutes before T; V = largest one-minute |dH| from T to T+12; Pc5 = H over T-120 to T+120 min, second-order Butterworth band-pass at 1.6-7 mHz run forward and backward, the maximum minus the lowest value before it; PA = the sum of one-minute |dH| over T-119 to T+120. R = South / North. Tests 1-5 Same pair over other events; every pair by local time and by season; every pair fitted to the same local time and day of the year; the trend. Test 6 Near-noon amplitude and speed, relative to event size, at every station with a dip latitude of 12-30. Test 7 Quiet days: PA and the day's range at SMS and SJG, with no sudden commencement. Test 8 SJG's one-second file on 10 May resampled to one minute five ways, to bound filtering (with the reproduction). Test 9 The step ratio dZ/dH at each station, a rough pointer to coastal and ground effects. Also The South Atlantic Anomaly's minimum and area from IGRF-14, the field at the stations, and field-line tracing, after a check of our synthesis against NOAA's WMM2025 test values. Inputs and their terms EMBRACE (INPE): one-minute files, https://embracedata.inpe.br/magnetometer/ . The programme's conditions: its mandatory acknowledgement, and no redistribution to third parties. If you publish results from these data, include the statement "The author(s) thank(s) the Embrace/INPE Program of the MCTI for providing the data for this work". INTERMAGNET, through BGS's GIN, https://imag-data.bgs.ac.uk/GIN_V1/ : CC BY-NC 4.0 unless noted, see https://intermagnet.org/data_conditions.html . Definitive data to 2020: https://doi.org/10.5880/INTERMAGNET.1991.2020 . BCMT's observatories (here IPM, PPT, PHU, KOU, CLF): https://doi.org/10.18715/BCMT.MAG.DEF , CC BY-NC 4.0 in its DataCite record. Observatori de l'Ebre: annual lists of sudden commencements, https://www.obsebre.es/en/variations/rapid , for ISGI (https://isgi.unistra.fr). Used to choose events and, through the mean of their five listed amplitudes, to size each one. IGRF-14 and WMM2025 (NOAA, public). GFZ's daily Ap (CC BY 4.0, doi:10.5880/Kp.0001), to choose quiet days. """ import concurrent.futures as cf import csv import datetime as dt import glob import hashlib import json import math import os import random import re import statistics as st import sys import threading import time import urllib.error import urllib.request import zipfile # ------------------------------------------------------------------ field synthesis and field lines A_REF = 6371.2 WGS84_A = 6378.137 WGS84_F = 1.0 / 298.257223563 WGS84_E2 = WGS84_F * (2.0 - WGS84_F) NMAX = 13 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 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 geocentric_to_geodetic(x, y, z): """Geodetic latitude (deg), longitude (deg), height (km) from geocentric Cartesian km.""" lon = math.degrees(math.atan2(y, x)) p = math.hypot(x, y) lat = math.atan2(z, p * (1 - WGS84_E2)) for _ in range(8): s = math.sin(lat) N = WGS84_A / math.sqrt(1 - WGS84_E2 * s * s) h = p / math.cos(lat) - N lat = math.atan2(z, p * (1 - WGS84_E2 * N / (N + h))) s = math.sin(lat) N = WGS84_A / math.sqrt(1 - WGS84_E2 * s * s) return math.degrees(lat), lon, p / math.cos(lat) - N 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_sph(g, h, r, theta, lam): """B_r, B_theta, B_phi (nT) at geocentric radius r (km), colatitude theta and longitude lam (radians).""" st = math.sin(theta) P, D = legendre(theta) 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] return br, bt, bp / st 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) br, bt, bp = field_sph(g, h, r, math.pi / 2.0 - phic, math.radians(lon)) xc, zc = -bt, -br cd, sd = math.cos(delta), math.sin(delta) return xc * cd + zc * sd, bp, zc * cd - xc * sd def elements(g, h, lat, lon, hkm=0.0): X, Y, Z = field(g, h, lat, lon, hkm) H = math.hypot(X, Y) I = math.degrees(math.atan2(Z, H)) D = math.degrees(math.atan2(Y, X)) F = math.sqrt(X * X + Y * Y + Z * Z) return {"X": X, "Y": Y, "Z": Z, "H": H, "D": D, "I": I, "F": F, "dip_lat": math.degrees(math.atan(math.tan(math.radians(I)) / 2.0))} def validate_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]) worst = {k: 0.0 for k in ("X", "Y", "Z", "F")} count = 0 for line in tests.splitlines(): if line.startswith("#") or not line.strip(): continue q = [float(x) for x in line.split()] t, hkm, lat, lon = q[0], q[1], q[2], q[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, q[4]), ("Y", Y, q[5]), ("Z", Z, q[6]), ("F", F, q[8])): worst[key] = max(worst[key], abs(mine - ref)) count += 1 return {"test_points": count, "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())} # ------------------------------------------------------------------ field-line tracing def bcart(g, h, x, y, z): r = math.sqrt(x * x + y * y + z * z) theta = math.acos(z / r) lam = math.atan2(y, x) br, bt, bp = field_sph(g, h, r, theta, lam) st, ct = math.sin(theta), math.cos(theta) cl, sl = math.cos(lam), math.sin(lam) bx = br * st * cl + bt * ct * cl - bp * sl by = br * st * sl + bt * ct * sl + bp * cl bz = br * ct - bt * st return bx, by, bz def trace(g, h, lat, lon, hkm=0.0, end_km=0.0, step_km=10.0, max_steps=40000): """Trace from (lat, lon, hkm) to the other footprint at geodetic height end_km. Returns the conjugate point's latitude and longitude, the apex radius in units of A_REF and the apex latitude.""" r, phic, _ = geodetic_to_geocentric(lat, hkm) lam = math.radians(lon) x, y, z = r * math.cos(phic) * math.cos(lam), r * math.cos(phic) * math.sin(lam), r * math.sin(phic) bx, by, bz = bcart(g, h, x, y, z) sign = 1.0 if (bx * x + by * y + bz * z) > 0 else -1.0 # head upward along the field line ds = step_km * sign rmax, n, apexlat = r, 0, lat def f(px, py, pz): bbx, bby, bbz = bcart(g, h, px, py, pz) b = math.sqrt(bbx * bbx + bby * bby + bbz * bbz) return bbx / b, bby / b, bbz / b while n < max_steps: k1 = f(x, y, z) k2 = f(x + 0.5 * ds * k1[0], y + 0.5 * ds * k1[1], z + 0.5 * ds * k1[2]) k3 = f(x + 0.5 * ds * k2[0], y + 0.5 * ds * k2[1], z + 0.5 * ds * k2[2]) k4 = f(x + ds * k3[0], y + ds * k3[1], z + ds * k3[2]) prev = (x, y, z) x += ds / 6.0 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) y += ds / 6.0 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) z += ds / 6.0 * (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2]) n += 1 rr = math.sqrt(x * x + y * y + z * z) if rr > rmax: rmax, apexlat = rr, math.degrees(math.asin(z / rr)) la, lo, hh = geocentric_to_geodetic(x, y, z) if n > 3 and hh < end_km: _, _, ph = geocentric_to_geodetic(*prev) w = (ph - end_km) / (ph - hh) la, lo, _ = geocentric_to_geodetic(prev[0] + w * (x - prev[0]), prev[1] + w * (y - prev[1]), prev[2] + w * (z - prev[2])) return {"lat": la, "lon": lo, "apex_re": rmax / A_REF, "apex_lat": apexlat, "steps": n} raise RuntimeError("field line did not return") def gc_km(lat1, lon1, lat2, lon2): a1, b1, a2, b2 = map(math.radians, (lat1, lon1, lat2, lon2)) s = math.sin((a2 - a1) / 2) ** 2 + math.cos(a1) * math.cos(a2) * math.sin((b2 - b1) / 2) ** 2 return 2 * 6371.2 * math.asin(min(1.0, math.sqrt(s))) def dipole_pole(g, h): g10, g11, h11 = g[1][0], g[1][1], h[1][1] b0 = math.sqrt(g10 ** 2 + g11 ** 2 + h11 ** 2) return 90 - math.degrees(math.acos(-g10 / b0)), math.degrees(math.atan2(-h11, -g11)), b0 def dipole_lat(g, h, lat, lon): plat, plon, _ = dipole_pole(g, h) r, phic, _ = geodetic_to_geocentric(lat, 0.0) a, b = phic, math.radians(lon) pa, pb = math.radians(plat), math.radians(plon) return math.degrees(math.asin(math.sin(a) * math.sin(pa) + math.cos(a) * math.cos(pa) * math.cos(b - pb))) # ------------------------------------------------------------------ fetching and reading the files UA = "Mozilla/5.0 (factchecked.science verification; non-commercial fact check)" GIN = ("https://imag-data.bgs.ac.uk/GIN_V1/GINServices?Request=GetData&format=IAGA2002&testObsys=0" "&observatoryIagaCode={code}&samplesPerDay={rate}&publicationState={state}" "&dataStartDate={date}&dataDuration=1&orientation=native&recordTermination=UNIX") EMBRACE = "https://embracedata.inpe.br/magnetometer/{code}/{year}/{lcode}{dd}{mon}.{yy}m" MONTHS = ["jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec"] EMBRACE_COORDS = { # EMBRACE's readme (https://embracedata.inpe.br/magnetometer/readme_magnetometer.html) "ALF": (-9.87033, -56.10403), "ARA": (-5.65130, -48.07352), "CAR": (-7.37916, -36.52972), "CBA": (-15.55472, -56.06945), "CHI": (-36.64122, -71.99485), "CXP": (-22.70207, -45.01443), "EUS": (-3.87998, -38.42467), "JAT": (-17.93185, -51.71840), "MED": (-25.29967, -54.11451), "MAN": (-3.10803, -59.97465), "PAL": (-10.17889, -48.35889), "PVE": (-8.8350, -63.94000), "RGA": (-53.78572, -67.75145), "SJC": (-23.20862, -45.96353), "SLZ": (-2.59417, -44.20972), "STM": (-2.68630, -54.53920), "SMS": (-29.44357, -53.82273), "TCM": (-26.79612, -65.31555), } FILL = 88000.0 _lock = threading.Lock() def log_fetch(root, relpath, url, data, status): line = "\t".join([relpath.replace("\\", "/"), url, dt.datetime.now(dt.timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ"), str(len(data)), hashlib.sha256(data).hexdigest(), str(status)]) with _lock: with open(os.path.join(root, "fetch-log.tsv"), "a", encoding="utf-8", newline="\n") as f: f.write(line + "\n") def fetch(root, folder, name, url, pause=0.1): """Return the bytes of url, cached at //. A 400 or 404 is cached as a '.missing' marker.""" d = os.path.join(root, folder) os.makedirs(d, exist_ok=True) path = os.path.join(d, name) missing = path + ".missing" if os.path.exists(path): with open(path, "rb") as f: return f.read() if os.path.exists(missing): return None req = urllib.request.Request(url, headers={"User-Agent": UA}) err = None for attempt in range(4): try: with urllib.request.urlopen(req, timeout=120) as r: data = r.read() time.sleep(pause) tmp = path + ".tmp" with open(tmp, "wb") as f: f.write(data) os.replace(tmp, path) log_fetch(root, folder + "/" + name, url, data, 200) return data except urllib.error.HTTPError as e: if e.code in (400, 404): with open(missing, "w", encoding="utf-8", newline="\n") as f: f.write(f"HTTP {e.code} on {url}\n") log_fetch(root, folder + "/" + name + ".missing", url, b"", e.code) time.sleep(pause) return None err = e except Exception as e: # network hiccup: retry politely err = e time.sleep(5 * (attempt + 1)) # a file the server keeps refusing (HTTP 5xx) counts as unavailable for that station and day, and is logged if isinstance(err, urllib.error.HTTPError) and err.code >= 500: with open(missing, "w", encoding="utf-8", newline="\n") as f: f.write(f"HTTP {err.code} on {url}, after 4 tries\n") log_fetch(root, folder + "/" + name + ".missing", url, b"", err.code) return None raise RuntimeError(f"failed {url}: {err}") def embrace_day(root, code, day): """{minute of day: (H, D_deg, Z)} from EMBRACE's daily file; None if the file does not exist.""" url = EMBRACE.format(code=code, year=day.year, lcode=code.lower(), dd=f"{day.day:02d}", mon=MONTHS[day.month - 1], yy=f"{day.year % 100:02d}") name = f"{code.lower()}{day.day:02d}{MONTHS[day.month - 1]}.{day.year % 100:02d}m" data = fetch(root, f"embrace-magnetometer/{code}/{day.year}", name, url) if data is None: return None out = {} for line in data.decode("latin-1").splitlines(): p = line.split() if len(p) < 10 or not p[0].isdigit(): continue try: if (int(p[0]), int(p[1]), int(p[2])) != (day.day, day.month, day.year): continue minute = int(p[3]) * 60 + int(p[4]) d, h, z = float(p[5]), float(p[6]), float(p[7]) except ValueError: continue if h <= 0 or h > FILL: continue out[minute] = (h, d, z if abs(z) <= FILL else None) return out def gin_day(root, code, day, state="adj-or-rep", rate="minute"): """{minute of day: (H, D_deg, Z)} and header dict from BGS GIN IAGA-2002; (None, None) if unavailable.""" url = GIN.format(code=code, rate=rate, state=state, date=day.isoformat()) data = fetch(root, f"gin-intermagnet-{rate}/{code}/{day.year}", f"{code.lower()}{day.isoformat()}{state}.txt", url) if data is None: return None, None out, hdr, cols = {}, {}, None for line in data.decode("latin-1").splitlines(): if line.startswith("DATE"): cols = line.split() continue if cols is None: if "|" in line: body = line.rstrip("|").rstrip() key, val = body[:24].strip(), body[24:].strip() if key: hdr[key] = val continue p = line.split() if len(p) < 7: continue hh, mm = p[1].split(":")[:2] minute = int(hh) * 60 + int(mm) try: a, b, c = float(p[3]), float(p[4]), float(p[5]) except ValueError: continue kinds = [x[-1] for x in cols[3:7]] vals = dict(zip(kinds, [a, b, c])) # A minute is dropped only if H itself can't be had: a fill in Z, or in D, leaves H good. if "X" in vals and "Y" in vals: if abs(vals["X"]) > FILL or abs(vals["Y"]) > FILL: continue h = math.hypot(vals["X"], vals["Y"]) d = math.degrees(math.atan2(vals["Y"], vals["X"])) elif "H" in vals and "D" in vals: if abs(vals["H"]) > FILL: continue h = vals["H"] d = vals["D"] / 60.0 if abs(vals["D"]) <= FILL else None # IAGA-2002 D is in minutes of arc else: continue z = vals.get("Z") out[minute] = (h, d, z if z is not None and abs(z) <= FILL else None) return out, hdr def series(root, code, t0, t1, source=None): """{datetime: (H, D, Z)} between t0 and t1 inclusive, one-minute, from whichever days the span touches.""" out = {} day = t0.date() while day <= t1.date(): if source == "embrace" or (source is None and code in EMBRACE_COORDS): rows = embrace_day(root, code, day) else: rows, _ = gin_day(root, code, day) if rows: for m, v in rows.items(): t = dt.datetime(day.year, day.month, day.day) + dt.timedelta(minutes=m) if t0 <= t <= t1: out[t] = v day += dt.timedelta(days=1) return out # ------------------------------------------------------------------ the sudden-commencement metrics PAPER_PAIRS = [("SMS", "SJG"), ("ASC", "GUI"), ("HBK", "PAG"), ("GAN", "ABG"), ("KDU", "KNY")] PAPER_STATIONS = [c for p in PAPER_PAIRS for c in p] EPOCH_2024 = 2024.357 def read_events(folder): """Every SC in Ebre's annual lists from 2013: dict with time (floored minute), size, type, list file.""" out = [] for f in sorted(glob.glob(os.path.join(folder, "ssc_*.txt"))): if int(os.path.basename(f)[4:8]) < 2013: # SMS's public record begins in 2013; older lists have other formats continue for line in open(f, encoding="latin-1"): m = re.match(r"\s*(\d{4}) (\d\d) (\d\d) (\d\d)\s+(\d\d(?:\.\d)?)\s+(.*)$", line) if not m: continue rest = m.group(6).split() # rise times (5), amplitudes (5), qualifications (5; a '-' can stand for one), observatories (5), type typ = rest[-1] obs = rest[-6:-1] amps = [float(x) for x in rest[5:10]] t = dt.datetime(int(m.group(1)), int(m.group(2)), int(m.group(3)), int(m.group(4)), int(float(m.group(5)))) out.append({"t": t, "minute_text": m.group(5), "amps": amps, "obs": obs, "type": typ, "size": st.mean(amps), "file": os.path.basename(f)}) out.sort(key=lambda e: e["t"]) return out def read_capabilities(path): txt = open(path, encoding="utf-8").read() obs = {} for m in re.finditer(r"(\w+)\s*(.*?)\s*([-\d.]+)\s*" r"([-\d.]+)", txt): lon = float(m.group(4)) obs[m.group(1)] = {"name": m.group(2), "lat": float(m.group(3)), "lon": lon if lon <= 180 else lon - 360} return obs def band_stations(obs, model): """Station code -> (lat, lon, dip latitude, F, network) for 12 <= |dip| <= 30 in IGRF-14 at 2024.36.""" g, h = igrf_at(model, EPOCH_2024) out = {} for code, o in obs.items(): e = elements(g, h, o["lat"], o["lon"]) if 12.0 <= abs(e["dip_lat"]) <= 30.0: out[code] = (o["lat"], o["lon"], e["dip_lat"], e["F"], "INTERMAGNET", o["name"]) for code, (la, lo) in EMBRACE_COORDS.items(): e = elements(g, h, la, lo) if 12.0 <= abs(e["dip_lat"]) <= 30.0: out[code] = (la, lo, e["dip_lat"], e["F"], "EMBRACE", code) return out # ------------------------------------------------------------------ the filter (test-rules.md) def butter_bp_sections(f1, f2, fs, order=2): """Digital Butterworth band-pass as cascaded biquads, by the bilinear transform of the analogue prototype.""" import cmath w1 = 2 * fs * math.tan(math.pi * f1 / fs) w2 = 2 * fs * math.tan(math.pi * f2 / fs) w0 = math.sqrt(w1 * w2) bw = w2 - w1 secs = [] for k in range(order): theta = math.pi * (2 * k + 1) / (2 * order) + math.pi / 2 p = cmath.exp(1j * theta) disc = cmath.sqrt((p * bw) ** 2 - 4 * w0 ** 2) for s in ((p * bw + disc) / 2, (p * bw - disc) / 2): if s.imag < 0: continue a1, a0 = -2 * s.real, abs(s) ** 2 K = 2 * fs A0 = K * K + a1 * K + a0 secs.append((bw * K / A0, 0.0, -bw * K / A0, (-2 * K * K + 2 * a0) / A0, (K * K - a1 * K + a0) / A0)) return secs SECS = butter_bp_sections(1.6e-3, 7e-3, 1 / 60.0, 2) def sos(x, secs): y = list(x) for b0, b1, b2, a1, a2 in secs: z1 = z2 = 0.0 out = [] for v in y: o = b0 * v + z1 z1 = b1 * v - a1 * o + z2 z2 = b2 * v - a2 * o out.append(o) y = out return y def filtfilt(x, pad=60): pre = [2 * x[0] - x[i] for i in range(pad, 0, -1)] post = [2 * x[-1] - x[-i - 1] for i in range(1, pad + 1)] y = sos(pre + list(x) + post, SECS) y = sos(y[::-1], SECS)[::-1] return y[pad:pad + len(x)] def metrics(s, T): """A, V, Pc5, PA and step ratios from {datetime: (H, D, Z)} around onset T, or a reason for exclusion.""" at = lambda k: s.get(T + dt.timedelta(minutes=k)) # noqa: E731 H = lambda k: (at(k)[0] if at(k) else None) # noqa: E731 win = [H(k) for k in range(-120, 121)] n241 = sum(v is not None for v in win) base = [H(k) for k in range(-60, 0) if H(k) is not None] w = [H(k) for k in range(0, 13) if H(k) is not None] row = {"n241": n241, "nbase": len(base), "nw": len(w)} if n241 < 0.95 * 241: return {**row, "excl": "n241<95%"} if len(base) < 54: return {**row, "excl": "baseline<90%"} if len(w) < 12: return {**row, "excl": "window<12/13"} diffs = [abs(H(k) - H(k - 1)) for k in range(-119, 121) if H(k) is not None and H(k - 1) is not None] if diffs and max(diffs) > 500: return {**row, "excl": "jump>500nT"} run = longest = 1 for k in range(-119, 121): if H(k) is not None and H(k - 1) is not None and H(k) == H(k - 1): run += 1 longest = max(longest, run) else: run = 1 if longest >= 15: return {**row, "excl": "frozen run>=15"} mean_base = sum(base) / len(base) A = max(w) - mean_base vs = [abs(H(k) - H(k - 1)) for k in range(0, 13) if H(k) is not None and H(k - 1) is not None] V = max(vs) if vs else None PA = sum(diffs) # gaps up to 5 minutes linearly interpolated for the band-pass xs = list(win) i = 0 ok = True while i < 241: if xs[i] is None: j = i while j < 241 and xs[j] is None: j += 1 if j - i > 5 or i == 0 or j == 241: if i == 0 or j == 241: fill = xs[j] if i == 0 and j < 241 else xs[i - 1] if j - i > 5: ok = False for k in range(i, j): xs[k] = fill else: ok = False for k in range(i, j): xs[k] = xs[i - 1] + (xs[j] - xs[i - 1]) * (k - i + 1) / (j - i + 1) else: for k in range(i, j): xs[k] = xs[i - 1] + (xs[j] - xs[i - 1]) * (k - i + 1) / (j - i + 1) i = j else: i += 1 Pc5 = None if ok: mean = sum(xs) / len(xs) y = filtfilt([v - mean for v in xs]) imax = max(range(len(y)), key=lambda q: y[q]) if imax > 0: Pc5 = y[imax] - min(y[:imax]) def step(idx): pre = [at(k)[idx] for k in range(-5, 0) if at(k) and at(k)[idx] is not None] post = [at(k)[idx] for k in range(3, 8) if at(k) and at(k)[idx] is not None] if len(pre) < 4 or len(post) < 4: return None return sum(post) / len(post) - sum(pre) / len(pre) dH, dZ = step(0), step(2) return {**row, "excl": "", "A": A, "V": V, "Pc5": Pc5, "PA": PA, "dH": dH, "dZ": dZ, "base": mean_base} def one(root, code, ev, network): T = ev["t"] s = series(root, code, T - dt.timedelta(minutes=120), T + dt.timedelta(minutes=120), source="embrace" if network == "EMBRACE" else "gin") if not s: return code, ev, {"n241": 0, "nbase": 0, "nw": 0, "excl": "no data"} return code, ev, metrics(s, T) def run_metrics(root, outcsv, workers=3, min_size=15.0): model = read_igrf(os.path.join(root, "igrf14", "igrf14coeffs.txt")) obs = read_capabilities(os.path.join(root, "gin-probe", "gin-capabilities.xml")) band = band_stations(obs, model) print("band stations:", len(band), sorted(band), flush=True) events = [e for e in read_events(os.path.join(root, "ebre-sc-lists")) if e["type"] in ("SSC", "SI")] print("positive events:", len(events), flush=True) sel = [e for e in events if e["size"] >= min_size] print("events with size >=", min_size, ":", len(sel), flush=True) stations = dict(band) for code in PAPER_STATIONS: if code not in stations: if code in EMBRACE_COORDS: la, lo = EMBRACE_COORDS[code] stations[code] = (la, lo, None, None, "EMBRACE", code) else: o = obs[code] stations[code] = (o["lat"], o["lon"], None, None, "INTERMAGNET", o["name"]) tasks = [] for e in sel: for code, info in stations.items(): # the band stations beyond the paper's are run only on the primary set (size >= 20 nT) if code not in PAPER_STATIONS and e["size"] < 20.0: continue tasks.append((code, e, info[4])) print("station-event tasks:", len(tasks), flush=True) rows = [] done = 0 fields = ["time", "code", "type", "size", "list", "excl", "n241", "nbase", "nw", "A", "V", "Pc5", "PA", "dH", "dZ"] with cf.ThreadPoolExecutor(max_workers=workers) as ex: futs = [ex.submit(one, root, c, e, n) for c, e, n in tasks] for fu in cf.as_completed(futs): code, ev, r = fu.result() rows.append({"time": ev["t"].strftime("%Y-%m-%dT%H:%M"), "code": code, "type": ev["type"], "size": round(ev["size"], 2), "list": ev["file"], **{k: (round(v, 3) if isinstance(v, float) else v) for k, v in r.items() if k != "base"}}) done += 1 if done % 200 == 0: print(done, "of", len(tasks), flush=True) rows.sort(key=lambda r: (r["time"], r["code"])) with open(outcsv, "w", encoding="utf-8", newline="\n") as f: w = csv.DictWriter(f, fieldnames=fields, lineterminator="\n", restval="") w.writeheader() for r in rows: w.writerow(r) print("wrote", outcsv, len(rows), flush=True) # ------------------------------------------------------------------ tests 1-7 and 9 of test-rules.md PAIRS = [("SMS", "SJG"), ("ASC", "GUI"), ("HBK", "PAG"), ("GAN", "ABG"), ("KDU", "KNY")] METRICS = ("A", "V", "Pc5", "PA") EVENT = "2024-05-10T17:05" SEED = 20261101 NBOOT = 2000 LT0 = 13.1 DOY0 = 131 NEIGHBOURS = ("PIL", "VSS", "JAT", "SJC", "CXP", "CHI", "MED", "TCM") def num(x): return float(x) if x not in ("", None) else None def load(csvpath): rows = {} excl = {} for r in csv.DictReader(open(csvpath, encoding="utf-8")): t = r["time"] code = r["code"] if r["excl"]: excl.setdefault(code, {}).setdefault(r["excl"], 0) excl[code][r["excl"]] += 1 rows.setdefault(t, {"size": float(r["size"]), "type": r["type"], "stations": {}}) rows[t]["stations"][code] = {"excl": r["excl"], "A": num(r["A"]), "V": num(r["V"]), "Pc5": num(r["Pc5"]), "PA": num(r["PA"]), "dH": num(r["dH"]), "dZ": num(r["dZ"])} return rows, excl def lt_of(t, lon): return (t.hour + t.minute / 60.0 + lon / 15.0) % 24.0 def circ_mean_h(a, b): x = math.cos(a * math.pi / 12) + math.cos(b * math.pi / 12) y = math.sin(a * math.pi / 12) + math.sin(b * math.pi / 12) return (math.atan2(y, x) * 12 / math.pi) % 24.0 def in_band(h, lo, hi): if lo <= hi: return lo <= h < hi return h >= lo or h < hi def boot_median(vals, rng, n=NBOOT): if len(vals) < 3: return (None, None) meds = [] for _ in range(n): s = [vals[rng.randrange(len(vals))] for _ in vals] meds.append(st.median(s)) meds.sort() return (meds[int(0.025 * n)], meds[int(0.975 * n)]) def solve(X, y): """Least squares by normal equations and Gaussian elimination with pivoting.""" p = len(X[0]) A = [[0.0] * (p + 1) for _ in range(p)] for row, v in zip(X, y): for i in range(p): ri = row[i] Ai = A[i] for j in range(p): Ai[j] += ri * row[j] Ai[p] += ri * v for i in range(p): piv = max(range(i, p), key=lambda q: abs(A[q][i])) A[i], A[piv] = A[piv], A[i] if abs(A[i][i]) < 1e-12: return None for k in range(i + 1, p): f = A[k][i] / A[i][i] if f: for j in range(i, p + 1): A[k][j] -= f * A[i][j] b = [0.0] * p for i in range(p - 1, -1, -1): b[i] = (A[i][p] - sum(A[i][j] * b[j] for j in range(i + 1, p))) / A[i][i] return b def design(lt, doy, year=None): a = 2 * math.pi * (lt - LT0) / 24.0 d = 2 * math.pi * (doy - DOY0) / 365.25 row = [1.0, math.cos(a), math.sin(a), math.cos(d), math.sin(d)] if year is not None: row.append(year - 2020.0) return row def fit_at_event(recs, rng, nboot=NBOOT): """recs: (ln R, lt, doy). exp(a + b1 + c1) with bootstrap interval.""" X = [design(l, d) for _, l, d in recs] y = [v for v, _, _ in recs] b = solve(X, y) if b is None: return None est = math.exp(b[0] + b[1] + b[3]) amp_lt = math.hypot(b[1], b[2]) amp_season = math.hypot(b[3], b[4]) vals = [] n = len(recs) for _ in range(nboot): idx = [rng.randrange(n) for _ in range(n)] bb = solve([X[i] for i in idx], [y[i] for i in idx]) if bb is not None: vals.append(math.exp(bb[0] + bb[1] + bb[3])) vals.sort() m = len(vals) return {"n": n, "R_at_LT13.1_day131": round(est, 3), "ci95": [round(vals[int(0.025 * m)], 3), round(vals[int(0.975 * m)], 3)], "b_cos_lt": round(b[1], 3), "b_sin_lt": round(b[2], 3), "c_cos_doy": round(b[3], 3), "c_sin_doy": round(b[4], 3), "amp_lt_lnR": round(amp_lt, 3), "amp_season_lnR": round(amp_season, 3)} def fit_trend(recs, rng, nboot=NBOOT): """recs: (ln R, lt, doy, year). The year coefficient, as percent per decade, with a bootstrap interval.""" X = [design(l, d, yr) for _, l, d, yr in recs] y = [v for v, _, _, _ in recs] b = solve(X, y) if b is None: return None vals = [] n = len(recs) for _ in range(nboot): idx = [rng.randrange(n) for _ in range(n)] bb = solve([X[i] for i in idx], [y[i] for i in idx]) if bb is not None: vals.append(bb[5]) vals.sort() m = len(vals) pct = lambda z: 100 * (math.exp(10 * z) - 1) # noqa: E731 return {"n": n, "pct_per_decade": round(pct(b[5]), 1), "ci95": [round(pct(vals[int(0.025 * m)]), 1), round(pct(vals[int(0.975 * m)]), 1)]} def spearman(xs, ys): def ranks(v): order = sorted(range(len(v)), key=lambda i: v[i]) r = [0.0] * len(v) i = 0 while i < len(v): j = i while j + 1 < len(v) and v[order[j + 1]] == v[order[i]]: j += 1 for k in range(i, j + 1): r[order[k]] = (i + j) / 2.0 + 1 i = j + 1 return r rx, ry = ranks(xs), ranks(ys) mx, my = st.mean(rx), st.mean(ry) sxy = sum((a - mx) * (b - my) for a, b in zip(rx, ry)) sxx = sum((a - mx) ** 2 for a in rx) syy = sum((b - my) ** 2 for b in ry) return sxy / math.sqrt(sxx * syy) def pair_records(rows, s, n, coords, metric, sel): """(time, R, LT of the pair, doy, year) for events passing sel(event).""" out = [] for tstr, ev in sorted(rows.items()): if not sel(ev, tstr): continue a, b = ev["stations"].get(s), ev["stations"].get(n) if not a or not b or a["excl"] or b["excl"]: continue x, y = a[metric], b[metric] if x is None or y is None: continue if metric == "A" and (x <= 3 or y <= 3): continue if x <= 0 or y <= 0: continue t = dt.datetime.strptime(tstr, "%Y-%m-%dT%H:%M") lt = circ_mean_h(lt_of(t, coords[s][1]), lt_of(t, coords[n][1])) doy = t.timetuple().tm_yday year = t.year + (doy - 1) / 365.25 out.append((tstr, x / y, lt, doy, year)) return out def summary(vals, rng): if not vals: return {"n": 0} lo, hi = boot_median(vals, rng) return {"n": len(vals), "median": round(st.median(vals), 3), "ci95": [round(lo, 3), round(hi, 3)] if lo else None, "share_above_1": round(sum(v > 1 for v in vals) / len(vals), 3)} def clean(rr, minutes, need): """The rules' exclusions for a stretch of one station-day {minute: (H, D, Z)}: enough steps, no jump, no frozen run.""" steps = run = 0 longest = 1 run = 1 for m in minutes: if m in rr and m - 1 in rr: steps += 1 d = abs(rr[m][0] - rr[m - 1][0]) if d > 500: return False run = run + 1 if d == 0 else 1 longest = max(longest, run) else: run = 1 return steps >= need and longest < 15 def run_tests(root, csvpath): rng = random.Random(SEED) rows, excl = load(csvpath) model = read_igrf(os.path.join(root, "igrf14", "igrf14coeffs.txt")) obs = read_capabilities(os.path.join(root, "gin-probe", "gin-capabilities.xml")) band = band_stations(obs, model) coords = {k: (v["lat"], v["lon"]) for k, v in obs.items()} coords.update(EMBRACE_COORDS) g, h = igrf_at(model, EPOCH_2024) res = {"rules": "test-rules.md", "seed": SEED, "n_boot": NBOOT, "exclusions_by_station": excl} primary = lambda ev, t: ev["size"] >= 20 # noqa: E731 sets = {"primary: size >= 20 nT, SSC and SI": primary, "size >= 15 nT": lambda ev, t: ev["size"] >= 15, "size >= 30 nT": lambda ev, t: ev["size"] >= 30, "SSC only, size >= 20 nT": lambda ev, t: ev["size"] >= 20 and ev["type"] == "SSC", "primary without 10 May 2024": lambda ev, t: ev["size"] >= 20 and t != EVENT, "SSC only 2013-2024, size >= 20 nT": lambda ev, t: ev["size"] >= 20 and ev["type"] == "SSC" and t < "2025", "all positive events size >= 15, 2013-2024": lambda ev, t: ev["size"] >= 15 and t < "2025"} print("Events in CSV:", len(rows), " primary:", sum(1 for e in rows.values() if e["size"] >= 20)) print("Exclusions:", json.dumps(excl)) # ---------------------------------------------------------------- test 1 print("\n=== Test 1: SMS/SJG across events ===") t1 = {} for name, sel in sets.items(): t1[name] = {} for m in METRICS: recs = pair_records(rows, "SMS", "SJG", coords, m, sel) vals = [r[1] for r in recs] sm = summary(vals, rng) ev = [r[1] for r in recs if r[0] == EVENT] if ev: sm["event_10May2024"] = round(ev[0], 3) sm["event_rank_from_top"] = 1 + sum(v > ev[0] for v in vals) t1[name][m] = sm print(name, {m: (t1[name][m].get("median"), t1[name][m].get("ci95"), t1[name][m].get("n"), t1[name][m].get("share_above_1"), t1[name][m].get("event_10May2024"), t1[name][m].get("event_rank_from_top")) for m in METRICS}) res["test1_same_pair"] = t1 # ---------------------------------------------------------------- test 2 print("\n=== Test 2: every pair, by local time ===") t2 = {} lt_on_event = {} tev = dt.datetime.strptime(EVENT, "%Y-%m-%dT%H:%M") for s, n in PAIRS: lt_on_event[f"{s}/{n}"] = round(circ_mean_h(lt_of(tev, coords[s][1]), lt_of(tev, coords[n][1])), 2) res["pair_lt_on_10May2024"] = lt_on_event print("Pair LT on 10 May 2024:", lt_on_event) for s, n in PAIRS: key = f"{s}/{n}" t2[key] = {} for m in METRICS: recs = pair_records(rows, s, n, coords, m, primary) vals_all = [r[1] for r in recs] noon = [r[1] for r in recs if in_band(r[2], 11, 15)] lt0 = lt_on_event[key] near_own = [r[1] for r in recs if abs(((r[2] - lt0 + 12) % 24) - 12) < 2] bins = {} for lo, hi, lab in ((9, 15, "09-15"), (15, 21, "15-21"), (21, 3, "21-03"), (3, 9, "03-09")): bins[lab] = summary([r[1] for r in recs if in_band(r[2], lo, hi)], rng) ev = [r[1] for r in recs if r[0] == EVENT] t2[key][m] = {"all": summary(vals_all, rng), "near_noon_11_15": summary(noon, rng), "at_own_10May_LT_pm2h": summary(near_own, rng), "lt_bins": bins, "event_10May2024": round(ev[0], 3) if ev else None} print(key, "R_A all", t2[key]["A"]["all"], "| near noon", t2[key]["A"]["near_noon_11_15"], "| own LT", t2[key]["A"]["at_own_10May_LT_pm2h"], "| 10 May", t2[key]["A"]["event_10May2024"]) print(" R_V all", t2[key]["V"]["all"].get("median"), "near noon", t2[key]["V"]["near_noon_11_15"].get("median"), t2[key]["V"]["near_noon_11_15"].get("n"), "| R_Pc5 all", t2[key]["Pc5"]["all"].get("median"), "noon", t2[key]["Pc5"]["near_noon_11_15"].get("median"), "| R_PA all", t2[key]["PA"]["all"].get("median"), "noon", t2[key]["PA"]["near_noon_11_15"].get("median")) print(" R_A by LT:", {k: (v.get("median"), v.get("n")) for k, v in t2[key]["A"]["lt_bins"].items()}) print(" R_PA by LT:", {k: (v.get("median"), v.get("n")) for k, v in t2[key]["PA"]["lt_bins"].items()}) res["test2_local_time"] = t2 # ---------------------------------------------------------------- test 3 print("\n=== Test 3: season, near noon and at all local times ===") t3 = {} def season(doy_month): return "May-Aug" if doy_month in (5, 6, 7, 8) else ("Nov-Feb" if doy_month in (11, 12, 1, 2) else "Mar-Apr, Sep-Oct") for s, n in PAIRS: key = f"{s}/{n}" t3[key] = {} for m in METRICS: recs = pair_records(rows, s, n, coords, m, primary) parts = {"near_noon": {}, "all_lt": {}} for lab in ("May-Aug", "Nov-Feb", "Mar-Apr, Sep-Oct"): parts["near_noon"][lab] = summary([r[1] for r in recs if in_band(r[2], 11, 15) and season(int(r[0][5:7])) == lab], rng) parts["all_lt"][lab] = summary([r[1] for r in recs if season(int(r[0][5:7])) == lab], rng) t3[key][m] = parts print(key, "R_A near noon:", {k: (v.get("median"), v.get("n")) for k, v in t3[key]["A"]["near_noon"].items()}, "| all LT:", {k: (v.get("median"), v.get("n")) for k, v in t3[key]["A"]["all_lt"].items()}) res["test3_season"] = t3 # ---------------------------------------------------------------- test 4 print("\n=== Test 4: fitted R at LT 13.1 and day 131, per pair ===") t4 = {} for s, n in PAIRS: key = f"{s}/{n}" t4[key] = {} for m in METRICS: recs = pair_records(rows, s, n, coords, m, primary) t4[key][m] = fit_at_event([(math.log(r[1]), r[2], r[3]) for r in recs], rng) print(key, {m: (t4[key][m]["R_at_LT13.1_day131"], t4[key][m]["ci95"], t4[key][m]["n"]) for m in METRICS}) print(" lnR amplitudes (LT, season) for A:", t4[key]["A"]["amp_lt_lnR"], t4[key]["A"]["amp_season_lnR"], " coefficients", {k: t4[key]["A"][k] for k in ("b_cos_lt", "b_sin_lt", "c_cos_doy", "c_sin_doy")}) res["test4_fit"] = t4 # ---------------------------------------------------------------- test 5 print("\n=== Test 5: trend ===") t5 = {} for s, n in PAIRS: key = f"{s}/{n}" t5[key] = {} for m in METRICS: recs = pair_records(rows, s, n, coords, m, primary) fit = fit_trend([(math.log(r[1]), r[2], r[3], r[4]) for r in recs], rng) per = {} for lab, a, b in (("2013-2017", "2013", "2018"), ("2018-2022", "2018", "2023"), ("2023-2026", "2023", "2027")): per[lab] = summary([r[1] for r in recs if a <= r[0] < b], rng) t5[key][m] = {"fit": fit, "periods": per} print(key, {m: (t5[key][m]["fit"]["pct_per_decade"], t5[key][m]["fit"]["ci95"]) for m in METRICS}) print(" R_A by period", {k: (v.get("median"), v.get("ci95"), v.get("n")) for k, v in t5[key]["A"]["periods"].items()}) print(" R_PA by period", {k: (v.get("median"), v.get("ci95"), v.get("n")) for k, v in t5[key]["PA"]["periods"].items()}) res["test5_trend"] = t5 # ---------------------------------------------------------------- test 6 print("\n=== Test 6: band stations near their own noon ===") stations = dict(band) t6 = {} for code, info in stations.items(): lat, lon, dip, F = info[0], info[1], info[2], info[3] recs = [] for tstr, ev in rows.items(): if ev["size"] < 20: continue r = ev["stations"].get(code) if not r or r["excl"] or r["A"] is None: continue t = dt.datetime.strptime(tstr, "%Y-%m-%dT%H:%M") lt = lt_of(t, lon) if not in_band(lt, 11, 15): continue recs.append((tstr, r["A"] / ev["size"], (r["V"] or 0) / ev["size"], int(tstr[5:7]))) allr = [ev["stations"][code]["A"] / ev["size"] for ev in rows.values() if ev["size"] >= 20 and code in ev["stations"] and not ev["stations"][code]["excl"] and ev["stations"][code]["A"] is not None] t6[code] = {"lat": round(lat, 2), "lon": round(lon, 2), "dip_lat": round(dip, 1), "F": round(F), "network": info[4], "n_noon": len(recs), "A_rel_noon": round(st.median([r[1] for r in recs]), 3) if recs else None, "V_rel_noon": round(st.median([r[2] for r in recs]), 3) if recs else None, "A_rel_noon_MayAug": round(st.median([r[1] for r in recs if r[3] in (5, 6, 7, 8)]), 3) if [r for r in recs if r[3] in (5, 6, 7, 8)] else None, "n_MayAug": len([r for r in recs if r[3] in (5, 6, 7, 8)]), "A_rel_noon_NovFeb": round(st.median([r[1] for r in recs if r[3] in (11, 12, 1, 2)]), 3) if [r for r in recs if r[3] in (11, 12, 1, 2)] else None, "n_NovFeb": len([r for r in recs if r[3] in (11, 12, 1, 2)]), "A_rel_all_lt": round(st.median(allr), 3) if allr else None, "n_all": len(allr)} ok = {k: v for k, v in t6.items() if v["n_noon"] >= 5} print("station lat lon dip F n | A_rel V_rel | MayAug (n) NovFeb (n) | all LT (n)") for k, v in sorted(ok.items(), key=lambda kv: kv[1]["F"]): print(f"{k:4s}{v['lat']:7.1f}{v['lon']:7.1f}{v['dip_lat']:6.1f}{v['F']:7d}{v['n_noon']:4d} | {v['A_rel_noon']:5.2f} {v['V_rel_noon']:5.2f} | {v['A_rel_noon_MayAug']} ({v['n_MayAug']}) {v['A_rel_noon_NovFeb']} ({v['n_NovFeb']}) | {v['A_rel_all_lt']} ({v['n_all']})") xs = [v["F"] for v in ok.values()] res["test6_band"] = {"stations": t6} if len(ok) >= 5: res["test6_band"]["spearman_A_rel_vs_F"] = round(spearman(xs, [v["A_rel_noon"] for v in ok.values()]), 3) res["test6_band"]["spearman_V_rel_vs_F"] = round(spearman(xs, [v["V_rel_noon"] for v in ok.values()]), 3) res["test6_band"]["spearman_A_rel_vs_absdip"] = round(spearman([abs(v["dip_lat"]) for v in ok.values()], [v["A_rel_noon"] for v in ok.values()]), 3) weak = [v for v in ok.values() if v["F"] < 25000] rest = [v for v in ok.values() if v["F"] >= 25000] res["test6_band"]["weak_field_stations"] = sorted(k for k, v in ok.items() if v["F"] < 25000) res["test6_band"]["median_A_rel_weak"] = round(st.median([v["A_rel_noon"] for v in weak]), 3) if weak else None res["test6_band"]["median_A_rel_rest"] = round(st.median([v["A_rel_noon"] for v in rest]), 3) if rest else None res["test6_band"]["median_V_rel_weak"] = round(st.median([v["V_rel_noon"] for v in weak]), 3) if weak else None res["test6_band"]["median_V_rel_rest"] = round(st.median([v["V_rel_noon"] for v in rest]), 3) if rest else None res["test6_band"]["n_stations"] = len(ok) print({k: v for k, v in res["test6_band"].items() if k != "stations"}) print("\nSMS against neighbours (median SMS/X over shared events) and X/SJG") nb = {} for x in NEIGHBOURS: rec_x = {} for lab, cond in (("near_noon_at_SMS", True), ("all_lt", False)): for metric in ("A", "V", "PA"): vals, vals2 = [], [] for tstr, ev in rows.items(): if ev["size"] < 20: continue a, b, c = ev["stations"].get("SMS"), ev["stations"].get(x), ev["stations"].get("SJG") if not a or not b or a["excl"] or b["excl"] or a[metric] is None or b[metric] is None or b[metric] <= 3: continue t = dt.datetime.strptime(tstr, "%Y-%m-%dT%H:%M") if cond and not in_band(lt_of(t, coords["SMS"][1]), 11, 15): continue vals.append(a[metric] / b[metric]) if c and not c["excl"] and c[metric] and c[metric] > 3: vals2.append(b[metric] / c[metric]) rec_x[f"{lab}_{metric}"] = {"SMS_over_X": summary(vals, rng), "X_over_SJG": summary(vals2, rng)} nb[x] = rec_x print(x, "near noon A: SMS/X", rec_x["near_noon_at_SMS_A"]["SMS_over_X"].get("median"), "n", rec_x["near_noon_at_SMS_A"]["SMS_over_X"].get("n"), "| X/SJG", rec_x["near_noon_at_SMS_A"]["X_over_SJG"].get("median"), "n", rec_x["near_noon_at_SMS_A"]["X_over_SJG"].get("n"), "|| all LT A: SMS/X", rec_x["all_lt_A"]["SMS_over_X"].get("median"), "n", rec_x["all_lt_A"]["SMS_over_X"].get("n"), "| X/SJG", rec_x["all_lt_A"]["X_over_SJG"].get("median"), "|| PA all LT SMS/X", rec_x["all_lt_PA"]["SMS_over_X"].get("median"), "X/SJG", rec_x["all_lt_PA"]["X_over_SJG"].get("median")) res["test6_neighbours"] = nb # Rule 6's last bullet, first run after the fresh review of 1 November 2026 found it missing: SJG against the other Northern band stations. print("\nSJG against the other Northern band stations (median SJG/X over shared events)") north = {} for x in sorted(k for k, v in ok.items() if v["dip_lat"] > 0 and k != "SJG"): rec_x = {} for lab, cond in (("near_noon_at_SJG", True), ("all_lt", False)): for metric in ("A", "V", "PA"): vals = [] for tstr, ev in rows.items(): if ev["size"] < 20: continue a, b = ev["stations"].get("SJG"), ev["stations"].get(x) if not a or not b or a["excl"] or b["excl"] or a[metric] is None or b[metric] is None or b[metric] <= 3: continue t = dt.datetime.strptime(tstr, "%Y-%m-%dT%H:%M") if cond and not in_band(lt_of(t, coords["SJG"][1]), 11, 15): continue vals.append(a[metric] / b[metric]) rec_x[f"{lab}_{metric}"] = summary(vals, rng) north[x] = rec_x print(x, "near noon A", rec_x["near_noon_at_SJG_A"].get("median"), "n", rec_x["near_noon_at_SJG_A"].get("n"), "| all LT A", rec_x["all_lt_A"].get("median"), "n", rec_x["all_lt_A"].get("n")) res["test6_SJG_vs_north"] = north # ---------------------------------------------------------------- test 7 print("\n=== Test 7: quiet days 2024 ===") url = "https://kp.gfz.de/app/json/?start=2024-01-01T00:00:00Z&end=2024-12-31T23:59:59Z&index=Ap" ap = json.loads(fetch(root, "gfz-kp", "ap_2024.json", url).decode("utf-8")) aps = {d[:10]: v for d, v in zip(ap["datetime"], ap["Ap"])} sc_days = {t[:10] for t in rows} all_events = read_events(os.path.join(root, "ebre-sc-lists")) sc_days |= {e["t"].date().isoformat() for e in all_events} q = [] for mon in ("01", "04", "07", "10"): days = sorted([d for d in aps if d[5:7] == mon and d not in sc_days], key=lambda d: aps[d])[:5] q += days t7 = [] for d in sorted(q): day = dt.date.fromisoformat(d) a = embrace_day(root, "SMS", day) b, _ = gin_day(root, "SJG", day) if not a or not b: t7.append({"day": d, "note": "missing"}) continue rec = {"day": d, "Ap": aps[d]} pa = {} ok_window = ok_day = True for name, rr in (("SMS", a), ("SJG", b)): mins = [m for m in range(15 * 60 + 5, 19 * 60 + 5) if m in rr and m - 1 in rr] steps = [abs(rr[m][0] - rr[m - 1][0]) for m in mins] pa[name] = sum(steps) rec[name + "_n"] = len(mins) ok_window = ok_window and clean(rr, range(15 * 60 + 5, 19 * 60 + 5), 228) ok_day = ok_day and clean(rr, range(1, 1440), 1368) vals = [rr[m][0] for m in rr] rec[name + "_range"] = round(max(vals) - min(vals), 1) rec[name + "_PA"] = round(pa[name], 1) rec["window_valid"] = ok_window rec["day_valid"] = ok_day if ok_window: rec["R_PA"] = round(pa["SMS"] / pa["SJG"], 2) if ok_day: rec["R_range"] = round(rec["SMS_range"] / rec["SJG_range"], 2) t7.append(rec) print(rec) good = [r for r in t7 if "R_PA" in r] gooday = [r for r in t7 if "R_range" in r] res["test7_quiet"] = {"days": t7, "n_chosen": len(t7), "n_window_valid": len(good), "n_day_valid": len(gooday), "median_R_PA": round(st.median([r["R_PA"] for r in good]), 2) if good else None, "median_R_range": round(st.median([r["R_range"] for r in gooday]), 2) if gooday else None, "median_PA_SMS": round(st.median([r["SMS_PA"] for r in good]), 1) if good else None, "median_PA_SJG": round(st.median([r["SJG_PA"] for r in good]), 1) if good else None} print({k: v for k, v in res["test7_quiet"].items() if k != "days"}) # ---------------------------------------------------------------- test 9 print("\n=== Test 9: step ratio dZ/dH ===") t9 = {} for code in [c for p in PAIRS for c in p] + ["PIL", "VSS", "JAT", "SJC", "CXP"]: vals = [] for ev in rows.values(): r = ev["stations"].get(code) if not r or r["excl"] or r["dH"] is None or r["dZ"] is None or abs(r["dH"]) < 10: continue vals.append(r["dZ"] / r["dH"]) if vals: t9[code] = {"n": len(vals), "median_dZ_over_dH": round(st.median(vals), 3), "median_abs": round(st.median([abs(v) for v in vals]), 3)} res["test9_induction"] = t9 print(t9) return res # ------------------------------------------------------------------ checks added after the results ISLANDS = {"SJG", "HON", "PPT", "API", "IPM", "ASC", "GUI", "CKI"} # small islands far from a continent's interior WEAK = {"SMS", "PIL", "SJC", "JAT", "CXP", "VSS"} QUIET_OTHERS = ("VSS", "PIL", "FRD", "BSL", "TUC", "HON", "GUI") OUTSIDE = ("GUI", "KNY", "ABG", "HON") # the stations of the "who moves" table: outside South America, with long public records SPANS = (("1999-2026", 1998, 2027), ("1999-2012", 1998, 2013), ("2005-2026", 2005, 2027), ("2013-2026", 2013, 2027)) def read_old_list(path): """Ebre's lists for 2000-2012: rise times (5), amplitudes (5), qualifications (5), observatories (5), maybe a type.""" out, bad = [], [] for line in open(path, encoding="latin-1"): m = re.match(r"\s*(\d{4}) (\d\d) (\d\d) (\d\d)\s+(\d\d(?:\.\d)?)\s+(.*)$", line) if not m: continue rest = m.group(6).split() try: t = dt.datetime(int(m.group(1)), int(m.group(2)), int(m.group(3)), int(m.group(4)), int(float(m.group(5)))) if len(rest) >= 12 and rest[7].isalpha(): # the older format: duration, one amplitude, 5 qualifications, 5 observatories amps = [float(rest[1])] else: amps = [float(x) for x in rest[5:10]] except ValueError: bad.append(line.strip()) continue typ = rest[-1] if rest[-1] in ("SSC", "SI", "SI-") else "" out.append({"t": t, "amps": amps, "size": st.mean(amps), "type": typ or "SC"}) return out, bad def fit_trend_generic(recs, rng, extra=False, nboot=2000): """recs: (ln R, lt, doy, year[, ln size]). Year coefficient as percent per decade with a bootstrap interval.""" X = [] for r in recs: row = design(r[1], r[2], r[3]) if extra: row = row + [r[4]] X.append(row) y = [r[0] for r in recs] b = solve(X, y) if b is None: return None vals = [] n = len(recs) for _ in range(nboot): idx = [rng.randrange(n) for _ in range(n)] bb = solve([X[i] for i in idx], [y[i] for i in idx]) if bb is not None: vals.append(bb[5]) vals.sort() m = len(vals) pct = lambda z: 100 * (math.exp(10 * z) - 1) # noqa: E731 return {"n": n, "pct_per_decade": round(pct(b[5]), 1), "ci95": [round(pct(vals[int(0.025 * m)]), 1), round(pct(vals[int(0.975 * m)]), 1)], "ci98.75": [round(pct(vals[int(0.00625 * m)]), 1), round(pct(vals[int(0.99375 * m)]), 1)]} def run_added(root, csvpath, t6): rng = random.Random(SEED + 1) rows, _ = load(csvpath) model = read_igrf(os.path.join(root, "igrf14", "igrf14coeffs.txt")) obs = read_capabilities(os.path.join(root, "gin-probe", "gin-capabilities.xml")) coords = {k: (v["lat"], v["lon"]) for k, v in obs.items()} coords.update(EMBRACE_COORDS) res = {"note": "All added after the results of test-rules.md were read."} primary = lambda ev, t: ev["size"] >= 20 # noqa: E731 # ---------------------------------------------------------------- A. long baseline print("=== A. VSS/SJG long baseline ===") ev_old, bad = [], [] for y in range(1998, 2013): url = f"https://www.obsebre.es/php/geomagnetisme/vrapides/ssc_{y}_d.txt" fetch(root, "ebre-sc-lists", f"ssc_{y}_d.txt", url) e, b = read_old_list(os.path.join(root, "ebre-sc-lists", f"ssc_{y}_d.txt")) ev_old += e bad += b print("old-list lines that could not be read:", len(bad), bad[:3]) res["A_old_list_lines_unread"] = len(bad) ev_old = [e for e in ev_old if e["size"] >= 20] print("events 1998-2012 with size >= 20:", len(ev_old)) long_rows = {} for e in ev_old: rec = {"size": e["size"], "type": e["type"], "stations": {}} for code in ("VSS", "SJG"): s = series(root, code, e["t"] - dt.timedelta(minutes=120), e["t"] + dt.timedelta(minutes=120), source="gin") mm = metrics(s, e["t"]) if s else {"excl": "no data"} rec["stations"][code] = {"excl": mm.get("excl", ""), "A": mm.get("A"), "V": mm.get("V"), "Pc5": mm.get("Pc5"), "PA": mm.get("PA"), "dH": mm.get("dH"), "dZ": mm.get("dZ")} long_rows[e["t"].strftime("%Y-%m-%dT%H:%M")] = rec merged = dict(long_rows) merged.update({t: v for t, v in rows.items() if v["size"] >= 20}) per_bins = [("1998-2002", 1998, 2003), ("2003-2007", 2003, 2008), ("2008-2012", 2008, 2013), ("2013-2017", 2013, 2018), ("2018-2022", 2018, 2023), ("2023-2026", 2023, 2027)] resA = {} for m in METRICS: recs = pair_records(merged, "VSS", "SJG", coords, m, primary) binsum = {lab: summary([r[1] for r in recs if a <= int(r[0][:4]) < b], rng) for lab, a, b in per_bins} fit = fit_trend_generic([(math.log(r[1]), r[2], r[3], r[4]) for r in recs], rng) resA[m] = {"n": len(recs), "first": recs[0][0] if recs else None, "last": recs[-1][0] if recs else None, "all": summary([r[1] for r in recs], rng), "by_period": binsum, "fit": fit} print(m, "n", len(recs), "first", resA[m]["first"], "all", resA[m]["all"].get("median"), resA[m]["all"].get("share_above_1"), "fit", fit) print(" by period", {k: (v.get("median"), v.get("ci95"), v.get("n")) for k, v in binsum.items()}) res["A_long_baseline_VSS_SJG"] = resA # F ratio over the same span F = lambda code, y: elements(*igrf_at(model, y), *coords[code])["F"] # noqa: E731 res["F_ratio_SJG_over_VSS"] = {str(y): round(F("SJG", float(y)) / F("VSS", float(y)), 4) for y in (2000, 2013, 2026)} print("F(SJG)/F(VSS):", res["F_ratio_SJG_over_VSS"]) # ---------------------------------------------------------------- B. the trend, more closely print("\n=== B. trends ===") resB = {} for s, n in (("SMS", "SJG"), ("VSS", "SJG"), ("PIL", "SJG"), ("JAT", "SJG"), ("CXP", "SJG"), ("SMS", "VSS"), ("SMS", "PIL"), ("VSS", "PIL")): key = f"{s}/{n}" resB[key] = {} for m in ("A", "V", "Pc5", "PA"): recs = pair_records(rows, s, n, coords, m, primary) sizes = {t: rows[t]["size"] for t in rows} if len(recs) < 12: continue plain = fit_trend_generic([(math.log(r[1]), r[2], r[3], r[4]) for r in recs], rng) cov = fit_trend_generic([(math.log(r[1]), r[2], r[3], r[4], math.log(sizes[r[0]])) for r in recs], rng, extra=True) resB[key][m] = {"plain": plain, "with_size": cov} print(key, {m: (v["plain"]["n"], v["plain"]["pct_per_decade"], v["plain"]["ci95"], "size-adj", v["with_size"]["pct_per_decade"], v["with_size"]["ci95"]) for m, v in resB[key].items()}) res["B_trends"] = resB byyear = {} for m in ("A", "V", "Pc5", "PA"): recs = pair_records(rows, "SMS", "SJG", coords, m, primary) byyear[m] = {str(y): [round(st.median([r[1] for r in recs if r[0][:4] == str(y)]), 2), sum(1 for r in recs if r[0][:4] == str(y))] for y in range(2013, 2027) if any(r[0][:4] == str(y) for r in recs)} print("SMS/SJG", m, "by year", byyear[m]) res["B_by_year_SMS_SJG"] = byyear # the one-minute product at SJG by year (headers of the day files) prod = {} for f in sorted(glob.glob(os.path.join(root, "gin-intermagnet-minute", "SJG", "*", "*.txt"))): yr = os.path.basename(os.path.dirname(f)) txt = open(f, encoding="latin-1").read(2500) mt = re.search(r"Data Interval Type\s+(.*?)\s*\|", txt) md = re.search(r"Data Type\s+(.*?)\s*\|", txt) prod.setdefault(yr, set()).add(((mt.group(1) if mt else "?"), (md.group(1) if md else "?"))) res["B_SJG_product_by_year"] = {k: sorted(map(list, v)) for k, v in prod.items()} # who moves: South American stations and SJG, each against stations outside the region print("Against stations outside South America (percent per decade, 95% interval):") resB2 = {} for s, n in (("VSS", "GUI"), ("VSS", "KNY"), ("VSS", "ABG"), ("VSS", "HON"), ("PIL", "GUI"), ("PIL", "KNY"), ("PIL", "ABG"), ("PIL", "HON"), ("SMS", "GUI"), ("SMS", "KNY"), ("SMS", "ABG"), ("SMS", "HON"), ("SJG", "GUI"), ("SJG", "KNY"), ("SJG", "ABG"), ("SJG", "HON"), ("HON", "GUI"), ("HON", "KNY"), ("GUI", "KNY"), ("ABG", "KNY")): key = f"{s}/{n}" resB2[key] = {} for m in ("A", "V"): recs = pair_records(rows, s, n, coords, m, primary) resB2[key][m] = fit_trend_generic([(math.log(r[1]), r[2], r[3], r[4]) for r in recs], rng) print(" ", key, "A", resB2[key]["A"]["pct_per_decade"], resB2[key]["A"]["ci95"], "| V", resB2[key]["V"]["pct_per_decade"], resB2[key]["V"]["ci95"], "| n", resB2[key]["V"]["n"]) res["B2_against_outside_stations"] = resB2 print("SJG one-minute product by year:", {k: sorted(v) for k, v in prod.items() if k >= "2013"}) # ---------------------------------------------------------------- C. EMBRACE against INTERMAGNET print("\n=== C. EMBRACE neighbours against VSS ===") resC = {} for s, n in (("CXP", "VSS"), ("SJC", "VSS")): key = f"{s}/{n}" resC[key] = {m: summary([r[1] for r in pair_records(rows, s, n, coords, m, primary)], rng) for m in ("A", "V", "Pc5", "PA")} print(key, {m: (v.get("median"), v.get("ci95"), v.get("n")) for m, v in resC[key].items()}) res["C_embrace_vs_intermagnet"] = resC # ---------------------------------------------------------------- D. all four above 1, by season print("\n=== D. all four ratios above 1, and all four below 1 ===") resD = {} season = lambda mon: "May-Aug" if mon in (5, 6, 7, 8) else ("Nov-Feb" if mon in (11, 12, 1, 2) else "Equinox") # noqa: E731 for s, n in PAIRS: key = f"{s}/{n}" byev = {} for m in METRICS: for r in pair_records(rows, s, n, coords, m, primary): byev.setdefault(r[0], {})[m] = r[1] full = {t: v for t, v in byev.items() if len(v) == 4} resD[key] = {} for lab in ("all", "May-Aug", "Nov-Feb", "Equinox"): sel = [v for t, v in full.items() if lab == "all" or season(int(t[5:7])) == lab] resD[key][lab] = {"n": len(sel), "all_above_1": sum(all(x > 1 for x in v.values()) for v in sel), "all_below_1": sum(all(x < 1 for x in v.values()) for v in sel)} print(key, resD[key]) res["D_all_four"] = resD # ---------------------------------------------------------------- E. by setting print("\n=== E. near-noon A_rel by setting (exploratory) ===") ok = {k: v for k, v in t6.items() if v["n_noon"] >= 5} groups = {"South American weak field": [k for k in ok if k in WEAK], "small islands": [k for k in ok if k in ISLANDS], "other": [k for k in ok if k not in WEAK and k not in ISLANDS]} resE = {} for lab, codes in groups.items(): resE[lab] = {"stations": sorted(codes), "median_A_rel": round(st.median([ok[c]["A_rel_noon"] for c in codes]), 3), "median_V_rel": round(st.median([ok[c]["V_rel_noon"] for c in codes]), 3)} print(lab, resE[lab]) rank = sorted(ok, key=lambda c: ok[c]["A_rel_noon"]) resE["lowest_A_rel"] = rank[:4] resE["highest_A_rel"] = rank[-4:] rankv = sorted(ok, key=lambda c: ok[c]["V_rel_noon"]) resE["lowest_V_rel"] = rankv[:4] resE["highest_V_rel"] = rankv[-4:] others = [ok[c]["A_rel_noon"] for c in ok if c not in WEAK and c != "SJG"] med = st.median(others) resE["median_A_rel_excl_weak_and_SJG"] = round(med, 3) a, b = ok["SMS"]["A_rel_noon"], ok["SJG"]["A_rel_noon"] resE["log_share"] = {"SMS_above_rest": round(math.log(a / med) / math.log(a / b), 3), "rest_above_SJG": round(math.log(med / b) / math.log(a / b), 3)} print({k: resE[k] for k in ("lowest_A_rel", "highest_A_rel", "lowest_V_rel", "highest_V_rel", "median_A_rel_excl_weak_and_SJG", "log_share")}) # permutation tests across stations (the table has 22 stations, so these are rough) codes = sorted(ok) F_all = [ok[c]["F"] for c in codes] A_all = [ok[c]["A_rel_noon"] for c in codes] V_all = [ok[c]["V_rel_noon"] for c in codes] def perm_p(xs, ys, nperm=20000, stream=None): obs = spearman(xs, ys) hits = 0 ys2 = list(ys) for _ in range(nperm): (stream or rng).shuffle(ys2) if abs(spearman(xs, ys2)) >= abs(obs) - 1e-12: hits += 1 return round(obs, 3), round((hits + 1) / (nperm + 1), 4) resE["spearman_A_rel_F"] = perm_p(F_all, A_all) resE["spearman_V_rel_F"] = perm_p(F_all, V_all) weak_v = [ok[c]["A_rel_noon"] for c in groups["South American weak field"]] other_v = [ok[c]["A_rel_noon"] for c in groups["other"]] u = sum((a > b) + 0.5 * (a == b) for a in weak_v for b in other_v) pool = weak_v + other_v k = len(weak_v) hits = total = 0 for _ in range(20000): rng.shuffle(pool) uu = sum((a > b) + 0.5 * (a == b) for a in pool[:k] for b in pool[k:]) total += 1 if abs(uu - k * len(other_v) / 2.0) >= abs(u - k * len(other_v) / 2.0) - 1e-9: hits += 1 resE["weak_vs_other_U_of"] = [u, k * len(other_v)] resE["weak_vs_other_perm_p"] = round((hits + 1) / (total + 1), 4) # the cleanest single comparison: Kakadu has SMS's hemisphere and dip latitude and twice its field resE["SMS_vs_KDU"] = {c: {k2: ok[c][k2] for k2 in ("dip_lat", "F", "n_noon", "A_rel_noon", "V_rel_noon", "A_rel_noon_MayAug", "n_MayAug", "A_rel_noon_NovFeb", "n_NovFeb")} for c in ("SMS", "KDU")} rest_codes = [c for c in ok if c not in WEAK] resE["season_excess"] = {} for key, lab in (("A_rel_noon_MayAug", "May-Aug"), ("A_rel_noon_NovFeb", "Nov-Feb")): wk = st.median([ok[c][key] for c in groups["South American weak field"] if ok[c][key] is not None]) rs = st.median([ok[c][key] for c in rest_codes if ok[c][key] is not None]) resE["season_excess"][lab] = {"weak_field_median": round(wk, 3), "rest_median": round(rs, 3), "excess_pct": round(100 * (wk / rs - 1), 1)} print({k: resE[k] for k in ("spearman_A_rel_F", "spearman_V_rel_F", "weak_vs_other_U_of", "weak_vs_other_perm_p", "SMS_vs_KDU", "season_excess")}) res["E_setting"] = resE # ---------------------------------------------------------------- F. quiet days, 2022-2024 print("\n=== F. quiet days 2022-2024 ===") url = "https://kp.gfz.de/app/json/?start=2022-01-01T00:00:00Z&end=2024-12-31T23:59:59Z&index=Ap" ap = json.loads(fetch(root, "gfz-kp", "ap_2022_2024.json", url).decode("utf-8")) aps = {d[:10]: v for d, v in zip(ap["datetime"], ap["Ap"])} sc_days = {e["t"].date().isoformat() for e in read_events(os.path.join(root, "ebre-sc-lists"))} q = [] for yr in ("2022", "2023", "2024"): for mon in ("01", "04", "07", "10"): days = sorted([d for d in aps if d[:4] == yr and d[5:7] == mon and d not in sc_days], key=lambda d: aps[d])[:5] q += days out = [] for d in sorted(q): day = dt.date.fromisoformat(d) a = embrace_day(root, "SMS", day) b, _ = gin_day(root, "SJG", day) if not a or not b: out.append({"day": d, "note": "missing"}) continue rec = {"day": d, "Ap": aps[d]} pa = {} okw = okd = True for name, rr in (("SMS", a), ("SJG", b)): mins = [m for m in range(15 * 60 + 5, 19 * 60 + 5) if m in rr and m - 1 in rr] pa[name] = sum(abs(rr[m][0] - rr[m - 1][0]) for m in mins) okw = okw and clean(rr, range(15 * 60 + 5, 19 * 60 + 5), 228) okd = okd and clean(rr, range(1, 1440), 1368) vals = [rr[m][0] for m in rr] rec[name + "_range"] = round(max(vals) - min(vals), 1) rec[name + "_PA"] = round(pa[name], 1) if okw: rec["R_PA"] = round(pa["SMS"] / pa["SJG"], 2) if okd: rec["R_range"] = round(rec["SMS_range"] / rec["SJG_range"], 2) # the same quiet window at two INTERMAGNET observatories in the region, to separate SMS's own noise from a regional signal, # and at observatories outside it (the last five added after the fresh review), to see how far the factor belongs to SJG for name in QUIET_OTHERS: rr, _ = gin_day(root, name, day) if rr and clean(rr, range(15 * 60 + 5, 19 * 60 + 5), 228): rec[name + "_PA"] = round(sum(abs(rr[m][0] - rr[m - 1][0]) for m in range(15 * 60 + 5, 19 * 60 + 5) if m in rr and m - 1 in rr), 1) if clean(b, range(15 * 60 + 5, 19 * 60 + 5), 228): rec[name + "_over_SJG"] = round(rec[name + "_PA"] / pa["SJG"], 2) out.append(rec) goodw = [r for r in out if "R_PA" in r] goodd = [r for r in out if "R_range" in r] bymon = {mon: {"n_window": sum(1 for r in goodw if r["day"][5:7] == mon), "median_R_PA": round(st.median([r["R_PA"] for r in goodw if r["day"][5:7] == mon]), 2) if any(r["day"][5:7] == mon for r in goodw) else None, "n_day": sum(1 for r in goodd if r["day"][5:7] == mon), "median_R_range": round(st.median([r["R_range"] for r in goodd if r["day"][5:7] == mon]), 2) if any(r["day"][5:7] == mon for r in goodd) else None} for mon in ("01", "04", "07", "10")} res["F_quiet"] = {"chosen": len(out), "n_window_valid": len(goodw), "n_day_valid": len(goodd), "median_R_PA": round(st.median([r["R_PA"] for r in goodw]), 2), "share_R_PA_above_1": round(sum(r["R_PA"] > 1 for r in goodw) / len(goodw), 3), "median_R_range": round(st.median([r["R_range"] for r in goodd]), 2), "median_PA_SMS": round(st.median([r["SMS_PA"] for r in goodw]), 1), "median_PA_SJG": round(st.median([r["SJG_PA"] for r in goodw]), 1), "by_month": bymon, "days": out} for name in QUIET_OTHERS: vals = [r[name + "_over_SJG"] for r in out if name + "_over_SJG" in r] res["F_quiet"][name + "_over_SJG_PA"] = {"n": len(vals), "median": round(st.median(vals), 2) if vals else None} vals = [r["R_PA"] / r["VSS_over_SJG"] for r in out if "R_PA" in r and "VSS_over_SJG" in r] res["F_quiet"]["SMS_over_VSS_PA"] = {"n": len(vals), "median": round(st.median(vals), 2) if vals else None} # PA does respond to a jolt: how far events stand above a quiet window, at each station (added after the fresh review) qa, qb = res["F_quiet"]["median_PA_SMS"], res["F_quiet"]["median_PA_SJG"] pa_rows = [] for tstr, ev in rows.items(): x, y = ev["stations"].get("SMS"), ev["stations"].get("SJG") if ev["size"] >= 20 and x and y and not x["excl"] and not y["excl"] and x["PA"] and y["PA"]: pa_rows.append((tstr, x["PA"], y["PA"])) ev10 = [r for r in pa_rows if r[0] == EVENT][0] res["F_quiet"]["PA_events_against_quiet"] = { "n_events": len(pa_rows), "quiet_median_SMS": qa, "quiet_median_SJG": qb, "10May2024": {"PA_SMS": round(ev10[1], 1), "PA_SJG": round(ev10[2], 1), "times_quiet_SMS": round(ev10[1] / qa, 1), "times_quiet_SJG": round(ev10[2] / qb, 1), "excess_ratio": round((ev10[1] - qa) / (ev10[2] - qb), 2)}, "median_times_quiet_SMS": round(st.median([r[1] / qa for r in pa_rows]), 1), "median_times_quiet_SJG": round(st.median([r[2] / qb for r in pa_rows]), 1), "median_excess_ratio": round(st.median([(r[1] - qa) / (r[2] - qb) for r in pa_rows]), 2)} print({k: v for k, v in res["F_quiet"].items() if k != "days"}) # ---------------------------------------------------------------- G. 10 May among May-August events print("\n=== G. 10 May 2024 among May-August events ===") resG = {} for m in METRICS: recs = [r for r in pair_records(rows, "SMS", "SJG", coords, m, primary) if season(int(r[0][5:7])) == "May-Aug"] vals = [r[1] for r in recs] ev = [r[1] for r in recs if r[0] == EVENT][0] resG[m] = {"n": len(vals), "median": round(st.median(vals), 3), "event": round(ev, 3), "rank_from_top": 1 + sum(v > ev for v in vals)} res["G_within_season"] = resG print(resG) # ---------------------------------------------------------------- H. the long baseline, span by span and against outside stations print("\n=== H. VSS/SJG by span, and against stations outside South America ===") rng2 = random.Random(SEED + 2) # a stream of its own, so that these additions leave the stream of the checks above alone resH = {"note": "Added after the fresh review of 1 November 2026.", "spans_VSS_SJG": {}, "outside_stations": list(OUTSIDE), "over_outside": {}, "SMS_SJG_leave_out": {}} log_recs = lambda recs: [(math.log(r[1]), r[2], r[3], r[4]) for r in recs] # noqa: E731 for m in METRICS: recs = pair_records(merged, "VSS", "SJG", coords, m, primary) known = {"1999-2026": resA[m]["fit"], "2013-2026": resB["VSS/SJG"][m]["plain"]} # the same fits as above, so that the page prints one interval for each resH["spans_VSS_SJG"][m] = {lab: known.get(lab) or fit_trend_generic(log_recs([r for r in recs if a <= int(r[0][:4]) < b]), rng2) for lab, a, b in SPANS} print("VSS/SJG", m, {lab: (v["n"], v["pct_per_decade"], v["ci95"]) for lab, v in resH["spans_VSS_SJG"][m].items()}) for tstr, rec in long_rows.items(): tt = dt.datetime.strptime(tstr, "%Y-%m-%dT%H:%M") for code in OUTSIDE: s = series(root, code, tt - dt.timedelta(minutes=120), tt + dt.timedelta(minutes=120), source="gin") mm = metrics(s, tt) if s else {"excl": "no data"} rec["stations"][code] = {"excl": mm.get("excl", ""), "A": mm.get("A"), "V": mm.get("V"), "Pc5": mm.get("Pc5"), "PA": mm.get("PA"), "dH": mm.get("dH"), "dZ": mm.get("dZ")} def usable(ev, code, m): r = ev["stations"].get(code) return bool(r) and not r["excl"] and r[m] is not None and r[m] > 3 def over_outside(tops, m, lon, need): """(time, R, LT, day of year, year): the mean of tops over the mean of the outside stations, where at least 3 of those are usable.""" out_ = [] for tstr, ev in sorted(merged.items()): ou = [ev["stations"][c][m] for c in OUTSIDE if usable(ev, c, m)] top = [ev["stations"][c][m] for c in tops if usable(ev, c, m)] if len(ou) < 3 or len(top) < need: continue tt = dt.datetime.strptime(tstr, "%Y-%m-%dT%H:%M") doy = tt.timetuple().tm_yday out_.append((tstr, st.mean(top) / st.mean(ou), lt_of(tt, lon), doy, tt.year + (doy - 1) / 365.25)) return out_ for lab0, tops, lon, need, spans in (("VSS", ("VSS",), coords["VSS"][1], 1, SPANS), ("SJG", ("SJG",), coords["SJG"][1], 1, SPANS), ("SMS, VSS and PIL", ("SMS", "VSS", "PIL"), -55.0, 2, SPANS[3:])): resH["over_outside"][lab0] = {} for m in METRICS: recs = over_outside(tops, m, lon, need) fits = {} for lab, a, b in spans: sel = [r for r in recs if a <= int(r[0][:4]) < b] if len(sel) >= 12: fits[lab] = fit_trend_generic(log_recs(sel), rng2) resH["over_outside"][lab0][m] = {"n": len(recs), "by_period": {lab: summary([r[1] for r in recs if a <= int(r[0][:4]) < b], rng2) for lab, a, b in per_bins}, "fits": fits} print(lab0, "over outside,", m, "n", len(recs), "periods", [(v.get("median"), v.get("n")) for v in resH["over_outside"][lab0][m]["by_period"].values()], "fits", {k: (v["n"], v["pct_per_decade"], v["ci95"]) for k, v in fits.items()}) for m in METRICS: recs = pair_records(rows, "SMS", "SJG", coords, m, primary) top = max(recs, key=lambda r: r[1]) resH["SMS_SJG_leave_out"][m] = {"without_2023": fit_trend_generic(log_recs([r for r in recs if r[0][:4] != "2023"]), rng2), "without_2026": fit_trend_generic(log_recs([r for r in recs if r[0][:4] != "2026"]), rng2), "n_2023": sum(r[0][:4] == "2023" for r in recs), "n_2026": sum(r[0][:4] == "2026" for r in recs), "largest_ratio": {"time": top[0], "R": round(top[1], 2)}} print("SMS/SJG", m, resH["SMS_SJG_leave_out"][m]) res["H_long_baseline_by_span"] = resH # ---------------------------------------------------------------- I. the season, pair by pair print("\n=== I. the season, pair by pair ===") resI = {"note": "Added after the fresh review of 1 November 2026.", "shift_lnR_A_May_Aug_to_Nov_Feb": {}, "only_pair_with_all_four_above_1": {}} full = {} for s, n in PAIRS: key = f"{s}/{n}" recs = pair_records(rows, s, n, coords, "A", primary) xa = [math.log(r[1]) for r in recs if season(int(r[0][5:7])) == "May-Aug"] xb = [math.log(r[1]) for r in recs if season(int(r[0][5:7])) == "Nov-Feb"] draws = sorted(st.median([xb[rng2.randrange(len(xb))] for _ in xb]) - st.median([xa[rng2.randrange(len(xa))] for _ in xa]) for _ in range(NBOOT)) resI["shift_lnR_A_May_Aug_to_Nov_Feb"][key] = {"n_May_Aug": len(xa), "n_Nov_Feb": len(xb), "shift": round(st.median(xb) - st.median(xa), 2), "ci95": [round(draws[int(0.025 * NBOOT)], 2), round(draws[int(0.975 * NBOOT)], 2)]} byev = {} for m in METRICS: for r in pair_records(rows, s, n, coords, m, primary): byev.setdefault(r[0], {})[m] = r[1] full[key] = {tstr: v for tstr, v in byev.items() if len(v) == 4} print(resI["shift_lnR_A_May_Aug_to_Nov_Feb"]) common = set.intersection(*[set(v) for v in full.values()]) above = lambda k, tstr: all(x > 1 for x in full[k][tstr].values()) # noqa: E731 for lab in ("May-Aug", "Nov-Feb", "Equinox"): sel = [tstr for tstr in common if season(int(tstr[5:7])) == lab] resI["only_pair_with_all_four_above_1"][lab] = {"n_events_all_five_pairs_usable": len(sel), "SMS_SJG_all_above_1": sum(above("SMS/SJG", tstr) for tstr in sel), "SMS_SJG_the_only_pair": sum(above("SMS/SJG", tstr) and not any(above(k, tstr) for k in full if k != "SMS/SJG") for tstr in sel)} print(resI["only_pair_with_all_four_above_1"]) res["I_season_by_pair"] = resI # ---------------------------------------------------------------- J. the ratio's two ends, by yardstick print("\n=== J. the two ends of SMS/SJG, by yardstick ===") bandok = sorted(ok) nonsa = [c for c in bandok if c not in WEAK] yards = {"all stations outside South America except SJG": [c for c in nonsa if c != "SJG"], "non-island stations": [c for c in nonsa if c not in ISLANDS], "small islands except SJG": [c for c in nonsa if c in ISLANDS and c != "SJG"], "southern stations outside South America": [c for c in nonsa if ok[c]["dip_lat"] < 0], "northern stations except SJG": [c for c in nonsa if ok[c]["dip_lat"] > 0 and c != "SJG"]} page_yard = "all stations outside South America except SJG" times = sorted(tstr for tstr, ev in rows.items() if ev["size"] >= 20) near = {c: {} for c in bandok} for tstr in times: tt = dt.datetime.strptime(tstr, "%Y-%m-%dT%H:%M") for c in bandok: r = rows[tstr]["stations"].get(c) if r and not r["excl"] and r["A"] is not None and in_band(lt_of(tt, coords[c][1]), 11, 15): near[c][tstr] = (r["A"], r["V"] or 0.0) rel = {c: {tstr: (v[0] / rows[tstr]["size"], v[1] / rows[tstr]["size"]) for tstr, v in near[c].items()} for c in bandok} def station_medians(pick, idx, table): res_ = {} for c in bandok: v = [table[c][tstr][idx] for tstr in pick if tstr in table[c]] if v: res_[c] = st.median(v) return res_ def two_ends(mm, group): y = st.median([mm[c] for c in group if c in mm]) wk = st.median([mm[c] for c in WEAK if c in mm]) return {"median": round(y, 3), "SMS_share": round(math.log(mm["SMS"] / y) / math.log(mm["SMS"] / mm["SJG"]), 2), "SJG_share": round(math.log(y / mm["SJG"]) / math.log(mm["SMS"] / mm["SJG"]), 2), "SMS_over_it_pct": round(100 * (mm["SMS"] / y - 1)), "weak_field_stations_over_it_pct": round(100 * (wk / y - 1))} ma = station_medians(times, 0, rel) mv = station_medians(times, 1, rel) resJ = {"note": "Added after the fresh review of 1 November 2026. Near-noon amplitude relative to event size, medians per station.", "SMS": round(ma["SMS"], 3), "SJG": round(ma["SJG"], 3), "weak_field_median": round(st.median([ma[c] for c in WEAK]), 3), "yardsticks": {lab: {"stations": g, **two_ends(ma, g)} for lab, g in yards.items()}} # the event sized from the page's stations themselves def own_size(tstr): v = [rows[tstr]["stations"][c]["A"] for c in yards[page_yard] if c in rows[tstr]["stations"] and not rows[tstr]["stations"][c]["excl"] and rows[tstr]["stations"][c]["A"]] return st.median(v) if len(v) >= 8 else None sizes2 = {tstr: own_size(tstr) for tstr in times} rel2 = {c: {tstr: (v[0] / sizes2[tstr], v[1] / sizes2[tstr]) for tstr, v in near[c].items() if sizes2[tstr]} for c in bandok} ma2 = station_medians(times, 0, rel2) low2 = sorted(ma2, key=lambda c: ma2[c]) resJ["event_sized_by_the_stations_themselves"] = {"SMS": round(ma2["SMS"], 3), "SJG": round(ma2["SJG"], 3), **two_ends(ma2, yards[page_yard]), "lowest": [[c, round(ma2[c], 2)] for c in low2[:3]]} # by season resJ["by_season"] = {} for lab, mon in (("May-Aug", (5, 6, 7, 8)), ("Nov-Feb", (11, 12, 1, 2)), ("Equinox", (3, 4, 9, 10))): pick = [tstr for tstr in times if int(tstr[5:7]) in mon] mm = station_medians(pick, 0, rel) wk = st.median([mm[c] for c in WEAK if c in mm]) resJ["by_season"][lab] = {"n_SMS": sum(tstr in rel["SMS"] for tstr in pick), "n_SJG": sum(tstr in rel["SJG"] for tstr in pick), "SMS": round(mm["SMS"], 2), "SJG": round(mm["SJG"], 2), **two_ends(mm, yards[page_yard]), "weak_median": round(wk, 3), "weak_over_southern_pct": round(100 * (wk / st.median([mm[c] for c in yards["southern stations outside South America"] if c in mm]) - 1)), "weak_over_northern_pct": round(100 * (wk / st.median([mm[c] for c in yards["northern stations except SJG"] if c in mm]) - 1)), "n_range_stations": [min(sum(tstr in rel[c] for tstr in pick) for c in bandok), max(sum(tstr in rel[c] for tstr in pick) for c in bandok)]} # bootstrap over events: every station keeps the events it shares with the others shares = {lab: [] for lab in yards} excess, ppt, kdu_a, kdu_v, sms_sjg = [], [], [], [], [] per_station = {c: [] for c in bandok} lowest = used = 0 for _ in range(NBOOT): pick = [times[rng2.randrange(len(times))] for _ in times] xa = station_medians(pick, 0, rel) xv = station_medians(pick, 1, rel) if len(xa) < len(bandok): continue # a station lost all its near-noon events in this draw used += 1 for lab, g in yards.items(): y = st.median([xa[c] for c in g]) shares[lab].append(math.log(xa["SMS"] / y) / math.log(xa["SMS"] / xa["SJG"])) excess.append(st.median([xa[c] for c in WEAK]) / st.median([xa[c] for c in yards[page_yard]])) ppt.append(xa["PPT"] - xa["SJG"]) kdu_a.append(xa["SMS"] / xa["KDU"]) kdu_v.append(xv["SMS"] / xv["KDU"]) sms_sjg.append(xa["SMS"] / xa["SJG"]) lowest += xa["SJG"] == min(xa.values()) for c in bandok: per_station[c].append(xa[c]) def limits(v, nd=2): v = sorted(v) return [round(v[int(0.025 * len(v))], nd), round(v[int(0.975 * len(v))], nd)] resJ["bootstrap"] = {"draws_used": used, "SMS_share_ci95": {lab: limits(v) for lab, v in shares.items()}, "weak_field_over_page_yardstick_ci95": limits(excess), "weak_field_over_page_yardstick": round(st.median([ma[c] for c in WEAK]) / st.median([ma[c] for c in yards[page_yard]]), 2), "SJG_lowest_in_share_of_draws": round(lowest / used, 2), "PPT_minus_SJG_ci95": limits(ppt), "SMS_over_KDU_A": round(ma["SMS"] / ma["KDU"], 2), "SMS_over_KDU_A_ci95": limits(kdu_a), "SMS_over_KDU_V": round(mv["SMS"] / mv["KDU"], 2), "SMS_over_KDU_V_ci95": limits(kdu_v), "SMS_over_SJG_ci95": limits(sms_sjg), "station_ci95": {c: limits(v) for c, v in per_station.items()}} # how the event size is built evs = {e["t"].strftime("%Y-%m-%dT%H:%M"): e for e in read_events(os.path.join(root, "ebre-sc-lists"))} def size_no_sjg(tstr): e = evs[tstr] v = [x for x, o in zip(e["amps"], e["obs"]) if o != "SJG"] return st.mean(v) if v else e["size"] rel3 = {c: {tstr: (v[0] / size_no_sjg(tstr), v[1] / size_no_sjg(tstr)) for tstr, v in near[c].items()} for c in bandok} ma3 = station_medians(times, 0, rel3) low3 = sorted(ma3, key=lambda c: ma3[c]) resJ["event_size"] = {"n_events": len(times), "events_with_station_among_the_five": {c: sum(c in evs[tstr]["obs"] for tstr in times) for c in ("SJG", "HON", "KNY", "ABG", "GUI")}, "without_SJG_amplitude": {"SMS": round(ma3["SMS"], 2), "SJG": round(ma3["SJG"], 2), **two_ends(ma3, yards[page_yard]), "lowest": [[c, round(ma3[c], 2)] for c in low3[:3]]}} # field strength without South America's stations resJ["spearman_with_F_without_South_America"] = {"n": len(nonsa), "A": perm_p([ok[c]["F"] for c in nonsa], [ok[c]["A_rel_noon"] for c in nonsa], stream=rng2), "V": perm_p([ok[c]["F"] for c in nonsa], [ok[c]["V_rel_noon"] for c in nonsa], stream=rng2)} res["J_two_ends"] = resJ print(json.dumps({k: v for k, v in resJ.items() if k != "note"}, ensure_ascii=False)) return res # ------------------------------------------------------------------ the reproduction, and test 8 PUB = { # the preprint's Tables 4 and 6 (PDF pp. 19, 21): A (nT), V (nT/min), Pc5 (nT), PA (nT) "SMS": (120.66, 65.50, 77.71, 1029.00), "SJG": (65.30, 35.34, 30.48, 588.36), "ASC": (78.10, 46.10, 32.03, 619.80), "GUI": (71.35, 53.77, 42.34, 668.71), "HBK": (94.08, 40.10, 48.63, 678.71), "PAG": (116.02, 67.90, 79.62, 975.60), "GAN": (90.85, 40.18, 49.10, 597.88), "ABG": (104.74, 53.16, 57.15, 708.06), "KDU": (83.89, 38.55, 33.13, 519.69), "KNY": (99.26, 56.80, 49.18, 678.71), "PVE": (162.16, 113.70, 76.78, 1434.70), "HUA": (248.25, 117.50, 141.30, 2318.56), } PUB_QC = {"SMS": 0.905, "SJG": 1.000, "ASC": 0.820, "GUI": 0.929, "HBK": 0.889, "PAG": 0.996, "GAN": 0.825, "ABG": 0.908, "KDU": 0.858, "KNY": 0.941} # Table 3, PDF p. 9 PUB_R = {("SMS", "SJG"): (1.85, 1.85, 2.55, 1.75), ("GAN", "ABG"): (0.87, 0.76, 0.86, 0.84), ("ASC", "GUI"): (1.09, 0.86, 0.76, 0.93), ("KDU", "KNY"): (0.85, 0.68, 0.67, 0.77), ("HBK", "PAG"): (0.81, 0.59, 0.61, 0.70)} # Table 5, PDF p. 20 T_PAPER = dt.datetime(2024, 5, 10, 17, 4) def run_repro(root): T = T_PAPER res = {"onset": "2024-05-10 17:04 UT", "stations": {}, "pairs": {}, "headers": {}} print("code | A ours pub | V ours pub | Pc5 ours pub | PA ours pub | QC-like ours pub | data type") for code, pub in PUB.items(): s = series(root, code, T - dt.timedelta(minutes=120), T + dt.timedelta(minutes=120), source="embrace" if code in ("SMS", "PVE") else "gin") m = metrics(s, T) vals = [s[T + dt.timedelta(minutes=k)][0] for k in range(-120, 121) if (T + dt.timedelta(minutes=k)) in s] C = len(vals) / 241.0 Pnorep = len(set(vals)) / float(len(vals)) qc = (C + 1 + 1 + Pnorep) / 4.0 # physical bounds and non-zero fractions taken as 1 (no value fails them here) hdr = {} if code not in ("SMS", "PVE"): _, hdr = gin_day(root, code, T.date()) res["stations"][code] = {"A": round(m["A"], 2), "V": round(m["V"], 2), "Pc5": round(m["Pc5"], 2), "PA": round(m["PA"], 2), "pub": pub, "n241": m["n241"], "Pnorep": round(Pnorep, 3), "qc_like": round(qc, 3), "pub_qc": PUB_QC.get(code), "data_type": hdr.get("Data Type", "EMBRACE reported")} print(f"{code} | {m['A']:7.2f} {pub[0]:7.2f} | {m['V']:6.2f} {pub[1]:6.2f} | {m['Pc5']:6.2f} {pub[2]:6.2f} | " f"{m['PA']:8.2f} {pub[3]:8.2f} | {qc:5.3f} {PUB_QC.get(code)} | {hdr.get('Data Type', 'EMBRACE reported')}") print() for (a, b), pub in PUB_R.items(): ra, rb = res["stations"][a], res["stations"][b] ours = [round(ra[k] / rb[k], 2) for k in ("A", "V", "Pc5", "PA")] res["pairs"][f"{a}/{b}"] = {"ours": ours, "pub": pub} print(f"{a}/{b} R ours {ours} published {pub}") # how close overall for i, k in enumerate(("A", "V", "Pc5", "PA")): err = [abs(res["stations"][c][k] / PUB[c][i] - 1) for c in PUB] res.setdefault("median_abs_pct_diff", {})[k] = round(100 * st.median(err), 2) print("median absolute difference from the published value, percent:", res["median_abs_pct_diff"]) # ---- filtering: SJG's one-second file on 10 May 2024, resampled to one minute three ways url = GIN.format(code="SJG", rate="second", state="adj-or-rep", date=T.date().isoformat()) data = fetch(root, "gin-intermagnet-second/SJG/2024", "sjg2024-05-10-second.txt", url) sec = {} for line in data.decode("latin-1").splitlines(): if not line.startswith("2024-05-10"): continue p = line.split() hh, mm, ss = p[1].split(":") x, y = float(p[3]), float(p[4]) if abs(x) > 88000 or abs(y) > 88000: continue sec[int(hh) * 3600 + int(mm) * 60 + int(float(ss))] = math.hypot(x, y) def gauss(tc, sigma=15.90, half=45): s = w = 0.0 for k in range(-half, half + 1): v = sec.get(tc + k) if v is not None: g = math.exp(-k * k / (2 * sigma * sigma)) s += g * v w += g return s / w def mean60(tc): vals = [sec[tc + k] for k in range(-30, 30) if tc + k in sec] return sum(vals) / len(vals) def mean60fwd(tc): vals = [sec[tc + k] for k in range(0, 60) if tc + k in sec] return sum(vals) / len(vals) def inst(tc): return sec[tc] def inst30(tc): return sec[tc + 30] t0m = 17 * 60 + 4 filt = {} s_sms = res["stations"]["SMS"] for name, f in (("INTERMAGNET-style Gaussian", gauss), ("60-s mean, centred", mean60), ("60-s mean, forward", mean60fwd), ("instantaneous at :00", inst), ("instantaneous at :30", inst30)): H = {m: f(m * 60) for m in range(t0m - 120, t0m + 121)} base = sum(H[m] for m in range(t0m - 60, t0m)) / 60.0 A = max(H[m] for m in range(t0m, t0m + 13)) - base V = max(abs(H[m] - H[m - 1]) for m in range(t0m, t0m + 13)) PA = sum(abs(H[m] - H[m - 1]) for m in range(t0m - 119, t0m + 121)) filt[name] = {"A": round(A, 2), "V": round(V, 2), "PA": round(PA, 2), "R_A": round(s_sms["A"] / A, 2), "R_V": round(s_sms["V"] / V, 2), "R_PA": round(s_sms["PA"] / PA, 2)} print(name, filt[name]) res["sjg_filtering"] = filt # Timing at onset, and one bad file: derived numbers for the page's limits, no series check = {} for code, src in (("SMS", "embrace"), ("SLZ", "embrace"), ("ARA", "embrace"), ("VSS", "embrace"), ("VSS", "gin")): s = series(root, code, T - dt.timedelta(minutes=10), T + dt.timedelta(minutes=20), source=src) at = lambda k: s[T + dt.timedelta(minutes=k)][0] # noqa: E731 base = sum(at(k) for k in range(-10, 0)) / 10.0 first = next((k for k in range(-1, 16) if at(k) - at(k - 1) >= 10), None) check[f"{code} ({src})"] = {"first_rise_of_10_nT_in_a_minute": (T + dt.timedelta(minutes=first)).strftime("%H:%M") if first is not None else None, "largest_rise_in_first_6_min": round(max(at(k) - base for k in range(0, 7)), 1)} res["onset_check"] = check print("onset check:", check) return res # ------------------------------------------------------------------ the anomaly, the stations and their field lines WGS84_B = WGS84_A * (1 - WGS84_F) THRESHOLDS = (32000, 26000, 25000, 24000) YEARS = [2000.0, 2013.0, 2014.0, 2015.0, 2016.0, 2017.0, 2018.0, 2019.0, 2020.0, 2021.0, 2022.0, 2023.0, 2024.0, 2024.357, 2025.0, 2026.0, 2030.0] STATIONS = {"SMS": (-29.44357, -53.82273), "SJG": (18.113, -66.151), "ASC": (-7.949, -14.378), "GUI": (28.321, -16.441), "HBK": (-25.882, 27.707), "PAG": (42.511, 24.177), "GAN": (-0.695, 73.153), "ABG": (18.638, 72.872), "KDU": (-12.690, 132.470), "KNY": (31.424, 130.880), "PIL": (-31.67, -63.88), "VSS": (-22.40, -43.65), "JAT": (-17.93185, -51.7184), "SJC": (-23.20862, -45.96353), "CXP": (-22.70207, -45.01443), "CHI": (-36.64122, -71.99485), "HUA": (-12.05, -75.33), "PVE": (-8.835, -63.94)} PAIRS = [("SMS", "SJG"), ("ASC", "GUI"), ("HBK", "PAG"), ("GAN", "ABG"), ("KDU", "KNY")] def band_area(lat1, lat2): e = math.sqrt(WGS84_E2) def q(phi): s = math.sin(phi) return s / (1 - WGS84_E2 * s * s) + (1 / (2 * e)) * math.log((1 + e * s) / (1 - e * s)) return 0.5 * WGS84_B ** 2 * (q(math.radians(lat2)) - q(math.radians(lat1))) EARTH = 2 * math.pi * band_area(-90.0, 90.0) def anomaly(g, h): below = {t: 0.0 for t in THRESHOLDS} best = (1e9, 0.0, 0.0) for i in range(180): la = -89.5 + i area = band_area(la - 0.5, la + 0.5) * math.radians(1.0) for j in range(360): lo = -179.5 + j X, Y, Z = field(g, h, la, lo) F = math.sqrt(X * X + Y * Y + Z * Z) for t in THRESHOLDS: if F < t: below[t] += area if F < best[0]: best = (F, la, lo) F0, la0, lo0 = best fine = (1e9, la0, lo0) k = int(round(1.5 / 0.05)) for a in range(-k, k + 1, 2): # coarse-to-fine: 0.1 degree, then 0.05 around the best for b in range(-k, k + 1, 2): la, lo = la0 + a * 0.05, lo0 + b * 0.05 X, Y, Z = field(g, h, la, lo) F = math.sqrt(X * X + Y * Y + Z * Z) if F < fine[0]: fine = (F, la, lo) F1, la1, lo1 = fine for a in (-2, -1, 0, 1, 2): for b in (-2, -1, 0, 1, 2): la, lo = la1 + a * 0.05, lo1 + b * 0.05 X, Y, Z = field(g, h, la, lo) F = math.sqrt(X * X + Y * Y + Z * Z) if F < fine[0]: fine = (F, la, lo) return fine, {str(t): round(100 * v / EARTH, 3) for t, v in below.items()} def run_anomaly(root): check = validate_wmm(os.path.join(root, "wmm2025", "WMM2025COF.zip")) print("WMM2025 test values:", check) model = read_igrf(os.path.join(root, "igrf14", "igrf14coeffs.txt")) res = {"wmm_check": check, "anomaly": [], "stations": {}, "pairs": {}} for y in YEARS: g, h = igrf_at(model, y) (F, la, lo), area = anomaly(g, h) row = {"epoch": y, "min_F_nT": round(F, 1), "lat": round(la, 2), "lon": round(lo, 2), "pct_below": area} res["anomaly"].append(row) print(row, flush=True) for code, (la, lo) in STATIONS.items(): res["stations"][code] = {"lat": la, "lon": lo} for y in (2000.0, 2013.0, 2015.0, 2024.357, 2026.0): g, h = igrf_at(model, y) e = elements(g, h, la, lo) res["stations"][code][str(y)] = {k: round(e[k], 2) for k in ("F", "H", "I", "dip_lat")} # field-line tracing at the event's epoch g, h = igrf_at(model, 2024.357) for s, n in PAIRS: ts = trace(g, h, *STATIONS[s]) tn = trace(g, h, *STATIONS[n]) res["pairs"][f"{s}-{n}"] = { "conj_of_south": [round(ts["lat"], 2), round(ts["lon"], 2)], "conj_of_south_to_north_km": round(gc_km(ts["lat"], ts["lon"], *STATIONS[n]), 1), "conj_of_north": [round(tn["lat"], 2), round(tn["lon"], 2)], "conj_of_north_to_south_km": round(gc_km(tn["lat"], tn["lon"], *STATIONS[s]), 1), "apex_re_south": round(ts["apex_re"], 3), "apex_re_north": round(tn["apex_re"], 3)} F = lambda code, y: elements(*igrf_at(model, y), *STATIONS[code])["F"] # noqa: E731 res["F_ratio_SJG_over_SMS"] = {str(y): round(F("SJG", y) / F("SMS", y), 4) for y in (1990.0, 2000.0, 2005.0, 2010.0, 2013.0, 2015.0, 2020.0, 2024.357, 2026.0, 2030.0)} res["F_rate_2013_2026_nT_per_yr"] = {c: round((F(c, 2026.0) - F(c, 2013.0)) / 13.0, 1) for c in ("SMS", "SJG")} res["F_pct_change_2013_2026"] = {c: round(100 * (F(c, 2026.0) / F(c, 2013.0) - 1), 2) for c in ("SMS", "SJG")} best = [r for r in res["anomaly"] if r["epoch"] == 2024.357][0] res["SMS_km_from_minimum_2024"] = round(gc_km(best["lat"], best["lon"], *STATIONS["SMS"]), 0) print(json.dumps({k: res[k] for k in ("pairs", "F_ratio_SJG_over_SMS", "F_rate_2013_2026_nT_per_yr", "F_pct_change_2013_2026", "SMS_km_from_minimum_2024")}, indent=1)) return res # ------------------------------------------------------------------ 10 May 2024 at every station def run_all_stations(root): T = T_PAPER model = read_igrf(os.path.join(root, "igrf14", "igrf14coeffs.txt")) g, h = igrf_at(model, EPOCH_2024) obs = read_capabilities(os.path.join(root, "gin-probe", "gin-capabilities.xml")) rows = [] t0, t1 = T - dt.timedelta(minutes=120), T + dt.timedelta(minutes=120) todo = [(c, o["lat"], o["lon"], "gin", o["name"]) for c, o in obs.items() if abs(o["lat"]) <= 62] todo += [(c, la, lo, "embrace", "EMBRACE") for c, (la, lo) in EMBRACE_COORDS.items()] for code, la, lo, src, name in todo: s = series(root, code, t0, t1, source=src) if not s: continue m = metrics(s, T) if m.get("excl"): rows.append({"code": code, "network": src, "excl": m["excl"]}) continue e = elements(g, h, la, lo) rows.append({"code": code, "network": src, "name": name[:30], "lat": round(la, 2), "lon": round(lo, 2), "LT": round((17 + 4 / 60 + lo / 15) % 24, 2), "F": round(e["F"]), "dip_lat": round(e["dip_lat"], 1), "A": round(m["A"], 1), "V": round(m["V"], 1), "Pc5": round(m["Pc5"], 1) if m["Pc5"] else None, "PA": round(m["PA"], 1)}) rows.sort(key=lambda r: r.get("LT", 99)) for r in rows: print(r) return rows # ------------------------------------------------------------------ putting it together FIXED_INPUTS = { "igrf14/igrf14coeffs.txt": ("https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt", "8f8d88403028fc4ee92c4f38d97b46e0a87e2cfc496045b43c9e26c1d6b0903c"), "wmm2025/WMM2025COF.zip": ("https://www.ncei.noaa.gov/sites/default/files/2025-01/WMM2025COF.zip", "999080a12c95122c2031aeac8be77229d97584530e87e0f1a6df7d946ecf90c5"), "gin-probe/gin-capabilities.xml": ("https://imag-data.bgs.ac.uk/GIN_V1/GINServices?Request=GetCapabilities", None), } def ensure_inputs(root): os.makedirs(root, exist_ok=True) for rel, (url, sha) in FIXED_INPUTS.items(): folder, name = rel.split("/") data = fetch(root, folder, name, url) if sha: same = hashlib.sha256(data).hexdigest() == sha print(f"{rel}: {'same file as ours' if same else 'differs from the file we used'}") for y in range(2013, 2023): fetch(root, "ebre-sc-lists", f"ssc_{y}_d.txt", f"https://www.obsebre.es/php/geomagnetisme/vrapides/ssc_{y}_d.txt") for y in range(2023, 2027): fetch(root, "ebre-sc-lists", f"ssc_{y}_p.txt", f"https://www.obsebre.es/php/geomagnetisme/vrapides/ssc_{y}_p.txt") def drop_days(d): """Per-day quiet lists stay out of the output: the medians are what the page uses.""" if isinstance(d, dict): return {k: drop_days(v) for k, v in d.items() if k != "days"} return d def main(): root = os.path.abspath(sys.argv[1] if len(sys.argv) > 1 else "sama_numbers_data") outdir = os.path.abspath(sys.argv[2] if len(sys.argv) > 2 else os.path.dirname(os.path.abspath(__file__))) sys.stdout.reconfigure(encoding="utf-8") os.makedirs(outdir, exist_ok=True) ensure_inputs(root) out = {"what": "Derived numbers for https://factchecked.science/sama-sudden-commencement/ , from sama_numbers.py. No magnetometer series.", "terms": "Numbers derived from EMBRACE/INPE, INTERMAGNET, Observatori de l'Ebre/ISGI, GFZ, IGRF-14 and WMM2025. " "The author(s) thank(s) the Embrace/INPE Program of the MCTI for providing the data for this work."} print("\n==== field and anomaly ====") out["field"] = run_anomaly(root) print("\n==== reproduction and filtering (test 8) ====") out["reproduction"] = run_repro(root) print("\n==== 10 May 2024 at every public station below 62 degrees of latitude ====") out["all_stations_10may2024"] = run_all_stations(root) print("\n==== event metrics ====") csvpath = os.path.join(outdir, "sama-event-metrics.csv") run_metrics(root, csvpath) print("\n==== tests 1-7 and 9 ====") tests = run_tests(root, csvpath) out["tests"] = drop_days(tests) print("\n==== checks added after the results ====") out["added_after_results"] = drop_days(run_added(root, csvpath, tests["test6_band"]["stations"])) with open(os.path.join(outdir, "sama-numbers.json"), "w", encoding="utf-8", newline="\n") as f: json.dump(out, f, indent=1, ensure_ascii=False) print("\nwrote", os.path.join(outdir, "sama-numbers.json"), "and", csvpath) if __name__ == "__main__": main()