"""Where Earth's core field accelerates, 1900-2025: a re-run of two public field models. Fact-Checked Science, October 2026, for the page "The core field accelerates most near the equator, in patterns that persist, not collapse" (/jerks-secular-acceleration/). It reads two model files published by their makers and prints the numbers the page takes from them (a few, such as percentage changes, ratios and distances, follow from the printed values by arithmetic); with --grids DIR it also writes the time-RMS grids from which the page's maps and latitude chart are drawn: CHAOS-8.5_spline-coefficients.dat CHAOS-8 geomagnetic field model, version 8.5 (Kloss et al., Zenodo, 2 February 2026, CC BY 4.0, DOI 10.5281/zenodo.18432111), distributed by DTU Space from the CHAOS-8.5 section of https://www.spacecenter.dk/files/magnetic-models/CHAOS-8/archive.html COV-OBS.x2-int COV-OBS.x2 mean model (Huder et al. 2020; the paper is CC BY 4.0, the download page states no licence), from https://www.spacecenter.dk/files/magnetic-models/COV-OBSx2/COV-OBS.x2-int Put both beside this script (or give their paths as arguments) and run python3 sa_cmb.py or python3 sa_cmb.py --grids grids It needs only Python 3 (no NumPy) and takes about half a minute. What it computes: the second time derivative (secular acceleration, SA) of the radial field at the core-mantle boundary, radius 3,485 km, from the models' B-spline Gauss coefficients, truncated at degree 9 (and 13 as a check) on 2-degree cells; its RMS over time in each period; where it is strongest; its share near the equator; its strength beneath named places and inside the Bermuda Triangle; its strength beneath the South Atlantic Anomaly however the anomaly is outlined; its global power by degree band; and the surface field strength along the equator and worldwide. """ import hashlib import math import os import sys A, C = 6371.2, 3485.0 # reference radius and core-mantle boundary radius, km STEP = 2.0 # grid cell size, degrees EXPECTED = { # SHA-256 of the files we used (fetched 25 October 2026) "CHAOS-8.5_spline-coefficients.dat": "ff1805332fc9050d0932d4c43d26e88bea269a88e558369955b25d857a9d89e9", "COV-OBS.x2-int": "6c588becf4e3896b07be80c3f24697d299432120a7f4ab1139a914395f9d2a17", } PLACES = [ # (name, latitude, longitude); east positive ("Bermuda Triangle, centroid of Miami-Bermuda-San Juan", 25.5, -70.4), ("San Juan, Puerto Rico", 18.47, -66.11), ("Dragon's Triangle, about 30N 141E (Izu-Ogasawara)", 30.0, 141.0), ("Old West Pacific Anomaly core, about 13N 115E", 13.0, 115.0), ] BERMUDA = [(25.76, -80.19), (32.30, -64.78), (18.47, -66.11)] # Miami, Bermuda, San Juan (Gaddis) LAT_BANDS = [(0, 10), (10, 20), (20, 30), (30, 45), (45, 60), (60, 90)] SAA_OUTLINES = (24000, 26000, 28000, 32000) # surface F below this, nT; 32,000 is # Campuzano et al. 2021's definition # ---------- the models ------------------------------------------------------------------- class Model: def __init__(self, name, knots, order, nmax, coefs): self.name, self.knots, self.p, self.nmax, self.coefs = name, knots, order - 1, nmax, coefs self.nb = len(coefs) def gauss(self, t, deriv=0): """Gauss coefficients g10, g11, h11, g20, ... (deriv 0), their SV (1) or SA (2) at year t.""" i = span(self.nb, 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 read_chaos(path): lines = open(path, encoding="ascii").read().split("\n") nmax, nspl, order = (int(x) for x in lines[1].split()) knots = [float(x) for x in lines[2:2 + nspl + order]] vals = [float(x) for line in lines[2 + nspl + order:] for x in line.split()] n = nmax * (nmax + 2) return Model("CHAOS-8.5", knots, order, nmax, [vals[j * n:(j + 1) * n] for j in range(nspl)]) def read_covobs(path): lines = open(path, encoding="ascii").read().split("\n") head = lines[1].split() nmax, nspl, order = int(head[0]), int(head[1]), int(head[2]) knots = [float(x) for x in head[3:]] vals = [float(x) for line in lines[2:] for x in line.split()] n = nmax * (nmax + 2) return Model("COV-OBS.x2", knots, order, nmax, [vals[j * n:(j + 1) * n] for j in range(nspl)]) 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, 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 # ---------- spherical harmonics ---------------------------------------------------------- def idx(n, m, kind): base = (n - 1) * (n + 1) return base if m == 0 else base + 2 * m - 1 + (kind == "h") def schmidt(nmax, theta): """Schmidt semi-normalized P_n^m(cos theta) and dP/dtheta.""" c, s = math.cos(theta), math.sin(theta) P, dP = {(0, 0): 1.0}, {(0, 0): 0.0} for m in range(nmax + 1): if m == 1: P[(1, 1)], dP[(1, 1)] = s, c elif m > 1: f = math.sqrt((2 * m - 1) / (2 * m)) P[(m, m)] = f * s * P[(m - 1, m - 1)] dP[(m, m)] = f * (c * P[(m - 1, m - 1)] + s * dP[(m - 1, m - 1)]) for n in range(max(m + 1, 1), nmax + 1): a1 = (2 * n - 1) / math.sqrt(n * n - m * m) a2 = math.sqrt(((n - 1) ** 2 - m * m) / (n * n - m * m)) if n - 1 >= m else 0.0 P[(n, m)] = a1 * c * P[(n - 1, m)] - a2 * P.get((n - 2, m), 0.0) dP[(n, m)] = a1 * (c * dP[(n - 1, m)] - s * P[(n - 1, m)]) - a2 * dP.get((n - 2, m), 0.0) return P, dP class Grid: def __init__(self, step, nmax): self.lats = [-90 + step / 2 + i * step for i in range(int(round(180 / step)))] self.lons = [-180 + step / 2 + j * step for j in range(int(round(360 / step)))] self.P = [schmidt(nmax, math.radians(90 - la))[0] for la in self.lats] self.cos = [[math.cos(m * math.radians(lo)) for lo in self.lons] for m in range(nmax + 1)] self.sin = [[math.sin(m * math.radians(lo)) for lo in self.lons] for m in range(nmax + 1)] self.w = [math.cos(math.radians(la)) for la in self.lats] def radial(self, g, r, nmax): out = [] for P in self.P: Am, Bm = [0.0] * (nmax + 1), [0.0] * (nmax + 1) for n in range(1, nmax + 1): f = (n + 1) * (A / r) ** (n + 2) for m in range(n + 1): Am[m] += g[idx(n, m, "g")] * P[(n, m)] * f if m: Bm[m] += g[idx(n, m, "h")] * P[(n, m)] * f out.append([sum(Am[m] * self.cos[m][j] + Bm[m] * self.sin[m][j] for m in range(nmax + 1)) for j in range(len(self.lons))]) return out def mean_square_radial(g, r, n0, n1): tot = 0.0 for n in range(n0, n1 + 1): s = sum(g[idx(n, m, "g")] ** 2 + (g[idx(n, m, "h")] ** 2 if m else 0.0) for m in range(n + 1)) tot += (n + 1) ** 2 / (2 * n + 1) * (A / r) ** (2 * n + 4) * s return tot def intensity_row(g, nmax, lat, lons, r=A): """Total field F along one latitude at radius r.""" theta = math.radians(90 - lat) P, dP = schmidt(nmax, theta) st = math.sin(theta) out = [] for lon in lons: phi = math.radians(lon) br = bt = bp = 0.0 for n in range(1, nmax + 1): f = (A / r) ** (n + 2) for m in range(n + 1): gg = g[idx(n, m, "g")] hh = g[idx(n, m, "h")] if m else 0.0 cm, sm = math.cos(m * phi), math.sin(m * phi) br += (n + 1) * f * (gg * cm + hh * sm) * P[(n, m)] bt -= f * (gg * cm + hh * sm) * dP[(n, m)] if st > 1e-12: bp += f * m * (gg * sm - hh * cm) * P[(n, m)] / st out.append(math.sqrt(br * br + bt * bt + bp * bp)) return out # ---------- the analysis ----------------------------------------------------------------- def epochs(t0, t1, dt): return [t0 + k * dt for k in range(int(round((t1 - t0) / dt)) + 1)] def time_rms(model, ts, grid, nmax): ss = [[0.0] * len(grid.lons) for _ in grid.lats] for t in ts: for i, row in enumerate(grid.radial(model.gauss(t, 2), C, nmax)): for j, v in enumerate(row): ss[i][j] += v * v return [[math.sqrt(x / len(ts)) for x in row] for row in ss] def area_rms(M, g, cells=None): num = den = 0.0 for i in range(len(M)): for j in range(len(M[i])): if cells is None or (i, j) in cells: num += g.w[i] * M[i][j] ** 2 den += g.w[i] return math.sqrt(num / den) def lat_shares(M, g, limits=(10, 20, 30)): tot = sum(g.w[i] * v * v for i in range(len(M)) for v in M[i]) return [(L, sum(g.w[i] * v * v for i in range(len(M)) for v in M[i] if abs(g.lats[i]) <= L) / tot, math.sin(math.radians(L))) for L in limits] def peaks(M, g, top): found = [] for i in range(len(M)): for j in range(len(M[0])): v = M[i][j] if all(v >= M[i + di][(j + dj) % len(M[0])] for di in (-1, 0, 1) for dj in (-1, 0, 1) if 0 <= i + di < len(M)): found.append((v, g.lats[i], g.lons[j])) return sorted(found, reverse=True)[:top] def percentile(M, g, lat, lon): i = min(range(len(g.lats)), key=lambda k: abs(g.lats[k] - lat)) j = min(range(len(g.lons)), key=lambda k: abs(((g.lons[k] - lon + 180) % 360) - 180)) v = M[i][j] below = sum(g.w[a] for a in range(len(M)) for b in range(len(M[0])) if M[a][b] < v) return v, 100 * below / sum(g.w[a] * len(M[a]) for a in range(len(M))) def coverage(M, g, T): return sum(g.w[i] for i in range(len(M)) for v in M[i] if v > T) / sum(g.w[i] * len(M[i]) for i in range(len(M))) def band_rms(model, ts, bands): vals = [[math.sqrt(mean_square_radial(model.gauss(t, 2), C, a, b)) for a, b in bands] for t in ts] return [math.sqrt(sum(v[k] ** 2 for v in vals) / len(vals)) for k in range(len(bands))] def band_shares(M, g, bands=LAT_BANDS): """Share of the SA power (area-weighted square of the time-RMS) and of the area in each |lat| band.""" tp = sum(g.w[i] * v * v for i in range(len(M)) for v in M[i]) ta = sum(g.w[i] * len(M[i]) for i in range(len(M))) return [(a, b, sum(g.w[i] * v * v for i in range(len(M)) for v in M[i] if a <= abs(g.lats[i]) < b) / tp, sum(g.w[i] * len(M[i]) for i in range(len(M)) if a <= abs(g.lats[i]) < b) / ta) for a, b in bands] def in_triangle(lat, lon, tri=BERMUDA): (y1, x1), (y2, x2), (y3, x3) = tri d = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3) a = ((y2 - y3) * (lon - x3) + (x3 - x2) * (lat - y3)) / d b = ((y3 - y1) * (lon - x3) + (x1 - x3) * (lat - y3)) / d return a >= 0 and b >= 0 and 1 - a - b >= 0 def bermuda(M, g): """Range of the time-RMS SA over the cells inside the Bermuda Triangle, and whether the core's strongest cell lies inside it.""" inside = [(M[i][j], g.lats[i], g.lons[j]) for i in range(len(M)) for j in range(len(M[0])) if in_triangle(g.lats[i], g.lons[j])] top = max((M[i][j], g.lats[i], g.lons[j]) for i in range(len(M)) for j in range(len(M[0]))) hi = max(inside) return (f"inside the Bermuda Triangle ({len(inside)} cells): {ns(min(inside)[0])} to {ns(hi[0])}, the highest at " f"{hi[1]:+.0f} {hi[2]:+.0f}, {hi[0] / area_rms(M, g):.2f} x the global level; the core's strongest cell " f"({ns(top[0])} at {top[1]:+.0f} {top[2]:+.0f}) is {'INSIDE' if in_triangle(top[1], top[2]) else 'outside'} it") def write_grid(M, g, folder, name): import csv os.makedirs(folder, exist_ok=True) with open(os.path.join(folder, name), "w", newline="", encoding="utf-8") as fh: w = csv.writer(fh) w.writerow(["lat\\lon"] + [f"{x:g}" for x in g.lons]) for la, row in zip(g.lats, M): w.writerow([f"{la:g}"] + [f"{v:.1f}" for v in row]) def gridname(model, a, b, n): return f"rms_{model.name.replace('.', '')}_{a}-{b}_N{n}.csv" def ns(x): return f"{x:,.0f}" def main(): here = os.path.dirname(os.path.abspath(__file__)) args = sys.argv[1:] grids = None if "--grids" in args: k = args.index("--grids") grids = args[k + 1] del args[k:k + 2] paths = args[:2] if len(args) >= 2 else [os.path.join(here, f) for f in EXPECTED] for p, (name, sha) in zip(paths, EXPECTED.items()): if not os.path.exists(p): sys.exit(f"Missing {name}: download it as the header of this script says and put it beside the script.") got = hashlib.sha256(open(p, "rb").read()).hexdigest() print(f"{name}: SHA-256 {'matches' if got == sha else 'DIFFERS from'} the file we used") chaos, cov = read_chaos(paths[0]), read_covobs(paths[1]) g9, g13 = Grid(STEP, 9), Grid(STEP, 13) print("\n1. CHAOS-8.5, 2000-2025: time-RMS of radial SA at the core-mantle boundary (nT/yr^2)") ch = epochs(2000.0, 2025.0, 0.25) M9, M13 = time_rms(chaos, ch, g9, 9), time_rms(chaos, ch, g13, 13) for label, M, g in (("n<=9", M9, g9), ("n<=13", M13, g13)): print(f" {label}: area-RMS {ns(area_rms(M, g))}; share of SA power within " + ", ".join(f"+-{L} deg {s:.0%} (area {a:.0%})" for L, s, a in lat_shares(M, g))) print(" strongest patches, n<=9:", "; ".join(f"{la:+.0f} {lo:+.0f} {ns(v)}" for v, la, lo in peaks(M9, g9, 10))) print(" strongest patches, n<=13:", "; ".join(f"{la:+.0f} {lo:+.0f} {ns(v)}" for v, la, lo in peaks(M13, g13, 10))) for name, la, lo in PLACES: v9, p9 = percentile(M9, g9, la, lo) v13, p13 = percentile(M13, g13, la, lo) print(f" beneath {name}: {ns(v9)}, stronger than {p9:.1f}% of the core surface; n<=13 {ns(v13)}, {p13:.1f}%") print(" n<=9 " + bermuda(M9, g9)) print(" n<=13 " + bermuda(M13, g13)) top9 = peaks(M9, g9, 1)[0][0] print(f" the strongest patch, n<=9: {top9 / area_rms(M9, g9):.2f} x the global level") print(" share of SA power (and of the area) by latitude band, n<=9: " + "; ".join(f"{a}-{b} deg {p:.1%} ({ar:.1%})" for a, b, p, ar in band_shares(M9, g9))) if grids: write_grid(M9, g9, grids, gridname(chaos, 2000, 2025, 9)) write_grid(M13, g13, grids, gridname(chaos, 2000, 2025, 13)) print("\n2. The South Atlantic Anomaly and the West Pacific in 2025 (CHAOS-8.5, surface, n<=13)") g25 = chaos.gauss(2025.0) lats1, lons1 = [-89.5 + k for k in range(180)], [-179.5 + k for k in range(360)] F = [intensity_row(g25, 13, la, lons1) for la in lats1] mins = sorted((F[i][j], lats1[i], lons1[j]) for i in range(1, 179) for j in range(360) if all(F[i][j] <= F[i + di][(j + dj) % 360] for di in (-1, 0, 1) for dj in (-1, 0, 1))) print(" surface intensity minima:", "; ".join(f"{ns(v)} nT at {la:+.1f} {lo:+.1f}" for v, la, lo in mins[:3])) wp = [F[i][j] for i, la in enumerate(lats1) for j, lo in enumerate(lons1) if -30 <= la <= 40 and 100 <= lo <= 180] print(f" West Pacific, 30S-40N 100E-180: F from {ns(min(wp))} to {ns(max(wp))} nT (no minimum there)") F2 = [intensity_row(g25, 13, la, g9.lons) for la in g9.lats] glob = area_rms(M9, g9) for thr in SAA_OUTLINES: cells = {(i, j) for i in range(len(g9.lats)) for j in range(len(g9.lons)) if F2[i][j] < thr} share = sum(g9.w[i] for i, j in cells) / sum(g9.w[i] * len(g9.lons) for i in range(len(g9.lats))) for label, M, g in (("n<=9", M9, g9), ("n<=13", M13, g13)): print(f" SA beneath F<{ns(thr)} nT ({share:.1%} of the area, {label}): {ns(area_rms(M, g, cells))}, " f"{area_rms(M, g, cells) / area_rms(M, g):.2f} x the global level") box = {(i, j) for i, la in enumerate(g9.lats) for j, lo in enumerate(g9.lons) if -70 <= la <= -20 and -90 <= lo <= 70} print(f" SA in Campuzano et al. 2021's SAA box, 20-70S 90W-70E (n<=9): {area_rms(M9, g9, box) / glob:.2f} x global") print("\n3. COV-OBS.x2 by period, and CHAOS-8.5 beside it (n<=9 unless stated)") bands = ((1, 4), (5, 9), (10, 14)) runs = [(cov, a, b, 0.5) for a, b in ((1900, 1925), (1925, 1950), (1950, 1975), (1975, 2000), (2000, 2020))] runs.append((chaos, 2000, 2020, 0.25)) for model, a, b, dt in runs: ts = epochs(float(a), float(b), dt) M = time_rms(model, ts, g9, 9) tb = [t for t in (epochs(float(a), float(b), dt)) if t < b or b in (2020, 2025)] r4, r9, r14 = band_rms(model, tb, bands) pk = peaks(M, g9, 2) print(f" {model.name} {a}-{b}: area-RMS {ns(area_rms(M, g9))}; bands n1-4 {ns(r4)}, n5-9 {ns(r9)}, " f"n10-14 {ns(r14)}; area above 250: {coverage(M, g9, 250):.0%}; within +-10 deg " f"{lat_shares(M, g9, (10,))[0][1]:.0%}; strongest " + ", ".join(f"{la:+.0f} {lo:+.0f} ({ns(v)})" for v, la, lo in pk)) print(" stronger than this share of the core surface: " + "; ".join(f"{name.split(',')[0]} {percentile(M, g9, la, lo)[1]:.1f}%" for name, la, lo in PLACES)) print(" " + bermuda(M, g9)) if grids: write_grid(M, g9, grids, gridname(model, a, b, 9)) print("\n4. Global SA power by scale (area-RMS at the CMB, nT/yr^2)") for a in range(1900, 2020, 10): r4, r9 = band_rms(cov, [t for t in epochs(float(a), float(a + 10), 0.5) if t < a + 10], bands[:2]) print(f" COV-OBS.x2 {a}s: n1-4 {ns(r4)}, n5-9 {ns(r9)}") for a in range(2000, 2025, 5): ts = [t for t in epochs(float(a), float(a + 5), 0.25) if t < a + 5 or a == 2020] r4, r9, r14 = band_rms(chaos, ts, bands) print(f" CHAOS-8.5 {a}-{a + 5}: n1-4 {ns(r4)}, n5-9 {ns(r9)}, n1-9 {ns(math.hypot(r4, r9))}") ser = [(t, math.sqrt(mean_square_radial(chaos.gauss(t, 2), C, 1, 9))) for t in ch] print(" CHAOS-8.5 pulses (local maxima of global SA power, n<=9):", ", ".join(f"{t:.2f}" for k, (t, v) in enumerate(ser[1:-1], 1) if v >= ser[k - 1][1] and v >= ser[k + 1][1])) print("\n5. Surface field strength (nT)") for model, t, nmax in ((cov, 1900.0, 14), (cov, 2020.0, 14), (chaos, 2000.0, 13), (chaos, 2025.0, 13)): g = model.gauss(t) eq = intensity_row(g, nmax, 0.0, list(range(-180, 180, 2))) num = den = 0.0 for la in range(-88, 90, 4): row = intensity_row(g, nmax, la, list(range(-180, 180, 4))) num += math.cos(math.radians(la)) * sum(row) den += math.cos(math.radians(la)) * len(row) print(f" {model.name} {t:.0f}: along the equator {ns(sum(eq) / len(eq))}, worldwide {ns(num / den)}, " f"dipole {math.sqrt(g[0] ** 2 + g[1] ** 2 + g[2] ** 2):,.1f}") if __name__ == "__main__": main()