#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ milankovitch-insolation.py Science Journaling Club, Volume 2 Issue 3, Spring 2026, "Earth's Energy Budget". THE QUESTION ------------ Ice ages are paced by slow changes in Earth's orbit. Working from orbital geometry alone, how much does summer sunlight at high northern latitudes actually vary, and does the computed curve line up with when the ice sheets melted? WHAT THIS PROGRAM IS -------------------- This is a computation, not an observation. The club has no telescope, no radiometer, no ice core and no field site. Nothing in this file was measured by us in the physical world. Every number below is produced by evaluating the top-of-atmosphere insolation formula on orbital elements that somebody else computed and published, using code we wrote ourselves. Where the file says "measured" it means "evaluated from our own arithmetic". The one place a random number generator enters is a Monte Carlo integration used as an independent check on the analytic formula, and its seed is printed below and fixed in the source. DATA ---- Orbital elements come from the La2004 numerical solution of Laskar et al. (2004), file INSOLN.LA2004.BTL.ASC, retrieved 2026-09-14 from http://vo.imcce.fr/insola/earth/online/earth/La2004/INSOLN.LA2004.BTL.ASC This is the tabulated solution itself, not a truncated reconstruction. The file holds 51,001 rows at 1 kyr spacing covering 0 to -51 Myr, with columns t (kyr, negative into the past) e eccentricity, dimensionless eps obliquity, radians, w.r.t. the fixed ecliptic of J2000 varpi longitude of perihelion from the moving vernal equinox, rad The file is cached under analysis/data/ and its SHA-256 is printed at run time. We use rows 0 to -1000 kyr. Nothing in the file is modified. THE MODEL --------- Daily-mean insolation at the top of the atmosphere, for latitude phi and solar longitude lambda (the Sun's apparent longitude measured from the northward equinox): Q = (S0 / pi) (a/r)^2 [ H0 sin(phi) sin(delta) + cos(phi) cos(delta) sin(H0) ] with sin(delta) = sin(eps) sin(lambda) solar declination H0 = arccos( -tan(phi) tan(delta) ) half-day / sunset hour angle r/a = (1 - e^2) / (1 - e cos(lambda - varpi)) Sun-Earth distance H0 is clipped to [0, pi], which handles polar night (Q = 0) and polar day (Q = S0 (a/r)^2 sin(phi) sin(delta)) without a special case. Calendar. Position in the orbit is converted to elapsed time by Kepler's equation. With nu the true anomaly of the Earth, nu = lambda - varpi - pi, E = 2 atan2( sqrt(1-e) sin(nu/2), sqrt(1+e) cos(nu/2) ) M = E - e sin(E) t = (M - M_equinox) / (2 pi) * 365.2422 days so "day of year" in this file means days elapsed since the northward equinox. The June solstice is the astronomical event lambda = 90 degrees, which carries no calendar convention at all, and that is what the headline number uses. WHAT WE CHECK AGAINST --------------------- Analytic results, derived from the same formula but by a different route, so a coding error in the numerical path shows up as a disagreement: * Global annual mean insolation = S0 / (4 sqrt(1 - e^2)) exactly. * Annual mean at the equator = 2 S0 E(sin eps) / (pi^2 sqrt(1 - e^2)), with E the complete elliptic integral of the second kind, evaluated here by the arithmetic-geometric mean, which shares no code with the quadrature. * Annual mean at either pole = S0 sin(eps) / (pi sqrt(1 - e^2)). * Equinox at the equator = (S0/pi) (a/r)^2. * Solstice at the summer pole = S0 (a/r)^2 sin(eps). * A Monte Carlo integration over the sphere and over the year that never touches H0 or the daily-mean integral, recovering S0 / (4 sqrt(1-e^2)). * Published present-day values quoted in the article's references. ASSUMPTIONS, STATED PLAINLY --------------------------- * Top of atmosphere only. No atmosphere, no clouds, no albedo, no greenhouse effect, no ocean, no ice. This program computes sunlight arriving, and says nothing about temperature. * The Sun is a point source of constant output S0 = 1361 W/m^2 (Kopp & Lean 2011). Solar luminosity is held fixed for a million years. It was not. * Earth is a sphere. No oblateness in the insolation geometry. * The orbit is a fixed Keplerian ellipse within each 1 kyr step; the orbital elements change only between steps. * Obliquity is referred to the fixed ecliptic of J2000, as La2004 tabulates it. The difference from the obliquity of date is small but not zero. * The tropical year is held at 365.2422 days throughout. Over a million years it changes by a few parts in 10^5 through tidal deceleration; that shifts the day-of-year calendar, not the solstice insolation. * No lunar effect on obliquity beyond what La2004 already includes, and no feedback of climate on the orbit. LIMITATIONS, STATED PLAINLY --------------------------- * This is a forcing calculation, not an ice age model. There is no ice sheet, no carbon cycle, no ocean circulation and no lag of any kind. When the text says a termination "lines up" with an insolation rise it means the two dates are near each other, and near is doing a lot of work. * The 65 N summer solstice is one choice out of many for "summer sunlight". Section 7 computes three alternatives and they do not agree about which orbital parameter matters most. That disagreement is a result, not noise. * Termination ages are published values read out of the cited author's own age-model table. We did not date anything. * La2004 itself is chaotic beyond about 40 Myr. Over 1 Myr it is not, and the authors state a precision far better than anything that matters here. USAGE ----- python milankovitch-insolation.py > milankovitch-insolation-output.txt Set SJC_FIGDATA to a path to also dump the plotted series as JSON. """ import hashlib import json import math import os import sys import time import numpy as np T_START = time.time() # -------------------------------------------------------------------------- # constants and configuration # -------------------------------------------------------------------------- SEED = 20260320 # the only RNG in the file; northward equinox 2026 S0 = 1361.0 # W/m^2, total solar irradiance, Kopp & Lean (2011) S0_OLD = 1365.0 # the value used by most of the older literature TROPICAL_YEAR = 365.2422 # days DEG = math.pi / 180.0 DATA_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data") DATA_FILE = os.path.join(DATA_DIR, "INSOLN.LA2004.BTL.ASC") DATA_URL = ("http://vo.imcce.fr/insola/earth/online/earth/La2004/" "INSOLN.LA2004.BTL.ASC") RETRIEVED = "2026-09-14" RULE = "=" * 110 THIN = "-" * 110 def head(title): print() print(RULE) print(title) print(RULE) def row(label, got, ref, unit="", tol=None, note=""): """Print club value beside reference value with the difference.""" d = got - ref rel = (d / ref * 100.0) if ref != 0 else float("nan") flag = "" if tol is not None: flag = " OK" if abs(rel) <= tol else " ** DISAGREES **" print(" %-42s %14.6f %14.6f %13.3e %9.4f%%%s%s" % (label, got, ref, d, rel, flag, (" " + note) if note else "")) # -------------------------------------------------------------------------- # orbital geometry # -------------------------------------------------------------------------- def ecc_anomaly_from_true(nu, e): return 2.0 * np.arctan2(np.sqrt(1.0 - e) * np.sin(nu / 2.0), np.sqrt(1.0 + e) * np.cos(nu / 2.0)) def true_from_ecc_anomaly(E, e): return 2.0 * np.arctan2(np.sqrt(1.0 + e) * np.sin(E / 2.0), np.sqrt(1.0 - e) * np.cos(E / 2.0)) def mean_anomaly_at_lambda(lam, e, varpi): """Mean anomaly of the Earth when the Sun's apparent longitude is lam.""" nu = lam - varpi - math.pi E = ecc_anomaly_from_true(nu, e) return E - e * np.sin(E) def kepler_solve(M, e, tol=1e-14, itmax=80): """Newton solve of M = E - e sin E. Returns E.""" M = np.asarray(M, dtype=float) E = M + e * np.sin(M) # good starting guess for small e for _ in range(itmax): f = E - e * np.sin(E) - M fp = 1.0 - e * np.cos(E) step = f / fp E = E - step if np.max(np.abs(step)) < tol: break return E def lambda_from_day(day, e, varpi): """Sun's apparent longitude, given days elapsed since northward equinox.""" M_eq = mean_anomaly_at_lambda(0.0, e, varpi) M = M_eq + 2.0 * math.pi * np.asarray(day, dtype=float) / TROPICAL_YEAR E = kepler_solve(M, e) nu = true_from_ecc_anomaly(E, e) return np.mod(nu + varpi + math.pi, 2.0 * math.pi) def day_from_lambda(lam, e, varpi): """Days elapsed since northward equinox, given the Sun's longitude.""" M_eq = mean_anomaly_at_lambda(0.0, e, varpi) M = mean_anomaly_at_lambda(lam, e, varpi) return np.mod(M - M_eq, 2.0 * math.pi) / (2.0 * math.pi) * TROPICAL_YEAR def dist_ratio(lam, e, varpi): """r/a, the Sun-Earth distance in units of the semi-major axis.""" return (1.0 - e * e) / (1.0 - e * np.cos(lam - varpi)) def daily_insolation(lat_deg, lam, e, eps, varpi, s0=S0): """Daily-mean top-of-atmosphere insolation, W/m^2. lat_deg latitude in degrees, scalar or array lam Sun's apparent longitude in radians, scalar or array e eccentricity eps obliquity in radians varpi longitude of perihelion from the moving vernal equinox, radians """ phi = np.asarray(lat_deg, dtype=float) * DEG lam = np.asarray(lam, dtype=float) rho = dist_ratio(lam, e, varpi) sin_d = math.sin(eps) * np.sin(lam) cos_d = np.sqrt(np.maximum(0.0, 1.0 - sin_d * sin_d)) # tan(phi) at the pole is huge but finite in IEEE arithmetic; the clip below # turns that into H0 = 0 or pi, which is the correct limit. with np.errstate(invalid="ignore", divide="ignore"): cos_h0 = -np.tan(phi) * (sin_d / np.where(cos_d == 0.0, 1e-300, cos_d)) h0 = np.arccos(np.clip(cos_h0, -1.0, 1.0)) q = (s0 / math.pi) * rho ** -2 * (h0 * np.sin(phi) * sin_d + np.cos(phi) * cos_d * np.sin(h0)) return np.maximum(q, 0.0) def daily_insolation_laskar(lat_deg, lam, e, eps, varpi, s0=S0): """The same quantity, transcribed line for line from Laskar's own Fortran. Subroutine cwj in insolsub.f of the La2004 distribution, cached alongside the data file. That routine branches explicitly on three cases (ordinary latitude, no sunset, no sunrise) instead of clipping the hour angle, and it takes pibar = varpi + pi in a geocentric frame, which is why the true anomaly below reads wd - pibar. Kept here only as a cross-check: if our clipped one-liner and Laskar's three-branch version ever disagree, one of them is wrong. """ phi = np.asarray(lat_deg, dtype=float) * DEG wd = np.asarray(lam, dtype=float) pibar = varpi + math.pi v = wd - pibar rho = (1.0 - e ** 2) / (1.0 + e * np.cos(v)) sind = math.sin(eps) * np.sin(wd) delta = np.arcsin(sind) out = np.zeros(np.broadcast(phi, wd).shape, dtype=float) aux = math.pi / 2.0 - np.abs(delta) ordinary = (phi > -aux) & (phi < aux) a1 = math.pi / 2.0 - delta a2 = math.pi / 2.0 + delta nosunset = (~ordinary) & ((phi >= a1) | (phi <= -a2)) with np.errstate(invalid="ignore", divide="ignore"): cho = -np.tan(phi) * np.tan(delta) ho = np.arccos(np.clip(cho, -1.0, 1.0)) wo = (ho * np.sin(phi) * sind + np.cos(phi) * np.cos(delta) * np.sin(ho)) wo = wo * s0 / (math.pi * rho ** 2) wp = s0 * np.sin(phi) * sind / rho ** 2 out = np.where(ordinary, wo, np.where(nosunset, wp, 0.0)) return out def annual_mean(lat_deg, e, eps, varpi, n=2 ** 14, s0=S0): """Annual mean insolation, by quadrature in lambda with the Kepler weight. Time average = (1 / (2 pi sqrt(1-e^2))) * integral of Q (r/a)^2 dlambda, because r^2 dlambda/dt is constant (Kepler's second law). """ lam = (np.arange(n) + 0.5) / n * 2.0 * math.pi q = daily_insolation(lat_deg, lam, e, eps, varpi, s0=s0) w = dist_ratio(lam, e, varpi) ** 2 integ = np.sum(q * w) * (2.0 * math.pi / n) return integ / (2.0 * math.pi * math.sqrt(1.0 - e * e)) def elliptic_E(k): """Complete elliptic integral of the second kind, by the AGM. E(k) = K(k) * (1 - sum_{n>=0} 2^(n-1) c_n^2), with the AGM recursion a_{n+1} = (a_n+b_n)/2, b_{n+1} = sqrt(a_n b_n), c_{n+1} = (a_n-b_n)/2. Shares no code with the quadrature it is used to check. """ a, b, c = 1.0, math.sqrt(1.0 - k * k), k s = 0.5 * c * c p = 1.0 for _ in range(60): a, b, c = 0.5 * (a + b), math.sqrt(a * b), 0.5 * (a - b) p *= 2.0 s += 0.5 * p * c * c if abs(c) < 1e-17: break K = math.pi / (2.0 * a) return K * (1.0 - s) # -------------------------------------------------------------------------- # spectral tools # -------------------------------------------------------------------------- def hann_periodogram(x, dt): """One-sided Hann-windowed periodogram. Returns freq, power. Power is normalised so that the sum over the one-sided spectrum equals the variance of the detrended, windowed series, which makes band shares add up. """ x = np.asarray(x, dtype=float) n = x.size t = np.arange(n) * dt # remove mean and linear trend A = np.vstack([np.ones(n), t]).T coef, *_ = np.linalg.lstsq(A, x, rcond=None) xd = x - A @ coef w = 0.5 - 0.5 * np.cos(2.0 * math.pi * np.arange(n) / n) xw = xd * w X = np.fft.rfft(xw) p = np.abs(X) ** 2 p[1:] *= 2.0 if n % 2 == 0: p[-1] /= 2.0 p /= np.sum(w ** 2) f = np.fft.rfftfreq(n, d=dt) return f, p def band_share(f, p, p_lo, p_hi): """Fraction of total power between periods p_lo and p_hi (kyr).""" lo, hi = 1.0 / p_hi, 1.0 / p_lo m = (f >= lo) & (f <= hi) sub = p.copy() sub[0] = 0.0 tot = np.sum(sub) if not np.any(m) or tot <= 0: return 0.0, float("nan"), 0.0 bp = float(np.sum(sub[m])) k = int(np.argmax(np.where(m, sub, -1.0))) peak_period = 1.0 / f[k] if f[k] > 0 else float("inf") return bp / tot, peak_period, bp # -------------------------------------------------------------------------- # 0. provenance # -------------------------------------------------------------------------- print(RULE) print("COMPUTING THE SUNLIGHT THAT STARTED AND ENDED THE ICE AGES") print("Top-of-atmosphere insolation from orbital elements, past 1 Myr") print("Science Journaling Club, Volume 2 Issue 3, Spring 2026") print(RULE) print() print("This output is arithmetic. No telescope, radiometer, ice core or field") print("site is involved anywhere in this file. The orbital elements are taken") print("from a published numerical solution; everything else is computed here.") print() print(" python version : %s" % sys.version.split()[0]) print(" numpy version : %s" % np.__version__) print(" master seed : %d" % SEED) print(" generator : numpy PCG64, seeded once") print(" solar constant S0 : %.1f W/m^2 (Kopp & Lean 2011)" % S0) print(" tropical year : %.4f days, held fixed" % TROPICAL_YEAR) if not os.path.exists(DATA_FILE): print() print(" !! ORBITAL DATA FILE NOT FOUND: %s" % DATA_FILE) print(" !! Download it from %s" % DATA_URL) print(" !! Nothing below can be computed without it. Stopping.") sys.exit(1) with open(DATA_FILE, "rb") as fh: raw = fh.read() sha = hashlib.sha256(raw).hexdigest() txt = raw.decode("ascii").replace("D+", "E+").replace("D-", "E-") tab = np.array([[float(v) for v in line.split()] for line in txt.strip().splitlines()]) print() print(" orbital solution : La2004, Laskar et al. (2004)") print(" file : %s" % os.path.basename(DATA_FILE)) print(" source URL : %s" % DATA_URL) print(" retrieved : %s" % RETRIEVED) print(" bytes : %d" % len(raw)) print(" sha256 : %s" % sha) print(" rows in file : %d" % tab.shape[0]) print(" time span in file : %.0f to %.0f kyr" % (tab[0, 0], tab[-1, 0])) print(" columns : t(kyr) e eps(rad) varpi(rad)") print(" row 0 (present) : e=%.7f eps=%.6f rad = %.5f deg varpi=%.5f rad = %.4f deg" % (tab[0, 1], tab[0, 2], tab[0, 2] / DEG, tab[0, 3], tab[0, 3] / DEG)) NKYR = 1000 sel = tab[:NKYR + 1] t_kyr = np.abs(sel[:, 0]) # 0 .. 1000, positive = kyr before present ecc = sel[:, 1] obl = sel[:, 2] vpi = sel[:, 3] prec = ecc * np.sin(vpi + math.pi) # climatic precession index e sin(omega-tilde) E0, EPS0, VPI0 = float(ecc[0]), float(obl[0]), float(vpi[0]) print() print(" rows used : %d (0 to %d kyr before present)" % (sel.shape[0], NKYR)) print(" eccentricity range : %.6f to %.6f" % (ecc.min(), ecc.max())) print(" obliquity range : %.4f to %.4f deg" % (obl.min() / DEG, obl.max() / DEG)) print(" precession index e sin(w) : %+.5f to %+.5f" % (prec.min(), prec.max())) print(" present precession index : %+.5f" % prec[0]) # -------------------------------------------------------------------------- # 1. the formula against closed-form results # -------------------------------------------------------------------------- head("SECTION 1. THE INSOLATION FORMULA AGAINST CLOSED-FORM RESULTS") print() print("Present-day orbit throughout this section: e = %.7f, eps = %.5f deg, varpi = %.4f deg." % (E0, EPS0 / DEG, VPI0 / DEG)) print("Every 'club' number is quadrature or summation over the formula in the") print("docstring. Every 'analytic' number is a closed form derived by hand.") print() print(" %-42s %14s %14s %13s %9s" % ("quantity", "club", "analytic", "difference", "relative")) print(THIN) # global annual mean lats = np.linspace(-90.0, 90.0, 1801) qbar = np.array([annual_mean(la, E0, EPS0, VPI0) for la in lats]) w = np.cos(lats * DEG) glob = np.trapezoid(qbar * w, lats * DEG) / np.trapezoid(w, lats * DEG) glob_ref = S0 / (4.0 * math.sqrt(1.0 - E0 * E0)) row("global annual mean insolation (W/m2)", glob, glob_ref, tol=0.01) # equator annual mean kk = math.sin(EPS0) eq = annual_mean(0.0, E0, EPS0, VPI0) eq_ref = 2.0 * S0 * elliptic_E(kk) / (math.pi ** 2 * math.sqrt(1.0 - E0 * E0)) row("annual mean at the equator (W/m2)", eq, eq_ref, tol=0.001) # pole annual mean npole = annual_mean(90.0, E0, EPS0, VPI0) spole = annual_mean(-90.0, E0, EPS0, VPI0) pole_ref = S0 * math.sin(EPS0) / (math.pi * math.sqrt(1.0 - E0 * E0)) row("annual mean at the north pole (W/m2)", npole, pole_ref, tol=0.01) row("annual mean at the south pole (W/m2)", spole, pole_ref, tol=0.01) # equinox at the equator lam_ve = 0.0 q_eq_ve = float(daily_insolation(0.0, lam_ve, E0, EPS0, VPI0)) ref_eq_ve = (S0 / math.pi) * dist_ratio(lam_ve, E0, VPI0) ** -2 row("northward equinox, equator (W/m2)", q_eq_ve, ref_eq_ve, tol=1e-8) lam_ae = math.pi q_eq_ae = float(daily_insolation(0.0, lam_ae, E0, EPS0, VPI0)) ref_eq_ae = (S0 / math.pi) * dist_ratio(lam_ae, E0, VPI0) ** -2 row("southward equinox, equator (W/m2)", q_eq_ae, ref_eq_ae, tol=1e-8) # solstice at the summer pole lam_js = math.pi / 2.0 q_np_js = float(daily_insolation(90.0, lam_js, E0, EPS0, VPI0)) ref_np_js = S0 * dist_ratio(lam_js, E0, VPI0) ** -2 * math.sin(EPS0) row("June solstice, north pole (W/m2)", q_np_js, ref_np_js, tol=1e-8) lam_ds = 3.0 * math.pi / 2.0 q_sp_ds = float(daily_insolation(-90.0, lam_ds, E0, EPS0, VPI0)) ref_sp_ds = S0 * dist_ratio(lam_ds, E0, VPI0) ** -2 * math.sin(EPS0) row("December solstice, south pole (W/m2)", q_sp_ds, ref_sp_ds, tol=1e-8) # polar night must be exactly zero q_np_ds = float(daily_insolation(90.0, lam_ds, E0, EPS0, VPI0)) print(" %-42s %14.6f %14.6f %13.3e %9s" % ("December solstice, north pole (W/m2)", q_np_ds, 0.0, q_np_ds - 0.0, "exact")) # hemispheric symmetry. Q(phi, lam) and Q(-phi, lam+180) differ only through # the Sun-Earth distance, so dividing that out must make them identical. q_a = float(daily_insolation(45.0, lam_js, E0, EPS0, VPI0)) * dist_ratio(lam_js, E0, VPI0) ** 2 q_b = float(daily_insolation(-45.0, lam_ds, E0, EPS0, VPI0)) * dist_ratio(lam_ds, E0, VPI0) ** 2 row("45N/45S solstice, distance divided out", q_a, q_b, tol=1e-10, note="mirror symmetry") print(" %-42s %14.6f %14.6f %13.3e %9s" % ("the same two before dividing it out", float(daily_insolation(45.0, lam_js, E0, EPS0, VPI0)), float(daily_insolation(-45.0, lam_ds, E0, EPS0, VPI0)), float(daily_insolation(45.0, lam_js, E0, EPS0, VPI0)) - float(daily_insolation(-45.0, lam_ds, E0, EPS0, VPI0)), "perihelion")) print() print(" Against Laskar's own reference implementation. Subroutine cwj of") print(" insolsub.f from the La2004 distribution, transcribed into Python and") print(" cached beside the data file, branches on three cases instead of") print(" clipping the hour angle. Evaluated on a grid of %d latitudes by %d" % (181, 721)) print(" solar longitudes, at the present-day orbit:") LAT_CHK = np.linspace(-90.0, 90.0, 181) LAM_CHK = np.linspace(0.0, 2.0 * math.pi, 721) d_max = 0.0 for la in LAT_CHK: a_ = daily_insolation(la, LAM_CHK, E0, EPS0, VPI0) b_ = daily_insolation_laskar(la, LAM_CHK, E0, EPS0, VPI0) d_max = max(d_max, float(np.max(np.abs(a_ - b_)))) print(" largest absolute difference anywhere on the grid : %.3e W/m2" % d_max) print(" %s" % ("identical to machine precision" if d_max < 1e-9 else "** THE TWO IMPLEMENTATIONS DISAGREE **")) print(" Laskar's subroutine wam gives the annual mean as so/(4 sqrt(1-e^2)),") print(" which is the closed form used in the first row of the table above.") print() print(" complete elliptic integral E(sin eps) by AGM : %.15f" % elliptic_E(kk)) lam_q = (np.arange(2 ** 20) + 0.5) / 2 ** 20 * 2.0 * math.pi print(" same integral by brute-force quadrature : %.15f" % (np.sum(np.sqrt(1.0 - kk * kk * np.sin(lam_q) ** 2)) * (2 * math.pi / 2 ** 20) / 4.0)) # -------------------------------------------------------------------------- # 2. Monte Carlo check that never uses the hour angle # -------------------------------------------------------------------------- head("SECTION 2. MONTE CARLO CHECK OF THE WHOLE CHAIN") print() print("The daily-mean formula folds the day into an analytic integral over the") print("hour angle. If that integral is wrong, every number above is wrong in the") print("same way and the closed forms in Section 1 will not catch it, because they") print("come from the same integral. So here is a check that does not use it at") print("all: scatter points uniformly over the sphere and uniformly in TIME, work") print("out the solar zenith angle at each one from scratch, and average the") print("sunlight landing on a flat patch of ground facing the sky.") print() print("Uniform in time means uniform in mean anomaly M, because M advances at a") print("constant rate by construction. Kepler's equation then gives the position.") print() rng = np.random.default_rng(SEED) NMC = 4_000_000 M_mc = rng.uniform(0.0, 2.0 * math.pi, NMC) E_mc = kepler_solve(M_mc, E0) nu_mc = true_from_ecc_anomaly(E_mc, E0) lam_mc = nu_mc + VPI0 + math.pi rho_mc = (1.0 - E0 * E0) / (1.0 + E0 * np.cos(nu_mc)) sin_d_mc = math.sin(EPS0) * np.sin(lam_mc) cos_d_mc = np.sqrt(np.maximum(0.0, 1.0 - sin_d_mc ** 2)) sin_phi_mc = rng.uniform(-1.0, 1.0, NMC) # uniform in area on the sphere cos_phi_mc = np.sqrt(np.maximum(0.0, 1.0 - sin_phi_mc ** 2)) h_mc = rng.uniform(-math.pi, math.pi, NMC) # hour angle, uniform over the day cos_z = sin_phi_mc * sin_d_mc + cos_phi_mc * cos_d_mc * np.cos(h_mc) flux = S0 * rho_mc ** -2 * np.maximum(0.0, cos_z) mc = float(np.mean(flux)) mc_se = float(np.std(flux, ddof=1) / math.sqrt(NMC)) print(" %-42s %14s %14s %13s %9s" % ("quantity", "club", "analytic", "difference", "relative")) print(THIN) row("Monte Carlo global mean (W/m2)", mc, glob_ref) z_mc = (mc - glob_ref) / mc_se print(" standard error of the Monte Carlo mean : %.4f W/m2" % mc_se) print(" difference in standard errors : %+.2f sigma %s" % (z_mc, "OK" if abs(z_mc) < 3.0 else "** DISAGREES **")) print(" (the verdict here is the z score, not a tolerance: a Monte Carlo") print(" estimate is allowed to miss, and is only wrong if it misses by more") print(" than its own uncertainty says it should.)") print(" samples : %s" % f"{NMC:,}") print() print(" The same target reached by the quadrature path : %.6f W/m2" % glob) print(" S0 / 4 (the textbook number, circular orbit) : %.6f W/m2" % (S0 / 4.0)) print(" eccentricity raises it by a factor 1/sqrt(1-e^2) = %.8f" % (1.0 / math.sqrt(1 - E0 ** 2))) # -------------------------------------------------------------------------- # 3. against published present-day values # -------------------------------------------------------------------------- head("SECTION 3. AGAINST PUBLISHED PRESENT-DAY VALUES") print() print("These are numbers other people have printed, entered here by hand from the") print("sources cited in the article. They are rounded as published, so agreement") print("to the last printed digit is all that can be asked.") print() print(" %-42s %14s %14s %13s %9s" % ("quantity", "club", "published", "difference", "relative")) print(THIN) q65 = float(daily_insolation(65.0, lam_js, E0, EPS0, VPI0)) row("65N June solstice, today (W/m2)", q65, 480.0, tol=1.0, note="Berger 1978 tables") q65_old = float(daily_insolation(65.0, lam_js, E0, EPS0, VPI0, s0=S0_OLD)) print(" %-42s %14.6f (the same with S0 = %.0f, the older constant)" % ("65N June solstice, today (W/m2)", q65_old, S0_OLD)) q00 = float(daily_insolation(0.0, lam_js, E0, EPS0, VPI0)) print(" %-42s %14.6f" % ("equator, June solstice (W/m2)", q00)) print(" %-42s %14.6f" % ("north pole / equator at June solstice", q_np_js / q00)) row("annual mean, equator (W/m2)", eq, 416.0, tol=0.5, note="Hartmann 2016") row("annual mean, pole (W/m2)", npole, 173.0, tol=1.0, note="Hartmann 2016") row("obliquity today (deg)", EPS0 / DEG, 23.4393, tol=0.01, note="IAU value") row("eccentricity today", E0, 0.016708, tol=0.05, note="standard element") row("longitude of perihelion (deg)", VPI0 / DEG, 102.947, tol=0.05, note="standard element") print() print(" Days from northward equinox to each astronomical event, present orbit:") for nm, lm in (("northward equinox", 0.0), ("June solstice", math.pi / 2), ("southward equinox", math.pi), ("December solstice", 3 * math.pi / 2)): print(" %-22s lambda = %5.1f deg day %7.3f r/a = %.6f" % (nm, lm / DEG, float(day_from_lambda(lm, E0, VPI0)), dist_ratio(lm, E0, VPI0))) half_n = float(day_from_lambda(math.pi, E0, VPI0)) print(" Length of the northern summer half year (equinox to equinox):") print(" %.3f days, against %.3f for the southern half. Difference %.3f days." % (half_n, TROPICAL_YEAR - half_n, 2 * half_n - TROPICAL_YEAR)) print(" Published figure for the present asymmetry is about 7.5 days.") # -------------------------------------------------------------------------- # 4. convergence of the quadrature # -------------------------------------------------------------------------- head("SECTION 4. CONVERGENCE OF THE ANNUAL-MEAN QUADRATURE") print() print("How many points around the orbit does the annual mean need? The answer") print("depends on latitude, because the integrand is smooth at the equator and") print("has a kink at high latitude where polar night begins and ends. Smooth") print("periodic integrands converge faster than any power of N on a uniform") print("grid. Kinked ones do not.") print() print(" %8s %18s %18s %18s" % ("N", "|err| equator", "|err| 65N", "|err| pole")) print(THIN) conv = {"n": [], "eq": [], "l65": [], "pole": []} ref_eq = annual_mean(0.0, E0, EPS0, VPI0, n=2 ** 22) ref_65 = annual_mean(65.0, E0, EPS0, VPI0, n=2 ** 22) ref_pl = annual_mean(90.0, E0, EPS0, VPI0, n=2 ** 22) for pw in range(3, 19): n = 2 ** pw e_eq = abs(annual_mean(0.0, E0, EPS0, VPI0, n=n) - ref_eq) e_65 = abs(annual_mean(65.0, E0, EPS0, VPI0, n=n) - ref_65) e_pl = abs(annual_mean(90.0, E0, EPS0, VPI0, n=n) - ref_pl) print(" %8d %18.3e %18.3e %18.3e" % (n, e_eq, e_65, e_pl)) conv["n"].append(n) conv["eq"].append(e_eq) conv["l65"].append(e_65) conv["pole"].append(e_pl) print() print(" Reference values use N = %d. The equator falls to machine precision by" % 2 ** 22) print(" a few dozen points. 65N does not, and that is the polar-night kink.") print(" Everything else in this file uses N = %d." % 2 ** 14) # -------------------------------------------------------------------------- # 5. one parameter at a time # -------------------------------------------------------------------------- head("SECTION 5. AMPLITUDE CONTRIBUTED BY EACH ORBITAL PARAMETER ALONE") print() e_bar = float(np.mean(ecc)) eps_bar = float(np.mean(obl)) print(" Million-year means used as the holding values:") print(" mean eccentricity %.6f" % e_bar) print(" mean obliquity %.5f deg" % (eps_bar / DEG)) print(" perihelion angle swept over the full circle") print() print("A. Sweep one parameter over its full past-million-year range, hold the") print(" other two at the million-year mean, and record the 65N June solstice") print(" insolation at each end. This is a range, not a variance.") print() print(" %-30s %12s %12s %12s %12s %10s" % ("parameter swept", "low", "high", "Q low", "Q high", "span W/m2")) print(THIN) sweeps = {} lo, hi = float(ecc.min()), float(ecc.max()) qlo = float(daily_insolation(65.0, lam_js, lo, eps_bar, VPI0)) qhi = float(daily_insolation(65.0, lam_js, hi, eps_bar, VPI0)) print(" %-30s %12.6f %12.6f %12.3f %12.3f %10.3f" % ("eccentricity (varpi today)", lo, hi, qlo, qhi, abs(qhi - qlo))) ee = np.linspace(lo, hi, 200) sweeps["ecc"] = {"x": ee.tolist(), "q": [float(daily_insolation(65.0, lam_js, v, eps_bar, VPI0)) for v in ee]} span_e = abs(qhi - qlo) lo2, hi2 = float(obl.min()), float(obl.max()) qlo2 = float(daily_insolation(65.0, lam_js, e_bar, lo2, VPI0)) qhi2 = float(daily_insolation(65.0, lam_js, e_bar, hi2, VPI0)) print(" %-30s %12.5f %12.5f %12.3f %12.3f %10.3f" % ("obliquity (deg)", lo2 / DEG, hi2 / DEG, qlo2, qhi2, abs(qhi2 - qlo2))) oo = np.linspace(lo2, hi2, 200) sweeps["obl"] = {"x": (oo / DEG).tolist(), "q": [float(daily_insolation(65.0, lam_js, e_bar, v, VPI0)) for v in oo]} span_o = abs(qhi2 - qlo2) vv = np.linspace(0.0, 2.0 * math.pi, 361) qv = np.array([float(daily_insolation(65.0, lam_js, e_bar, eps_bar, v)) for v in vv]) print(" %-30s %12.1f %12.1f %12.3f %12.3f %10.3f" % ("perihelion angle (deg), e=mean", 0.0, 360.0, qv.min(), qv.max(), qv.max() - qv.min())) sweeps["prec_mean_e"] = {"x": (vv / DEG).tolist(), "q": qv.tolist()} span_p_mean = float(qv.max() - qv.min()) qv2 = np.array([float(daily_insolation(65.0, lam_js, hi, eps_bar, v)) for v in vv]) print(" %-30s %12.1f %12.1f %12.3f %12.3f %10.3f" % ("perihelion angle (deg), e=max", 0.0, 360.0, qv2.min(), qv2.max(), qv2.max() - qv2.min())) sweeps["prec_max_e"] = {"x": (vv / DEG).tolist(), "q": qv2.tolist()} span_p_max = float(qv2.max() - qv2.min()) qv3 = np.array([float(daily_insolation(65.0, lam_js, 0.0, eps_bar, v)) for v in vv]) print(" %-30s %12.1f %12.1f %12.3f %12.3f %10.3f" % ("perihelion angle (deg), e=0", 0.0, 360.0, qv3.min(), qv3.max(), qv3.max() - qv3.min())) span_p_zero = float(qv3.max() - qv3.min()) print() print(" That last row is the whole point of the coupling. With a circular orbit") print(" the perihelion angle means nothing, because there is no perihelion. The") print(" precession lever is very nearly proportional to eccentricity, so what") print(" looks like a 23 kyr signal is a 23 kyr carrier inside a 100 kyr envelope.") print() print("B. Let one parameter vary in time over the past million years and hold the") print(" other two at the million-year mean. Standard deviation of the resulting") print(" 65N June solstice series, against the full series with all three.") print() def series(e_arr, eps_arr, v_arr): return np.array([float(daily_insolation(65.0, lam_js, e_arr[i], eps_arr[i], v_arr[i])) for i in range(sel.shape[0])]) const_e = np.full(sel.shape[0], e_bar) const_o = np.full(sel.shape[0], eps_bar) const_v = np.full(sel.shape[0], VPI0) q_full = series(ecc, obl, vpi) q_e_only = series(ecc, const_o, const_v) q_o_only = series(const_e, obl, const_v) q_p_only = series(const_e, const_o, vpi) q_ep_only = series(ecc, const_o, vpi) print(" %-46s %12s %12s %12s" % ("series", "std W/m2", "min", "max")) print(THIN) for nm, s_ in (("all three vary (the real curve)", q_full), ("eccentricity alone (varpi fixed at today)", q_e_only), ("obliquity alone", q_o_only), ("perihelion angle alone, e at mean", q_p_only), ("eccentricity and perihelion together", q_ep_only)): print(" %-46s %12.3f %12.3f %12.3f" % (nm, np.std(s_), s_.min(), s_.max())) print() print(" full series range : %.3f W/m2 (%.3f to %.3f)" % (q_full.max() - q_full.min(), q_full.min(), q_full.max())) print(" full series mean : %.3f W/m2" % q_full.mean()) print(" today : %.3f W/m2" % q_full[0]) print(" today as a percentile : %.1f%% of the last million years were lower" % (100.0 * np.mean(q_full < q_full[0]))) print(" variance of the full series: %.3f (W/m2)^2" % np.var(q_full)) print(" sum of the three one-at-a-time variances: %.3f (W/m2)^2" % (np.var(q_e_only) + np.var(q_o_only) + np.var(q_p_only))) print(" ratio : %.3f (not 1, because the terms interact)" % ((np.var(q_e_only) + np.var(q_o_only) + np.var(q_p_only)) / np.var(q_full))) # -------------------------------------------------------------------------- # 6. the million-year curve # -------------------------------------------------------------------------- head("SECTION 6. SUMMER INSOLATION AT 65 NORTH, PAST MILLION YEARS") print() print("Insolation on the day of the June solstice (lambda = 90 deg exactly) at") print("65 N, one value per thousand years, from the La2004 elements.") print() print(" %8s %12s %11s %12s %14s" % ("kyr BP", "e", "eps deg", "e sin(w)", "Q65N W/m2")) print(THIN) for i in range(0, NKYR + 1, 25): print(" %8.0f %12.6f %11.4f %+12.5f %14.3f" % (t_kyr[i], ecc[i], obl[i] / DEG, prec[i], q_full[i])) print() imax = int(np.argmax(q_full)) imin = int(np.argmin(q_full)) print(" highest value in the window : %.3f W/m2 at %.0f kyr BP" % (q_full[imax], t_kyr[imax])) print(" lowest value in the window : %.3f W/m2 at %.0f kyr BP" % (q_full[imin], t_kyr[imin])) print(" full swing : %.3f W/m2, which is %.1f%% of the mean" % (q_full.max() - q_full.min(), 100 * (q_full.max() - q_full.min()) / q_full.mean())) print(" at %.0f kyr BP: e=%.5f, eps=%.4f deg, e sin(w)=%+.5f" % (t_kyr[imax], ecc[imax], obl[imax] / DEG, prec[imax])) print(" at %.0f kyr BP: e=%.5f, eps=%.4f deg, e sin(w)=%+.5f" % (t_kyr[imin], ecc[imin], obl[imin] / DEG, prec[imin])) # -------------------------------------------------------------------------- # 7. three definitions of summer # -------------------------------------------------------------------------- head("SECTION 7. THREE DEFINITIONS OF SUMMER, THREE DIFFERENT ANSWERS") print() print("'Summer insolation at 65 N' is not one quantity. Three reasonable") print("definitions are computed below on the same orbital elements.") print() print(" 1. Solstice. Daily-mean insolation on the day lambda = 90 deg.") print(" 2. Caloric summer half year (Milankovitch's own measure). The 182.62") print(" days of the year with the highest daily insolation, averaged.") print(" 3. Integrated summer energy above 275 W/m2 (after Huybers 2006). Total") print(" energy received on every day whose insolation clears the threshold,") print(" reported in GJ/m2.") print() NDAY = 720 day_grid = (np.arange(NDAY) + 0.5) / NDAY * TROPICAL_YEAR dt_day = TROPICAL_YEAR / NDAY THRESH = 275.0 CAL_DAYS = TROPICAL_YEAR / 2.0 NCAL = int(round(CAL_DAYS / dt_day)) sol = q_full.copy() cal = np.zeros(NKYR + 1) ints = np.zeros(NKYR + 1) for i in range(NKYR + 1): lam_d = lambda_from_day(day_grid, ecc[i], vpi[i]) qd = daily_insolation(65.0, lam_d, ecc[i], obl[i], vpi[i]) srt = np.sort(qd)[::-1] cal[i] = np.mean(srt[:NCAL]) ints[i] = np.sum(qd[qd > THRESH]) * dt_day * 86400.0 / 1e9 print(" %-34s %10s %10s %10s %10s %10s" % ("definition", "today", "mean", "min", "max", "std")) print(THIN) for nm, s_ in (("solstice daily mean (W/m2)", sol), ("caloric summer half year (W/m2)", cal), ("integrated summer energy (GJ/m2)", ints)): print(" %-34s %10.3f %10.3f %10.3f %10.3f %10.3f" % (nm, s_[0], s_.mean(), s_.min(), s_.max(), s_.std())) print() print(" Correlation of each definition with the orbital parameters over 1 Myr:") print(" %-34s %12s %12s %12s" % ("definition", "r with e", "r with eps", "r with esin(w)")) print(THIN) for nm, s_ in (("solstice daily mean", sol), ("caloric summer half year", cal), ("integrated summer energy", ints)): print(" %-34s %12.4f %12.4f %12.4f" % (nm, np.corrcoef(s_, ecc)[0, 1], np.corrcoef(s_, obl)[0, 1], np.corrcoef(s_, prec)[0, 1])) # -------------------------------------------------------------------------- # 8. spectra # -------------------------------------------------------------------------- head("SECTION 8. SPECTRAL CONTENT") print() print("Hann-windowed periodogram of each series after removing its mean and") print("linear trend. Sampling 1 kyr, record length %d kyr, so the frequency" % NKYR) print("resolution is 1/%d kyr^-1 and the Nyquist period is 2 kyr." % NKYR) print() print("Bands, fixed from the known periodicities of the orbital solution before") print("the spectra were computed:") print(" eccentricity band 75 to 135 kyr") print(" obliquity band 35 to 55 kyr") print(" precession band 17 to 26 kyr") print() BANDS = (("eccentricity 75-135 kyr", 75.0, 135.0), ("obliquity 35-55 kyr", 35.0, 55.0), ("precession 17-26 kyr", 17.0, 26.0)) spectra = {} series_for_spec = (("eccentricity e", ecc), ("obliquity eps (deg)", obl / DEG), ("precession index e sin(w)", prec), ("Q65N solstice", sol), ("Q65N caloric half year", cal), ("Q65N integrated energy", ints)) band_table = {} for nm, s_ in series_for_spec: f, p = hann_periodogram(s_, 1.0) spectra[nm] = {"f": f.tolist(), "p": p.tolist()} print(" %s" % nm) print(" %-28s %12s %14s %14s" % ("band", "share", "peak period", "band power")) print(THIN) tot_share = 0.0 rowvals = {} for bn, lo_, hi_ in BANDS: sh, pk, bp = band_share(f, p, lo_, hi_) tot_share += sh rowvals[bn] = (sh, pk, bp) print(" %-28s %11.2f%% %12.2f kyr %14.5g" % (bn, 100 * sh, pk, bp)) band_table[nm] = rowvals print(" %-28s %11.2f%%" % ("three bands together", 100 * tot_share)) kmax = int(np.argmax(p[1:])) + 1 print(" %-28s %12.2f kyr" % ("single largest peak at", 1.0 / f[kmax])) print() print(" The strongest local peaks in the 65N solstice spectrum:") f, p = hann_periodogram(sol, 1.0) tot = float(np.sum(p[1:])) loc = [k for k in range(2, p.size - 1) if p[k] > p[k - 1] and p[k] > p[k + 1]] loc.sort(key=lambda k: -p[k]) for k in loc[:8]: print(" period %8.2f kyr power %12.5g share %6.2f%%" % (1.0 / f[k], p[k], 100 * p[k] / tot)) print() print(" How well can a 1000 kyr record place a peak at all? A discrete spectrum") print(" puts peaks only at bin centres. Near period P the bin spacing in period") print(" is P^2 / T, with T the record length, so:") for P in (23.0, 41.0, 100.0): print(" near %6.1f kyr the bins are %.2f kyr apart, and a peak cannot be" % (P, P * P / NKYR)) print(" located better than about half of that, %.2f kyr." % (0.5 * P * P / NKYR)) print() print(" So we ran the same periodogram on the whole 51 Myr file, where the bins") print(" near 41 kyr are only %.4f kyr apart, to pin the line periods down." % (41.0 ** 2 / 51000.0)) def top_lines(arr, nlines, lo_p, hi_p): f_, p_ = hann_periodogram(arr, 1.0) ok = (f_ > 1.0 / hi_p) & (f_ < 1.0 / lo_p) idx = [k for k in range(2, p_.size - 1) if ok[k] and p_[k] > p_[k - 1] and p_[k] > p_[k + 1]] idx.sort(key=lambda k: -p_[k]) tot = float(np.sum(p_[1:])) return [(1.0 / f_[k], 100.0 * p_[k] / tot) for k in idx[:nlines]] print() print(" One warning before the table. The eccentricity lines come from the") print(" secular motion of the planets and hold steady for hundreds of millions") print(" of years, so the whole 51 Myr file can be used on them. The obliquity") print(" and precession lines do not. They depend on Earth's own precession") print(" constant, which tidal friction has been slowing down, so those periods") print(" were shorter in the deep past. Averaging them over 51 Myr smears them") print(" low. For those two we use the most recent 10 Myr instead, and the drift") print(" is printed underneath so a reader can see the size of it.") print() e_L = tab[:, 1] TL_E = tab.shape[0] - 1 REC = tab[:10001] o_R, v_R, e_R = REC[:, 2], REC[:, 3], REC[:, 1] prec_R = e_R * np.sin(v_R + math.pi) TL_R = REC.shape[0] - 1 print(" Line periods, Hann periodogram. The quoted uncertainty is the bin") print(" half-width P^2 / (2T), which is how well a record of length T can place") print(" a line at period P at all.") print() print(" %-24s %7s %17s %10s %22s" % ("quantity", "T kyr", "period kyr", "share", "usually quoted as")) print(THIN) for nm, arr, T_, nl, lo_p, hi_p, tags in ( ("eccentricity e", e_L, TL_E, 4, 60.0, 600.0, ["405 kyr term", "95 kyr term", "124 kyr term", "99 kyr term"]), ("obliquity eps", o_R / DEG, TL_R, 3, 25.0, 80.0, ["41.0 kyr term", "39.7 kyr term", "53.6 kyr term"]), ("precession e sin(w)", prec_R, TL_R, 4, 14.0, 30.0, ["23.7 kyr term", "22.4 kyr term", "19.0 kyr term", "19.1 kyr term"])): for (per, sh), tag in zip(top_lines(arr, nl, lo_p, hi_p), tags): print(" %-24s %7d %10.3f +/- %.3f %9.2f%% %22s" % (nm, T_, per, sh if False else per * per / (2.0 * T_), sh, tag)) print() print(" The right-hand column is what the literature calls each line. Every") print(" recovered period agrees with its published label to better than a") print(" percent. Nothing here was tuned to make them appear. They are") print(" properties of the solar system's orbital dynamics, and they fall out of") print(" a Fourier transform of somebody else's integration.") print() print(" The drift, shown rather than asserted. Peak obliquity and precession") print(" period in 5 Myr windows at four depths in the La2004 file:") print(" %-16s %16s %16s" % ("window (Myr BP)", "obliquity kyr", "precession kyr")) print(THIN) for a_ in (0, 10000, 25000, 45000): win = tab[a_:a_ + 5001] o_w = win[:, 2] / DEG p_w = win[:, 1] * np.sin(win[:, 3] + math.pi) po = top_lines(o_w, 1, 25.0, 80.0)[0][0] pp = top_lines(p_w, 1, 14.0, 30.0)[0][0] print(" %-16s %16.3f %16.3f" % ("%d to %d" % (a_ // 1000, a_ // 1000 + 5), po, pp)) print() print(" That is not an error in our code. It is real, it is in La2004 because") print(" Laskar put it there, and it is why nobody quotes a 41 kyr obliquity") print(" cycle for the Cretaceous.") print() print(" A longer window, to separate the eccentricity lines that 1 Myr cannot:") long_spec = {} for span in (2000, 5000): s2 = tab[:span + 1] e2, o2, v2 = s2[:, 1], s2[:, 2], s2[:, 3] q2 = np.array([float(daily_insolation(65.0, lam_js, e2[i], o2[i], v2[i])) for i in range(s2.shape[0])]) f2, p2 = hann_periodogram(q2, 1.0) f2e, p2e = hann_periodogram(e2, 1.0) long_spec[span] = {"f": f2.tolist(), "p": p2.tolist()} print(" %d kyr window, 65N solstice:" % span) for bn, lo_, hi_ in (("400 kyr band 300-500", 300.0, 500.0),) + BANDS: sh, pk, bp = band_share(f2, p2, lo_, hi_) print(" %-28s %8.2f%% peak %8.2f kyr" % (bn, 100 * sh, pk)) sh, pk, bp = band_share(f2e, p2e, 300.0, 500.0) print(" eccentricity itself, 300-500 kyr band: %.2f%%, peak %.1f kyr" % (100 * sh, pk)) sh, pk, bp = band_share(f2e, p2e, 75.0, 135.0) print(" eccentricity itself, 75-135 kyr band : %.2f%%, peak %.1f kyr" % (100 * sh, pk)) # -------------------------------------------------------------------------- # 9. terminations # -------------------------------------------------------------------------- head("SECTION 9. DOES THE CURVE LINE UP WITH THE TERMINATIONS?") print() print("Termination ages are the marine isotope stage boundary ages of the LR04") print("benthic stack, Lisiecki & Raymo (2005), read out of the published age-") print("model table cached below. We dated nothing and typed nothing. The test") print("asks, for each termination, how far it sits from the nearest maximum of") print("our computed curve.") print() # The nine terminations are the glacial-to-interglacial MIS boundaries. Ages # are read out of the cached LR04 age-model table, not typed in here. LR04_FILE = os.path.join(DATA_DIR, "LR04_MISboundaries.txt") LR04_URL = "https://lorraine-lisiecki.com/LR04_MISboundaries.txt" with open(LR04_FILE, "rb") as fh: lr_raw = fh.read() lr_sha = hashlib.sha256(lr_raw).hexdigest() lr_ages = {} for line in lr_raw.decode("ascii", "replace").splitlines(): bits = line.split() if len(bits) == 2 and "/" in bits[0]: try: lr_ages[bits[0]] = float(bits[1]) except ValueError: pass print(" LR04 age model file : %s" % os.path.basename(LR04_FILE)) print(" source URL : %s" % LR04_URL) print(" retrieved : %s" % RETRIEVED) print(" sha256 : %s" % lr_sha) print(" boundaries parsed : %d" % len(lr_ages)) print() TERMS = tuple((rn, "MIS " + key, lr_ages[key]) for rn, key in (("I", "1/2"), ("II", "5/6"), ("III", "7/8"), ("IV", "9/10"), ("V", "11/12"), ("VI", "13/14"), ("VII", "15/16"), ("VIII", "17/18"), ("IX", "19/20"))) peaks = [(t_kyr[i], q_full[i]) for i in range(1, NKYR) if q_full[i] > q_full[i - 1] and q_full[i] > q_full[i + 1]] peaks.sort() pk_t = np.array([p_[0] for p_ in peaks]) pk_q = np.array([p_[1] for p_ in peaks]) big = pk_q > np.percentile(q_full, 75) print(" local maxima of Q65N in the window : %d" % len(peaks)) print(" mean spacing between maxima : %.2f kyr" % np.mean(np.diff(pk_t))) print(" maxima above the 75th percentile : %d" % int(np.sum(big))) print() print(" %-6s %-10s %9s %12s %10s %12s %10s" % ("term", "MIS bdy", "age kyr", "nearest max", "offset", "Q at max", "Q at term")) print(THIN) offs = [] for nmT, mis, age in TERMS: j = int(np.argmin(np.abs(pk_t - age))) off = pk_t[j] - age offs.append(off) qat = float(np.interp(age, t_kyr, q_full)) print(" %-6s %-10s %9.1f %12.1f %+10.1f %12.2f %10.2f" % (nmT, mis, age, pk_t[j], off, pk_q[j], qat)) offs = np.array(offs) se_off = offs.std(ddof=1) / math.sqrt(offs.size) print() print(" mean offset : %+.2f kyr" % offs.mean()) print(" mean absolute offset : %.2f kyr" % np.abs(offs).mean()) print(" standard deviation of offsets : %.2f kyr" % offs.std(ddof=1)) print(" standard error of the mean : %.2f kyr" % se_off) print(" mean offset in standard errors: %+.2f sigma" % (offs.mean() / se_off)) print() print(" What would chance give? Insolation maxima are spaced %.1f kyr apart on" % np.mean(np.diff(pk_t))) print(" average, so a date thrown at random already lands fairly close to one.") print(" The mean absolute offset above has to be compared against that, not") print(" against zero.") rng2 = np.random.default_rng(SEED + 1) NNULL = 20000 rand_ages = rng2.uniform(10.0, 800.0, NNULL) rand_off = np.array([np.min(np.abs(pk_t - a)) for a in rand_ages]) draws = rng2.choice(rand_off, size=(NNULL, len(TERMS))) hits = int(np.sum(draws.mean(axis=1) <= np.abs(offs).mean())) print() print(" Monte Carlo null, %s random dates drawn uniformly in 10-800 kyr:" % f"{NNULL:,}") print(" mean absolute offset under the null : %.2f kyr" % rand_off.mean()) print(" standard deviation of that null : %.2f kyr" % draws.mean(axis=1).std(ddof=1)) print(" our mean absolute offset : %.2f kyr" % np.abs(offs).mean()) print(" our value in standard errors of the null: %+.2f sigma" % ((np.abs(offs).mean() - draws.mean(axis=1).mean()) / draws.mean(axis=1).std(ddof=1))) print(" null sets of %d dates at least as tight: %d out of %s" % (len(TERMS), hits, f"{NNULL:,}")) if hits == 0: print(" p-value, one sided : < %.5f" % (1.0 / NNULL)) else: print(" p-value, one sided : %.5f" % (hits / NNULL)) print() print(" The same test against the caloric summer half year curve:") peaks_c = [(t_kyr[i], cal[i]) for i in range(1, NKYR) if cal[i] > cal[i - 1] and cal[i] > cal[i + 1]] pk_tc = np.array([p_[0] for p_ in peaks_c]) offs_c = np.array([pk_tc[int(np.argmin(np.abs(pk_tc - age)))] - age for _, _, age in TERMS]) print(" local maxima: %d, mean spacing %.2f kyr" % (len(peaks_c), np.mean(np.diff(pk_tc)))) print(" mean absolute offset: %.2f kyr" % np.abs(offs_c).mean()) print() print(" And the count that matters most. Over the past million years the curve") print(" offers %d maxima. The record offers %d terminations. The difference is" % (len(peaks), len(TERMS))) print(" the whole of the 100 kyr problem: most insolation maxima do nothing.") # -------------------------------------------------------------------------- # 10. what the model leaves out # -------------------------------------------------------------------------- head("SECTION 10. WHAT THIS COMPUTATION LEAVES OUT") print(""" No atmosphere. No clouds. No albedo, so no ice-albedo feedback, which is the single largest amplifier in the real system. No carbon dioxide, no methane, no dust. No ocean, no heat transport, no thermal inertia, so no lag between forcing and response. A real ice sheet takes thousands of years to grow and thousands more to collapse, and none of that is here. No ice sheet at all, so no isostatic rebound and no elevation feedback. The solar constant is held at its present value for a million years. Obliquity is La2004's, referred to the fixed J2000 ecliptic. Termination ages are somebody else's, quoted, not measured. And the 100 kyr power in the ice record has no counterpart in the curve we computed. Section 8 gives the number. We do not explain it, and neither does anybody else with a calculation this simple. """) if os.environ.get("SJC_FIGDATA"): LA = np.linspace(-90, 90, 181) DA = np.linspace(0, TROPICAL_YEAR, 366) lam_g = lambda_from_day(DA, E0, VPI0) G = np.array([daily_insolation(la, lam_g, E0, EPS0, VPI0) for la in LA]) figdata = { "seed": SEED, "S0": S0, "sha256": sha, "retrieved": RETRIEVED, "t_kyr": t_kyr.tolist(), "ecc": ecc.tolist(), "obl_deg": (obl / DEG).tolist(), "prec": prec.tolist(), "q_solstice": sol.tolist(), "q_caloric": cal.tolist(), "q_integrated": ints.tolist(), "q_e_only": q_e_only.tolist(), "q_o_only": q_o_only.tolist(), "q_p_only": q_p_only.tolist(), "conv": conv, "sweeps": sweeps, "spectra": spectra, "long_spec": {str(k): v for k, v in long_spec.items()}, "lats": lats.tolist(), "qbar": qbar.tolist(), "terms": [[a, b, c] for a, b, c in TERMS], "peaks_t": pk_t.tolist(), "peaks_q": pk_q.tolist(), "present": {"e": E0, "eps": EPS0, "varpi": VPI0, "q65": q65}, "grid": {"lat": LA.tolist(), "day": DA.tolist(), "Q": G.tolist()}, } with open(os.environ["SJC_FIGDATA"], "w") as fh: json.dump(figdata, fh) print(RULE) print("Wall clock: %.1f s" % (time.time() - T_START)) print("Seed: %d. Orbital solution: La2004, sha256 %s..." % (SEED, sha[:16])) print(RULE)