#!/usr/bin/env python3 """Claim 6, "More vulnerable to solar storms": how much the field's weakening since 1900 changes Earth's exposure, next to what storms and the solar cycle do. Fact-Checked Science, "Earth's magnetic field" guide, section 7 ("Our analysis"). Standard library only. Run it python3 claim6_scalings.py [folder] [--grid] The folder (default: the script's own) must hold, or will receive: igrf14coeffs.txt https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt four CSV files from the Oulu neutron monitor's online query (URLs in OULU below); data of the Oulu Cosmic Ray Station, https://cosmicrays.oulu.fi/, Sodankyla Geophysical Observatory, for research, education and non-commercial services Missing files are downloaded and their SHA-256 compared with the copies we used on October 26, 2026 (the Oulu query can return slightly revised data). It prints everything and writes claim6-scalings.json. --grid adds the global grid of cutoffs, 1900 and 2025 (2.5 x 5 degrees; slow, tens of minutes in pure Python). What it computes A. IGRF-14 dipole moment, pole and tilt, 1900-2030. B. For 35 places: centred-dipole (CD) and apex latitude, and the Stormer vertical cutoff rigidity at 20 km, Rc = mu0 c M cos^4(lat) / (16 pi r^2), with the CD latitude and with the apex latitude (the latitude of the IGRF field line's apex mapped to the ground with cos^2(lat) = a / r_apex, a proxy that includes the non-dipole field); and NERC TPL-007-4's latitude factor alpha = 0.001 exp(0.115 L). C. The proxy against the trajectory-traced cutoffs of Herbst, Kopp & Heber (2013, Table 1). D. The auroral oval's edge under the three scalings in Glassmeier et al. (2004, p. 3669 and Table 1), cos(L) ~ M^(alpha - 1/2). E. The magnetopause stand-off distance, R ~ M^(1/3) P^(-1/6). F. Global radiocarbon production, re-implementing Kovaltsov, Mishev & Usoskin (2012), eqs. 4-7 and Table A.1, with eq. (5)'s geomagnetic transmission: the share of the globe a particle of rigidity P reaches is 1 - f(E) = 1 - sqrt(1 - sqrt(P / 1.9M)) for P <= 1.9M, and 1 above (f is sin of the latitude where the dipole cutoff 1.9M cos^4 equals P). Checked against the paper: within about +-15% for M = 6 to 10, and a factor of two between solar maximum and minimum over 1951-2010 (p. 5). Corrected on October 26, 2026, after review: the first version used P / 1.9M for the transmission, which overstated the field's effect. G. The Oulu neutron monitor's solar-cycle swings and two Forbush decreases (October 2003, May 2024). H, I. What a change of latitude does to TPL-007's benchmark field, and each place's net change in distance from a storm-time oval edge near 50 degrees, both effects combined. J. 2015-2020 rates at four places, beside Maffei, Livermore & Mound (2025). K. (--grid) The area-weighted mean cutoff, and the share of Earth's area below 1 GV. Geocentric latitudes are used for the field (a = 6371.2 km); the geodetic correction changes latitudes by at most about 0.2 degrees and does not affect differences between epochs. """ import csv import hashlib import json import math import os import sys import urllib.request A_KM = 6371.2 MU0 = 4e-7 * math.pi OULU_Q = "https://cosmicrays.oulu.fi/webform/onlinequery.cgi?station=OULU&" FILES = { "igrf14coeffs.txt": ("https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt", "8f8d88403028fc4ee92c4f38d97b46e0a87e2cfc496045b43c9e26c1d6b0903c"), "oulu-nm-yearly-1965-2025.csv": (OULU_Q + "startday=01&startmonth=01&startyear=1965&starttime=00%3A00&endday=01&endmonth=01" "&endyear=2026&endtime=00%3A00&resolution=525680&outputmode=csv", "9961245d8eb5d074794d70d5ff8091ab14a1c1065780e397f74d07023a9eb0ba"), "oulu-nm-monthly-1965-2026.csv": (OULU_Q + "startday=01&startmonth=01&startyear=1965&starttime=00%3A00&endday=01&endmonth=10" "&endyear=2026&endtime=00%3A00&resolution=44640&outputmode=csv", "6bf012e5aaa9058621a9cd9202aada3eadbfdae2b7531c19bb38a5972c20c5d3"), "oulu-nm-hourly-2003-10-15_2003-11-30.csv": (OULU_Q + "startday=15&startmonth=10&startyear=2003&starttime=00%3A00&endday=30" "&endmonth=11&endyear=2003&endtime=00%3A00&resolution=60&outputmode=csv", "959fd64c9ddf75aebaa4b10b69dbcc7f23246700332a014b404e12f658ccfd08"), "oulu-nm-hourly-2024-05-01_2024-05-31.csv": (OULU_Q + "startday=01&startmonth=05&startyear=2024&starttime=00%3A00&endday=31" "&endmonth=05&endyear=2024&endtime=00%3A00&resolution=60&outputmode=csv", "5b2775ca9611be6ef850ef473b645baddef5b1124fc1587898430766e41e95c4"), } OUT = {} LOG = [] def say(*a): s = " ".join(str(x) for x in a) LOG.append(s) print(s) 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 claim6_scalings.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 digest != want: say(f"Note: {name} differs from the copy we used; results may differ slightly.") return path # ---------------------------------------------------------------- IGRF-14 tools def read_igrf(path): """{epoch: {(g|h, n, m): nT}} for 1900-2025, plus 2030.0 = 2025 + 5 years of the secular variation.""" header, rows = None, [] with open(path, encoding="utf-8", errors="replace") as f: for line in f: p = line.split() if not p or line.startswith("#"): continue if p[0] == "g/h": header = p elif p[0] in ("g", "h") and header is not None: rows.append(p) cols = header[3:] epochs = {float(c): {} for c in cols if "-" not in c} sv = {} for p in rows: key = (p[0], int(p[1]), int(p[2])) for c, v in zip(cols, p[3:]): if "-" in c: sv[key] = float(v) else: epochs[float(c)][key] = float(v) last = max(epochs) epochs[last + 5.0] = {k: v + 5.0 * sv.get(k, 0.0) for k, v in epochs[last].items()} return epochs def dipole(c): g10, g11, h11 = c[("g", 1, 0)], c[("g", 1, 1)], c[("h", 1, 1)] b0 = math.sqrt(g10 ** 2 + g11 ** 2 + h11 ** 2) colat = math.degrees(math.acos(-g10 / b0)) return {"B0_nT": b0, "M_Am2": 4 * math.pi * (A_KM * 1e3) ** 3 / MU0 * b0 * 1e-9, "pole_lat": 90.0 - colat, "pole_lon": math.degrees(math.atan2(-h11, -g11)), "tilt_deg": colat} def cd_latitude(lat, lon, plat, plon): la, lo, pa, po = map(math.radians, (lat, lon, plat, plon)) s = math.sin(la) * math.sin(pa) + math.cos(la) * math.cos(pa) * math.cos(lo - po) return math.degrees(math.asin(max(-1.0, min(1.0, s)))) def legendre(nmax, theta): """Schmidt semi-normalised P_n^m(cos theta) and dP/dtheta.""" ct, st = math.cos(theta), math.sin(theta) P = [[0.0] * (nmax + 2) for _ in range(nmax + 2)] dP = [[0.0] * (nmax + 2) for _ in range(nmax + 2)] P[0][0] = 1.0 for n in range(1, nmax + 1): if n == 1: P[1][1], dP[1][1] = st, ct else: k = math.sqrt((2 * n - 1) / (2 * n)) P[n][n] = k * st * P[n - 1][n - 1] dP[n][n] = k * (ct * P[n - 1][n - 1] + st * dP[n - 1][n - 1]) for m in range(0, n): a, b, d = 2 * n - 1, math.sqrt(max(0.0, (n - 1) ** 2 - m ** 2)), math.sqrt(n ** 2 - m ** 2) p2 = P[n - 2][m] if n >= 2 else 0.0 dp2 = dP[n - 2][m] if n >= 2 else 0.0 P[n][m] = (a * ct * P[n - 1][m] - b * p2) / d dP[n][m] = (a * (ct * dP[n - 1][m] - st * P[n - 1][m]) - b * dp2) / d return P, dP def nmax_of(c): return max(n for (_, n, _) in c) def field(c, r_km, colat_deg, lon_deg, nmax): theta, phi = math.radians(colat_deg), math.radians(lon_deg) P, dP = legendre(nmax, theta) st = math.sin(theta) or 1e-10 br = bt = bp = 0.0 ar = A_KM / r_km for n in range(1, nmax + 1): f = ar ** (n + 2) for m in range(0, n + 1): g = c.get(("g", n, m), 0.0) h = c.get(("h", n, m), 0.0) if m > 0 else 0.0 cm, sm = math.cos(m * phi), math.sin(m * phi) gh = g * cm + h * sm br += (n + 1) * f * gh * P[n][m] bt += -f * gh * dP[n][m] bp += f * m * (g * sm - h * cm) * P[n][m] / st return br, bt, bp def bdir(c, x, y, z, nmax): r = math.sqrt(x * x + y * y + z * z) th = math.acos(max(-1.0, min(1.0, z / r))) ph = math.atan2(y, x) br, bt, bp = field(c, r, math.degrees(th), math.degrees(ph), nmax) st, ct, sp, cp = math.sin(th), math.cos(th), math.sin(ph), math.cos(ph) bx = br * st * cp + bt * ct * cp - bp * sp by = br * st * sp + bt * ct * sp + bp * cp bz = br * ct - bt * st b = math.sqrt(bx * bx + by * by + bz * bz) return bx / b, by / b, bz / b, br def apex_latitude(c, lat, lon, step_km=20.0, max_steps=200000): """Trace the field line up from the ground (RK4) to its apex; return (apex latitude, apex radius in km).""" nmax = nmax_of(c) th, ph = math.radians(90.0 - lat), math.radians(lon) x, y, z = A_KM * math.sin(th) * math.cos(ph), A_KM * math.sin(th) * math.sin(ph), A_KM * math.cos(th) sgn = 1.0 if bdir(c, x, y, z, nmax)[3] > 0 else -1.0 rmax = prev_r = A_KM h = step_km for _ in range(max_steps): k1 = bdir(c, x, y, z, nmax)[:3] k2 = bdir(c, x + 0.5 * h * sgn * k1[0], y + 0.5 * h * sgn * k1[1], z + 0.5 * h * sgn * k1[2], nmax)[:3] k3 = bdir(c, x + 0.5 * h * sgn * k2[0], y + 0.5 * h * sgn * k2[1], z + 0.5 * h * sgn * k2[2], nmax)[:3] k4 = bdir(c, x + h * sgn * k3[0], y + h * sgn * k3[1], z + h * sgn * k3[2], nmax)[:3] dx = [h * sgn * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6.0 for i in range(3)] x, y, z = x + dx[0], y + dx[1], z + dx[2] r = math.sqrt(x * x + y * y + z * z) rmax = max(rmax, r) if r < prev_r: if h > 0.5: # passed the apex: step back and refine x, y, z = x - dx[0], y - dx[1], z - dx[2] h /= 4.0 prev_r = math.sqrt(x * x + y * y + z * z) continue break prev_r = r h = max(step_km, 0.01 * r) if h >= step_km else h if r > 60 * A_KM: break lam = math.degrees(math.acos(math.sqrt(min(1.0, A_KM / rmax)))) return (lam if lat >= 0 else -lam), rmax def geodetic_to_geocentric(lat_deg): a, f = 6378.137, 1 / 298.257223563 la = math.radians(lat_deg) n = a / math.sqrt(1 - (2 * f - f * f) * math.sin(la) ** 2) return math.degrees(math.atan2(n * (1 - f) ** 2 * math.sin(la), n * math.cos(la))) def stormer_GV(moment, maglat, alt_km=20.0): r = (A_KM + alt_km) * 1e3 return MU0 * 299792458.0 * moment / (16 * math.pi * r * r) * math.cos(math.radians(maglat)) ** 4 / 1e9 def alpha_tpl007(L): """NERC TPL-007-4, Attachment 1, eq. (3): alpha = 0.001 exp(0.115 L), between 0.1 and 1.""" return min(1.0, max(0.1, 0.001 * math.exp(0.115 * abs(L)))) PLACES = [ # name, geodetic latitude, longitude (east positive) ("Anchorage", 61.22, -149.90), ("Edmonton", 53.55, -113.49), ("Winnipeg", 49.90, -97.14), ("Minneapolis", 44.98, -93.27), ("Chicago", 41.88, -87.63), ("Ottawa", 45.40, -75.55), ("Montreal", 45.50, -73.57), ("Quebec City", 46.81, -71.21), ("New York", 40.71, -74.01), ("Salem NJ", 39.46, -75.54), ("Washington DC", 38.91, -77.04), ("Atlanta", 33.75, -84.39), ("Seattle", 47.61, -122.33), ("Denver", 39.74, -104.99), ("San Antonio", 29.42, -98.49), ("Mexico City", 19.43, -99.13), ("British Virgin Isl.", 18.00, -64.00), ("Reykjavik", 64.15, -21.94), ("Oslo", 59.91, 10.75), ("Stockholm", 59.33, 18.07), ("Helsinki", 60.17, 24.94), ("Malmo", 55.60, 13.00), ("Edinburgh", 55.95, -3.19), ("London", 51.51, -0.13), ("Kiel", 54.34, 10.12), ("Berlin", 52.52, 13.40), ("Paris", 48.86, 2.35), ("Rome", 41.90, 12.50), ("Madrid", 40.42, -3.70), ("Moscow", 55.76, 37.62), ("Tokyo", 35.68, 139.69), ("Johannesburg", -26.20, 28.05), ("Buenos Aires", -34.60, -58.38), ("Sao Paulo", -23.55, -46.63), ("Sydney", -33.87, 151.21), ] HERBST = { # trajectory-traced vertical cutoffs, GV, for 1913, 1955, 1976 and 1997 (Herbst, Kopp & Heber 2013, Table 1) "Kiel": (54.34, 10.12, (2.64, 2.49, 2.53, 2.52)), "Moscow": (55.76, 37.62, (2.70, 2.38, 2.37, 2.29)), "Mexico City": (19.43, -99.13, (11.38, 9.79, 9.08, 8.21)), "British Virgin Isl.": (18.00, -64.00, (7.48, 7.84, 8.18, 9.74)), } R_FINLAY = 0.91 # the axial dipole about 9% weaker than in 1840 (Finlay, Aubert & Gillet 2016, p. 1) # Kovaltsov, Mishev & Usoskin (2012): Table A.1 yield functions (atoms sr per nucleon, divided by pi), eqs. 4-7 E_TAB = [0.1, 0.3, 0.5, 0.7, 1, 3, 7, 10, 19, 49, 99, 499, 999] YP = [0.025, 0.26, 0.72, 1.29, 2.07, 5.19, 8.32, 9.72, 12.40, 17.45, 23.24, 48.30, 72.73] YA4 = [0.036, 0.38, 0.89, 1.55, 2.16, 4.18, 7.17, 8.67, 12.40, 17.45, 23.24, 48.30, 72.73] ER = 0.938 def loginterp(x, xs, ys): i = 0 if x <= xs[0] else len(xs) - 2 if x >= xs[-1] else max(j for j in range(len(xs) - 1) if xs[j] <= x) x0, x1, y0, y1 = math.log(xs[i]), math.log(xs[i + 1]), math.log(ys[i]), math.log(ys[i + 1]) return math.exp(y0 + (y1 - y0) * (math.log(x) - x0) / (x1 - x0)) def j_lis(e): p = math.sqrt(e * (e + 2 * ER)) return 1.9e4 * p ** -2.78 / (1 + 0.4866 * p ** -2.51) def q_global(phi_MV, M22, emin=0.01, emax=1000.0, n=4000): """Global 14C production, atoms cm^-2 s^-1, for modulation potential phi (MV) and dipole moment M (1e22 A m^2).""" tot = 0.0 lo, hi = math.log(emin), math.log(emax) for k in range(n): e = math.exp(lo + (k + 0.5) * (hi - lo) / n) de = e * (hi - lo) / n for AZ, Y, scale in ((1.0, YP, 1.0), (2.0, YA4, 0.3)): Phi = phi_MV / 1000.0 / AZ J = scale * j_lis(e + Phi) * e * (e + 2 * ER) / ((e + Phi) * (e + Phi + 2 * ER)) P = AZ * math.sqrt(e * (e + 2 * ER)) Pc = 1.9 * M22 # eq. (5): f = sqrt(1 - sqrt(P / Pc)) for P <= Pc tot += math.pi * loginterp(e, E_TAB, Y) * J * (1.0 if P > Pc else 1.0 - math.sqrt(1.0 - math.sqrt(P / Pc))) * de return tot / 1e4 def read_oulu(path): out = [] with open(path, encoding="utf-8", errors="replace") as f: rd = csv.reader(f) next(rd, None) for r in rd: try: out.append((r[0], float(r[3]))) # timestamp, corrected count rate (counts per minute) except (ValueError, IndexError): pass return out def main(): sys.stdout.reconfigure(encoding="utf-8") args = [a for a in sys.argv[1:] if not a.startswith("--")] folder = args[0] if args else os.path.dirname(os.path.abspath(__file__)) paths = {n: fetch(folder, n) for n in FILES} E = read_igrf(paths["igrf14coeffs.txt"]) epochs = sorted(E) say("A. IGRF-14 dipole: moment, north geomagnetic pole, tilt") dip = {ep: dipole(E[ep]) for ep in epochs} for ep in epochs: d = dip[ep] say(f" {ep:6.1f} M = {d['M_Am2'] / 1e22:.4f} e22 A m2 pole {d['pole_lat']:.2f} N {d['pole_lon']:.2f} E tilt {d['tilt_deg']:.2f} deg" + (" (2030: forecast)" if ep == 2030.0 else "")) M1900, M2025 = dip[1900.0]["M_Am2"], dip[2025.0]["M_Am2"] ratio = M2025 / M1900 say(f" M(2025)/M(1900) = {ratio:.4f} ({(ratio - 1) * 100:.2f}%, {(ratio - 1) * 100 / 1.25:.2f}% a century)") OUT["dipole"] = {"M_1900_1e22": round(M1900 / 1e22, 4), "M_2025_1e22": round(M2025 / 1e22, 4), "ratio_2025_1900": round(ratio, 4), "tilt_1900": round(dip[1900.0]["tilt_deg"], 2), "tilt_2025": round(dip[2025.0]["tilt_deg"], 2)} say("\nB. Places, 1900 -> 2025: centred-dipole (CD) and apex latitude, Stormer cutoff (GV), TPL-007 alpha ratio") places = [] for name, lat, lon in PLACES: glat = geodetic_to_geocentric(lat) res = {} for ep in (1900.0, 2025.0): d = dip[ep] cd = cd_latitude(glat, lon, d["pole_lat"], d["pole_lon"]) ap = apex_latitude(E[ep], glat, lon)[0] res[ep] = (cd, ap, stormer_GV(d["M_Am2"], cd), stormer_GV(d["M_Am2"], ap)) (cd0, ap0, rc0, ra0), (cd1, ap1, rc1, ra1) = res[1900.0], res[2025.0] row = {"place": name, "lat": lat, "lon": lon, "cd_1900": round(cd0, 2), "cd_2025": round(cd1, 2), "d_cd": round(cd1 - cd0, 2), "apex_1900": round(ap0, 2), "apex_2025": round(ap1, 2), "d_apex": round(ap1 - ap0, 2), "rc_cd_1900": round(rc0, 2), "rc_cd_2025": round(rc1, 2), "d_rc_cd_pct": round((rc1 / rc0 - 1) * 100, 1), "rc_apex_1900": round(ra0, 2), "rc_apex_2025": round(ra1, 2), "d_rc_apex_pct": round((ra1 / ra0 - 1) * 100, 1), "alpha_ratio_cd": round(alpha_tpl007(cd1) / alpha_tpl007(cd0), 3), "alpha_ratio_apex": round(alpha_tpl007(ap1) / alpha_tpl007(ap0), 3)} places.append(row) say(f" {name:20s} CD {cd0:6.2f} -> {cd1:6.2f} ({cd1 - cd0:+5.2f}) apex {ap0:6.2f} -> {ap1:6.2f} ({ap1 - ap0:+5.2f}) " f"Rc CD {rc0:5.2f} -> {rc1:5.2f} ({row['d_rc_cd_pct']:+5.1f}%) Rc apex {ra0:5.2f} -> {ra1:5.2f} ({row['d_rc_apex_pct']:+5.1f}%)") OUT["places"] = places say("\nC. The apex proxy against Herbst, Kopp & Heber (2013, Table 1), GV in 1913, 1955, 1976, 1997") def coeffs_at(year): lo = max(e for e in epochs if e <= year) hi = min(e for e in epochs if e >= year) if lo == hi: return E[lo] w = (year - lo) / (hi - lo) return {k: (1 - w) * E[lo].get(k, 0.0) + w * E[hi].get(k, 0.0) for k in set(E[lo]) | set(E[hi])} OUT["herbst_check"] = {} for name, (lat, lon, vals) in HERBST.items(): glat = geodetic_to_geocentric(lat) ours = [] for yr in (1913.0, 1955.0, 1976.0, 1997.0): c = coeffs_at(yr) d = dipole(c) ours.append(stormer_GV(d["M_Am2"], apex_latitude(c, glat, lon)[0])) OUT["herbst_check"][name] = {"herbst": vals, "ours_apex": [round(v, 2) for v in ours]} say(f" {name:20s} Herbst {' '.join(f'{v:5.2f}' for v in vals)} ({(vals[3] / vals[0] - 1) * 100:+.0f}%) | " f"ours {' '.join(f'{v:5.2f}' for v in ours)} ({(ours[3] / ours[0] - 1) * 100:+.0f}%)") say("\nD. Auroral oval edge: shift toward the equator, degrees (Glassmeier et al. 2004: cos L ~ M^(alpha - 1/2))") scal = (("Siscoe & Chen (alpha=1/3)", -1.0 / 6.0), ("Siebert (alpha=0)", -0.5), ("Vogt & Glassmeier (alpha=1/2)", 0.0)) OUT["oval_shift_deg"] = {} for label, rr in (("2025 vs 1900", ratio), ("9% weaker", R_FINLAY)): OUT["oval_shift_deg"][label] = {} for edge in (65.0, 60.0, 55.0, 50.0, 45.0, 40.0, 30.0): shifts = {s: round(edge - math.degrees(math.acos(min(1.0, math.cos(math.radians(edge)) * rr ** ex))), 2) for s, ex in scal} OUT["oval_shift_deg"][label][str(edge)] = shifts say(f" {label}, edge at {edge:4.1f}: " + "; ".join(f"{s} {v:+.2f}" for s, v in shifts.items())) say("\nE. Magnetopause stand-off distance, R ~ M^(1/3) P^(-1/6), from 10 Earth radii") OUT["magnetopause"] = {"2025_vs_1900_RE": round(10 * ratio ** (1 / 3), 2), "9pct_weaker_RE": round(10 * R_FINLAY ** (1 / 3), 2), "pressure_2_to_x_nPa_RE": {p: round(10 * (p / 2.0) ** (-1 / 6), 2) for p in (10.0, 20.0, 50.0, 100.0)}} say(" ", OUT["magnetopause"]) say("\nF. Global 14C production Q (atoms cm^-2 s^-1), Kovaltsov, Mishev & Usoskin (2012) model") m19, m25 = M1900 / 1e22, M2025 / 1e22 OUT["production"] = {} for phi in (250, 300, 400, 500, 650, 800, 1000, 1300, 1500): q0, q1 = q_global(phi, m19), q_global(phi, m25) OUT["production"][str(phi)] = {"Q_1900": round(q0, 3), "Q_2025": round(q1, 3), "ratio": round(q1 / q0, 4)} say(f" phi {phi:5d} MV: Q(1900) {q0:.3f} Q(2025) {q1:.3f} ratio {q1 / q0:.4f}") OUT["production_check_650MV_M7.8"] = round(q_global(650, 7.8), 2) OUT["production_solar_cycle_ratio_300_vs_1300MV"] = round(q_global(300, m25) / q_global(1300, m25), 2) q78 = q_global(650, 7.8) OUT["production_check_M6_M10_vs_M7.8_at_650MV"] = [round(q_global(650, 6.0) / q78, 3), round(q_global(650, 10.0) / q78, 3)] OUT["production_check_300_1300MV_at_M7.8"] = [round(q_global(300, 7.8), 2), round(q_global(1300, 7.8), 2)] say(f" check: Q(650 MV, M = 7.8) = {OUT['production_check_650MV_M7.8']} (paper: 1951-2010 mean 1.64)") say(f" check: M = 6 and 10 against 7.8 at 650 MV: x{OUT['production_check_M6_M10_vs_M7.8_at_650MV'][0]} and " f"x{OUT['production_check_M6_M10_vs_M7.8_at_650MV'][1]} (paper: Holocene range within about +-15%)") say(f" check: Q(300 MV) and Q(1300 MV) at M = 7.8: {OUT['production_check_300_1300MV_at_M7.8']} (paper: 2.2 at the 2010 minimum, 1.1 at the 1990 maximum)") say(f" solar-cycle range at today's dipole, Q(300 MV) / Q(1300 MV) = {OUT['production_solar_cycle_ratio_300_vs_1300MV']}") say("\nG. Oulu neutron monitor, corrected count rates (Oulu Cosmic Ray Station, Sodankyla Geophysical Observatory)") ymap = {int(t[:4]): v for t, v in read_oulu(paths["oulu-nm-yearly-1965-2025.csv"])} cycles = [("cycle 20", range(1964, 1967), range(1968, 1972)), ("cycle 21", range(1975, 1979), range(1980, 1984)), ("cycle 22", range(1985, 1988), range(1989, 1992)), ("cycle 23", range(1995, 1999), range(2000, 2005)), ("cycle 24", range(2008, 2011), range(2012, 2016)), ("cycle 25", range(2019, 2021), range(2022, 2026))] OUT["oulu_solar_cycle_drop_pct"] = {} for label, ymin, ymax in cycles: hy = max((y for y in ymin if y in ymap), key=lambda y: ymap[y]) ly = min((y for y in ymax if y in ymap), key=lambda y: ymap[y]) drop = (ymap[hy] - ymap[ly]) / ymap[hy] * 100 OUT["oulu_solar_cycle_drop_pct"][label] = round(drop, 1) say(f" {label}: yearly {ymap[hy]:.0f} ({hy}) -> {ymap[ly]:.0f} ({ly}): {drop:.1f}% lower at solar maximum") OUT["oulu_forbush_pct"] = {} for name, b0, b1, label in (("oulu-nm-hourly-2003-10-15_2003-11-30.csv", "2003-10-15T00", "2003-10-27T23", "October 2003"), ("oulu-nm-hourly-2024-05-01_2024-05-31.csv", "2024-05-01T00", "2024-05-09T23", "May 2024")): h = read_oulu(paths[name]) base = [v for t, v in h if b0 <= t[:13] <= b1] b = sum(base) / len(base) tmin, vmin = min(h, key=lambda r: r[1]) OUT["oulu_forbush_pct"][label] = round((vmin / b - 1) * 100, 1) say(f" {label} Forbush decrease: baseline {b:.0f}, hourly minimum {vmin:.0f} at {tmin}: {(vmin / b - 1) * 100:.1f}%") say("\nH. TPL-007: what a change of latitude dL does to the benchmark field, exp(0.115 dL)") OUT["tpl007_factor"] = {str(dl): round(math.exp(0.115 * dl), 3) for dl in (-3, -2, -1, -0.5, 0.5, 1, 2, 3)} say(" ", OUT["tpl007_factor"]) say("\nI. Net change, 1900 -> 2025, in distance from a storm-time oval edge near 50 degrees (+ = farther)") delta50 = 50.0 - math.degrees(math.acos(math.cos(math.radians(50.0)) * ratio ** (-1.0 / 6.0))) say(f" the edge itself moves {delta50:.2f} degrees toward the equator (Siscoe & Chen scaling)") OUT["oval_edge_50_shift_deg"] = round(delta50, 2) OUT["net_distance_deg"] = {} for r in places: if r["place"] in ("Edmonton", "Winnipeg", "Minneapolis", "Chicago", "Montreal", "New York", "Washington DC", "Seattle", "Oslo", "Stockholm", "Helsinki", "Malmo", "Edinburgh", "London", "Berlin", "Paris", "Moscow"): net_cd, net_ap = -delta50 - r["d_cd"], -delta50 - r["d_apex"] OUT["net_distance_deg"][r["place"]] = {"cd": round(net_cd, 2), "apex": round(net_ap, 2)} say(f" {r['place']:14s} CD {net_cd:+5.2f} apex {net_ap:+5.2f}") say("\nJ. 2015-2020 change in |latitude| at four places, degrees a year (a fall = moving away from the auroral zone)") OUT["rates_2015_2020"] = {} for name, lat, lon in (("Edmonton", 53.55, -113.49), ("Leeds", 53.80, -1.55), ("Salekhard", 66.53, 66.60), ("Dunedin", -45.87, 170.50)): glat = geodetic_to_geocentric(lat) v = {ep: (cd_latitude(glat, lon, dip[ep]["pole_lat"], dip[ep]["pole_lon"]), apex_latitude(E[ep], glat, lon)[0]) for ep in (2015.0, 2020.0)} rcd = (abs(v[2020.0][0]) - abs(v[2015.0][0])) / 5.0 rap = (abs(v[2020.0][1]) - abs(v[2015.0][1])) / 5.0 OUT["rates_2015_2020"][name] = {"cd": round(rcd, 4), "apex": round(rap, 4)} say(f" {name:10s} CD {rcd:+.4f} apex {rap:+.4f}") if "--grid" in sys.argv: say("\nK. Global grid, 2.5 x 5 degrees, 20 km: area-weighted mean cutoff and the share of area below 1 GV") OUT["grid"] = {} for ep in (1900.0, 2025.0): c, d = E[ep], dip[ep] ws = rs = rcs = low = lowcd = 0.0 lat = -90 + 1.25 while lat < 90: w = math.cos(math.radians(lat)) lon = -180 + 2.5 while lon < 180: rc = stormer_GV(d["M_Am2"], apex_latitude(c, lat, lon, step_km=40.0)[0]) rc2 = stormer_GV(d["M_Am2"], cd_latitude(lat, lon, d["pole_lat"], d["pole_lon"])) ws, rs, rcs = ws + w, rs + w * rc, rcs + w * rc2 low, lowcd = low + w * (rc < 1.0), lowcd + w * (rc2 < 1.0) lon += 5.0 lat += 2.5 OUT["grid"][str(ep)] = {"mean_proxy_GV": round(rs / ws, 3), "mean_cd_GV": round(rcs / ws, 3), "area_below_1GV_proxy_pct": round(low / ws * 100, 1), "area_below_1GV_cd_pct": round(lowcd / ws * 100, 1)} say(f" {ep}: {OUT['grid'][str(ep)]}") with open(os.path.join(folder, "claim6-scalings.json"), "w", encoding="utf-8", newline="\n") as f: json.dump(OUT, f, indent=1) print("Wrote", os.path.join(folder, "claim6-scalings.json")) if __name__ == "__main__": main()