#!/usr/bin/env python3 """Today's dipole against the long record: PADM2M (axial dipole moment, 0-2 million years; Ziegler, Constable, Johnson and Tauxe 2011) and GGF100k (dipole and axial dipole moment, 0-100,000 years; Panovska, Constable and Korte 2018). Fact-Checked Science, "Earth's magnetic field" guide, section 1 ("Our analysis"). Standard library only. Run it python3 long_record.py [folder] The folder (default: the script's own) must hold, or will receive: padm2m.model.tar https://earthref.org/ERDA/download:1138/ (record: https://earthref.org/ERDA/1138/) ggf100k.model.zip https://earthref.org/ERDA/download:2382/ (record: https://earthref.org/ERDA/2382/) igrf14coeffs.txt https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt (today's dipole) Missing files are downloaded and their SHA-256 compared with the copies we used on October 26, 2026. Only the data files are read from the two archives, in memory (padm2m.adms.1k; DM_GGF100k.dat and ADM_GGF100k.dat); nothing in them is extracted to disk or run. It prints the results and writes long-record.json, with the axial dipole series charted in section 1 (PADM2M every 2,000 years, GGF100k every 250). The model files give no uncertainties. Panovska et al. (2018, section 4.2) list the 100,000-year averages of four reconstructions, 6.04 to 7.13 x 10^22 A m^2, which is a fair measure of how far they disagree. """ import hashlib import io import json import os import sys import tarfile import urllib.request import zipfile FILES = { "padm2m.model.tar": ("https://earthref.org/ERDA/download:1138/", "7085bcceb1b20159908ecc7f0fa29125e52a9c8c9daab2d65d20320b2b88d303"), "ggf100k.model.zip": ("https://earthref.org/ERDA/download:2382/", "b33fe63b20abf01ced3b0106e279abfb9223add86c5c5e020a383e84dd3f1293"), "igrf14coeffs.txt": ("https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt", "8f8d88403028fc4ee92c4f38d97b46e0a87e2cfc496045b43c9e26c1d6b0903c"), } A_REF = 6371.2e3 def fetch(folder, name): url, want = FILES[name] 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 long_record.py"}) with urllib.request.urlopen(req, timeout=120) as r, open(path, "wb") as f: f.write(r.read()) with open(path, "rb") as f: digest = hashlib.sha256(f.read()).hexdigest() if want and digest != want: print("Note: %s differs from the copy we used; results may differ." % name) return path, digest def two_columns(text): rows = [] for line in text.splitlines(): p = line.split() try: rows.append((float(p[0]), float(p[1]))) except (ValueError, IndexError): continue return sorted(rows) def stats(series, t0, t1): v = sorted(m for t, m in series if t0 <= t <= t1) return {"from_ka": t0, "to_ka": t1, "n": len(v), "mean": round(sum(v) / len(v), 3), "median": round(v[len(v) // 2], 3), "min": round(v[0], 3), "max": round(v[-1], 3)} def share_at_or_above(series, t0, t1, level): v = [m for t, m in series if t0 <= t <= t1] return round(100.0 * sum(1 for m in v if m >= level) / len(v), 1) def lows(series, below): out, run = [], [] for t, m in series: if m < below: run.append((m, t)) elif run: out.append(min(run)) run = [] if run: out.append(min(run)) return [{"ka": round(t, 2), "moment": round(m, 3)} for m, t in out] def today(igrf_path): with open(igrf_path, encoding="ascii") as f: coef = {} for line in f: if line[:1] in "gh" and not line.startswith("g/h"): p = line.split() coef[(p[0], int(p[1]), int(p[2]))] = float(p[-2]) # the 2025.0 column (the last is the SV) g10, g11, h11 = coef[("g", 1, 0)], coef[("g", 1, 1)], coef[("h", 1, 1)] k = A_REF ** 3 * 1e-9 / 1e-7 / 1e22 return {"M_2025": round(k * (g10 ** 2 + g11 ** 2 + h11 ** 2) ** 0.5, 3), "ADM_2025": round(k * abs(g10), 3)} 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 = {n: fetch(folder, n) for n in FILES} with tarfile.open(paths["padm2m.model.tar"][0]) as t: member = next(m for m in t.getmembers() if m.name.endswith("padm2m.adms.1k")) padm = [(age_myr * 1000.0, adm) for age_myr, adm in two_columns(t.extractfile(member).read().decode("ascii", "replace"))] with zipfile.ZipFile(paths["ggf100k.model.zip"][0]) as z: names = {os.path.basename(n): n for n in z.namelist()} dm = two_columns(z.read(names["DM_GGF100k.dat"]).decode("ascii", "replace")) adm = two_columns(z.read(names["ADM_GGF100k.dat"]).decode("ascii", "replace")) now = today(paths["igrf14coeffs.txt"][0]) out = {"what": "Today's dipole against PADM2M and GGF100k (factchecked.science/magnetic-field-guide/, Our analysis)", "inputs_sha256": {n: p[1] for n, p in paths.items()}, "today_igrf14": now, "padm2m": {"present_value": round(padm[0][1], 3), "last_100ka": stats(padm, 0, 100), "last_780ka": stats(padm, 0, 780), "last_2Ma": stats(padm, 0, 2000), "from_780_to_2000ka": stats(padm, 780, 2000), "pct_of_last_2Ma_at_or_above_today": share_at_or_above(padm, 0, 2000, now["ADM_2025"]), "pct_of_last_780ka_at_or_above_today": share_at_or_above(padm, 0, 780, now["ADM_2025"]), "lows_below_3": lows(padm, 3.0)}, "ggf100k_DM": {"value_1950": round(dm[0][1], 3), "last_100ka": stats(dm, 0, 100), "last_12ka": stats(dm, 0, 12), "pct_of_last_100ka_at_or_above_today": share_at_or_above(dm, 0, 100, now["M_2025"]), "lows_below_4": lows(dm, 4.0)}, "ggf100k_ADM": {"value_1950": round(adm[0][1], 3), "last_100ka": stats(adm, 0, 100), "pct_of_last_100ka_at_or_above_today": share_at_or_above(adm, 0, 100, now["ADM_2025"]), "lows_below_4": lows(adm, 4.0)}} out["today_over_mean"] = {"ADM_vs_PADM2M_2Ma": round(now["ADM_2025"] / out["padm2m"]["last_2Ma"]["mean"], 3), "ADM_vs_PADM2M_780ka": round(now["ADM_2025"] / out["padm2m"]["last_780ka"]["mean"], 3), "ADM_vs_PADM2M_100ka": round(now["ADM_2025"] / out["padm2m"]["last_100ka"]["mean"], 3), "ADM_vs_GGF100k_100ka": round(now["ADM_2025"] / out["ggf100k_ADM"]["last_100ka"]["mean"], 3), "M_vs_GGF100k_DM_100ka": round(now["M_2025"] / out["ggf100k_DM"]["last_100ka"]["mean"], 3)} out["chart_padm2m_every_2ka"] = [[round(t), round(m, 3)] for t, m in padm if round(t) % 2 == 0] out["chart_ggf100k_ADM_every_250yr"] = [[round(t, 2), round(m, 3)] for t, m in adm if abs(t * 4 - round(t * 4)) < 1e-6] with open(os.path.join(folder, "long-record.json"), "w", encoding="utf-8", newline="\n") as f: json.dump(out, f, indent=1) for k in ("today_igrf14", "today_over_mean"): print(k, out[k]) for k in ("padm2m", "ggf100k_DM", "ggf100k_ADM"): print(k, {kk: vv for kk, vv in out[k].items() if not kk.startswith("lows")}) print("Wrote", os.path.join(folder, "long-record.json")) if __name__ == "__main__": main()