#!/usr/bin/env python3 """How fast is Earth's magnetic field being lost? The presenter's figures against IGRF-14 and CHAOS-8. Fact-Checked Science, October 2026, for the Evidence review "field-loss-projections" ("Our analysis"). Standard library only. Run it python3 field_loss.py [folder] The folder (default: the script's own) must hold, or will receive, the inputs below. A missing file is downloaded from its maker; every file is then checked against the SHA-256 of the copy we used on October 27, 2026 (a different hash means the maker changed the file, so numbers may differ). The script prints a summary and writes field-loss-numbers.json into the folder. It takes a few minutes. Inputs igrf14coeffs.txt IGRF-14 Gauss coefficients 1900-2025 and secular variation 2025-30 (IAGA V-MOD), from https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt CHAOS-8.7_spline-coefficients.dat CHAOS-8.7 core-field model (Kloss et al., DTU Space, 30 September 2026, CC BY 4.0, DOI 10.5281/zenodo.22959079), from the "Latest release" page of https://www.spacecenter.dk/files/magnetic-models/CHAOS-8/ CHAOS-8.5_spline-coefficients.dat (optional) the release our page on geomagnetic jerks re-ran; used as a check. WMM2025COF.zip World Magnetic Model 2025 with NOAA's test values, to check the field evaluation and for WMM2025's own dipole forecast: https://www.ncei.noaa.gov/sites/default/files/2025-01/WMM2025COF.zip COV-OBS.x2-int COV-OBS.x2 mean model, 1840-2020 (Huder et al. 2020), for the field since 1840, worldwide and at the Antarctic core site behind the presenter's "55%": https://www.spacecenter.dk/files/magnetic-models/COV-OBSx2/COV-OBS.x2-int 2025-106_Sagnotti-et-al_data-csv.zip the data behind Sagnotti et al. (2026), GFZ Data Services, CC BY 4.0, doi:10.5880/fidgeo.2025.106: the HOLOANTA relative-paleointensity stack and core AV43 sample by sample, read as data from inside the zip: https://datapub.gfz.de/download/10.5880.FIDGEO.2025.106-CEani/2025-106_Sagnotti-et-al_data-csv.zip What it computes A. Arithmetic on the presenter's own stated figures (no model): the rates per decade and per century they imply, the field they imply today, and the date each rate would reach half or none of the field. B. The dipole moment, M = (4 pi a^3 / mu0) sqrt(g10^2 + g11^2 + h11^2), and the axial dipole |g10|, over the Swarm era (2014.0-2026.0) in IGRF-14 and CHAOS-8.7 (and 8.5), against what 5% per decade would have done. C. The field strength F at the WGS84 surface on a 1-degree grid in 2014.0 and 2025.0 (IGRF-14) and 2014.0 and 2026.0 (CHAOS-8.7, degree 13): the global mean, the largest local losses and gains in nT and in per cent, the share of the surface that lost 5% per decade or more, the Americas, and the four regions the presenter named (Canada's strong-field maximum, Siberia's, the South Atlantic Anomaly's minimum, the southern maximum). D. The fastest local rates of change: (dF/dt)/F worldwide in mid-2014 (IGRF-14's 2010-15 secular variation and CHAOS-8.7's derivative) and in IGRF-14's 2025-30 forecast; and the change over Swarm's first six months (1 January to 30 June 2014) in CHAOS-8.7, the period of the 2014 news stories. E. Forecasts: IGRF-14's and WMM2025's dipole for 2030; the pace needed to halve the axial dipole by the 2040s. F. Also: the South Atlantic Anomaly's area below 26,000 to 23,000 nT in 2014 and 2025-26; CHAOS-8.7's dipole weakening year by year, 1999.5-2026.5 (is it accelerating?); and the field since 1840 at the Ross Sea core site (74.141 S, 166.083 E) of the sediment record that the presenter read as "55% since the 1600s". G. The series behind the page's three charts: the dipole moment as a share of its 1850 value (COV-OBS.x2 for 1840-1895, IGRF-14 every five years from 1900 to 2025 and its forecast for 2030) beside the presenter's figures and their straight-line paces; CHAOS-8.7's yearly pace, with IGRF-14's five-year paces since 1900 and the 1840-2010 average for comparison; and core AV43's top samples against the modelled field at the core site, both as a share of their values in 1840. H. The 7 July 2026 claim that as much was lost in the last 80 years as in the 350 before, tested on the axial-dipole series of /next-excursion/ (IGRF-14 from 1900, 16 nT a year back to 1840, 2 nT a year from 1590, the trend of one reconstruction); and the 1590 field if the 1840-2010 pace had held since 1590, against a 55% loss since the 1600s. Conventions: Schmidt semi-normalised Legendre functions; IGRF coefficients interpolated linearly between epochs and extrapolated after 2025 with the 2025-30 secular variation, as IGRF is defined; CHAOS from its order-6 B-splines; X north, Y east, Z down; reference radius a = 6371.2 km. Per-decade rates are the change over the span divided by its length, as a share of the starting value. """ import hashlib import io import json import math import os import sys import urllib.request import zipfile URLS = { "igrf14coeffs.txt": "https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt", "CHAOS-8.7_spline-coefficients.dat": "https://www.spacecenter.dk/files/magnetic-models/CHAOS-8/_downloads/" "d51e04c49142b02b19f7544b6d32bec8/CHAOS-8.7_spline-coefficients.dat", "CHAOS-8.5_spline-coefficients.dat": "https://www.spacecenter.dk/files/magnetic-models/CHAOS-8/_downloads/" "bbcd3b0503d5a959b42c4aba3f7a4ba0/CHAOS-8.5_spline-coefficients.dat", "WMM2025COF.zip": "https://www.ncei.noaa.gov/sites/default/files/2025-01/WMM2025COF.zip", "COV-OBS.x2-int": "https://www.spacecenter.dk/files/magnetic-models/COV-OBSx2/COV-OBS.x2-int", "2025-106_Sagnotti-et-al_data-csv.zip": "https://datapub.gfz.de/download/10.5880.FIDGEO.2025.106-CEani/" "2025-106_Sagnotti-et-al_data-csv.zip", } HASHES = { "2025-106_Sagnotti-et-al_data-csv.zip": "78e3a9b0672eab3c7f48ecc061488a7728b62ee2ad79ccd7a6d023711fd14038", "COV-OBS.x2-int": "6c588becf4e3896b07be80c3f24697d299432120a7f4ab1139a914395f9d2a17", "igrf14coeffs.txt": "8f8d88403028fc4ee92c4f38d97b46e0a87e2cfc496045b43c9e26c1d6b0903c", "CHAOS-8.7_spline-coefficients.dat": "f759e725b763dd300cb97be07f44cda655226babde0cfe23f9e37ed7a425d203", "CHAOS-8.5_spline-coefficients.dat": "ff1805332fc9050d0932d4c43d26e88bea269a88e558369955b25d857a9d89e9", "WMM2025COF.zip": "999080a12c95122c2031aeac8be77229d97584530e87e0f1a6df7d946ecf90c5", } OPTIONAL = {"CHAOS-8.5_spline-coefficients.dat"} A_REF = 6371.2 # km WGS84_A = 6378.137 # km WGS84_F = 1.0 / 298.257223563 WGS84_E2 = WGS84_F * (2.0 - WGS84_F) WGS84_B = WGS84_A * (1.0 - WGS84_F) MU0_OVER_4PI = 1.0e-7 NMAX = 13 # IGRF-14 runs to 13 from 2000; CHAOS is truncated to 13 for like-for-like maps JUNE30_2014 = 2014.0 + 181.0 / 365.0 # 2014-07-01 00:00, the end of 30 June # ------------------------------------------------------------------ inputs def fetch(folder, name): path = os.path.join(folder, name) if not os.path.exists(path): print("downloading", URLS[name]) req = urllib.request.Request(URLS[name], headers={"User-Agent": "factchecked.science field_loss.py"}) try: with urllib.request.urlopen(req, timeout=120) as r, open(path, "wb") as f: f.write(r.read()) except Exception as e: # optional inputs may be skipped if name in OPTIONAL: print(" (skipped: %s)" % e) return None, None raise with open(path, "rb") as f: digest = hashlib.sha256(f.read()).hexdigest() if digest != HASHES[name]: print("Note: %s differs from the copy we used (SHA-256 %s...); results may differ." % (name, HASHES[name][:16])) return path, digest def zeros(): return [[0.0] * (NMAX + 1) for _ in range(NMAX + 1)] def read_igrf(path): epochs, rows = None, {} with open(path, encoding="ascii") as f: for line in f: if line.startswith("g/h"): epochs = [float(x) for x in line.split()[3:-1]] elif line[:1] in "gh" and not line.startswith("g/h"): p = line.split() rows[(p[0], int(p[1]), int(p[2]))] = [float(x) for x in p[3:]] return epochs, {k: v[:-1] for k, v in rows.items()}, {k: v[-1] for k, v in rows.items()} def igrf_at(model, t): epochs, main, sv = model g, h = zeros(), zeros() if t >= epochs[-1]: i0 = len(epochs) - 1 for (gh, n, m), v in main.items(): (g if gh == "g" else h)[n][m] = v[i0] + (t - epochs[-1]) * sv[(gh, n, m)] else: i0 = max(i for i, e in enumerate(epochs) if e <= t) w = (t - epochs[i0]) / (epochs[i0 + 1] - epochs[i0]) for (gh, n, m), v in main.items(): (g if gh == "g" else h)[n][m] = v[i0] + w * (v[i0 + 1] - v[i0]) return g, h def igrf_sv(model, t): """Secular variation in force at t: the slope of the linear segment (or the 2025-30 forecast).""" epochs, main, sv = model g, h = zeros(), zeros() if t >= epochs[-1]: for (gh, n, m), v in sv.items(): (g if gh == "g" else h)[n][m] = v else: i0 = max(i for i, e in enumerate(epochs) if e <= t) for (gh, n, m), v in main.items(): (g if gh == "g" else h)[n][m] = (v[i0 + 1] - v[i0]) / (epochs[i0 + 1] - epochs[i0]) return g, h class Chaos: """CHAOS core field from its B-spline file: the Gauss coefficients and their time derivatives.""" def __init__(self, path, name): nums, header = [], None with open(path, encoding="ascii") as f: for line in f: s = line.strip() if not s or s.startswith("#"): continue parts = s.split() if header is None: try: vals = [int(x) for x in parts] except ValueError: continue # a text title line, e.g. "CHAOS-8.5, extracted on ..." if len(vals) == 3: header = vals continue nums.extend(float(x) for x in parts) self.name = name self.nmax, self.nspl, order = header self.p = order - 1 self.knots = nums[:self.nspl + order] vals = nums[self.nspl + order:] n = self.nmax * (self.nmax + 2) assert len(vals) == self.nspl * n, "unexpected CHAOS file layout" self.coefs = [vals[j * n:(j + 1) * n] for j in range(self.nspl)] self.t_first, self.t_last = self.knots[self.p], self.knots[self.nspl] @classmethod def covobs(cls, path, name="COV-OBS.x2"): """COV-OBS.x2 mean model: title line, then nmax, nspl, order and the knots on one line, then values.""" self = cls.__new__(cls) lines = open(path, encoding="ascii").read().split("\n") head = lines[1].split() self.nmax, self.nspl, order = int(head[0]), int(head[1]), int(head[2]) self.knots = [float(x) for x in head[3:]] vals = [float(x) for line in lines[2:] for x in line.split()] n = self.nmax * (self.nmax + 2) self.coefs = [vals[j * n:(j + 1) * n] for j in range(self.nspl)] self.name, self.p = name, order - 1 self.t_first, self.t_last = self.knots[self.p], self.knots[self.nspl] return self def flat(self, t, deriv=0): i = span(self.nspl, self.p, t, self.knots) w = basis_derivs(i, t, self.p, deriv, self.knots)[deriv] out = [0.0] * len(self.coefs[0]) for j in range(self.p + 1): if w[j]: row = self.coefs[i - self.p + j] for q in range(len(out)): out[q] += w[j] * row[q] return out def gh(self, t, deriv=0, nmax=NMAX): c = self.flat(t, deriv) g, h = zeros(), zeros() for n in range(1, min(nmax, self.nmax) + 1): base = (n - 1) * (n + 1) g[n][0] = c[base] for m in range(1, n + 1): g[n][m] = c[base + 2 * m - 1] h[n][m] = c[base + 2 * m] return g, h def span(nb, p, u, U): if u >= U[nb]: return nb - 1 if u <= U[p]: return p lo, hi = p, nb mid = (lo + hi) // 2 while u < U[mid] or u >= U[mid + 1]: lo, hi = (lo, mid) if u < U[mid] else (mid, hi) mid = (lo + hi) // 2 return mid def basis_derivs(i, u, p, n, U): """Non-zero B-splines and their derivatives (Piegl and Tiller, The NURBS Book, algorithm A2.3).""" ndu = [[0.0] * (p + 1) for _ in range(p + 1)] left, right = [0.0] * (p + 1), [0.0] * (p + 1) ndu[0][0] = 1.0 for j in range(1, p + 1): left[j], right[j] = u - U[i + 1 - j], U[i + j] - u saved = 0.0 for r in range(j): ndu[j][r] = right[r + 1] + left[j - r] temp = ndu[r][j - 1] / ndu[j][r] ndu[r][j] = saved + right[r + 1] * temp saved = left[j - r] * temp ndu[j][j] = saved ders = [[0.0] * (p + 1) for _ in range(n + 1)] for j in range(p + 1): ders[0][j] = ndu[j][p] a = [[0.0] * (p + 1) for _ in range(2)] for r in range(p + 1): s1, s2 = 0, 1 a[0][0] = 1.0 for k in range(1, n + 1): d, rk, pk = 0.0, r - k, p - k if r >= k: a[s2][0] = a[s1][0] / ndu[pk + 1][rk] d = a[s2][0] * ndu[rk][pk] j1 = 1 if rk >= -1 else -rk j2 = k - 1 if (r - 1) <= pk else p - r for j in range(j1, j2 + 1): a[s2][j] = (a[s1][j] - a[s1][j - 1]) / ndu[pk + 1][rk + j] d += a[s2][j] * ndu[rk + j][pk] if r <= pk: a[s2][k] = -a[s1][k - 1] / ndu[pk + 1][r] d += a[s2][k] * ndu[r][pk] ders[k][r] = d s1, s2 = s2, s1 r = p for k in range(1, n + 1): for j in range(p + 1): ders[k][j] *= r r *= (p - k) return ders 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]) if n <= NMAX: 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 guide'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 class Row: def __init__(self, lat_deg, h_km=0.0): self.r, phic, self.delta = geodetic_to_geocentric(lat_deg, h_km) theta = math.pi / 2.0 - phic self.st = math.sin(theta) self.P, self.D = legendre(theta) self.rr = [(A_REF / self.r) ** (n + 2) for n in range(NMAX + 1)] def sums(self, g, h): out = [] for m in range(NMAX + 1): s = [0.0] * 6 for n in range(max(1, m), NMAX + 1): f, P, D = self.rr[n], self.P[n][m], self.D[n][m] s[0] += (n + 1) * f * g[n][m] * P s[1] += (n + 1) * f * h[n][m] * P s[2] += f * g[n][m] * D s[3] += f * h[n][m] * D s[4] += f * g[n][m] * P s[5] += f * h[n][m] * P out.append(s) return out def xyz(self, sums, cosm, sinm): br = bt = bp = 0.0 for m, (gr, hr, gt, ht, gp, hp) in enumerate(sums): c, s = cosm[m], sinm[m] br += gr * c + hr * s bt -= gt * c + ht * s bp -= m * (-gp * s + hp * c) bp /= self.st xc, zc = -bt, -br cd, sd = math.cos(self.delta), math.sin(self.delta) return xc * cd + zc * sd, bp, zc * cd - xc * sd def trig(lon_deg): lam = math.radians(lon_deg) return [math.cos(m * lam) for m in range(NMAX + 1)], [math.sin(m * lam) for m in range(NMAX + 1)] def field(g, h, lat, lon, hkm=0.0): row = Row(lat, hkm) c, s = trig(lon) return row.xyz(row.sums(g, h), c, s) def F_of(v): return math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) def validate(epoch, g, h, gd, hd, tests): worst, n = 0.0, 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) Xd, Yd, Zd = field(gd, hd, lat, lon, hkm) Fd = (X * Xd + Y * Yd + Z * Zd) / F for mine, ref in ((X, p[4]), (Y, p[5]), (Z, p[6]), (F, p[8]), (Xd, p[12]), (Yd, p[13]), (Zd, p[14]), (Fd, p[16])): worst = max(worst, abs(mine - ref)) n += 1 return {"test_points": n, "largest_abs_difference_nT_or_nT_per_yr": round(worst, 3), "pass_within_0.1": worst <= 0.1} # ------------------------------------------------------------------ the dipole def dipole(g, h): g10, g11, h11 = g[1][0], g[1][1], h[1][1] b0 = math.sqrt(g10 * g10 + g11 * g11 + h11 * h11) return {"M_1e22_Am2": (A_REF * 1e3) ** 3 * (b0 * 1e-9) / MU0_OVER_4PI / 1e22, "absg10_nT": abs(g10)} def pct(a, b): return 100.0 * (b - a) / a def per_decade(a, b, t0, t1): return pct(a, b) / (t1 - t0) * 10.0 # ------------------------------------------------------------------ the surface grid 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_AREA = 2 * math.pi * band_area(-90.0, 90.0) class Grid: def __init__(self, step=1.0): n_lat, n_lon = int(round(180 / step)), int(round(360 / step)) self.lat = [-90 + step / 2 + i * step for i in range(n_lat)] self.lon = [-180 + step / 2 + j * step for j in range(n_lon)] self.rows = [Row(la) for la in self.lat] self.trigs = [trig(lo) for lo in self.lon] self.area = [band_area(la - step / 2, la + step / 2) * math.radians(step) for la in self.lat] def vectors(self, g, h): out = [] for row in self.rows: s = row.sums(g, h) out.append([row.xyz(s, c, sn) for (c, sn) in self.trigs]) return out def Fgrid(V): return [[F_of(v) for v in row] for row in V] def rate_grid(V, Vd): """(dF/dt) in nT/yr from field vectors and their time derivatives.""" out = [] for row, rowd in zip(V, Vd): out.append([(v[0] * d[0] + v[1] * d[1] + v[2] * d[2]) / F_of(v) for v, d in zip(row, rowd)]) return out def mean_F(G, F): tot = w = 0.0 for i, row in enumerate(F): tot += sum(row) * G.area[i] w += len(row) * G.area[i] return tot / w def in_box(la, lo, box): la0, la1, lo0, lo1 = box return la0 <= la <= la1 and lo0 <= lo <= lo1 def extreme(G, F, box, want_max): best = None for i, la in enumerate(G.lat): for j, lo in enumerate(G.lon): if in_box(la, lo, box): v = F[i][j] if best is None or (v > best[0] if want_max else v < best[0]): best = (v, la, lo) return best def refine(g, h, lat0, lon0, want_max, half=1.0, step=0.05): best = None k = int(round(half / step)) lons = [lon0 + j * step for j in range(-k, k + 1)] trigs = [trig(lo) for lo in lons] for i in range(-k, k + 1): la = lat0 + i * step if abs(la) > 89.9: continue row = Row(la) s = row.sums(g, h) for lo, (c, sn) in zip(lons, trigs): v = F_of(row.xyz(s, c, sn)) if best is None or (v > best[0] if want_max else v < best[0]): best = (v, la, lo) return best REGIONS = { # (lat0, lat1, lon0, lon1) search boxes, east positive "Canada strong-field maximum": ((40.0, 85.0, -150.0, -50.0), True), "Siberia strong-field maximum": ((40.0, 85.0, 60.0, 160.0), True), "South Atlantic Anomaly minimum": ((-60.0, 10.0, -100.0, 30.0), False), "Southern strong-field maximum": ((-89.0, -40.0, -180.0, 180.0), True), } AMERICAS = (-60.0, 75.0, -170.0, -30.0) # a box over the Americas and adjacent oceans (not a land mask) def surface_change(G, gh0, gh1, t0, t1, label): """Field strength at the surface at two dates: global mean, local extremes, thresholds, regions.""" F0, F1 = Fgrid(G.vectors(*gh0)), Fgrid(G.vectors(*gh1)) yrs = t1 - t0 m0, m1 = mean_F(G, F0), mean_F(G, F1) rel = [[100.0 * (b - a) / a for a, b in zip(r0, r1)] for r0, r1 in zip(F0, F1)] d_nt = [[b - a for a, b in zip(r0, r1)] for r0, r1 in zip(F0, F1)] worst = min(((rel[i][j], G.lat[i], G.lon[j]) for i in range(len(G.lat)) for j in range(len(G.lon)))) best = max(((rel[i][j], G.lat[i], G.lon[j]) for i in range(len(G.lat)) for j in range(len(G.lon)))) worst_nt = min(((d_nt[i][j], G.lat[i], G.lon[j]) for i in range(len(G.lat)) for j in range(len(G.lon)))) thresholds = {} for per_dec in (5.0, 2.5, 1.0): lim = -per_dec * yrs / 10.0 a = sum(G.area[i] for i in range(len(G.lat)) for j in range(len(G.lon)) if rel[i][j] <= lim) thresholds["lost at least %.1f%% per decade" % per_dec] = round(100.0 * a / EARTH_AREA, 4) gain = sum(G.area[i] for i in range(len(G.lat)) for j in range(len(G.lon)) if rel[i][j] > 0) am_w = am_rel = 0.0 am_worst = None for i, la in enumerate(G.lat): for j, lo in enumerate(G.lon): if in_box(la, lo, AMERICAS): am_w += G.area[i] am_rel += rel[i][j] * G.area[i] if am_worst is None or rel[i][j] < am_worst[0]: am_worst = (rel[i][j], la, lo) regions = {} for name, (box, want_max) in REGIONS.items(): e0 = extreme(G, F0, box, want_max) e1 = extreme(G, F1, box, want_max) r0 = refine(*gh0, e0[1], e0[2], want_max) r1 = refine(*gh1, e1[1], e1[2], want_max) regions[name] = {"F_%s" % t0: round(r0[0], 1), "at_%s" % t0: [round(r0[1], 2), round(r0[2], 2)], "F_%s" % t1: round(r1[0], 1), "at_%s" % t1: [round(r1[1], 2), round(r1[2], 2)], "change_nT": round(r1[0] - r0[0], 1), "change_pct": round(pct(r0[0], r1[0]), 3)} return {"label": label, "from": t0, "to": t1, "global_mean_F_nT": [round(m0, 1), round(m1, 1)], "global_mean_change_pct": round(pct(m0, m1), 3), "global_mean_change_pct_per_decade": round(pct(m0, m1) / yrs * 10.0, 3), "largest_local_loss_pct": [round(worst[0], 3), worst[1], worst[2]], "largest_local_loss_pct_per_decade": round(worst[0] / yrs * 10.0, 3), "largest_local_gain_pct": [round(best[0], 3), best[1], best[2]], "largest_local_loss_nT": [round(worst_nt[0], 1), worst_nt[1], worst_nt[2]], "share_of_surface_pct": thresholds, "share_of_surface_gaining_pct": round(100.0 * gain / EARTH_AREA, 2), "americas_box_mean_change_pct": round(am_rel / am_w, 3), "americas_box_largest_loss_pct": [round(am_worst[0], 3), am_worst[1], am_worst[2]], "regions": regions} def rate_stats(G, gh, ghd, label): """Fastest local relative rate of change, (dF/dt)/F, worldwide.""" V, Vd = G.vectors(*gh), G.vectors(*ghd) F = Fgrid(V) R = rate_grid(V, Vd) rel = [[100.0 * r / f for r, f in zip(rr, ff)] for rr, ff in zip(R, F)] # % per year worst = min(((rel[i][j], G.lat[i], G.lon[j], R[i][j], F[i][j]) for i in range(len(G.lat)) for j in range(len(G.lon)))) worst_nt = min(((R[i][j], G.lat[i], G.lon[j], F[i][j]) for i in range(len(G.lat)) for j in range(len(G.lon)))) share = {} for per_dec in (5.0, 3.0, 2.0): a = sum(G.area[i] for i in range(len(G.lat)) for j in range(len(G.lon)) if rel[i][j] <= -per_dec / 10.0) share["falling at least %.0f%% per decade" % per_dec] = round(100.0 * a / EARTH_AREA, 4) meanrate = sum(sum(rr) * G.area[i] for i, rr in enumerate(R)) / sum(len(rr) * G.area[i] for i, rr in enumerate(R)) return {"label": label, "fastest_relative_fall_pct_per_decade": round(10 * worst[0], 3), "at": [worst[1], worst[2]], "there_dFdt_nT_per_yr": round(worst[3], 1), "there_F_nT": round(worst[4], 1), "fastest_fall_nT_per_yr": round(worst_nt[0], 1), "fastest_fall_at": [worst_nt[1], worst_nt[2]], "fastest_fall_there_pct_per_decade": round(1000 * worst_nt[0] / worst_nt[3], 3), "share_of_surface_pct": share, "global_mean_dFdt_nT_per_yr": round(meanrate, 2)} def six_months(G, ch): V0, V1 = G.vectors(*ch.gh(2014.0)), G.vectors(*ch.gh(JUNE30_2014)) F0, F1 = Fgrid(V0), Fgrid(V1) d = [[b - a for a, b in zip(r0, r1)] for r0, r1 in zip(F0, F1)] lo = min(((d[i][j], G.lat[i], G.lon[j], F0[i][j]) for i in range(len(G.lat)) for j in range(len(G.lon)))) hi = max(((d[i][j], G.lat[i], G.lon[j], F0[i][j]) for i in range(len(G.lat)) for j in range(len(G.lon)))) rel = min(((100 * d[i][j] / F0[i][j], G.lat[i], G.lon[j]) for i in range(len(G.lat)) for j in range(len(G.lon)))) yr = JUNE30_2014 - 2014.0 return {"span": [2014.0, round(JUNE30_2014, 4)], "largest_fall_nT": round(lo[0], 1), "fall_at": [lo[1], lo[2]], "fall_there_F_nT": round(lo[3], 1), "largest_rise_nT": round(hi[0], 1), "rise_at": [hi[1], hi[2]], "largest_relative_fall_pct": round(rel[0], 4), "rel_fall_at": [rel[1], rel[2]], "largest_relative_fall_as_pct_per_decade": round(rel[0] / yr * 10, 3), "largest_fall_as_nT_per_yr": round(lo[0] / yr, 1)} def saa_areas(G, gh, thresholds=(26000, 25000, 24000, 23000)): """Area of the weak region joined to the South Atlantic minimum below each level, % of Earth's surface.""" F = Fgrid(G.vectors(*gh)) nlat, nlon = len(G.lat), len(G.lon) i0, j0 = min(((i, j) for i in range(nlat) for j in range(nlon)), key=lambda p: F[p[0]][p[1]]) out = {} for t in thresholds: seen, stack, area = {(i0, j0)}, [(i0, j0)], 0.0 while stack: i, j = stack.pop() area += G.area[i] for di in (-1, 0, 1): for dj in (-1, 0, 1): ii, jj = i + di, (j + dj) % nlon if 0 <= ii < nlat and (ii, jj) not in seen and F[ii][jj] < t: seen.add((ii, jj)) stack.append((ii, jj)) out[str(t)] = round(100.0 * area / EARTH_AREA, 3) return out def yearly_dipole_paces(ch, t0=1999.5, t1=2026.5): """dM/dt from the spline's derivative, as % of M per century (positive = weakening).""" rows = [] t = t0 while t <= t1 + 1e-9: g, h = ch.gh(t, 0, 1) gd, hd = ch.gh(t, 1, 1) g10, g11, h11 = g[1][0], g[1][1], h[1][1] b0 = math.sqrt(g10 ** 2 + g11 ** 2 + h11 ** 2) db0 = (g10 * gd[1][0] + g11 * gd[1][1] + h11 * hd[1][1]) / b0 rows.append({"year": t, "M_pct_per_century_weakening": round(-100.0 * 100.0 * db0 / b0, 2), "g10_nT_per_yr_weakening": round(gd[1][0] if g10 < 0 else -gd[1][0], 2)}) # |g10| falling = positive t += 1.0 return rows def read_holoanta(zpath): """The HOLOANTA RPI stack (NRM/ARM) and core AV43 sample by sample (depth, RPI scaled to its maximum, age BP). Ages count back from AD 2000, as in Sagnotti et al. (2026).""" import csv z = zipfile.ZipFile(zpath) names = {os.path.basename(n): n for n in z.namelist()} def rows(name): return csv.reader(io.StringIO(z.read(names[name]).decode("utf-8-sig")), delimiter=";") stack = [] for r in rows("2025-106_Sagnotti-et-al_data_HoloAntaRPI.csv"): try: stack.append({"age_bp": float(r[1]), "rpi": float(r[2]), "sd": float(r[3]), "n": int(r[4])}) except (ValueError, IndexError): continue av43 = [] for r in rows("2025-106_Sagnotti-et-al_data_data.csv"): try: av43.append({"depth_cm": float(r[0]), "rpi": float(r[7]), "age_bp": float(r[9])}) except (ValueError, IndexError): continue return stack, av43 def site_series(models, lat, lon, times): """Field strength at one place through time, from several models (each within its span).""" out = [] for t in times: row = {"year": t} for name, fn in models: try: gh = fn(t) except ValueError: continue if gh is not None: row[name] = round(F_of(field(*gh, lat, lon)), 1) out.append(row) return out # ------------------------------------------------------------------ A. the presenter's own figures def presenter_arithmetic(): out = {} # 1. "10% ... in only about 150 years" (Disaster Cycle 00:38:12-00:38:18; mid-1800s to 2000 in the 22 Sep 2025 # video; 1859 to 2000 in the 10 Apr 2023 video) out["10% over ~150 years (to 2000)"] = {"per_decade_pct": round(10 / 15, 3), "per_century_pct": round(100 / 15, 2), "1859_to_2000_per_decade_pct": round(10 / 14.1, 3)} # 2. "5% loss per century" out["5% per century"] = {"per_decade_pct": 0.5} # 3. "5% per decade" and "15% total field loss" by 2010 = 10% + 5% in 2000-2010 out["5% per decade (from 2006 or 2010)"] = {"per_century_pct_linear": 50.0, "ratio_to_5_per_century": 10.0} # 4. "to lose that same 10% would take 20 years" (00:39:12-00:39:18): 0.5% a year out["10% in 20 years"] = {"per_decade_pct": 5.0} # 5. 22 Sep 2025: "We've lost 8% in the last 11 years" and "we should be down by about 30% now" r = 8 / 11 * 10 own = {"15% (2010) + 5%/decade to 2014 + 8% (2014-2025)": 15 + 5 * 0.4 + 8, "15% (2010) + 8%/11 years for 2010-2025": round(15 + 8 / 11 * 15, 2)} out["8% in 11 years (2014-2025)"] = {"per_decade_pct": round(r, 3), "per_century_pct": round(r * 10, 2), "his_inputs_give_loss_2025_pct": own, "he_says_2025_pct": 30} # 6. 15 Dec 2025: "about 30% since the 1800s, half in the last two decades" and "55% since the 1600s" out["30% since the 1800s, half in the last two decades"] = { "last_two_decades_per_decade_pct": 7.5, "earlier_15pct_over_1850_2005_per_decade_pct": round(15 / 15.5, 3)} out["55% since the 1600s (1 to 0.45)"] = { "per_century_pct_if_since_1625": round(55 / 4.0, 2), "per_century_pct_if_since_1675": round(55 / 3.5, 2), "field_then_relative_to_now": round(1 / 0.45, 3)} # 7. What each rate implies for the dates of half and none of the field, counting losses as shares of the # 1850s field, as the presenter does ("down by about 30%") dates = {} for label, start_year, start_loss, per_dec in ( ("his 2010 figures: 15% lost by 2010, then 5% per decade", 2010, 15.0, 5.0), ("his 2025 figures: 30% lost by 2025, then 8% per 11 years", 2025, 30.0, 80 / 11), ("his Dec 2025 figures: 30% by 2025, then 7.5% per decade", 2025, 30.0, 7.5)): dates[label] = {"loss_in_2026_pct": round(start_loss + per_dec * (2026 - start_year) / 10, 1), "loss_in_2045_pct": round(start_loss + per_dec * (2045 - start_year) / 10, 1), "half_gone_year": round(start_year + (50 - start_loss) / per_dec * 10, 1), "all_gone_year": round(start_year + (100 - start_loss) / per_dec * 10, 1)} out["projections on his own rates (straight lines)"] = dates # 8. What a field gone by the 2040s would need from his 2025 starting point (30% lost) need = {} for yr in (2040, 2045, 2049): need[str(yr)] = {"to_zero_pct_points_per_decade": round(70 / (yr - 2025) * 10, 1), "to_half_of_1850s_pct_points_per_decade": round(20 / (yr - 2025) * 10, 1)} out["what the 2040s would need, from his 30% in 2025"] = need return out # ------------------------------------------------------------------ H. the axial dipole back to 1590, as on /next-excursion/ def next_excursion_series(igrf): """The axial dipole |g10| as /next-excursion/ builds it: IGRF-14 from 1900 (linear between its epochs, as IGRF is defined); 1840 at the 2010 value plus 170 years at 16 nT a year (Finlay, Aubert and Gillet 2016's 1840-2010 mean), joined to 1900 by a straight line; and 1590-1840 at 2 nT a year, the trend of one reconstruction (Gubbins, Jones and Finlay 2006: 2.28 +/- 2.72 nT a year, rounded). Returns a function of the year.""" g2010 = dipole(*igrf_at(igrf, 2010.0))["absg10_nT"] g1900 = dipole(*igrf_at(igrf, 1900.0))["absg10_nT"] g1840 = g2010 + 170 * 16.0 def A(t): if t >= 1900.0: return dipole(*igrf_at(igrf, t))["absg10_nT"] if t >= 1840.0: return g1840 + (g1900 - g1840) * (t - 1840.0) / 60.0 return g1840 + 2.0 * (1840.0 - t) return A, g1840 def july_2026_check(igrf): """The 7 July 2026 video: on a chart he extended before 1900, 'the same amount lost in 80 years' as in the 350 years before, 'maybe even more' (00:06:27-00:06:46). Tested on /next-excursion/'s series, and what a 55% loss since the 1600s would need against the 1590 field at the post-1840 pace.""" A, g1840 = next_excursion_series(igrf) now = A(2025.0) l80, l350 = A(1945.0) - now, A(1595.0) - A(1945.0) l100, l335 = A(1925.0) - now, A(1590.0) - A(1925.0) out = {"basis": "axial dipole |g10|, nT: IGRF-14 from 1900; 1840 = 2010 + 170 yr x 16 nT/yr; 1840-1900 straight; " "1590-1840 at 2 nT/yr (Gubbins et al. 2006's trend, one reconstruction), as on /next-excursion/", "values_nT": {str(y): round(A(float(y)), 1) for y in (1590, 1595, 1840, 1900, 1925, 1945, 2025)}, "loss_1945_2025_nT": round(l80, 1), "loss_1595_1945_nT": round(l350, 1), "last_80_over_350_before": round(l80 / l350, 2), "loss_1925_2025_nT": round(l100, 1), "loss_1590_1925_nT": round(l335, 1), "last_100_over_335_before": round(l100 / l335, 2), "1590_over_2025_on_this_series": round(A(1590.0) / now, 3)} # if the decline had run at the 1840-2010 mean of 16 nT/yr since 1590 (a decline well before 1840), the earlier # spans lose more, and the 1590 field is still nowhere near the 2.2 times today's that a 55% loss would need g1590_16 = g1840 + 250 * 16.0 l350_16 = (g1840 + 245 * 16.0) - A(1945.0) out["if_16nT_per_yr_since_1590"] = {"1590_over_2025": round(g1590_16 / now, 3), "loss_1595_1945_nT": round(l350_16, 1), "last_80_over_350_before": round(l80 / l350_16, 2)} out["field_1590_needed_for_55pct_loss"] = round(1 / 0.45, 2) return out # ------------------------------------------------------------------ main def main(): sys.stdout.reconfigure(encoding="utf-8") folder = sys.argv[1] if len(sys.argv) > 1 else os.path.dirname(os.path.abspath(__file__)) paths, hashes = {}, {} for name in URLS: paths[name], hashes[name] = fetch(folder, name) out = {"what": "The presenter's field-loss figures against IGRF-14 and CHAOS-8 (factchecked.science, Our analysis)", "inputs_sha256": hashes} out["A_presenter_arithmetic"] = presenter_arithmetic() ew, gw, hw, gdw, hdw, tests = read_wmm(paths["WMM2025COF.zip"]) out["check_against_noaa_wmm2025_test_values"] = validate(ew, gw, hw, gdw, hdw, tests) print("Check against NOAA's WMM2025 test values:", out["check_against_noaa_wmm2025_test_values"]) igrf = read_igrf(paths["igrf14coeffs.txt"]) ch = Chaos(paths["CHAOS-8.7_spline-coefficients.dat"], "CHAOS-8.7") ch85 = Chaos(paths["CHAOS-8.5_spline-coefficients.dat"], "CHAOS-8.5") if paths.get("CHAOS-8.5_spline-coefficients.dat") else None out["chaos_valid_span"] = {m.name: [round(m.t_first, 3), round(m.t_last, 3)] for m in (ch, ch85) if m} # B. the dipole over the Swarm era dip = {} times = [1900.0, 2000.0, 2005.0, 2010.0, 2014.0, JUNE30_2014, 2015.0, 2020.0, 2024.0, 2025.0, 2026.0, 2026.5, 2030.0] for t in times: row = {"IGRF-14": dipole(*igrf_at(igrf, t))} if ch.t_first <= t <= min(ch.t_last, 2026.5): row["CHAOS-8.7"] = dipole(*ch.gh(t, 0, 1)) if ch85 and ch85.t_first <= t <= min(ch85.t_last, 2026.0): row["CHAOS-8.5"] = dipole(*ch85.gh(t, 0, 1)) dip[str(round(t, 4))] = {k: {kk: round(vv, 5) for kk, vv in v.items()} for k, v in row.items()} out["dipole"] = dip def span(model, t0, t1, key="M_1e22_Am2"): a, b = dip[str(round(t0, 4))][model][key], dip[str(round(t1, 4))][model][key] five = a * (1 - 0.05 * (t1 - t0) / 10) # 5% of the starting value each decade five_c = a * 0.95 ** ((t1 - t0) / 10) # 5% compounding each decade return {"from": t0, "to": t1, "start": a, "end": b, "change_pct": round(pct(a, b), 3), "per_decade_pct": round(per_decade(a, b, t0, t1), 3), "at_5pct_per_decade_end_would_be": round(five, 4), "compounding": round(five_c, 4), "at_5pct_per_decade_change_pct": round(pct(a, five), 2), "measured_as_share_of_5pct_per_decade": round(pct(a, b) / pct(a, five), 4)} swarm = {} for model in ("IGRF-14", "CHAOS-8.7", "CHAOS-8.5"): for t0, t1 in ((2014.0, 2025.0), (2014.0, 2026.0), (2014.0, 2024.0), (2015.0, 2025.0), (2000.0, 2010.0), (2005.0, 2025.0), (2025.0, 2030.0)): k0, k1 = str(round(t0, 4)), str(round(t1, 4)) if model in dip.get(k0, {}) and model in dip.get(k1, {}): swarm["%s %s-%s M" % (model, t0, t1)] = span(model, t0, t1) swarm["%s %s-%s |g10|" % (model, t0, t1)] = span(model, t0, t1, "absg10_nT") out["B_swarm_era_dipole"] = swarm w25, w30 = dipole(gw, hw), dipole([[gw[i][j] + 5 * gdw[i][j] for j in range(NMAX + 1)] for i in range(NMAX + 1)], [[hw[i][j] + 5 * hdw[i][j] for j in range(NMAX + 1)] for i in range(NMAX + 1)]) out["E_forecasts"] = { "IGRF-14 M 2025 -> 2030": [dip["2025.0"]["IGRF-14"]["M_1e22_Am2"], dip["2030.0"]["IGRF-14"]["M_1e22_Am2"]], "IGRF-14 M 2025-2030 per decade pct": swarm["IGRF-14 2025.0-2030.0 M"]["per_decade_pct"], "IGRF-14 |g10| 2025-2030 nT per yr": round((dip["2025.0"]["IGRF-14"]["absg10_nT"] - dip["2030.0"]["IGRF-14"]["absg10_nT"]) / 5, 2), "WMM2025 M 2025 -> 2030": [round(w25["M_1e22_Am2"], 4), round(w30["M_1e22_Am2"], 4)], "WMM2025 M per decade pct": round(per_decade(w25["M_1e22_Am2"], w30["M_1e22_Am2"], 2025, 2030), 3), } g10_2025 = dip["2025.0"]["IGRF-14"]["absg10_nT"] sv = out["E_forecasts"]["IGRF-14 |g10| 2025-2030 nT per yr"] halve = {} for yr in (2040, 2045, 2049): need = g10_2025 / 2 / (yr - 2025) halve[str(yr)] = {"nT_per_yr": round(need, 1), "times_IGRF14_forecast": round(need / sv, 1), "pct_of_2025_per_decade": round(100 * need * 10 / g10_2025, 1)} out["E_forecasts"]["halving the axial dipole from 2025 by"] = halve print("Dipole moment, 10^22 A m^2:") for t in ("2014.0", "2025.0", "2026.0"): print(" ", t, {k: v["M_1e22_Am2"] for k, v in dip[t].items()}) for k in ("IGRF-14 2014.0-2025.0 M", "IGRF-14 2014.0-2026.0 M", "CHAOS-8.7 2014.0-2026.0 M", "CHAOS-8.7 2014.0-2025.0 M"): s = swarm[k] print(" %-26s %+.3f%% (%.3f%% per decade); 5%% per decade would give %+.2f%%; measured = %.3f of that" % (k, s["change_pct"], s["per_decade_pct"], s["at_5pct_per_decade_change_pct"], s["measured_as_share_of_5pct_per_decade"])) # C and D. the surface print("Building the 1-degree grid ...") G = Grid(1.0) out["C_surface_IGRF14_2014_2025"] = surface_change(G, igrf_at(igrf, 2014.0), igrf_at(igrf, 2025.0), 2014.0, 2025.0, "IGRF-14") print("IGRF-14 2014-2025:", json.dumps({k: out["C_surface_IGRF14_2014_2025"][k] for k in ( "global_mean_change_pct", "largest_local_loss_pct", "largest_local_loss_pct_per_decade", "share_of_surface_pct", "americas_box_mean_change_pct")})) out["C_surface_CHAOS87_2014_2026"] = surface_change(G, ch.gh(2014.0), ch.gh(2026.0), 2014.0, 2026.0, "CHAOS-8.7, degree 13") out["C_surface_CHAOS87_2014_2025"] = surface_change(G, ch.gh(2014.0), ch.gh(2025.0), 2014.0, 2025.0, "CHAOS-8.7, degree 13") print("CHAOS-8.7 2014-2026:", json.dumps({k: out["C_surface_CHAOS87_2014_2026"][k] for k in ( "global_mean_change_pct", "largest_local_loss_pct", "largest_local_loss_pct_per_decade", "share_of_surface_pct")})) out["D_rates_IGRF14_2014.5"] = rate_stats(G, igrf_at(igrf, 2014.5), igrf_sv(igrf, 2014.5), "IGRF-14, 2010-15 secular variation, at 2014.5") out["D_rates_CHAOS87_2014.5"] = rate_stats(G, ch.gh(2014.5), ch.gh(2014.5, 1), "CHAOS-8.7 at 2014.5") out["D_rates_CHAOS87_2026.0"] = rate_stats(G, ch.gh(2026.0), ch.gh(2026.0, 1), "CHAOS-8.7 at 2026.0") out["D_rates_IGRF14_2025_2030"] = rate_stats(G, igrf_at(igrf, 2025.0), igrf_sv(igrf, 2025.0), "IGRF-14 forecast for 2025-30, at 2025.0") for k in ("D_rates_IGRF14_2014.5", "D_rates_CHAOS87_2014.5", "D_rates_CHAOS87_2026.0", "D_rates_IGRF14_2025_2030"): r = out[k] print("%-28s fastest relative fall %.2f%% per decade at %s; fastest fall %.1f nT/yr at %s (%.2f%% per decade); share of surface: %s" % (k, r["fastest_relative_fall_pct_per_decade"], r["at"], r["fastest_fall_nT_per_yr"], r["fastest_fall_at"], r["fastest_fall_there_pct_per_decade"], r["share_of_surface_pct"])) out["D_six_months_2014_CHAOS87"] = six_months(G, ch) print("CHAOS-8.7, 1 Jan - 30 Jun 2014:", out["D_six_months_2014_CHAOS87"]) # the South Atlantic Anomaly's area below several levels (cf. the presenter's "about 50%" and "about 15%") areas = {"IGRF-14 2014.0": saa_areas(G, igrf_at(igrf, 2014.0)), "IGRF-14 2025.0": saa_areas(G, igrf_at(igrf, 2025.0)), "CHAOS-8.7 2014.0": saa_areas(G, ch.gh(2014.0)), "CHAOS-8.7 2026.0": saa_areas(G, ch.gh(2026.0))} growth = {} for a, b in (("IGRF-14 2014.0", "IGRF-14 2025.0"), ("CHAOS-8.7 2014.0", "CHAOS-8.7 2026.0")): growth["%s -> %s" % (a, b)] = {t: {"pct_of_Earth": [areas[a][t], areas[b][t]], "relative_growth_pct": round(pct(areas[a][t], areas[b][t]), 1) if areas[a][t] else None} for t in areas[a]} out["C_saa_areas"] = growth print("South Atlantic Anomaly areas:", json.dumps(growth)) # is the dipole's decline accelerating in the Swarm era? CHAOS-8.7, year by year out["D_yearly_dipole_paces_CHAOS87"] = yearly_dipole_paces(ch) print("CHAOS-8.7 dipole weakening, % of M per century:", ", ".join("%g %.2f" % (r["year"], r["M_pct_per_century_weakening"]) for r in out["D_yearly_dipole_paces_CHAOS87"])) # the Antarctic core site of the paleointensity record behind the presenter's 55% (ANTA02-AV43, Wood Bay) cov = Chaos.covobs(paths["COV-OBS.x2-int"]) site = (-74.141, 166.083) def cov_gh(t): if not (max(cov.t_first, 1840.0) <= t <= min(cov.t_last, 2020.0)): return None return cov.gh(t) def igrf_gh(t): return igrf_at(igrf, t) if 1900.0 <= t <= 2030.0 else None ser = site_series([("COV-OBS.x2", cov_gh), ("IGRF-14", igrf_gh)], site[0], site[1], [1840.0 + 10.0 * k for k in range(19)] + [2014.0, 2025.0]) by = {r["year"]: r for r in ser} out["site_AV43"] = {"lat_lon": site, "series": ser, "COV-OBS.x2 1840->2000 pct": round(pct(by[1840.0]["COV-OBS.x2"], by[2000.0]["COV-OBS.x2"]), 2), "COV-OBS.x2 1840->2020 pct": round(pct(by[1840.0]["COV-OBS.x2"], by[2020.0]["COV-OBS.x2"]), 2), "IGRF-14 1900->2025 pct": round(pct(by[1900.0]["IGRF-14"], by[2025.0]["IGRF-14"]), 2), "COV-OBS.x2 dipole M 1840, 2000": [round(dipole(*cov.gh(1840.0))["M_1e22_Am2"], 4), round(dipole(*cov.gh(2000.0))["M_1e22_Am2"], 4)], "COV-OBS.x2 |g10| 1840 -> 2000 pct": round(pct(dipole(*cov.gh(1840.0))["absg10_nT"], dipole(*cov.gh(2000.0))["absg10_nT"]), 2), "COV-OBS.x2 valid span": [round(cov.t_first, 2), round(cov.t_last, 2)]} print("AV43 site:", json.dumps(out["site_AV43"])) # the measured loss since the presenter's start dates, from COV-OBS.x2 (1840-2020) joined to IGRF-14 (2025) since = {} end_igrf = dip["2025.0"]["IGRF-14"] for t0 in (1840.0, 1850.0, 1859.0): d0 = dipole(*cov.gh(t0)) d2000 = dipole(*cov.gh(2000.0)) d2010 = dipole(*cov.gh(2010.0)) since[str(t0)] = { "M_start": round(d0["M_1e22_Am2"], 4), "absg10_start_nT": round(d0["absg10_nT"], 1), "M_change_to_2000_pct": round(pct(d0["M_1e22_Am2"], d2000["M_1e22_Am2"]), 2), "M_change_to_2010_pct": round(pct(d0["M_1e22_Am2"], d2010["M_1e22_Am2"]), 2), "M_change_to_2025_pct (IGRF-14 end)": round(pct(d0["M_1e22_Am2"], end_igrf["M_1e22_Am2"]), 2), "absg10_change_to_2000_pct": round(pct(d0["absg10_nT"], d2000["absg10_nT"]), 2), "absg10_change_to_2010_pct": round(pct(d0["absg10_nT"], d2010["absg10_nT"]), 2), "absg10_change_to_2025_pct (IGRF-14 end)": round(pct(d0["absg10_nT"], end_igrf["absg10_nT"]), 2)} out["measured_loss_since_start_dates_COVOBSx2"] = since print("Measured loss since 1840/1850/1859 (COV-OBS.x2, then IGRF-14 for 2025):", json.dumps(since)) # the average field at the surface since then (COV-OBS.x2 to degree 13, then IGRF-14) meanF = {str(t): round(mean_F(G, Fgrid(G.vectors(*cov.gh(t)))), 1) for t in (1840.0, 1850.0, 1859.0, 2000.0, 2020.0)} meanF["2025.0 (IGRF-14)"] = out["C_surface_IGRF14_2014_2025"]["global_mean_F_nT"][1] meanF["change_1840_2025_pct"] = round(pct(meanF["1840.0"], meanF["2025.0 (IGRF-14)"]), 2) meanF["change_1859_2025_pct"] = round(pct(meanF["1859.0"], meanF["2025.0 (IGRF-14)"]), 2) meanF["change_2000_2025_pct"] = round(pct(meanF["2000.0"], meanF["2025.0 (IGRF-14)"]), 2) out["average_surface_field_since_1840"] = meanF print("Average surface field:", json.dumps(meanF)) # IGRF-14's five-year paces of the dipole moment, 1900-2030 (as on the site's guide), and the 1840-2010 average paces = [] for k in range(26): t0 = 1900.0 + 5 * k a, b = dipole(*igrf_at(igrf, t0))["M_1e22_Am2"], dipole(*igrf_at(igrf, t0 + 5))["M_1e22_Am2"] paces.append({"span": "%d-%d" % (t0, t0 + 5), "M_pct_per_century_weakening": round(100.0 * (a - b) / a * 20.0, 2)}) since45 = [p for p in paces if int(p["span"][:4]) >= 1945 and int(p["span"][:4]) < 2025] fastest = max(since45, key=lambda p: p["M_pct_per_century_weakening"]) m1840, m2010 = dipole(*cov.gh(1840.0))["M_1e22_Am2"], dipole(*cov.gh(2010.0))["M_1e22_Am2"] out["paces_for_comparison"] = { "IGRF-14 five-year paces of M, % per century": paces, "fastest five-year pace 1945-2025": fastest, "COV-OBS.x2 M 1840-2010 average, % per century": round(100.0 * (m1840 - m2010) / m1840 / 1.7, 2), "CHAOS-8.7 mean yearly pace 2000.5-2013.5": round(sum(r["M_pct_per_century_weakening"] for r in out["D_yearly_dipole_paces_CHAOS87"] if 2000.0 < r["year"] < 2014.0) / 14, 2), "CHAOS-8.7 mean yearly pace 2014.5-2025.5": round(sum(r["M_pct_per_century_weakening"] for r in out["D_yearly_dipole_paces_CHAOS87"] if 2014.0 < r["year"] < 2026.0) / 12, 2)} print("Paces for comparison:", json.dumps({k: v for k, v in out["paces_for_comparison"].items() if not k.startswith("IGRF-14 five")})) # H. the July 2026 claim (80 against 350 years), and the 1590 field against a 55% loss since the 1600s out["H_july_2026_and_1590"] = july_2026_check(igrf) print("July 2026 check:", json.dumps(out["H_july_2026_and_1590"])) # G. the charts' series: COV-OBS.x2 for 1840-1895, then IGRF-14 every five years (as on the guide) and its 2030 forecast m1850 = dipole(*cov.gh(1850.0))["M_1e22_Am2"] measured = [{"year": 1840.0 + 5 * k, "pct_of_1850": round(100 * dipole(*cov.gh(1840.0 + 5 * k))["M_1e22_Am2"] / m1850, 2), "model": "COV-OBS.x2"} for k in range(12)] measured += [{"year": 1900.0 + 5 * k, "pct_of_1850": round(100 * dipole(*igrf_at(igrf, 1900.0 + 5 * k))["M_1e22_Am2"] / m1850, 2), "model": "IGRF-14" + (" forecast" if k > 25 else "")} for k in range(27)] presenter_points = [{"year": 2000, "pct": 90, "said": "10% lost by 2000"}, {"year": 2010, "pct": 85, "said": "15% lost by 2010"}, {"year": 2025, "pct": 70, "said": "about 30% lost by 2025"}] lines = {"his 2010 pace: 5 points a decade from 85% in 2010": [[2010, 85.0], [2060, 85.0 - 5.0 * 5]], "his 2025 pace: 8 points per 11 years from 70% in 2025": [[2025, 70.0], [2060, round(70.0 - 80 / 11 * 3.5, 2)]]} out["chart_A_share_of_1850_field"] = {"measured_dipole_moment": measured, "presenter_points": presenter_points, "presenter_lines": lines, "M_1850_1e22_Am2": round(m1850, 4), "half_of_2025_dipole_pct_of_1850": round(50 * dipole(*igrf_at(igrf, 2025.0))["M_1e22_Am2"] / m1850, 2)} stack, av43 = read_holoanta(paths["2025-106_Sagnotti-et-al_data-csv.zip"]) top = [s for s in av43 if s["age_bp"] <= 300.0] ref = min(top, key=lambda s: abs(s["age_bp"] - 160.0)) # the sample nearest 1840 (160 years before 2000) site_1840 = by[1840.0]["COV-OBS.x2"] out["chart_C_antarctic_record"] = { "av43_top_samples": [{"year_ad": round(2000.0 - s["age_bp"], 1), "depth_cm": s["depth_cm"], "rpi_scaled_to_max": s["rpi"], "share_of_1840_sample": round(s["rpi"] / ref["rpi"], 3)} for s in top], "av43_reference_sample": {"year_ad": round(2000.0 - ref["age_bp"], 1), "rpi": ref["rpi"]}, "holoanta_stack_0_to_600_bp": [s for s in stack if s["age_bp"] <= 600.0], "site_field_share_of_1840": [{"year": r["year"], "model": "COV-OBS.x2" if "COV-OBS.x2" in r else "IGRF-14", "share": round((r.get("COV-OBS.x2") or r.get("IGRF-14")) / site_1840, 4)} for r in sorted(ser, key=lambda r: r["year"]) if (r.get("COV-OBS.x2") or (r["year"] > 2020 and r.get("IGRF-14")))]} print("AV43 top samples:", json.dumps(out["chart_C_antarctic_record"]["av43_top_samples"])) print("HOLOANTA stack, 0-600 BP:", json.dumps(out["chart_C_antarctic_record"]["holoanta_stack_0_to_600_bp"])) with open(os.path.join(folder, "field-loss-numbers.json"), "w", encoding="utf-8", newline="\n") as f: json.dump(out, f, indent=1) print("Wrote", os.path.join(folder, "field-loss-numbers.json")) if __name__ == "__main__": main()