#!/usr/bin/env python3 """Earth's dipole and the South Atlantic Anomaly in IGRF-14, 1900-2030. Fact-Checked Science, "Earth's magnetic field" guide, https://factchecked.science/magnetic-field-guide/ ("Our analysis" in sections 2 and 3). Standard library only. Run it python3 field_decline.py [folder] The folder (default: the script's own) must hold, or will receive, the two inputs below. If a file is missing the script downloads it from NOAA, then checks its SHA-256 against the copy we used on October 26, 2026 (a different hash means NOAA changed the file, so numbers may differ). It prints a summary and writes field-decline.json into the folder. Inputs igrf14coeffs.txt IGRF-14 Gauss coefficients 1900-2025 and secular variation 2025-30 (IAGA V-MOD): https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt WMM2025COF.zip World Magnetic Model 2025 coefficients and NOAA's official test values; used to check this script's field evaluation and for WMM2025's own dipole: https://www.ncei.noaa.gov/sites/default/files/2025-01/WMM2025COF.zip What it computes 1. A check of our spherical-harmonic code against NOAA's 12 WMM2025 test values (X, Y, Z, F, rates). 2. The dipole at every IGRF-14 epoch 1900-2025 and the 2030 forecast: moment M = (4 pi a^3 / mu0) sqrt(g10^2 + g11^2 + h11^2), a = 6371.2 km; the axial dipole |g10|; the tilt of the dipole axis; the geomagnetic north pole (geocentric and WGS84 geodetic latitude). 3. Rates of change per century over spans and five-year steps, and whether the decline is speeding up (least-squares trend of the five-year paces; quadratic fits to M). 4. The South Atlantic Anomaly at the WGS84 surface (1-degree grid, minima refined to 0.05 degrees): the minimum total intensity and its place, the area below 32,000 / 30,000 / 28,000 / 26,000 / 25,000 / 24,000 nT (exact ellipsoidal cell areas), local minima below 32,000 nT, the global mean field; at every epoch, 2014.0, 2026.0 and 2030.0; and, year by year, when the second minimum south-west of Africa appears. 5. How deep that second low is: its depth below the saddle joining it to the main minimum (the level at which the two lows' basins first connect, on a 0.25-degree grid over 60-10 S, 80 W-20 E). Conventions: Schmidt semi-normalised Legendre functions; coefficients interpolated linearly between epochs and extrapolated after 2025 with the 2025-30 secular variation, as IGRF is defined; X north, Y east, Z down. IGRF-14 runs to degree 10 before 2000 and degree 13 from 2000. """ import hashlib import io import json import math import os import sys import urllib.request import zipfile IGRF_URL = "https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt" WMM_URL = "https://www.ncei.noaa.gov/sites/default/files/2025-01/WMM2025COF.zip" HASHES = {"igrf14coeffs.txt": "8f8d88403028fc4ee92c4f38d97b46e0a87e2cfc496045b43c9e26c1d6b0903c", "WMM2025COF.zip": "999080a12c95122c2031aeac8be77229d97584530e87e0f1a6df7d946ecf90c5"} A_REF = 6371.2 # km, IGRF/WMM reference radius 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 # T m / A THRESHOLDS = (32000, 30000, 28000, 26000, 25000, 24000) NMAX = 13 # ------------------------------------------------------------------ inputs def fetch(folder, name, url): path = os.path.join(folder, name) if not os.path.exists(path): print("downloading", url) req = urllib.request.Request(url, headers={"User-Agent": "factchecked.science field_decline.py"}) with urllib.request.urlopen(req, timeout=60) as r, open(path, "wb") as f: f.write(r.read()) with open(path, "rb") as f: digest = hashlib.sha256(f.read()).hexdigest() if digest != 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 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(epochs, main, sv, t): g, h = zeros(), zeros() if t >= epochs[-1]: i0 = len(epochs) - 1 for (gh, n, m), v in main.items(): (g if gh == "g" else h)[n][m] = v[i0] + (t - epochs[-1]) * sv[(gh, n, m)] else: i0 = max(i for i, e in enumerate(epochs) if e <= t) w = (t - epochs[i0]) / (epochs[i0 + 1] - epochs[i0]) for (gh, n, m), v in main.items(): (g if gh == "g" else h)[n][m] = v[i0] + w * (v[i0 + 1] - v[i0]) return g, h def read_wmm(zpath): """WMM.COF and the test values, read from the zip in memory (data only).""" z = zipfile.ZipFile(zpath) names = {os.path.basename(n): n for n in z.namelist()} cof = z.read(names["WMM.COF"]).decode("ascii") tests = z.read(names["WMM2025_TEST_VALUES.txt"]).decode("ascii") g, h, gd, hd = zeros(), zeros(), zeros(), zeros() lines = cof.splitlines() epoch = float(lines[0].split()[0]) for line in lines[1:]: p = line.split() if len(p) < 6 or p[0].startswith("9999"): continue n, m = int(p[0]), int(p[1]) g[n][m], h[n][m], gd[n][m], hd[n][m] = map(float, p[2:6]) return epoch, g, h, gd, hd, tests # ------------------------------------------------------------------ field synthesis 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: """Everything that depends only on geodetic latitude and height.""" 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(xyz): return math.sqrt(sum(v * v for v in xyz)) # ------------------------------------------------------------------ 1. the check against NOAA's test values def validate(epoch, g, h, gd, hd, tests): worst = {k: 0.0 for k in ("X", "Y", "Z", "F", "Xd", "Yd", "Zd", "Fd")} n = 0 for line in tests.splitlines(): if line.startswith("#") or not line.strip(): continue p = [float(x) for x in line.split()] t, hkm, lat, lon = p[0], p[1], p[2], p[3] gg = [[g[i][j] + (t - epoch) * gd[i][j] for j in range(NMAX + 1)] for i in range(NMAX + 1)] hh = [[h[i][j] + (t - epoch) * hd[i][j] for j in range(NMAX + 1)] for i in range(NMAX + 1)] X, Y, Z = field(gg, hh, lat, lon, hkm) F = math.sqrt(X * X + Y * Y + Z * Z) Xd, Yd, Zd = field(gd, hd, lat, lon, hkm) Fd = (X * Xd + Y * Yd + Z * Zd) / F for key, mine, ref in (("X", X, p[4]), ("Y", Y, p[5]), ("Z", Z, p[6]), ("F", F, p[8]), ("Xd", Xd, p[12]), ("Yd", Yd, p[13]), ("Zd", Zd, p[14]), ("Fd", Fd, p[16])): worst[key] = max(worst[key], abs(mine - ref)) n += 1 return {"test_points": n, "largest_abs_difference": {k: round(v, 3) for k, v in worst.items()}, "pass_within_0.1": all(v <= 0.1 for v in worst.values())} # ------------------------------------------------------------------ 2-3. the dipole def dipole(g10, g11, h11): b0 = math.sqrt(g10 * g10 + g11 * g11 + h11 * h11) M = (A_REF * 1e3) ** 3 * (b0 * 1e-9) / MU0_OVER_4PI / 1e22 tilt = math.degrees(math.acos(-g10 / b0)) plat = 90.0 - tilt return {"B0_nT": b0, "M_1e22_Am2": M, "absg10_nT": abs(g10), "tilt_deg": tilt, "pole_lat": plat, "pole_lat_geodetic": math.degrees(math.atan(math.tan(math.radians(plat)) / (1 - WGS84_E2))), "pole_lon": math.degrees(math.atan2(-h11, -g11)), "ADM_1e22_Am2": (A_REF * 1e3) ** 3 * (abs(g10) * 1e-9) / MU0_OVER_4PI / 1e22} def gc_km(lat1, lon1, lat2, lon2, R=6371.0): p1, p2 = math.radians(lat1), math.radians(lat2) c = math.sin(p1) * math.sin(p2) + math.cos(p1) * math.cos(p2) * math.cos(math.radians(lon2 - lon1)) return R * math.acos(max(-1.0, min(1.0, c))) def span(series, t0, t1, key): a, b = series[t0][key], series[t1][key] pc = (b - a) / (t1 - t0) * 100.0 return {"from": t0, "to": t1, "start": round(a, 4), "end": round(b, 4), "change_pct": round(100 * (b - a) / a, 2), "per_century_abs": round(pc, 4), "per_century_pct_of_start": round(100 * pc / a, 2)} def lstsq(X, y): n = len(X[0]) A = [[sum(r[i] * r[j] for r in X) for j in range(n)] for i in range(n)] b = [sum(r[i] * yy for r, yy in zip(X, y)) for i in range(n)] for i in range(n): p = max(range(i, n), key=lambda k: abs(A[k][i])) A[i], A[p], b[i], b[p] = A[p], A[i], b[p], b[i] for k in range(n): if k != i: f = A[k][i] / A[i][i] A[k] = [a - f * c for a, c in zip(A[k], A[i])] b[k] -= f * b[i] return [b[i] / A[i][i] for i in range(n)] # ------------------------------------------------------------------ 4. the anomaly def band_area(lat1, lat2): """Area of the WGS84 ellipsoid between two geodetic latitudes, km^2 per radian of longitude.""" 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) def refine_min(g, h, lat0, lon0, half=1.5, step=0.05): best = (1e9, lat0, lon0) 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): row = Row(lat0 + i * step) s = row.sums(g, h) for lo, (c, sn) in zip(lons, trigs): F = F_of(row.xyz(s, c, sn)) if F < best[0]: best = (F, lat0 + i * step, lo) return best class Grid: def __init__(self): self.lat = [-89.5 + i for i in range(180)] self.lon = [-179.5 + j for j in range(360)] self.rows = [Row(la) for la in self.lat] self.trigs = [trig(lo) for lo in self.lon] self.area = [band_area(la - 0.5, la + 0.5) * math.radians(1.0) for la in self.lat] def anomaly(g, h, G): Fg = [] for row in G.rows: s = row.sums(g, h) Fg.append([F_of(row.xyz(s, c, sn)) for (c, sn) in G.trigs]) nlat, nlon = len(G.lat), len(G.lon) below = {t: 0.0 for t in THRESHOLDS} tot = wsum = 0.0 gmin = (1e9, 0, 0) for i in range(nlat): for j in range(nlon): F = Fg[i][j] tot += F * G.area[i] wsum += G.area[i] for t in THRESHOLDS: if F < t: below[t] += G.area[i] if F < gmin[0]: gmin = (F, i, j) minima = [] for i in range(1, nlat - 1): for j in range(nlon): F = Fg[i][j] if F < 32000 and all(Fg[i + di][(j + dj) % nlon] > F for di in (-1, 0, 1) for dj in (-1, 0, 1) if di or dj): Fr, la, lo = refine_min(g, h, G.lat[i], G.lon[j], half=1.0) minima.append({"F_nT": round(Fr), "lat": round(la, 2), "lon": round(lo, 2)}) Fr, la, lo = refine_min(g, h, G.lat[gmin[1]], G.lon[gmin[2]]) connected = {} for t in THRESHOLDS: # share of each area joined to the main minimum seen = {(gmin[1], gmin[2])} stack = [(gmin[1], gmin[2])] area = 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 Fg[ii][jj] < t: seen.add((ii, jj)) stack.append((ii, jj)) connected[str(t)] = round(area / below[t], 4) if below[t] else None return {"min_F_nT": round(Fr, 1), "min_lat": round(la, 2), "min_lon": round(lo, 2), "share_connected_to_minimum": connected, "local_minima_below_32000": sorted(minima, key=lambda d: d["F_nT"]), "area_below": {str(t): {"million_km2": round(below[t] / 1e6, 2), "pct_of_surface": round(100 * below[t] / EARTH_AREA, 3)} for t in THRESHOLDS}, "global_mean_F_nT": round(tot / wsum, 1)} def second_minimum(g, h): """A strict local minimum south-west of Africa (60-20 S, 25 W-20 E), on a 0.25-degree grid.""" lats = [-60.0 + 0.25 * i for i in range(161)] lons = [-25.0 + 0.25 * j for j in range(181)] tr = [trig(lo) for lo in lons] grid = [] for la in lats: row = Row(la) s = row.sums(g, h) grid.append([F_of(row.xyz(s, c, sn)) for (c, sn) in tr]) found = None for i in range(1, len(lats) - 1): for j in range(1, len(lons) - 1): F = grid[i][j] if all(grid[i + di][j + dj] > F for di in (-1, 0, 1) for dj in (-1, 0, 1) if di or dj): if found is None or F < found[0]: found = (F, lats[i], lons[j]) return None if found is None else {"F_nT": round(found[0]), "lat": found[1], "lon": found[2]} def second_low_depth(g, h, sec): """Depth of the second low below the saddle that joins it to the main minimum (0.25-degree grid).""" lats = [-60.0 + 0.25 * i for i in range(201)] lons = [-80.0 + 0.25 * j for j in range(401)] tr = [trig(lo) for lo in lons] F = [] for la in lats: row = Row(la) s = row.sums(g, h) F.append([F_of(row.xyz(s, c, sn)) for (c, sn) in tr]) nlat, nlon = len(lats), len(lons) main_ij = min(((i, j) for i in range(nlat) for j in range(nlon)), key=lambda p: F[p[0]][p[1]]) sec_ij = (int(round((sec["lat"] + 60.0) / 0.25)), int(round((sec["lon"] + 80.0) / 0.25))) parent = {} def find(p): while parent[p] != p: parent[p] = parent[parent[p]] p = parent[p] return p for v, i, j in sorted((F[i][j], i, j) for i in range(nlat) for j in range(nlon)): # fill from the bottom up parent[(i, j)] = (i, j) for di in (-1, 0, 1): for dj in (-1, 0, 1): q = (i + di, j + dj) if q in parent: a, b = find((i, j)), find(q) if a != b: parent[a] = b if main_ij in parent and sec_ij in parent and find(main_ij) == find(sec_ij): return {"second_F_nT": round(F[sec_ij[0]][sec_ij[1]], 1), "saddle_F_nT": round(v, 1), "depth_nT": round(v - F[sec_ij[0]][sec_ij[1]], 1)} return None # ------------------------------------------------------------------ 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__)) igrf, h1 = fetch(folder, "igrf14coeffs.txt", IGRF_URL) wmmz, h2 = fetch(folder, "WMM2025COF.zip", WMM_URL) out = {"what": "Dipole and South Atlantic Anomaly from IGRF-14 (factchecked.science/magnetic-field-guide/, Our analysis)", "inputs_sha256": {"igrf14coeffs.txt": h1, "WMM2025COF.zip": h2}} ew, gw, hw, gdw, hdw, tests = read_wmm(wmmz) 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"]) w25 = dipole(gw[1][0], gw[1][1], hw[1][1]) w30 = dipole(gw[1][0] + 5 * gdw[1][0], gw[1][1] + 5 * gdw[1][1], hw[1][1] + 5 * hdw[1][1]) out["wmm2025_dipole"] = {"M_2025": round(w25["M_1e22_Am2"], 4), "M_2030": round(w30["M_1e22_Am2"], 4), "pct_per_century": round(100 * (w30["M_1e22_Am2"] - w25["M_1e22_Am2"]) / 5 * 100 / w25["M_1e22_Am2"], 2), "tilt_2025": round(w25["tilt_deg"], 3), "pole_2025_geocentric": [round(w25["pole_lat"], 3), round(w25["pole_lon"], 3)]} epochs, main_c, sv = read_igrf(igrf) times = epochs + [2030.0] ser = {} for t in times: g, h = igrf_at(epochs, main_c, sv, t) ser[t] = dipole(g[1][0], g[1][1], h[1][1]) out["dipole_series"] = [{"epoch": t, "forecast": t > 2025, **{k: round(v, 4) for k, v in d.items()}} for t, d in ser.items()] spans = [(1900.0, 1950.0), (1950.0, 2000.0), (2000.0, 2025.0), (1900.0, 2025.0), (2025.0, 2030.0), (1900.0, 1925.0), (1925.0, 1950.0), (1950.0, 1975.0), (1975.0, 2000.0)] out["rates_M"] = [span(ser, a, b, "M_1e22_Am2") for a, b in spans] out["rates_absg10"] = [span(ser, a, b, "absg10_nT") for a, b in spans] ts = sorted(ser) five = [] for a, b in zip(ts[:-1], ts[1:]): five.append({"from": a, "to": b, "forecast": b > 2025, "M_pct_per_century": -span(ser, a, b, "M_1e22_Am2")["per_century_pct_of_start"], "g10_nT_per_yr": round((ser[a]["absg10_nT"] - ser[b]["absg10_nT"]) / (b - a), 2), "pole_km_per_yr": round(gc_km(ser[a]["pole_lat"], ser[a]["pole_lon"], ser[b]["pole_lat"], ser[b]["pole_lon"]) / (b - a), 1)}) out["five_year_paces"] = five # positive = weakening, % of the start value per century # is it speeding up? meas = [s for s in five if not s["forecast"]] tm = [(s["from"] + s["to"]) / 2 - 1962.5 for s in meas] b0, b1 = lstsq([[1, t] for t in tm], [s["M_pct_per_century"] for s in meas]) acc = {"mean_five_year_pace": round(sum(s["M_pct_per_century"] for s in meas) / len(meas), 2), "trend_of_paces_pct_per_century_per_century": round(100 * b1, 2), "quadratic_fits": []} pts_all = [(t, ser[t]["M_1e22_Am2"]) for t in ts if t <= 2025] for t0 in (1900.0, 1950.0, 1975.0): pts = [p for p in pts_all if p[0] >= t0] tc = (pts[0][0] + pts[-1][0]) / 2 a, b, c = lstsq([[1, t - tc, (t - tc) ** 2] for t, _ in pts], [m for _, m in pts]) rate = lambda t: -(b + 2 * c * (t - tc)) * 100 # 1e22 A m^2 per century, positive = weakening acc["quadratic_fits"].append({"from": t0, "to": 2025.0, "curvature": c, "pace_at_start_abs": round(rate(pts[0][0]), 4), "pace_at_2025_abs": round(rate(2025.0), 4), "pace_at_start_pct": round(100 * rate(pts[0][0]) / pts[0][1], 2), "pace_at_2025_pct": round(100 * rate(2025.0) / pts[-1][1], 2)}) out["acceleration"] = acc path = sum(gc_km(ser[a]["pole_lat"], ser[a]["pole_lon"], ser[b]["pole_lat"], ser[b]["pole_lon"]) for a, b in zip(ts[:-2], ts[1:-1])) out["geomagnetic_pole_1900_2025"] = {"straight_km": round(gc_km(ser[1900.0]["pole_lat"], ser[1900.0]["pole_lon"], ser[2025.0]["pole_lat"], ser[2025.0]["pole_lon"])), "path_km": round(path)} G = Grid() saa = [] for t in sorted(set(times + [2014.0, 2026.0])): g, h = igrf_at(epochs, main_c, sv, t) a = anomaly(g, h, G) a["epoch"] = t saa.append(a) print("%7.1f minimum %8.1f nT at %6.2f %7.2f below 26,000 nT %5.2f%% below 32,000 nT %5.2f%%" % ( t, a["min_F_nT"], a["min_lat"], a["min_lon"], a["area_below"]["26000"]["pct_of_surface"], a["area_below"]["32000"]["pct_of_surface"])) out["anomaly"] = saa out["earth_area_km2"] = round(EARTH_AREA) by = {a["epoch"]: a for a in saa} for t0, t1 in ((1900.0, 2025.0), (1970.0, 2020.0), (2014.0, 2025.0), (2020.0, 2025.0), (2025.0, 2026.0)): a, b = by[t0], by[t1] d = gc_km(a["min_lat"], a["min_lon"], b["min_lat"], b["min_lon"]) out["minimum_change_%d_%d" % (t0, t1)] = {"F_from": a["min_F_nT"], "F_to": b["min_F_nT"], "nT": round(b["min_F_nT"] - a["min_F_nT"], 1), "pct": round(100 * (b["min_F_nT"] - a["min_F_nT"]) / a["min_F_nT"], 2), "moved_km": round(d), "km_per_yr": round(d / (t1 - t0), 1)} second = [] for yr in range(1990, 2031): g, h = igrf_at(epochs, main_c, sv, float(yr)) second.append({"year": yr, "second_minimum": second_minimum(g, h)}) out["second_minimum_yearly"] = second out["second_minimum_first_year"] = next((s["year"] for s in second if s["second_minimum"]), None) depth = [] for s in second: if s["second_minimum"] and s["year"] in (2006, 2007, 2008, 2010, 2015, 2020, 2025, 2030): g, h = igrf_at(epochs, main_c, sv, float(s["year"])) depth.append({"year": s["year"], **second_low_depth(g, h, s["second_minimum"])}) out["second_low_depth"] = depth with open(os.path.join(folder, "field-decline.json"), "w", encoding="utf-8", newline="\n") as f: json.dump(out, f, indent=1) print("\nDipole moment (10^22 A m^2), tilt, geomagnetic north pole (geodetic):") for t in ts: d = ser[t] print(" %6.1f %.4f tilt %.2f %.2f N %.2f E%s" % (t, d["M_1e22_Am2"], d["tilt_deg"], d["pole_lat_geodetic"], d["pole_lon"], " (forecast)" if t > 2025 else "")) print("Per century:", ", ".join("%d-%d %+.2f%%" % (r["from"], r["to"], r["per_century_pct_of_start"]) for r in out["rates_M"][:5])) print("Five-year paces, %% a century: %s" % ", ".join("%d-%d %.1f" % (s["from"], s["to"], s["M_pct_per_century"]) for s in five)) print("Speeding up? mean pace %.2f%%; trend %+.2f%% a century per century; quadratic fits:" % (acc["mean_five_year_pace"], acc["trend_of_paces_pct_per_century_per_century"]), "; ".join("%d-2025 pace %.2f -> %.2f e22/century" % (q["from"], q["pace_at_start_abs"], q["pace_at_2025_abs"]) for q in acc["quadratic_fits"])) print("Second minimum first appears:", out["second_minimum_first_year"]) print("Its depth below the saddle to the main minimum:", ", ".join("%d %.1f nT" % (d["year"], d["depth_nT"]) for d in depth)) print("Wrote", os.path.join(folder, "field-decline.json")) if __name__ == "__main__": main()