#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ daisyworld.py -- Science Journaling Club, Volume 1 Issue 4, Summer 2025. QUESTION -------- Lovelock and Watson built a toy world of black and white daisies to argue that life can stabilise a planet's climate without anything resembling foresight. Does the model really hold temperature steady, and how hard do you have to push before it fails? We replicate the Watson & Lovelock (1983) Daisyworld model, sweep the solar luminosity from 0.6 to 1.6 times present, and measure three things the original paper asserts qualitatively but does not tabulate: (a) the luminosity interval over which regulation actually holds, (b) the residual temperature drift inside that interval, (c) how both change when the growth window is narrowed, the death rate is raised, the heat transport is altered, or the seeding rule is changed. THE COMPUTATION IS THE EXPERIMENT. The club has no laboratory, no telescope and no planet. Every number below is printed output from this file. Nothing here was observed. Daisyworld is not a description of Earth and we do not treat it as one; it is an existence proof about a class of feedback, and the point of re-running it is to find out how wide that class is. THE MODEL --------- State variables are three areas that must sum to one: a_b fraction of the planet covered by black daisies (albedo A_b = 0.25) a_w fraction covered by white daisies (albedo A_w = 0.75) a_g fraction of bare ground (albedo A_g = 0.50) Planetary albedo is the area-weighted mean A = a_b A_b + a_w A_w + a_g A_g. The planet is a zero-dimensional black body in radiative equilibrium with a star of luminosity L: sigma Te^4 = S L (1 - A), S = 917 W m^-2, sigma = 5.67032e-8. Heat moves between patches imperfectly. A patch of albedo A_i sits at a local temperature given by the Watson & Lovelock parameterisation T_i^4 = q (A - A_i) + Te^4, q = 2.06e9 K^4. q = 0 would be a perfectly conducting planet with no local temperature contrast and therefore no daisy feedback at all; q = S L / sigma would be a perfectly insulating one. 2.06e9 is the published intermediate value. Daisies grow into bare ground and die at a constant rate: d a_b / dt = a_b ( a_g beta(T_b) - gamma ) d a_w / dt = a_w ( a_g beta(T_w) - gamma ) d a_g / dt = - ( d a_b / dt + d a_w / dt ) with a parabolic growth curve, zero outside a window of half-width h about the optimum T_opt: beta(T) = max(0, 1 - (T - T_opt)^2 / h^2), T_opt = 295.5 K, h = 17.5 K. gamma = 0.3 per unit time. Time carries no physical units: only the ratio of growth to death matters. The third equation is redundant with the first two, and that is deliberate. We integrate all three and then check that their sum stays at one, which is the cheapest available test that the integrator is not silently manufacturing or destroying planet surface. ASSUMPTIONS, STATED PLAINLY --------------------------- 1. Zero spatial dimensions. There is no equator, no pole, no latitude, no ocean, no season, no day. "Local temperature" is an accounting device, not a place. 2. No atmosphere, no greenhouse effect, no clouds, no water vapour. Albedo is the only radiative property of the surface and it is set entirely by what is growing on it. 3. Two species, fixed albedos, one shared optimum, identical growth curves. No mutation, no evolution of albedo, no herbivores, no competition for anything except bare ground. 4. Constant heat-transport coefficient q, independent of the daisy cover whose heat it is supposed to be moving. 5. Seeds are always present. Before each luminosity step we floor each daisy area at SEED_FLOOR = 0.01 so a dead population can recolonise. This is the standard convention and it is load-bearing: section [10] runs the same sweep without it and gets a materially different answer. 6. Quasi-static luminosity. Each L is integrated to steady state before L moves, and the state is carried forward to the next L. A real star brightens on a timescale set by nuclear physics and a real biosphere responds on a timescale set by biology; the model assumes the second is infinitely faster than the first. LIMITATIONS, i.e. WHAT THIS CANNOT SHOW --------------------------------------- Daisyworld was built by people who already believed its conclusion, in order to demonstrate that the conclusion is possible. It succeeds at that and it cannot do more. A regulating result here is evidence that albedo-mediated feedback CAN stabilise a planetary temperature in a system with these assumptions. It is not evidence that Earth's biosphere does so, that any planet does so, or that regulation is a generic outcome of putting life on a rock. Nothing in this file constrains the real climate system in any way. REPRODUCING ----------- python daisyworld.py > daisyworld-output.txt Master seed 20250621, hard coded below. numpy PCG64 via SeedSequence.spawn. The only stochastic part of the study is the Monte Carlo in section [9]; every other number is deterministic. """ import sys import time import math import numpy as np # ---------------------------------------------------------------------------- # Constants and configuration # ---------------------------------------------------------------------------- MASTER_SEED = 20250621 SIGMA = 5.67032e-8 # Stefan-Boltzmann constant, W m^-2 K^-4 S_FLUX = 917.0 # Daisyworld solar constant at L = 1, W m^-2 A_BLACK = 0.25 # albedo of black daisies A_WHITE = 0.75 # albedo of white daisies A_GROUND = 0.50 # albedo of bare ground Q_HEAT = 2.06e9 # heat transport coefficient, K^4 GAMMA = 0.3 # death rate, per unit time T_OPT = 295.5 # growth optimum, K HALFWIDTH = 17.5 # growth window half-width, K SEED_FLOOR = 0.01 # daisy area restored before each luminosity step ALIVE = 0.01 # total daisy area above which we call the biota alive DT = 0.05 # RK4 step for the scalar runs T_MAX = 200.0 # integration ceiling per luminosity step STEADY_TOL = 1e-9 # |d a/dt| below this counts as steady L_LO, L_HI = 0.6, 1.6 L_N_MAIN = 501 # main sweep resolution, dL = 0.002 L_N_SCAN = 201 # parameter scans, dL = 0.005 L_N_MC = 101 # Monte Carlo, dL = 0.010 KELVIN = 273.15 # ---------------------------------------------------------------------------- # Scalar engine. Deliberately free of numpy: on a two-variable ODE the numpy # call overhead is about fifty times the cost of the arithmetic, and the whole # study is a few million RK4 steps. # ---------------------------------------------------------------------------- def sweep(Lgrid, Ab=A_BLACK, Aw=A_WHITE, Ag=A_GROUND, q=Q_HEAT, gamma=GAMMA, Topt=T_OPT, h=HALFWIDTH, S=S_FLUX, seed_floor=SEED_FLOOR, ab0=SEED_FLOOR, aw0=SEED_FLOOR, dt=DT, tmax=T_MAX, tol=STEADY_TOL): """Quasi-static luminosity sweep with continuation. Each L starts from the previous steady state, floored at seed_floor when that is not None. Returns a dict of arrays plus diagnostics.""" nmax = int(round(tmax / dt)) k = 1.0 / (h * h) qAb = q * Ab qAw = q * Aw ab, aw = float(ab0), float(aw0) ag = 1.0 - ab - aw B, W, G, TE, TB, TW, AL = [], [], [], [], [], [], [] total_steps = 0 worst_area = 0.0 unconverged = 0 for L in Lgrid: if seed_floor is not None: if ab < seed_floor: ab = seed_floor if aw < seed_floor: aw = seed_floor ag = 1.0 - ab - aw C = S * float(L) / SIGMA # so that Te^4 = C (1 - A) def dv(a, b, c): A = a * Ab + b * Aw + c * Ag te4 = C * (1.0 - A) qA = q * A t4 = qA - qAb + te4 Tb = t4 ** 0.25 if t4 > 0.0 else 0.0 u4 = qA - qAw + te4 Tw = u4 ** 0.25 if u4 > 0.0 else 0.0 x = 1.0 - k * (Tb - Topt) ** 2 y = 1.0 - k * (Tw - Topt) ** 2 d1 = a * (c * (x if x > 0.0 else 0.0) - gamma) d2 = b * (c * (y if y > 0.0 else 0.0) - gamma) return d1, d2, -(d1 + d2) done = False for _ in range(nmax): a1, b1, c1 = dv(ab, aw, ag) a2, b2, c2 = dv(ab + 0.5 * dt * a1, aw + 0.5 * dt * b1, ag + 0.5 * dt * c1) a3, b3, c3 = dv(ab + 0.5 * dt * a2, aw + 0.5 * dt * b2, ag + 0.5 * dt * c2) a4, b4, c4 = dv(ab + dt * a3, aw + dt * b3, ag + dt * c3) da = (a1 + 2.0 * a2 + 2.0 * a3 + a4) / 6.0 db = (b1 + 2.0 * b2 + 2.0 * b3 + b4) / 6.0 dc = (c1 + 2.0 * c2 + 2.0 * c3 + c4) / 6.0 ab += dt * da aw += dt * db ag += dt * dc total_steps += 1 e = abs(ab + aw + ag - 1.0) if e > worst_area: worst_area = e if abs(da) < tol and abs(db) < tol: done = True break if not done: unconverged += 1 A = ab * Ab + aw * Aw + ag * Ag te4 = C * (1.0 - A) qA = q * A t4 = qA - qAb + te4 u4 = qA - qAw + te4 B.append(ab) W.append(aw) G.append(ag) AL.append(A) TE.append(te4 ** 0.25 if te4 > 0.0 else 0.0) TB.append(t4 ** 0.25 if t4 > 0.0 else 0.0) TW.append(u4 ** 0.25 if u4 > 0.0 else 0.0) return { "L": np.asarray(Lgrid, dtype=float), "ab": np.array(B), "aw": np.array(W), "ag": np.array(G), "A": np.array(AL), "Te": np.array(TE), "Tb": np.array(TB), "Tw": np.array(TW), "steps": total_steps, "worst_area": worst_area, "unconverged": unconverged, } def relax(ab, aw, ag, L, dt=DT, tmax=T_MAX, tol=STEADY_TOL, Ab=A_BLACK, Aw=A_WHITE, Ag=A_GROUND, q=Q_HEAT, gamma=GAMMA, Topt=T_OPT, h=HALFWIDTH, S=S_FLUX): """Integrate a single state at a single luminosity to steady state.""" nmax = int(round(tmax / dt)) k = 1.0 / (h * h) qAb, qAw = q * Ab, q * Aw C = S * float(L) / SIGMA def dv(a, b, c): A = a * Ab + b * Aw + c * Ag te4 = C * (1.0 - A) qA = q * A t4 = qA - qAb + te4 Tb = t4 ** 0.25 if t4 > 0.0 else 0.0 u4 = qA - qAw + te4 Tw = u4 ** 0.25 if u4 > 0.0 else 0.0 x = 1.0 - k * (Tb - Topt) ** 2 y = 1.0 - k * (Tw - Topt) ** 2 d1 = a * (c * (x if x > 0.0 else 0.0) - gamma) d2 = b * (c * (y if y > 0.0 else 0.0) - gamma) return d1, d2, -(d1 + d2) steps = 0 for _ in range(nmax): a1, b1, c1 = dv(ab, aw, ag) a2, b2, c2 = dv(ab + 0.5 * dt * a1, aw + 0.5 * dt * b1, ag + 0.5 * dt * c1) a3, b3, c3 = dv(ab + 0.5 * dt * a2, aw + 0.5 * dt * b2, ag + 0.5 * dt * c2) a4, b4, c4 = dv(ab + dt * a3, aw + dt * b3, ag + dt * c3) da = (a1 + 2 * a2 + 2 * a3 + a4) / 6.0 db = (b1 + 2 * b2 + 2 * b3 + b4) / 6.0 dc = (c1 + 2 * c2 + 2 * c3 + c4) / 6.0 ab += dt * da aw += dt * db ag += dt * dc steps += 1 if abs(da) < tol and abs(db) < tol: break A = ab * Ab + aw * Aw + ag * Ag Te = (C * (1.0 - A)) ** 0.25 return ab, aw, ag, Te, steps # ---------------------------------------------------------------------------- # Vector engine, used only for the Monte Carlo: many planets, one code path. # ---------------------------------------------------------------------------- def sweep_vec(Lgrid, Ab, Aw, q, gamma, Topt, h, Ag=A_GROUND, S=S_FLUX, seed_floor=SEED_FLOOR, dt=0.10, tmax=150.0, tol=STEADY_TOL): """Identical dynamics, evaluated on numpy arrays so that N independent parameter draws are swept together.""" nmax = int(round(tmax / dt)) k = 1.0 / (h * h) qAb, qAw = q * Ab, q * Aw shape = np.broadcast(Ab, Aw, q, gamma, Topt, h).shape ab = np.full(shape, seed_floor) aw = np.full(shape, seed_floor) ag = 1.0 - ab - aw B, W, TE = [], [], [] total = 0 for L in Lgrid: np.maximum(ab, seed_floor, out=ab) np.maximum(aw, seed_floor, out=aw) ag = 1.0 - ab - aw C = S * float(L) / SIGMA def dv(a, b, c): A = a * Ab + b * Aw + c * Ag te4 = C * (1.0 - A) qA = q * A Tb = np.maximum(0.0, qA - qAb + te4) ** 0.25 Tw = np.maximum(0.0, qA - qAw + te4) ** 0.25 gb = np.maximum(0.0, 1.0 - k * (Tb - Topt) ** 2) gw = np.maximum(0.0, 1.0 - k * (Tw - Topt) ** 2) d1 = a * (c * gb - gamma) d2 = b * (c * gw - gamma) return d1, d2, -(d1 + d2) for _ in range(nmax): a1, b1, c1 = dv(ab, aw, ag) a2, b2, c2 = dv(ab + 0.5 * dt * a1, aw + 0.5 * dt * b1, ag + 0.5 * dt * c1) a3, b3, c3 = dv(ab + 0.5 * dt * a2, aw + 0.5 * dt * b2, ag + 0.5 * dt * c2) a4, b4, c4 = dv(ab + dt * a3, aw + dt * b3, ag + dt * c3) da = (a1 + 2 * a2 + 2 * a3 + a4) / 6.0 db = (b1 + 2 * b2 + 2 * b3 + b4) / 6.0 dc = (c1 + 2 * c2 + 2 * c3 + c4) / 6.0 ab = ab + dt * da aw = aw + dt * db ag = ag + dt * dc total += 1 if np.max(np.abs(da)) < tol and np.max(np.abs(db)) < tol: break A = ab * Ab + aw * Aw + ag * Ag B.append(ab.copy()) W.append(aw.copy()) TE.append((C * (1.0 - A)) ** 0.25) return np.array(B), np.array(W), np.array(TE), total # ---------------------------------------------------------------------------- # Small helpers # ---------------------------------------------------------------------------- def bare_planet(Lgrid, S=S_FLUX, A=A_GROUND): """Closed-form Stefan-Boltzmann temperature of a daisy-free planet.""" return (S * np.asarray(Lgrid, dtype=float) * (1.0 - A) / SIGMA) ** 0.25 def longest_run(L, mask): """Longest contiguous run of True. Returns (L_first, L_last, indices).""" m = np.asarray(mask, dtype=bool) if not m.any(): return None idx = np.where(m)[0] runs = np.split(idx, np.where(np.diff(idx) > 1)[0] + 1) run = max(runs, key=len) return float(L[run[0]]), float(L[run[-1]]), run def span(a): a = np.asarray(a, dtype=float) return float(a.max() - a.min()) def erode(mask, k=2): """Drop the k grid points nearest either end of every True run. Steady state is reached exactly in the interior of a run and only approached at a boundary, where the population is in the act of appearing or dying and the integration window is finite. This separates the two cases instead of quietly averaging them.""" m = np.asarray(mask, dtype=bool).copy() for _ in range(k): left = np.r_[False, m[:-1]] right = np.r_[m[1:], False] m = m & left & right return m def rule(ch="-", n=78): return ch * n def header(title): print() print(rule("=")) print(title) print(rule("=")) # ============================================================================ def main(): t_start = time.time() ss = np.random.SeedSequence(MASTER_SEED) streams = ss.spawn(4) print(rule("=")) print("DAISYWORLD: A PLANET THAT REGULATES ITS OWN TEMPERATURE") print("Science Journaling Club, Volume 1 Issue 4, Summer 2025 -- replication") print(rule("=")) print("This is a simulation. No planet, no daisy and no measurement of any") print("physical object is involved anywhere in this file. Every number below") print("is arithmetic performed by this program.") print() print("python : %s" % sys.version.split()[0]) print("numpy : %s" % np.__version__) print("master seed : %d" % MASTER_SEED) # ---------------------------------------------------------------- [0] header("[0] CONFIGURATION") print(" Stefan-Boltzmann sigma = %.6e W m^-2 K^-4" % SIGMA) print(" solar constant S at L = 1 = %.1f W m^-2" % S_FLUX) print(" albedo black/ground/white = %.2f / %.2f / %.2f" % (A_BLACK, A_GROUND, A_WHITE)) print(" heat transport q = %.3e K^4" % Q_HEAT) print(" death rate gamma = %.2f per unit time" % GAMMA) print(" growth optimum T_opt = %.2f K (%.2f C)" % (T_OPT, T_OPT - KELVIN)) print(" growth half-width h = %.2f K -> window %.2f to %.2f K" % (HALFWIDTH, T_OPT - HALFWIDTH, T_OPT + HALFWIDTH)) print(" seed floor per species = %.3f" % SEED_FLOOR) print(" RK4 step dt = %.3f, ceiling t = %.0f" % (DT, T_MAX)) print(" steady-state tolerance = %.0e" % STEADY_TOL) print(" main grid %.2f to %.2f in %d points (dL = %.4f)" % (L_LO, L_HI, L_N_MAIN, (L_HI - L_LO) / (L_N_MAIN - 1))) print(" scan grid %.2f to %.2f in %d points (dL = %.4f)" % (L_LO, L_HI, L_N_SCAN, (L_HI - L_LO) / (L_N_SCAN - 1))) print(" MC grid %.2f to %.2f in %d points (dL = %.4f)" % (L_LO, L_HI, L_N_MC, (L_HI - L_LO) / (L_N_MC - 1))) Lmain = np.linspace(L_LO, L_HI, L_N_MAIN) Lscan = np.linspace(L_LO, L_HI, L_N_SCAN) analytic = bare_planet(Lmain) print() print(" Where q = 2.06e9 comes from. Watson and Lovelock write the local") print(" temperature in a linearised form, T_i = T_e + q'(A - A_i), with") print(" q' = 20 K per unit albedo. Differentiating T^4 gives the quartic") print(" version used here, with q = 4 T_opt^3 q'. So:") q_prime = 20.0 q_derived = 4.0 * T_OPT ** 3 * q_prime print(" club 4 * %.2f^3 * %.1f = %.6e K^4" % (T_OPT, q_prime, q_derived)) print(" published = %.6e K^4" % Q_HEAT) print(" difference = %.3e K^4 (%.3f%%)" % (q_derived - Q_HEAT, 100.0 * (q_derived - Q_HEAT) / Q_HEAT)) print(" The published constant is that arithmetic rounded to three figures.") print(" We keep 2.06e9 so that the replication uses the published number.") # ---------------------------------------------------------------- [1] header("[1] VALIDATION: AREA CONSERVATION") print("The three areas are integrated as three independent state variables.") print("Nothing in the integrator forces them to sum to one. If RK4 and the") print("right-hand side are both correct, the sum is a conserved quantity and") print("only floating-point rounding can move it.") print() t0 = time.time() run = sweep(Lmain) print(" main sweep: %d luminosities, %d RK4 steps, %.1f s" % (L_N_MAIN, run["steps"], time.time() - t0)) print(" luminosity steps that hit the integration ceiling : %d of %d" % (run["unconverged"], L_N_MAIN)) worst = run["worst_area"] sums = run["ab"] + run["aw"] + run["ag"] eps = float(np.finfo(float).eps) print() print(" worst |a_b + a_w + a_g - 1| at ANY of the %d steps : %.3e" % (run["steps"], worst)) print(" worst at the %d recorded steady states : %.3e" % (L_N_MAIN, float(np.max(np.abs(sums - 1.0))))) print(" float64 machine epsilon : %.3e" % eps) print() print(" club value %.3e required value 0 difference %.3e" % (worst, worst)) print(" that is %.1f machine epsilons, i.e. rounding and nothing else" % (worst / eps)) print(" VERDICT: %s" % ("PASS, area conserved to machine precision" if worst < 1e-13 else "FAIL")) # step-size check print() print(" Step-size check. The same sweep on a coarser grid at three values") print(" of dt. If the answer moves with dt, the integrator is the result.") print(" dt regulated L range drift(C) mean T(C) steps") for dtv in (0.20, 0.10, 0.05, 0.025): r = sweep(Lscan, dt=dtv) al = (r["ab"] + r["aw"]) > ALIVE lo, hi, rr = longest_run(Lscan, al) tt = r["Te"][rr] - KELVIN print(" %6.3f %.3f to %.3f %8.5f %9.5f %7d" % (dtv, lo, hi, span(tt), float(tt.mean()), r["steps"])) Ttot = run["ab"] + run["aw"] # ---------------------------------------------------------------- [2] header("[2] VALIDATION: BARE-PLANET CONTROL vs ANALYTIC BLACK BODY") print("Run the same integrator with no daisies seeded at all. The answer") print("must be the Stefan-Boltzmann temperature of a grey sphere of albedo") print("0.5, which we can write down in closed form:") print() print(" Te = ( S L (1 - A) / sigma )^(1/4)") print() dead = sweep(Lmain, seed_floor=None, ab0=0.0, aw0=0.0) d = np.abs(dead["Te"] - analytic) print(" L club Te (K) analytic (K) difference (K)") for L in (0.6, 0.8, 1.0, 1.2, 1.4, 1.6): i = int(np.argmin(np.abs(Lmain - L))) print(" %5.2f %10.5f %10.5f %10.3e" % (Lmain[i], dead["Te"][i], analytic[i], dead["Te"][i] - analytic[i])) print() print(" max |club - analytic| over all %d luminosities : %.3e K" % (L_N_MAIN, float(d.max()))) print(" VERDICT: %s" % ("PASS" if d.max() < 1e-9 else "FAIL")) print() print(" bare planet warming across the full sweep:") print(" L = 0.60 -> %.2f C ; L = 1.60 -> %.2f C ; span %.2f C" % (analytic[0] - KELVIN, analytic[-1] - KELVIN, analytic[-1] - analytic[0])) lin = np.polyfit(Lmain, analytic - KELVIN, 1) resid = analytic - KELVIN - np.polyval(lin, Lmain) print(" best straight line: %.3f C per unit L, worst residual %.3f C" % (lin[0], float(np.max(np.abs(resid))))) print(" the bare curve is L^(1/4), so 'linear' is an approximation good to") print(" %.2f C over this range. We say near-linear and mean exactly that." % float(np.max(np.abs(resid)))) # ---------------------------------------------------------------- [3] header("[3] THE MAIN RESULT: DOES IT REGULATE?") alive = Ttot > ALIVE Llo, Lhi, reg = longest_run(Lmain, alive) Tin = run["Te"][reg] - KELVIN Tbare_in = analytic[reg] - KELVIN drift0 = span(Tin) bare0 = span(Tbare_in) supp0 = bare0 / drift0 print(" regulated interval (daisy cover > %.2f, longest contiguous run)" % ALIVE) print(" L from %.3f to %.3f, width %.3f" % (Llo, Lhi, Lhi - Llo)) print(" that is %.1f%% of the swept range" % (100 * (Lhi - Llo) / (L_HI - L_LO))) print() print(" planetary temperature INSIDE the regulated interval") print(" min %.3f C at L = %.3f" % (Tin.min(), Lmain[reg][int(np.argmin(Tin))])) print(" max %.3f C at L = %.3f" % (Tin.max(), Lmain[reg][int(np.argmax(Tin))])) print(" mean %.3f C" % Tin.mean()) print(" residual drift (max - min) = %.3f C" % drift0) print(" daisy growth optimum for comparison = %.2f C" % (T_OPT - KELVIN)) print(" mean |regulated T - optimum| = %.3f C" % float(np.mean(np.abs(Tin - (T_OPT - KELVIN))))) print(" standard deviation of regulated T = %.3f C" % float(Tin.std(ddof=1))) print() print(" bare planet OVER THE SAME luminosity interval") print(" min %.3f C, max %.3f C, span %.3f C" % (Tbare_in.min(), Tbare_in.max(), bare0)) print() print(" club value vs control: regulated drift %.3f C against bare span %.3f C" % (drift0, bare0)) print(" suppression factor = %.2f (bare span divided by regulated drift)" % supp0) print() print(" Watson & Lovelock 1983 state this qualitatively: temperature is held") print(" close to the optimum for growth over a wide range of luminosity, and") print(" their Figure 3 shows the two-species curve flat from roughly L = 0.7") print(" to roughly L = 1.6. Our quantitative version: over L in [%.3f, %.3f]" % (Llo, Lhi)) print(" the planet stays within %.2f C of %.2f C, while an identical bare" % (max(abs(Tin.max() - (T_OPT - KELVIN)), abs(Tin.min() - (T_OPT - KELVIN))), T_OPT - KELVIN)) print(" planet moves through %.1f C." % bare0) print() print(" The drift above is dominated by the two edges, where regulation is") print(" in the act of failing. Two narrower statements of the same result:") tight = np.abs(run["Te"] - T_OPT) < 2.0 if tight.any(): tl, th, trr = longest_run(Lmain, tight) print(" |T - T_opt| < 2.0 C over L = %.3f to %.3f, width %.3f" % (tl, th, th - tl)) print(" drift inside that band = %.3f C" % span(run["Te"][trr])) print(" bare planet span over the same L = %.3f C" % span(analytic[trr])) print(" suppression factor there = %.2f" % (span(analytic[trr]) / span(run["Te"][trr]))) inner = reg[int(0.05 * len(reg)):int(0.95 * len(reg))] print(" middle 90%% of the regulated interval (L = %.3f to %.3f):" % (Lmain[inner][0], Lmain[inner][-1])) print(" drift %.3f C against bare span %.3f C, suppression %.2f" % (span(run["Te"][inner]), span(analytic[inner]), span(analytic[inner]) / span(run["Te"][inner]))) print() print(" Both edges are jumps, not ramps. One grid step is dL = %.4f." % (Lmain[1] - Lmain[0])) dT = np.diff(run["Te"]) i_cold = reg[0] - 1 i_hot = reg[-1] print(" edge L before L after T before(C) T after(C) jump(C) cover before -> after") print(" cold %6.3f %6.3f %10.3f %10.3f %+8.3f %.5f -> %.5f" % (Lmain[i_cold], Lmain[i_cold + 1], run["Te"][i_cold] - KELVIN, run["Te"][i_cold + 1] - KELVIN, dT[i_cold], Ttot[i_cold], Ttot[i_cold + 1])) print(" hot %6.3f %6.3f %10.3f %10.3f %+8.3f %.5f -> %.5f" % (Lmain[i_hot], Lmain[i_hot + 1], run["Te"][i_hot] - KELVIN, run["Te"][i_hot + 1] - KELVIN, dT[i_hot], Ttot[i_hot], Ttot[i_hot + 1])) jmp = int(np.argmax(np.abs(dT))) print(" largest single-step change anywhere: %+.3f C between L = %.3f and %.3f" % (dT[jmp], Lmain[jmp], Lmain[jmp + 1])) print(" over that same step the bare planet changes by %+.4f C" % (analytic[jmp + 1] - analytic[jmp])) print(" The hot edge is the catastrophic collapse of life under gradual") print(" forcing described by Ackland, Clark and Lenton; the cold edge is the") print(" same discontinuity run the other way, colonisation rather than death.") print() print(" Phases of the sweep:") lb = run["ab"] > ALIVE lw = run["aw"] > ALIVE both = lb & lw only_b = lb & ~lw only_w = lw & ~lb if only_b.any(): print(" black only : L %.3f to %.3f" % (Lmain[only_b][0], Lmain[only_b][-1])) if both.any(): print(" coexistence : L %.3f to %.3f" % (Lmain[both][0], Lmain[both][-1])) if only_w.any(): print(" white only : L %.3f to %.3f" % (Lmain[only_w][0], Lmain[only_w][-1])) print(" peak total cover %.4f at L = %.3f" % (float(Ttot.max()), Lmain[int(np.argmax(Ttot))])) print(" bare ground never falls below %.4f, because steady state requires" % float(run["ag"][reg].min())) print(" a_g = gamma / beta >= gamma = %.2f." % GAMMA) print() print(" FIGURE DATA 1: L, planetary T (C), bare T (C), black, white, bare ground") for i in list(range(0, L_N_MAIN, 10)) + [L_N_MAIN - 1]: print(" %.3f %9.4f %9.4f %8.5f %8.5f %8.5f" % (Lmain[i], run["Te"][i] - KELVIN, analytic[i] - KELVIN, run["ab"][i], run["aw"][i], run["ag"][i])) print() print(" FIGURE DATA 2: L, planetary albedo A, local T_black (C), local T_white (C)") for i in list(range(0, L_N_MAIN, 20)) + [L_N_MAIN - 1]: print(" %.3f %8.5f %9.4f %9.4f" % (Lmain[i], run["A"][i], run["Tb"][i] - KELVIN, run["Tw"][i] - KELVIN)) # ---------------------------------------------------------------- [4] header("[4] VALIDATION: STEADY-STATE ALGEBRA") print("At a steady state with a living population, growth must exactly") print("balance death: a_g * beta(T_i) = gamma for every surviving type i.") print("That is an algebraic identity we never imposed on the integrator.") print() kk = 1.0 / HALFWIDTH ** 2 bb = np.maximum(0.0, 1.0 - kk * (run["Tb"] - T_OPT) ** 2) bw = np.maximum(0.0, 1.0 - kk * (run["Tw"] - T_OPT) ** 2) resb = np.abs(run["ag"] * bb - GAMMA) resw = np.abs(run["ag"] * bw - GAMMA) rb, rw = resb[lb], resw[lw] ib, iw = erode(lb, 2), erode(lw, 2) print(" black daisies alive at %d of %d luminosities" % (int(lb.sum()), L_N_MAIN)) print(" club |a_g beta(T_b) - gamma| max = %.3e required 0" % float(rb.max())) print(" median = %.3e" % float(np.median(rb))) print(" max away from a live/dead edge = %.3e (%d points)" % (float(resb[ib].max()), int(ib.sum()))) print(" white daisies alive at %d of %d luminosities" % (int(lw.sum()), L_N_MAIN)) print(" club |a_g beta(T_w) - gamma| max = %.3e required 0" % float(rw.max())) print(" median = %.3e" % float(np.median(rw))) print(" max away from a live/dead edge = %.3e (%d points)" % (float(resw[iw].max()), int(iw.sum()))) interior_worst = max(float(resb[ib].max()), float(resw[iw].max())) bad = ((lb & (resb > 1e-6)) | (lw & (resw > 1e-6))) print() print(" %d of the %d live luminosities carry a residual above 1e-6." % (int(bad.sum()), int((lb | lw).sum()))) ib_w = int(np.argmax(np.where(lw, resw, -1.0))) print(" The worst is at L = %.3f, where white cover is %.6f and falling." % (Lmain[ib_w], run["aw"][ib_w])) print() print(" Our reading is that these are unfinished transients, not errors: a") print(" population on its way to zero approaches steady state slowly, and our") print(" integration window is t = %.0f. That is a testable claim. If it is" % T_MAX) print(" right, the residuals shrink when the window is lengthened and the") print(" physical answer does not move. So we lengthened it fourfold.") print() t0 = time.time() longrun = sweep(Lmain, tmax=4 * T_MAX) lb2 = longrun["ab"] > ALIVE lw2 = longrun["aw"] > ALIVE bb2_ = np.maximum(0.0, 1.0 - kk * (longrun["Tb"] - T_OPT) ** 2) bw2_ = np.maximum(0.0, 1.0 - kk * (longrun["Tw"] - T_OPT) ** 2) resb2 = np.abs(longrun["ag"] * bb2_ - GAMMA) resw2 = np.abs(longrun["ag"] * bw2_ - GAMMA) worst2 = max(float(resb2[lb2].max()), float(resw2[lw2].max())) bad2 = ((lb2 & (resb2 > 1e-6)) | (lw2 & (resw2 > 1e-6))) al2 = (longrun["ab"] + longrun["aw"]) > ALIVE lo2_, hi2_, rr2_ = longest_run(Lmain, al2) print(" t = %-5.0f t = %-5.0f change" % (T_MAX, 4 * T_MAX)) print(" worst |a_g beta - gamma| %10.3e %10.3e %6.2fx" % (float(max(rb.max(), rw.max())), worst2, float(max(rb.max(), rw.max())) / worst2)) print(" luminosities above 1e-6 %10d %10d" % (int(bad.sum()), int(bad2.sum()))) print(" regulated interval low %10.3f %10.3f" % (Llo, lo2_)) print(" regulated interval high %10.3f %10.3f" % (Lhi, hi2_)) print(" residual drift (C) %10.4f %10.4f" % (drift0, span(longrun["Te"][rr2_]))) print(" mean regulated T (C) %10.4f %10.4f" % (Tin.mean(), float(np.mean(longrun["Te"][rr2_] - KELVIN)))) print(" seconds %10.1f" % (time.time() - t0)) print() print(" The residual falls and the physics does not move. We report the") print(" t = %.0f numbers throughout, because that is the setting the whole" % T_MAX) print(" study was run at, and note that the slow points are all dying") print(" populations near an edge.") print(" VERDICT: PASS at %d of %d live luminosities to 1e-6; the rest are" % (int((lb | lw).sum() - bad.sum()), int((lb | lw).sum()))) print(" transients that shrink by %.1fx when the window is quadrupled." % (float(max(rb.max(), rw.max())) / worst2)) # ---------------------------------------------------------------- [5] header("[5] VALIDATION: THE COEXISTENCE SYMMETRY, AN ANALYTIC PREDICTION") print("Where both colours coexist, both satisfy a_g beta(T) = gamma with the") print("same a_g and the same gamma, so beta(T_b) = beta(T_w). The growth") print("curve is a parabola about T_opt, and the two local temperatures are") print("not equal, so they must sit symmetrically either side of the optimum:") print() print(" T_black + T_white = 2 T_opt = %.2f K, exactly." % (2 * T_OPT)) print() print("Nobody told the integrator this. It is a consequence of the model and") print("it is the sharpest single check available anywhere in the sweep.") print() sym = run["Tb"][both] + run["Tw"][both] ico = erode(both, 2) symi = run["Tb"][ico] + run["Tw"][ico] print(" coexistence at %d of %d luminosities, L = %.3f to %.3f" % (int(both.sum()), L_N_MAIN, Lmain[both][0], Lmain[both][-1])) print(" club mean (T_b + T_w) = %.9f K" % float(sym.mean())) print(" analytic 2 T_opt = %.9f K" % (2 * T_OPT)) print(" difference of means = %.3e K" % float(sym.mean() - 2 * T_OPT)) print(" worst single luminosity = %.3e K (at a boundary)" % float(np.max(np.abs(sym - 2 * T_OPT)))) print(" worst away from an edge = %.3e K (%d points)" % (float(np.max(np.abs(symi - 2 * T_OPT))), int(ico.sum()))) print(" median over the band = %.3e K" % float(np.median(np.abs(sym - 2 * T_OPT)))) both2 = lb2 & lw2 sym2 = longrun["Tb"][both2] + longrun["Tw"][both2] print(" same quantity at t = %.0f = %.3e K worst, %.3e K median" % (4 * T_MAX, float(np.max(np.abs(sym2 - 2 * T_OPT))), float(np.median(np.abs(sym2 - 2 * T_OPT))))) print(" VERDICT: PASS; the median point satisfies the identity to %.0e K and" % float(np.median(np.abs(sym - 2 * T_OPT)))) print(" the worst points are the same edge transients as section [4].") print() print(" The same algebra goes further, and predicts the numbers rather than") print(" a relation between them. Write T_b = T_opt + d and T_w = T_opt - d.") print(" Subtracting the two local-temperature equations kills T_e entirely:") print() print(" (T_opt + d)^4 - (T_opt - d)^4 = q (A_w - A_b)") print(" 8 T_opt^3 d + 8 T_opt d^3 = q (A_w - A_b)") print() print(" which is one cubic in d with no reference to luminosity, albedo, or") print(" the size of either population. Solve it and everything else follows:") print(" beta = 1 - (d/h)^2 for both colours, and a_g = gamma / beta.") print() co = np.roots([8.0 * T_OPT, 0.0, 8.0 * T_OPT ** 3, -Q_HEAT * (A_WHITE - A_BLACK)]) d_pred = float(min(r.real for r in co if abs(r.imag) < 1e-8 and r.real > 0)) beta_pred = 1.0 - (d_pred / HALFWIDTH) ** 2 ag_pred = GAMMA / beta_pred print(" predicted d = %.9f K" % d_pred) print(" predicted beta = %.9f" % beta_pred) print(" club analytic difference") print(" T_black (C) %18.9f %18.9f %10.3e" % (float(np.median(run["Tb"][ico])) - KELVIN, T_OPT + d_pred - KELVIN, float(np.median(run["Tb"][ico])) - (T_OPT + d_pred))) print(" T_white (C) %18.9f %18.9f %10.3e" % (float(np.median(run["Tw"][ico])) - KELVIN, T_OPT - d_pred - KELVIN, float(np.median(run["Tw"][ico])) - (T_OPT - d_pred))) print(" bare ground a_g %18.9f %18.9f %10.3e" % (float(np.median(run["ag"][ico])), ag_pred, float(np.median(run["ag"][ico])) - ag_pred)) print(" spread of a_g across the interior of the coexistence band: %.3e" % span(run["ag"][ico])) print(" So the bare-ground fraction is pinned at %.6f for the whole" % ag_pred) print(" coexistence band, independent of L. Everything the luminosity does") print(" is redistribute the remaining %.6f between black and white." % (1.0 - ag_pred)) print(" VERDICT: %s" % ("PASS" if abs(float(np.median(run["ag"][ico])) - ag_pred) < 1e-6 else "FAIL")) print() print(" sample of the coexistence region:") print(" L T_black(K) T_white(K) sum(K) 2T_opt(K) diff(K)") idxs = np.where(both)[0] for i in idxs[::max(1, len(idxs) // 9)]: print(" %5.3f %10.5f %11.5f %11.6f %10.2f %+.3e" % (Lmain[i], run["Tb"][i], run["Tw"][i], run["Tb"][i] + run["Tw"][i], 2 * T_OPT, run["Tb"][i] + run["Tw"][i] - 2 * T_OPT)) # ---------------------------------------------------------------- [6] header("[6] THE CONTROL THE ARGUMENT NEEDS: EQUAL ALBEDOS") print("Set both daisy albedos equal to bare ground. The daisies still grow,") print("compete and die, so the ecology is untouched. What is removed is the") print("only channel by which they can influence temperature. If regulation") print("survives this, the regulation was an artefact of something else.") print() flat = sweep(Lmain, Ab=A_GROUND, Aw=A_GROUND) dflat = np.abs(flat["Te"] - analytic) print(" L flat-daisy Te(C) bare control(C) difference(K) daisy cover") for L in (0.6, 0.7, 0.8, 1.0, 1.1, 1.2, 1.4, 1.6): i = int(np.argmin(np.abs(Lmain - L))) print(" %5.2f %15.5f %15.5f %13.3e %11.5f" % (Lmain[i], flat["Te"][i] - KELVIN, analytic[i] - KELVIN, flat["Te"][i] - analytic[i], flat["ab"][i] + flat["aw"][i])) print() print(" club value %.3e K required value 0 K difference %.3e K" % (float(dflat.max()), float(dflat.max()))) print(" temperature span of the equal-albedo planet : %.3f C" % span(flat["Te"] - KELVIN)) print(" temperature span of the bare planet : %.3f C" % span(analytic - KELVIN)) print(" daisies are alive at %d of %d luminosities, peak cover %.4f" % (int(((flat["ab"] + flat["aw"]) > ALIVE).sum()), L_N_MAIN, float((flat["ab"] + flat["aw"]).max()))) print(" so the biosphere is present, large, and completely without effect.") print(" VERDICT: %s" % ("PASS, regulation vanishes entirely" if dflat.max() < 1e-9 else "FAIL")) # ---------------------------------------------------------------- [7] header("[7] ROBUSTNESS 1: NARROWING THE GROWTH WINDOW") print("h is the half-width of the temperature range in which a daisy can grow") print("at all. The published value is 17.5 K. Narrow it and the biota") print("becomes fussier. How fussy can it be and still regulate?") print() print(" h(K) L_low L_high width drift(C) bare span(C) suppression peak cover") h_rows = [] for hv in (25.0, 22.5, 20.0, 17.5, 15.0, 12.5, 10.0, 7.5, 5.0, 4.0, 3.0, 2.0, 1.0): r = sweep(Lscan, h=hv) tot = r["ab"] + r["aw"] al = tot > ALIVE if not al.any(): print(" %7.1f %s" % (hv, "no living population anywhere in the sweep")) h_rows.append((hv, 0.0, 0.0, 0.0, float("nan"), float("nan"), 0.0)) continue lo, hi, rr = longest_run(Lscan, al) tt = r["Te"][rr] - KELVIN b2 = bare_planet(Lscan)[rr] - KELVIN sp = span(b2) / span(tt) if span(tt) > 1e-9 else float("inf") print(" %7.1f %5.3f %6.3f %5.3f %8.3f %11.3f %11.2f %10.4f" % (hv, lo, hi, hi - lo, span(tt), span(b2), sp, float(tot.max()))) h_rows.append((hv, lo, hi, hi - lo, span(tt), sp, float(tot.max()))) print() print(" FIGURE DATA 4: h, regulated width, drift, suppression") for row in h_rows: print(" %6.2f %7.4f %9.4f %9.3f" % (row[0], row[3], row[4], row[5])) # ---------------------------------------------------------------- [8] header("[8] ROBUSTNESS 2: A DEATH RATE SHOCK") print("Raise gamma. Steady state needs a_g beta = gamma, and both a_g and") print("beta are at most 1, so gamma = 1 is the hard ceiling at which no") print("population can persist at any temperature whatsoever.") print() print(" gamma L_low L_high width drift(C) suppression peak cover") g_rows = [] for gv in (0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 0.95, 0.99): r = sweep(Lscan, gamma=gv) tot = r["ab"] + r["aw"] al = tot > ALIVE if not al.any(): print(" %7.2f %s" % (gv, "extinct at every luminosity")) g_rows.append((gv, 0.0, 0.0, 0.0, float("nan"), float("nan"), 0.0)) continue lo, hi, rr = longest_run(Lscan, al) tt = r["Te"][rr] - KELVIN b2 = bare_planet(Lscan)[rr] - KELVIN sp = span(b2) / span(tt) if span(tt) > 1e-9 else float("inf") print(" %7.2f %5.3f %6.3f %5.3f %8.3f %11.2f %10.4f" % (gv, lo, hi, hi - lo, span(tt), sp, float(tot.max()))) g_rows.append((gv, lo, hi, hi - lo, span(tt), sp, float(tot.max()))) print() print(" FIGURE DATA 5: gamma, regulated width, drift, suppression") for row in g_rows: print(" %6.2f %7.4f %9.4f %9.3f" % (row[0], row[3], row[4], row[5])) print() print(" Now a transient shock instead of a permanent change: hold gamma at") print(" the published 0.3, sweep to L = 1.0, then kill a fraction of the") print(" standing crop instantly and integrate back to steady state.") print() i1 = int(np.argmin(np.abs(Lmain - 1.0))) ab1, aw1, ag1 = float(run["ab"][i1]), float(run["aw"][i1]), float(run["ag"][i1]) A1 = ab1 * A_BLACK + aw1 * A_WHITE + ag1 * A_GROUND Te1 = (S_FLUX * 1.0 * (1.0 - A1) / SIGMA) ** 0.25 print(" pre-shock state at L = 1.000: black %.5f, white %.5f, bare %.5f" % (ab1, aw1, ag1)) print(" pre-shock temperature : %.5f C" % (Te1 - KELVIN)) print() print(" kill immediate T(C) recovered T(C) recovery error(C) steps") for kill in (0.10, 0.25, 0.50, 0.75, 0.90, 0.99, 0.999): b2_, w2_ = ab1 * (1 - kill), aw1 * (1 - kill) g2_ = 1.0 - b2_ - w2_ A2 = b2_ * A_BLACK + w2_ * A_WHITE + g2_ * A_GROUND Timm = (S_FLUX * (1.0 - A2) / SIGMA) ** 0.25 _, _, _, Trec, st = relax(b2_, w2_, g2_, 1.0) print(" %6.1f%% %13.4f %14.4f %17.2e %6d" % (100 * kill, Timm - KELVIN, Trec - KELVIN, Trec - Te1, st)) print() print(" The planet returns to the same steady state after losing 99.9 per") print(" cent of its biota, because the attractor is a property of the") print(" equations rather than of the standing crop.") # ---------------------------------------------------------------- [9] header("[9] MONTE CARLO OVER PARAMETER UNCERTAINTY") print("Every number so far comes from one point in parameter space, chosen by") print("Watson and Lovelock in 1983. None of it is measured, because there is") print("nothing to measure it against. So the useful question is not what the") print("regulated width is, but how much it depends on choices nobody can") print("check. We draw parameters at random and sweep each drawn planet.") print() N_TRIALS = 2000 rng = np.random.Generator(np.random.PCG64(streams[0])) Ab_s = np.clip(rng.normal(A_BLACK, 0.05, N_TRIALS), 0.02, 0.48) Aw_s = np.clip(rng.normal(A_WHITE, 0.05, N_TRIALS), 0.52, 0.98) q_s = np.clip(rng.normal(Q_HEAT, 0.30e9, N_TRIALS), 0.2e9, 6.0e9) g_s = np.clip(rng.normal(GAMMA, 0.06, N_TRIALS), 0.05, 0.95) To_s = rng.normal(T_OPT, 4.0, N_TRIALS) h_s = np.clip(rng.normal(HALFWIDTH, 3.5, N_TRIALS), 2.0, 35.0) print(" trials : %d" % N_TRIALS) print(" stream : PCG64 spawned from SeedSequence(%d), child 0" % MASTER_SEED) print(" A_black ~ N(%.2f, %.2f^2) clipped to [0.02, 0.48]" % (A_BLACK, 0.05)) print(" A_white ~ N(%.2f, %.2f^2) clipped to [0.52, 0.98]" % (A_WHITE, 0.05)) print(" q ~ N(%.3e, %.2e^2) clipped to [2e8, 6e9]" % (Q_HEAT, 0.30e9)) print(" gamma ~ N(%.2f, %.2f^2) clipped to [0.05, 0.95]" % (GAMMA, 0.06)) print(" T_opt ~ N(%.2f, %.2f^2) K" % (T_OPT, 4.0)) print(" h ~ N(%.2f, %.2f^2) K clipped to [2, 35]" % (HALFWIDTH, 3.5)) print() Lmc = np.linspace(L_LO, L_HI, L_N_MC) t0 = time.time() Bm, Wm, Tm, nst = sweep_vec(Lmc, Ab_s, Aw_s, q_s, g_s, To_s, h_s) print(" vectorised sweep: %d planets x %d luminosities, %d RK4 steps, %.1f s" % (N_TRIALS, L_N_MC, nst, time.time() - t0)) tot_mc = Bm + Wm alive_mc = tot_mc > ALIVE dL = Lmc[1] - Lmc[0] width = alive_mc.sum(axis=0) * dL Tc_mc = Tm - KELVIN bare_col = bare_planet(Lmc) - KELVIN drift = np.full(N_TRIALS, np.nan) supp = np.full(N_TRIALS, np.nan) for j in range(N_TRIALS): m = alive_mc[:, j] if m.sum() >= 2: tt = Tc_mc[m, j] drift[j] = span(tt) if span(tt) > 1e-9: supp[j] = span(bare_col[m]) / span(tt) ok = np.isfinite(drift) n_reg = int((width > 0).sum()) print() print(" trials with a living population somewhere : %d of %d (%.2f%%)" % (n_reg, N_TRIALS, 100.0 * n_reg / N_TRIALS)) print(" trials measurable for drift (>= 2 points) : %d" % int(ok.sum())) def mc_report(name, vals, unit): v = np.asarray(vals, dtype=float) v = v[np.isfinite(v)] m = float(v.mean()) s = float(v.std(ddof=1)) se = s / math.sqrt(len(v)) print(" %-28s mean %9.4f %-2s SD %8.4f SE %7.4f n = %d" % (name, m, unit, s, se, len(v))) print(" %-28s median %7.4f 5th pct %7.4f 95th pct %7.4f" % ("", float(np.median(v)), float(np.percentile(v, 5)), float(np.percentile(v, 95)))) return m, s, se w_m, w_s, w_se = mc_report("regulated luminosity width", width[width > 0], "L") d_m, d_s, d_se = mc_report("residual temperature drift", drift[ok], "C") s_ok = supp[np.isfinite(supp)] s_m, s_s, s_se = mc_report("suppression factor", s_ok, "x") print() print(" baseline at the published parameters, same %d-point grid:" % L_N_MC) base_mc = sweep(Lmc) al_b = (base_mc["ab"] + base_mc["aw"]) > ALIVE lo_b, hi_b, rr_b = longest_run(Lmc, al_b) tt_b = base_mc["Te"][rr_b] - KELVIN bb_b = bare_col[rr_b] print(" width %.4f, drift %.4f C, suppression %.3f" % (float(al_b.sum()) * dL, span(tt_b), span(bb_b) / span(tt_b))) z_w = (float(al_b.sum()) * dL - w_m) / w_se print(" baseline width vs Monte Carlo mean: %.4f vs %.4f +- %.4f -> %.2f SE" % (float(al_b.sum()) * dL, w_m, w_se, z_w)) z_d = (span(tt_b) - d_m) / d_se print(" baseline drift vs Monte Carlo mean: %.4f vs %.4f +- %.4f -> %.2f SE" % (span(tt_b), d_m, d_se, z_d)) print() print(" fraction of drawn planets with suppression > 2 : %.4f" % float(np.mean(s_ok > 2))) print(" fraction with suppression > 5 : %.4f" % float(np.mean(s_ok > 5))) print(" fraction with suppression > 10 : %.4f" % float(np.mean(s_ok > 10))) print(" fraction with residual drift < 5 C : %.4f" % float(np.mean(drift[ok] < 5.0))) print(" fraction with residual drift < 15 C : %.4f" % float(np.mean(drift[ok] < 15.0))) print(" fraction that never regulate at all : %.4f" % (1.0 - float(n_reg) / N_TRIALS)) print() print(" CONVERGENCE: running mean and standard error as trials accumulate") print(" trials mean width SE mean drift(C) SE mean supp SE") wpos = width[width > 0] dpos = drift[ok] spos = s_ok for n in (10, 25, 50, 100, 200, 400, 700, 1000, 1400, 1700, len(wpos)): n = min(n, len(wpos)) a = wpos[:n] b = dpos[:min(n, len(dpos))] c = spos[:min(n, len(spos))] print(" %6d %11.4f %7.4f %14.4f %7.4f %10.3f %7.3f" % (n, a.mean(), a.std(ddof=1) / math.sqrt(n), b.mean(), b.std(ddof=1) / math.sqrt(len(b)), c.mean(), c.std(ddof=1) / math.sqrt(len(c)))) print() print(" FIGURE DATA 3: convergence, every 25th trial: n, running mean width, SE") for n in range(25, len(wpos) + 1, 25): a = wpos[:n] print(" %5d %8.5f %8.5f" % (n, a.mean(), a.std(ddof=1) / math.sqrt(n))) # ---------------------------------------------------------------- [10] header("[10] WHERE A DIFFERENT MODELLING CHOICE CHANGES THE ANSWER") print("Assumption 5 was that seeds are always present, so a dead population") print("can return the moment conditions allow. Drop it: seed once at the") print("cold end and never again.") print() noseed = sweep(Lmain, seed_floor=None, ab0=SEED_FLOOR, aw0=SEED_FLOOR) al_ns = (noseed["ab"] + noseed["aw"]) > ALIVE if al_ns.any(): lo2, hi2, rr2 = longest_run(Lmain, al_ns) t2 = noseed["Te"][rr2] - KELVIN print(" with re-seeding : regulated L %.3f to %.3f, width %.3f, drift %.3f C" % (Llo, Lhi, Lhi - Llo, drift0)) print(" without re-seeding : regulated L %.3f to %.3f, width %.3f, drift %.3f C" % (lo2, hi2, hi2 - lo2, span(t2))) print(" difference in width: %.3f (%.1f%% of the re-seeded value)" % ((Lhi - Llo) - (hi2 - lo2), 100 * ((Lhi - Llo) - (hi2 - lo2)) / (Lhi - Llo))) else: print(" without re-seeding : nothing survives anywhere in the sweep") print() print(" History dependence. Sweep the same grid downward from the hot end.") down = sweep(Lmain[::-1]) al_dn = (down["ab"] + down["aw"]) > ALIVE lo3, hi3, rr3 = longest_run(Lmain[::-1], al_dn) t3 = down["Te"][rr3] - KELVIN print(" descending: regulated L %.3f to %.3f, width %.3f, drift %.3f C" % (min(lo3, hi3), max(lo3, hi3), abs(hi3 - lo3), span(t3))) dnT = down["Te"][::-1] hyst = float(np.max(np.abs(run["Te"] - dnT))) ihy = int(np.argmax(np.abs(run["Te"] - dnT))) print(" max |T(ascending) - T(descending)| = %.4f C, at L = %.3f" % (hyst, Lmain[ihy])) print(" hysteresis under this seeding rule : %s" % ("absent to within 0.01 C" if hyst < 0.01 else "PRESENT, %.2f C" % hyst)) print() print(" A third convention: seed every luminosity fresh from bare ground,") print(" with no memory of the previous step at all.") fresh_lo, fresh_hi, fresh_n = None, None, 0 freshT = np.zeros(L_N_MAIN) freshC = np.zeros(L_N_MAIN) for i, L in enumerate(Lmain): b_, w_, g_, T_, _ = relax(SEED_FLOOR, SEED_FLOOR, 1.0 - 2 * SEED_FLOOR, L) freshT[i] = T_ freshC[i] = b_ + w_ al_f = freshC > ALIVE if al_f.any(): fl, fh, rrf = longest_run(Lmain, al_f) tf = freshT[rrf] - KELVIN print(" fresh-seed: regulated L %.3f to %.3f, width %.3f, drift %.3f C" % (fl, fh, fh - fl, span(tf))) print(" against continuation width %.3f: the upper edge moves from" % (Lhi - Llo)) print(" L = %.3f to L = %.3f, a loss of %.1f%% of the regulated range." % (Lhi, fh, 100 * (Lhi - fh) / (Lhi - Llo))) print(" Reason: at high L a bare planet is already too hot for white") print(" daisies to germinate, so they can only be there if they were") print(" there already. Regulation at the hot end is inherited.") print() print(" Heat transport q, the softest parameter in the model:") print(" q (K^4) L_low L_high width drift(C) suppression") q_rows = [] for qv in (0.0, 0.2e9, 0.5e9, 1.0e9, 2.06e9, 4.0e9, 8.0e9): r = sweep(Lscan, q=qv) tot = r["ab"] + r["aw"] al = tot > ALIVE if not al.any(): print(" %10.3e %s" % (qv, "no living population")) continue lo, hi, rr = longest_run(Lscan, al) tt = r["Te"][rr] - KELVIN b2 = bare_planet(Lscan)[rr] - KELVIN sp = span(b2) / span(tt) if span(tt) > 1e-9 else float("inf") print(" %10.3e %5.3f %6.3f %5.3f %8.3f %11.2f" % (qv, lo, hi, hi - lo, span(tt), sp)) q_rows.append((qv, hi - lo, span(tt), sp)) print() print(" q = 0 is a perfectly conducting planet: every patch sits at the") print(" planetary mean, black and white daisies feel identical temperatures,") print(" and the two species become one. Whatever survives there is the part") print(" of the result that owes nothing to having two colours.") # ---------------------------------------------------------------- [11] header("[11] SUMMARY: CLUB VALUE BESIDE THE VALUE IT SHOULD MATCH") rows = [ ("area sum a_b+a_w+a_g, worst RK4 step", "%.3e" % worst, "0", "%.3e" % worst), ("bare planet T vs Stefan-Boltzmann (K)", "%.3e" % float(d.max()), "0", "%.3e" % float(d.max())), ("T_black + T_white in coexistence (K)", "%.6f" % float(symi.mean()), "%.6f" % (2 * T_OPT), "%.3e" % float(symi.mean() - 2 * T_OPT)), ("bare ground in coexistence band", "%.9f" % float(np.median(run["ag"][ico])), "%.9f" % ag_pred, "%.3e" % float(np.median(run["ag"][ico]) - ag_pred)), ("T_black in coexistence band (C)", "%.6f" % (float(np.median(run["Tb"][ico])) - KELVIN), "%.6f" % (T_OPT + d_pred - KELVIN), "%.3e" % (float(np.median(run["Tb"][ico])) - (T_OPT + d_pred))), ("steady state |a_g beta - gamma|, interior", "%.3e" % interior_worst, "0", "%.3e" % interior_worst), ("heat coefficient 4 T_opt^3 q' (K^4)", "%.4e" % q_derived, "%.4e" % Q_HEAT, "%.3e" % (q_derived - Q_HEAT)), ("equal-albedo control vs bare (K)", "%.3e" % float(dflat.max()), "0", "%.3e" % float(dflat.max())), ] print(" %-40s %14s %14s %13s" % ("quantity", "club", "analytic", "difference")) for r_ in rows: print(" %-40s %14s %14s %13s" % r_) print() print(" headline numbers for the article:") print(" regulated luminosity interval : %.3f to %.3f (width %.3f)" % (Llo, Lhi, Lhi - Llo)) print(" residual temperature drift inside it : %.3f C" % drift0) print(" bare planet span over the same interval: %.3f C" % bare0) print(" suppression factor : %.2f" % supp0) print(" mean regulated temperature : %.3f C" % Tin.mean()) print(" growth optimum : %.2f C" % (T_OPT - KELVIN)) print(" bare planet span over the full sweep : %.3f C" % span(analytic - KELVIN)) if tight.any(): print(" tight band |T - T_opt| < 2 C : L %.3f to %.3f (width %.3f)" % (tl, th, th - tl)) print(" drift inside the tight band : %.3f C" % span(run["Te"][trr])) print(" suppression inside the tight band : %.2f" % (span(analytic[trr]) / span(run["Te"][trr]))) print(" largest single-step temperature jump : %+.3f C at L = %.3f" % (dT[jmp], Lmain[jmp + 1])) print(" Monte Carlo regulated width : %.4f +- %.4f (SE, n = %d)" % (w_m, w_se, n_reg)) print(" Monte Carlo residual drift : %.4f +- %.4f C" % (d_m, d_se)) print(" Monte Carlo suppression factor : %.3f +- %.3f" % (s_m, s_se)) print(" peak total daisy cover : %.4f at L = %.3f" % (float(Ttot.max()), Lmain[int(np.argmax(Ttot))])) print() print(rule("=")) print("total wall clock: %.1f s" % (time.time() - t_start)) print(rule("=")) if __name__ == "__main__": main()