#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
snowball-flicker.py
Science Journaling Club -- Field Notes -- analysis for
"Snowball Earth Was Flickering: A Planet Stuck in a 56-Million-Year Loop"

WHAT THIS IS
------------
A deliberately small, transparent, zero-dimensional (globally averaged) energy
balance model (EBM) coupled to a one-box carbonate-silicate carbon cycle, built
by the club to get an *intuition* for the result reported by Minsky, Wordsworth,
Johnston & Knoll (2026, PNAS, doi:10.1073/pnas.2525919123).

THIS IS NOT THEIR MODEL. They use a coupled box model of Neoproterozoic climate,
carbon, and oxygen. We use two ordinary differential equations. Every number this
script prints should be read as "a toy model that behaves the same way", not as a
reproduction of the paper.

PART A -- ICE-ALBEDO BISTABILITY (the Budyko-Sellers hysteresis)
----------------------------------------------------------------
Energy balance for a single global-mean surface temperature T (kelvin):

    C_T dT/dt = ASR(T) - OLR(T, pCO2)

    ASR(T)         = (S/4) * (1 - alpha(T))              absorbed shortwave
    alpha(T)       = a_ice + (a_warm - a_ice)/2 * (1 + tanh((T - T_a)/dT_a))
    OLR(T, pCO2)   = OLR_0 + B*(T - T_0) - F(pCO2)       outgoing longwave
    F(pCO2)        = k_CO2 * ln(pCO2 / pCO2_ref)         greenhouse forcing

alpha(T) is the classic smooth ice-albedo feedback: cold planet -> bright planet.
The tanh is a stand-in for "how much of the globe is ice-covered", not a fit to data.

Because ASR(T) is S-shaped and OLR(T) is a straight line, the two can cross at
one, two, or three temperatures. Three crossings = BISTABILITY: an ice-free
branch and a frozen branch, separated by an unstable branch that nothing can sit
on. Sweeping the solar constant S traces the hysteresis loop, and the two fold
(saddle-node) points are where a branch simply ceases to exist and the planet
jumps.

PART B -- THE CARBON CYCLE AND THE LIMIT CYCLE
-----------------------------------------------
Add a second, much slower equation for the exchangeable carbon inventory,
expressed directly as pCO2 (bar):

    d(pCO2)/dt = (V - W(T, pCO2)) / (beta * N_bar)

    V  = constant volcanic/metamorphic CO2 outgassing (mol C / yr)
    W  = silicate weathering sink, WHAK-type (Walker, Hays & Kasting 1981):

    W(T,pCO2) = W_0 * fW * (pCO2/pCO2_ref)^n * exp((T - T_0)/T_e) * g(T)

    g(T) = 0.5*(1 + tanh((T - T_a)/dT_a))   ice-free land fraction: weathering
                                            SHUTS OFF when the planet freezes.
    fW   = "weatherability" multiplier -- the Franklin basalt knob. Fresh flood
           basalt weathers far faster than old cratonic granite, so emplacing a
           large igneous province at the tropics multiplies W at fixed T and pCO2.

The shut-off is the whole story. On a warm planet, weathering is a thermostat:
too hot -> more weathering -> less CO2 -> cooler. On a frozen planet the
thermostat is unplugged, CO2 piles up at the outgassing rate, and the planet
eventually thaws. If the weathering-balanced steady state falls *inside* the
bistable gap -- i.e. on the unstable branch, where nothing can rest -- the system
cannot settle anywhere and instead relaxation-oscillates: a LIMIT CYCLE.

ASSUMPTIONS WE ARE MAKING (all of them debatable)
-------------------------------------------------
 1. Zero dimensions. No latitudes, no ocean circulation, no ice sheets, no
    seasons, no clouds. A "Jormungand"/slushball state with a thin open tropical
    belt (Abbot, Voigt & Koll 2011) cannot exist in this model by construction.
 2. Linear OLR in T (Budyko). True to about +/-10 K, extrapolated here over ~70 K.
 3. Logarithmic CO2 forcing with k_CO2 = 5.35 W/m^2 per e-fold, calibrated for
    modern Earth near 280 ppm and extrapolated to ~0.03 bar. This is the single
    largest source of error in the deglaciation threshold we report.
 4. Linear partitioning of carbon between atmosphere and ocean (buffer factor
    beta), held constant even when the ocean is sealed under ice.
 5. Constant outgassing V. No seafloor weathering, no organic carbon burial,
    no carbonate compensation, no continental drift.
 6. Faint young Sun at 717 Ma via the standard Gough parameterisation.

Python 3.12. numpy + scipy. RNG seeded (only used for a tiny robustness check).

Run:  python snowball-flicker.py
"""

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy.integrate import solve_ivp

RNG = np.random.default_rng(20260427)  # anchor paper publication date, why not

# ----------------------------------------------------------------------------
# PARAMETERS
# ----------------------------------------------------------------------------

S_NOW = 1361.0          # W/m^2, present-day total solar irradiance
AGE_GA = 0.717          # Ga before present -- Sturtian onset
# Gough (1981) faint young Sun:  S(t) = S_now / (1 + 0.4*(1 - t/t_now))
S_STURT = S_NOW / (1.0 + 0.4 * (AGE_GA / 4.57))

A_ICE = 0.50            # planetary albedo, fully glaciated (dust-darkened sea glacier)
A_WARM = 0.30           # planetary albedo, ice-free
T_A = 262.0             # K, centre of the ice-albedo transition
DT_A = 9.0              # K, width of the ice-albedo transition

T_0 = 288.0             # K, modern global mean surface temperature
OLR_0 = 238.2           # W/m^2, modern OLR at T_0 and pCO2_ref
B_OLR = 1.30            # W/m^2/K, Budyko OLR slope
K_CO2 = 5.35            # W/m^2 per e-fold of CO2
PCO2_REF = 280e-6       # bar, pre-industrial reference

C_T = 12.67             # W yr /m^2/K  (~100 m ocean mixed layer, 4.0e8 J/m^2/K)

V_OUT = 7.5e12          # mol C / yr, volcanic + metamorphic outgassing
W_0 = 7.5e12            # mol C / yr, modern silicate weathering sink (balances V)
N_EXP = 0.30            # weathering CO2 exponent (WHAK)
T_E = 13.7              # K, weathering temperature e-folding scale
BETA = 2.5              # atmosphere+ocean exchangeable carbon buffer factor
N_BAR = 1.182e20        # mol CO2 per bar of partial pressure over the whole Earth

MYR = 1.0e6

# ----------------------------------------------------------------------------
# MODEL FUNCTIONS
# ----------------------------------------------------------------------------

def albedo(T):
    return A_ICE + 0.5 * (A_WARM - A_ICE) * (1.0 + np.tanh((T - T_A) / DT_A))

def ice_free_fraction(T):
    return 0.5 * (1.0 + np.tanh((T - T_A) / DT_A))

def asr(T, S):
    return 0.25 * S * (1.0 - albedo(T))

def forcing(pco2):
    return K_CO2 * np.log(pco2 / PCO2_REF)

def olr(T, pco2):
    return OLR_0 + B_OLR * (T - T_0) - forcing(pco2)

def net_flux(T, S, pco2):
    return asr(T, S) - olr(T, pco2)

def weathering(T, pco2, fW):
    return W_0 * fW * (pco2 / PCO2_REF) ** N_EXP * np.exp((T - T_0) / T_E) * ice_free_fraction(T)

# ----------------------------------------------------------------------------
# ROOT FINDING / BRANCH TRACING
# ----------------------------------------------------------------------------

def equilibria(S, pco2, lo=150.0, hi=430.0, n=12000):
    """All roots of net_flux(T)=0, with stability (dNet/dT<0 => stable)."""
    T = np.linspace(lo, hi, n)
    f = net_flux(T, S, pco2)
    roots = []
    sign_change = np.where(np.sign(f[:-1]) != np.sign(f[1:]))[0]
    for i in sign_change:
        a, b = T[i], T[i + 1]
        for _ in range(80):
            m = 0.5 * (a + b)
            if np.sign(net_flux(a, S, pco2)) == np.sign(net_flux(m, S, pco2)):
                a = m
            else:
                b = m
        r = 0.5 * (a + b)
        h = 1e-4
        slope = (net_flux(r + h, S, pco2) - net_flux(r - h, S, pco2)) / (2 * h)
        roots.append((r, "stable" if slope < 0 else "unstable"))
    return roots

def folds_in_T():
    """Analytic fold temperatures: where d(Net)/dT = 0, i.e. (S/4)|dalpha/dT| = B.
    Independent of pCO2 because CO2 forcing shifts OLR but not its slope."""
    amp = 0.5 * (A_ICE - A_WARM)          # 0.10
    # dalpha/dT = -amp/dT_a * sech^2(x)
    k = 0.25 * S_STURT * amp / DT_A
    if k < B_OLR:
        return None
    sech2 = B_OLR / k
    x = np.arccosh(1.0 / np.sqrt(sech2))
    return T_A - DT_A * x, T_A + DT_A * x   # cold fold (thaw), warm fold (freeze)

def pco2_at_fold(T_fold):
    """CO2 level at which the given fold temperature is an equilibrium."""
    # asr = OLR_0 + B(T-T0) - k ln(p/pref)
    need_F = OLR_0 + B_OLR * (T_fold - T_0) - asr(T_fold, S_STURT)
    return PCO2_REF * np.exp(need_F / K_CO2)

# ----------------------------------------------------------------------------
# TWO-EQUATION COUPLED SYSTEM
# ----------------------------------------------------------------------------

def rhs(t, y, fW, V=V_OUT):
    T, p = y
    p = max(p, 1e-8)
    dT = net_flux(T, S_STURT, p) / C_T
    dp = (V - weathering(T, p, fW)) / (BETA * N_BAR)
    return [dT, dp]

def integrate(fW, t_end_myr=140.0, T0=300.0, p0=4.0e-4, V=V_OUT, max_step_yr=None):
    t_end = t_end_myr * MYR
    kw = dict(method="Radau", rtol=1e-9, atol=[1e-6, 1e-14], dense_output=True)
    if max_step_yr:
        kw["max_step"] = max_step_yr
    sol = solve_ivp(rhs, [0.0, t_end], [T0, p0], args=(fW, V), **kw)
    return sol

def cycle_stats(sol, discard_frac=0.45):
    """Measure period and amplitude from upward crossings of the albedo midpoint."""
    t_end = sol.t[-1]
    t = np.linspace(discard_frac * t_end, t_end, 400000)
    T, p = sol.sol(t)
    crossings = []
    above = T > T_A
    for i in range(1, len(t)):
        if above[i] and not above[i - 1]:
            # linear interpolation of the crossing time
            f = (T_A - T[i - 1]) / (T[i] - T[i - 1])
            crossings.append(t[i - 1] + f * (t[i] - t[i - 1]))
    if len(crossings) < 3:
        return dict(cycling=False, n_cycles=len(crossings), Tmin=T.min(), Tmax=T.max(),
                    pmin=p.min(), pmax=p.max(), period=np.nan, amp=np.nan,
                    frozen_frac=np.mean(T < T_A))
    per = np.diff(np.array(crossings))
    # analyse the last complete cycle
    a, b = crossings[-2], crossings[-1]
    m = (t >= a) & (t <= b)
    Tc, pc = T[m], p[m]
    frozen = np.mean(Tc < T_A)
    return dict(cycling=True, n_cycles=len(crossings), period=float(np.mean(per[-5:])),
                period_sd=float(np.std(per[-5:])), Tmin=float(Tc.min()), Tmax=float(Tc.max()),
                amp=float(Tc.max() - Tc.min()), pmin=float(pc.min()), pmax=float(pc.max()),
                frozen_frac=float(frozen))

# ----------------------------------------------------------------------------
# OUTPUT
# ----------------------------------------------------------------------------

OUT = []
def say(s=""):
    OUT.append(s)
    print(s)

def rule(title):
    say()
    say("=" * 74)
    say(title)
    say("=" * 74)

# ============================================================================
rule("0.  SETUP")
say(f"Present-day solar constant           S_now    = {S_NOW:.1f} W/m^2")
say(f"Solar constant at 717 Ma (Gough)     S_sturt  = {S_STURT:.1f} W/m^2"
    f"   ({100*S_STURT/S_NOW:.1f}% of today)")
say(f"Insolation deficit vs today                   = {S_NOW - S_STURT:.1f} W/m^2 "
    f"(= {0.25*(S_NOW-S_STURT)*0.7:.1f} W/m^2 of absorbed flux at albedo 0.30)")
say(f"Albedo:  ice-free {A_WARM:.2f}  ->  fully glaciated {A_ICE:.2f}"
    f"   (transition centred {T_A:.0f} K, width {DT_A:.0f} K)")
say(f"OLR slope B = {B_OLR:.2f} W/m^2/K  =>  equilibrium sensitivity to 2xCO2 "
    f"= {K_CO2*np.log(2)/B_OLR:.2f} K")
say(f"Outgassing V = {V_OUT:.2e} mol C/yr; modern weathering W_0 = {W_0:.2e} mol C/yr")

# ============================================================================
rule("1.  ICE-ALBEDO HYSTERESIS: SWEEPING THE SOLAR CONSTANT (pCO2 fixed at 280 ppm)")

S_grid = np.arange(1000.0, 1801.0, 2.0)
warm_S, warm_T, cold_S, cold_T, unst_S, unst_T = [], [], [], [], [], []
n_states = []
for S in S_grid:
    eq = equilibria(S, PCO2_REF)
    n_states.append(len(eq))
    st = sorted([r for r, k in eq if k == "stable"])
    un = sorted([r for r, k in eq if k == "unstable"])
    if len(st) == 2:
        cold_S.append(S); cold_T.append(st[0])
        warm_S.append(S); warm_T.append(st[1])
    elif len(st) == 1:
        if st[0] > T_A:
            warm_S.append(S); warm_T.append(st[0])
        else:
            cold_S.append(S); cold_T.append(st[0])
    for u in un:
        unst_S.append(S); unst_T.append(u)

n_states = np.array(n_states)
bistable = S_grid[n_states == 3]
S_freeze = bistable.min()    # below this, only the frozen branch survives
S_thaw = bistable.max()      # above this, only the ice-free branch survives
say(f"Bistable window in solar constant:  {S_freeze:.0f}  to  {S_thaw:.0f} W/m^2")
say(f"Width of the bistable window:       {S_thaw - S_freeze:.0f} W/m^2 "
    f"= {100*(S_thaw-S_freeze)/S_freeze:.1f}% of the lower edge")
say(f"Sturtian solar constant {S_STURT:.0f} W/m^2 sits "
    f"{'INSIDE' if S_freeze <= S_STURT <= S_thaw else 'OUTSIDE'} that window "
    f"-> both a frozen and an ice-free Earth are possible at 717 Ma.")
say(f"Today's {S_NOW:.0f} W/m^2 also sits "
    f"{'INSIDE' if S_freeze <= S_NOW <= S_thaw else 'OUTSIDE'} it.")
say()
say("Equilibrium temperatures at selected solar constants (K):")
say(f"{'S (W/m^2)':>10} {'cold branch':>13} {'unstable':>11} {'warm branch':>13} {'#states':>8}")
for S in [1100, 1200, 1250, S_STURT, 1300, 1361, 1450, 1600]:
    eq = equilibria(float(S), PCO2_REF)
    st = sorted([r for r, k in eq if k == "stable"])
    un = sorted([r for r, k in eq if k == "unstable"])
    c = f"{st[0]:.1f}" if st and st[0] < T_A else "-"
    w = f"{st[-1]:.1f}" if st and st[-1] > T_A else "-"
    u = f"{un[0]:.1f}" if un else "-"
    say(f"{S:>10.0f} {c:>13} {u:>11} {w:>13} {len(eq):>8}")
say()
say("HYSTERESIS DEMO -- ramp S down then back up (quasi-static, 280 ppm CO2):")
T = 300.0
down_jump = None
for S in np.arange(1400.0, 1099.0, -1.0):
    eq = [r for r, k in equilibria(S, PCO2_REF) if k == "stable"]
    T = min(eq, key=lambda r: abs(r - T))
    if T < T_A and down_jump is None:
        down_jump = (S, T)
up_jump = None
for S in np.arange(1100.0, 1801.0, 1.0):
    eq = [r for r, k in equilibria(S, PCO2_REF) if k == "stable"]
    T = min(eq, key=lambda r: abs(r - T))
    if T > T_A and up_jump is None:
        up_jump = (S, T)
say(f"  cooling: Earth snaps into the snowball at S = {down_jump[0]:.0f} W/m^2 "
    f"(lands at T = {down_jump[1]:.1f} K = {down_jump[1]-273.15:.1f} C)")
say(f"  warming: Earth escapes the snowball at   S = {up_jump[0]:.0f} W/m^2 "
    f"(lands at T = {up_jump[1]:.1f} K = {up_jump[1]-273.15:.1f} C)")
say(f"  the two thresholds differ by {up_jump[0]-down_jump[0]:.0f} W/m^2 -- "
    f"THAT GAP IS THE HYSTERESIS.")

# ============================================================================
rule("2.  THE SAME LOOP IN CO2 SPACE (solar constant fixed at the Sturtian value)")

T_thaw_fold, T_freeze_fold = folds_in_T()
p_thaw = pco2_at_fold(T_thaw_fold)
p_freeze = pco2_at_fold(T_freeze_fold)
say("Fold (saddle-node) temperatures are set by the condition (S/4)|d(alpha)/dT| = B.")
say(f"  cold-branch fold (deglaciation): T = {T_thaw_fold:.1f} K = {T_thaw_fold-273.15:.1f} C")
say(f"  warm-branch fold (freeze-in)   : T = {T_freeze_fold:.1f} K = {T_freeze_fold-273.15:.1f} C")
say()
say(f"FREEZE THRESHOLD    pCO2 = {p_freeze:.3e} bar = {p_freeze*1e6:.0f} ppmv "
    f"({p_freeze/PCO2_REF:.2f}x pre-industrial)")
say(f"DEGLACIATION THRESHOLD pCO2 = {p_thaw:.4f} bar = {p_thaw*1e6:.0f} ppmv "
    f"({p_thaw/PCO2_REF:.0f}x pre-industrial)")
say(f"Ratio thaw/freeze = {p_thaw/p_freeze:.0f}x  "
    f"({np.log10(p_thaw/p_freeze):.2f} orders of magnitude of CO2)")
say("  [literature deglaciation estimates span ~0.01-0.1 bar; we land inside that range,")
say("   but the log-forcing extrapolation makes this an order-of-magnitude claim only]")
say()
say("Branch temperatures at the two thresholds:")
for label, p in [("just before freeze-in", p_freeze), ("just after deglaciation", p_thaw)]:
    eq = equilibria(S_STURT, p)
    st = sorted([r for r, k in eq if k == "stable"])
    say(f"  {label:<26} pCO2={p:.3e} bar  stable T = "
        + ", ".join(f"{r:.1f} K ({r-273.15:+.1f} C)" for r in st))
say()
say("Peak-to-peak temperature swing implied by the two jumps:")
Tcold_at_freeze = sorted([r for r, k in equilibria(S_STURT, p_freeze) if k == "stable"])[0]
Twarm_at_thaw = sorted([r for r, k in equilibria(S_STURT, p_thaw) if k == "stable"])[-1]
say(f"  hothouse {Twarm_at_thaw-273.15:+.1f} C  ->  snowball {Tcold_at_freeze-273.15:+.1f} C "
    f"= {Twarm_at_thaw - Tcold_at_freeze:.1f} K amplitude")

# ============================================================================
rule("3.  WHERE DOES THE CARBON CYCLE WANT TO SIT? (the bifurcation in weatherability)")

say("A steady state needs W(T,pCO2) = V with T on a STABLE branch.")
say("If the only solution sits on the unstable branch, no steady state exists.")
say()
say(f"{'fW':>6} {'warm-branch steady state':>30} {'exists?':>10}")
fw_scan = [0.5, 1.0, 2.0, 3.0, 3.5, 3.75, 4.0, 5.0, 8.0, 12.0]
for fW in fw_scan:
    # solve W(T_eq(p), p) = V along the warm branch
    ps = np.logspace(np.log10(p_freeze * 0.2), np.log10(0.2), 3000)
    Ts, Ws, keep = [], [], []
    for p in ps:
        st = sorted([r for r, k in equilibria(S_STURT, p) if k == "stable"])
        warm = [r for r in st if r > T_freeze_fold]
        if warm:
            Ts.append(warm[0]); Ws.append(weathering(warm[0], p, fW)); keep.append(p)
    Ws = np.array(Ws); keep = np.array(keep); Ts = np.array(Ts)
    idx = np.where(np.sign(Ws[:-1] - V_OUT) != np.sign(Ws[1:] - V_OUT))[0]
    if len(idx):
        i = idx[0]
        say(f"{fW:>6.2f} {f'pCO2={keep[i]:.2e} bar, T={Ts[i]-273.15:+.1f} C':>30} {'YES':>10}")
    else:
        say(f"{fW:>6.2f} {'none on the warm branch':>30} {'NO':>10}")

# bisect the critical weatherability
_WARM_BRANCH = None

def warm_branch_grid():
    """(pCO2, warm-branch T) on a fixed grid. The branch geometry does not depend
    on fW, so build it once and reuse it for every step of the bisection below."""
    global _WARM_BRANCH
    if _WARM_BRANCH is None:
        ps = np.logspace(np.log10(p_freeze * 0.2), np.log10(0.2), 2000)
        keep, Ts = [], []
        for p in ps:
            st = sorted([r for r, k in equilibria(S_STURT, p) if k == "stable"])
            warm = [r for r in st if r > T_freeze_fold]
            if warm:
                keep.append(p); Ts.append(warm[0])
        _WARM_BRANCH = (np.array(keep), np.array(Ts))
    return _WARM_BRANCH

def warm_ss_exists(fW):
    """True iff W(T,pCO2) = V has a solution *on the warm branch*.

    NOTE (fixed): the earlier version returned True as soon as any grid point had
    W >= V. Weathering is enormous at the high-pCO2 end of the warm branch, so that
    test was satisfied for every fW, the bisection collapsed onto its lower bound,
    and fW* printed as 1.00 -- contradicting both the table above and the forward
    integrations below. What matters is whether the weathering curve CROSSES V while
    still on the branch, which is exactly the sign-change test used in the table."""
    keep, Ts = warm_branch_grid()
    if len(keep) < 2:
        return False
    Ws = weathering(Ts, keep, fW)
    return bool(np.any(np.sign(Ws[:-1] - V_OUT) != np.sign(Ws[1:] - V_OUT)))

# warm_ss_exists is TRUE for small fW and FALSE for large fW, so the bracket has
# to close downward on a False and upward on a True. The previous version had the
# two assignments the other way round, which walked `lo` all the way to the upper
# bound and printed fW* = 20.00 -- contradicting the table above (a warm steady
# state exists at 3.50 and not at 3.75) and the forward integrations below.
lo, hi = 1.0, 20.0
for _ in range(45):
    mid = 0.5 * (lo + hi)
    if warm_ss_exists(mid):
        lo = mid
    else:
        hi = mid
FW_CRIT = 0.5 * (lo + hi)
say()
say(f"CRITICAL WEATHERABILITY  fW* = {FW_CRIT:.2f}")
say(f"  Below {FW_CRIT:.2f}x, the planet finds a warm steady state and stays there forever.")
say(f"  Above {FW_CRIT:.2f}x, no steady state exists on either stable branch -> LIMIT CYCLE.")
say("  Basalt weathers roughly 5-10x faster than granite at the same temperature,")
say("  so emplacing the Franklin LIP in the wet tropics plausibly crosses this line.")

# ============================================================================
rule("4.  FORWARD INTEGRATION: DOES IT ACTUALLY CYCLE?")

say("Two coupled ODEs, stiff (climate relaxes in ~13 yr, carbon in ~10^6 yr).")
say("Integrated with an implicit Radau solver for 140 Myr from a warm start.")
say()
say(f"{'fW':>5} {'regime':>13} {'period (Myr)':>13} {'amp (K)':>9} "
    f"{'T range (C)':>20} {'pCO2 range (bar)':>24} {'% frozen':>9} {'cycles/56 Myr':>14}")

results = {}
for fW in [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 16.0]:
    sol = integrate(fW)
    st = cycle_stats(sol)
    results[fW] = st
    if st["cycling"]:
        n56 = 56.0 / (st["period"] / MYR)
        trange = "%+.1f to %+.1f" % (st["Tmin"] - 273.15, st["Tmax"] - 273.15)
        prange = "%.2e - %.2e" % (st["pmin"], st["pmax"])
        say(f"{fW:>5.1f} {'LIMIT CYCLE':>13} {st['period']/MYR:>13.3f} {st['amp']:>9.1f} "
            f"{trange:>20} {prange:>24} {100*st['frozen_frac']:>8.1f}% {n56:>14.0f}")
    else:
        trange = "%+.1f" % (st["Tmax"] - 273.15)
        prange = "%.2e" % st["pmax"]
        say(f"{fW:>5.1f} {'fixed point':>13} {'-':>13} {'-':>9} "
            f"{trange:>20} {prange:>24} {100*st['frozen_frac']:>8.1f}% {'-':>14}")

# ============================================================================
rule("5.  THE HEADLINE RUN  (fW = 5, a Franklin-scale weatherability boost)")

FW_NOM = 5.0
sol = integrate(FW_NOM, t_end_myr=140.0)
st = results[FW_NOM]
say(f"Period                    P = {st['period']/MYR:.3f} Myr "
    f"(std over last 5 cycles {st['period_sd']/MYR:.4f} Myr)")
say(f"Temperature amplitude     dT = {st['amp']:.1f} K")
say(f"  hothouse peak  {st['Tmax']-273.15:+.1f} C     snowball floor {st['Tmin']-273.15:+.1f} C")
say(f"pCO2 range                {st['pmin']*1e6:.0f} ppmv -> {st['pmax']*1e6:.0f} ppmv "
    f"({st['pmax']/st['pmin']:.0f}x, {np.log10(st['pmax']/st['pmin']):.2f} orders of magnitude)")
say(f"Fraction of each cycle spent frozen      {100*st['frozen_frac']:.1f}%")
say(f"Fraction of each cycle spent ice-free    {100*(1-st['frozen_frac']):.1f}%")
say(f"Cycles fitted into the 56 Myr Sturtian   {56.0/(st['period']/MYR):.1f}")
say(f"Mean duration of one frozen interval     {st['frozen_frac']*st['period']/MYR:.3f} Myr")
say(f"Mean duration of one ice-free interval   {(1-st['frozen_frac'])*st['period']/MYR:.3f} Myr")
say()
say("Timescale accounting -- why the period is what it is:")
dM = BETA * N_BAR * (p_thaw - p_freeze)
say(f"  carbon that must be added to go freeze->thaw:  {dM:.3e} mol C")
say(f"  at V = {V_OUT:.2e} mol/yr that recharge takes  {dM/V_OUT/MYR:.3f} Myr")
say(f"  -> the FROZEN half of the cycle is outgassing-limited and CANNOT be shortened")
say(f"     by weatherability. Only the warm half responds to fW.")

# ============================================================================
rule("6.  HOW THE PERIOD RESPONDS TO WEATHERABILITY (the Franklin knob)")

say(f"{'fW':>6} {'period (Myr)':>13} {'frozen (Myr)':>13} {'ice-free (Myr)':>15} "
    f"{'cycles/56 Myr':>14} {'% frozen':>9}")
sweep = [4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 16.0]
sweep_rows = []
for fW in sweep:
    s = results[fW]
    if not s["cycling"]:
        continue
    P = s["period"] / MYR
    fr = s["frozen_frac"] * P
    ic = P - fr
    sweep_rows.append((fW, P, fr, ic, 56.0 / P, 100 * s["frozen_frac"]))
    say(f"{fW:>6.1f} {P:>13.3f} {fr:>13.3f} {ic:>15.3f} {56.0/P:>14.0f} "
        f"{100*s['frozen_frac']:>8.1f}%")
say()
p4, p16 = results[4.0]["period"] / MYR, results[16.0]["period"] / MYR
say(f"Quadrupling weatherability from 4x to 16x shortens the period by only "
    f"{100*(1-p16/p4):.0f}% ({p4:.2f} -> {p16:.2f} Myr)")
say(f"and drives the frozen fraction from {100*results[4.0]['frozen_frac']:.0f}% up to "
    f"{100*results[16.0]['frozen_frac']:.0f}% of each cycle.")
say("The cycle has a FLOOR: you cannot flicker faster than volcanoes can refill the sky.")
say(f"Asymptotic floor (pure recharge time) = {dM/V_OUT/MYR:.3f} Myr.")

# ============================================================================
rule("7.  SENSITIVITY: WHAT IF OUTGASSING WAS LOWER? (Dutkiewicz et al. 2024)")

say("Independent tectonic reconstructions argue mid-ocean-ridge outgassing at 717 Ma")
say("was exceptionally low. In this model that lengthens the frozen half directly.")
say(f"{'V / V_modern':>13} {'period (Myr)':>13} {'frozen (Myr)':>13} {'cycles/56 Myr':>14}")
for frac in [0.4, 0.6, 0.8, 1.0, 1.3]:
    s2 = cycle_stats(integrate(FW_NOM, t_end_myr=200.0, V=V_OUT * frac))
    if s2["cycling"]:
        P = s2["period"] / MYR
        say(f"{frac:>13.1f} {P:>13.3f} {s2['frozen_frac']*P:>13.3f} {56.0/P:>14.0f}")
    else:
        say(f"{frac:>13.1f} {'fixed point':>13} {'-':>13} {'-':>14}")

# ============================================================================
rule("8.  ROBUSTNESS: JITTER THE PARAMETERS")

say("200 draws, each parameter perturbed by a uniform +/-10% (seeded RNG),")
say("to check that cycling is not an artefact of one lucky parameter set.")
_B, _TA, _DTA, _AI, _TE = B_OLR, T_A, DT_A, A_ICE, T_E
n_cyc, periods = 0, []
for _ in range(200):
    j = RNG.uniform(0.9, 1.1, 5)
    B_OLR, T_A, DT_A, A_ICE, T_E = _B * j[0], _TA * j[1], _DTA * j[2], min(_AI * j[3], 0.68), _TE * j[4]
    try:
        s3 = cycle_stats(integrate(FW_NOM, t_end_myr=120.0))
        if s3["cycling"]:
            n_cyc += 1
            periods.append(s3["period"] / MYR)
    except Exception:
        pass
B_OLR, T_A, DT_A, A_ICE, T_E = _B, _TA, _DTA, _AI, _TE
periods = np.array(periods)
say(f"Limit cycles found in {n_cyc}/200 draws = {100*n_cyc/200:.0f}%")
if len(periods):
    say(f"Period across those draws: median {np.median(periods):.2f} Myr, "
        f"10th-90th percentile {np.percentile(periods,10):.2f}-{np.percentile(periods,90):.2f} Myr")
say("So the OSCILLATION is robust; the PERIOD is not tightly pinned. That is the honest summary.")

# ============================================================================
rule("9.  FIGURE DATA")

say("--- FIG 1: hysteresis in (S, T), 280 ppm CO2 ---")
say("warm_branch S,T (every 40th point):")
say("  " + "; ".join(f"{s:.0f},{t:.1f}" for s, t in list(zip(warm_S, warm_T))[::40]))
say("cold_branch S,T (every 40th point):")
say("  " + "; ".join(f"{s:.0f},{t:.1f}" for s, t in list(zip(cold_S, cold_T))[::40]))
say("unstable_branch S,T (every 20th point):")
say("  " + "; ".join(f"{s:.0f},{t:.1f}" for s, t in list(zip(unst_S, unst_T))[::20]))

say()
say("--- FIG 2: the S-curve in (log10 pCO2, T) at S_sturt, plus the limit cycle ---")
p_grid = np.logspace(-5.0, -0.7, 260)
rows = []
for p in p_grid:
    eq = equilibria(S_STURT, p)
    st_ = sorted([r for r, k in eq if k == "stable"])
    un_ = sorted([r for r, k in eq if k == "unstable"])
    rows.append((p, st_, un_))
say("warm branch (log10 p, T):")
say("  " + "; ".join(f"{np.log10(p):.3f},{s[-1]:.1f}" for p, s, u in rows if s and s[-1] > T_freeze_fold)[:2000])
say("cold branch (log10 p, T):")
say("  " + "; ".join(f"{np.log10(p):.3f},{s[0]:.1f}" for p, s, u in rows if s and s[0] < T_thaw_fold)[:2000])
say("unstable branch (log10 p, T):")
say("  " + "; ".join(f"{np.log10(p):.3f},{u[0]:.1f}" for p, s, u in rows if u)[:2000])
say(f"fold points: thaw (log10 p, T) = ({np.log10(p_thaw):.3f}, {T_thaw_fold:.1f}) ; "
    f"freeze = ({np.log10(p_freeze):.3f}, {T_freeze_fold:.1f})")

say()
say("--- FIG 3: time series of the headline run (one cycle, fW=5) ---")
t_end = sol.t[-1]
tt = np.linspace(0.60 * t_end, 0.60 * t_end + 3.2 * st["period"], 700)
TT, PP = sol.sol(tt)
tt = (tt - tt[0]) / MYR
say("t_Myr,T_C,log10pCO2 (every 14th sample):")
say("  " + "; ".join(f"{a:.3f},{b-273.15:.1f},{np.log10(c):.3f}"
                      for a, b, c in list(zip(tt, TT, PP))[::14]))

say()
say("--- FIG 4: period vs weatherability ---")
say("fW,period_Myr,frozen_Myr,icefree_Myr:")
for r in sweep_rows:
    say(f"  {r[0]:.1f},{r[1]:.4f},{r[2]:.4f},{r[3]:.4f}")
say(f"fW_crit = {FW_CRIT:.2f}; recharge floor = {dM/V_OUT/MYR:.3f} Myr")

say()
say("--- SPINE MARKER: 11 points equally spaced by ARC LENGTH around one cycle ---")
tt2 = np.linspace(0.70 * t_end, 0.70 * t_end + st["period"], 40000)
T2, P2 = sol.sol(tt2)
x2 = np.log10(P2); y2 = T2
xn = (x2 - np.log10(p_freeze)) / (np.log10(p_thaw) - np.log10(p_freeze))
yn = (y2 - st["Tmin"]) / (st["Tmax"] - st["Tmin"])
ds = np.hypot(np.diff(xn), np.diff(yn))
s_cum = np.concatenate([[0], np.cumsum(ds)])
targets = np.linspace(0, s_cum[-1] * 0.999, 11)
say("section,log10pCO2,T_C,phase_frac")
for i, tg in enumerate(targets, start=1):
    k = int(np.searchsorted(s_cum, tg))
    k = min(k, len(tt2) - 1)
    say(f"  {i},{x2[k]:.3f},{y2[k]-273.15:.1f},{(tt2[k]-tt2[0])/st['period']:.4f}")

# ============================================================================
rule("10.  KEY NUMBERS FOR THE ARTICLE")
say(f"S_sturt                        {S_STURT:.0f} W/m^2 ({100*S_STURT/S_NOW:.1f}% of modern)")
say(f"bistable solar window          {S_freeze:.0f}-{S_thaw:.0f} W/m^2")
say(f"freeze threshold               {p_freeze*1e6:.0f} ppmv CO2")
say(f"deglaciation threshold         {p_thaw:.3f} bar = {p_thaw*1e6:.0f} ppmv CO2")
say(f"CO2 swing                      {p_thaw/p_freeze:.0f}x")
say(f"critical weatherability fW*    {FW_CRIT:.2f}x")
say(f"period at fW=5                 {results[5.0]['period']/MYR:.2f} Myr")
say(f"temperature amplitude          {results[5.0]['amp']:.0f} K "
    f"({results[5.0]['Tmin']-273.15:+.0f} C to {results[5.0]['Tmax']-273.15:+.0f} C)")
say(f"cycles in 56 Myr at fW=5       {56.0/(results[5.0]['period']/MYR):.0f}")
say(f"frozen fraction at fW=5        {100*results[5.0]['frozen_frac']:.0f}%")
say(f"period floor (recharge time)   {dM/V_OUT/MYR:.2f} Myr")
say(f"robustness                     cycles in {100*n_cyc/200:.0f}% of jittered draws")
say()
say("REMINDER: club toy model. Not the Minsky et al. box model. Numbers are")
say("illustrative of the MECHANISM, not predictions about the real Cryogenian.")

with open(__file__.replace(".py", "-output.txt"), "w", encoding="utf-8") as f:
    f.write("\n".join(OUT) + "\n")