#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
layered-greenhouse.py
Science Journaling Club, Volume 2 Issue 3, Spring 2026, "Earth's Energy Budget".

QUESTION
--------
The N-layer grey atmosphere is the standard classroom picture of the greenhouse
effect. How many layers does it take to reach Earth's actual surface temperature
of about 288 K, and at what point does the model stop telling the truth about the
real atmosphere?

WHAT THIS FILE IS
-----------------
This is a computation, not an observation. Nobody in the club measured a flux, a
temperature or a spectrum. Every physical input below is a published number taken
from the literature and cited in the accompanying article; everything else is
arithmetic performed by this file. The computation is the experiment.

MODEL
-----
(1) N-layer grey atmosphere in radiative equilibrium.
    Shortwave: the atmosphere is perfectly transparent, the surface absorbs
    F = S(1-alpha)/4 averaged over the globe.
    Longwave: each layer is a perfect absorber and a perfect emitter (emissivity
    1), isothermal, and radiates sigma*T^4 upward and sigma*T^4 downward.
    Unknowns are x_j = sigma*T_j^4 for j = 0 (surface) .. N (top layer).
      surface:   F + x_1 = x_0
      layer i:   x_{i-1} + x_{i+1} = 2 x_i,  with x_{N+1} = 0 (space)
    This is solved as a dense linear system by numpy.linalg.solve for each N and
    checked against the exact analytic solution T_s = T_e (N+1)^(1/4).

(2) Single layer with tunable longwave emissivity eps (a grey slab that is not a
    perfect absorber). Layer absorbs eps*sigma*T_s^4 and emits eps*sigma*T_a^4
    each way. Result: sigma*T_s^4 = F / (1 - eps/2). We invert this for the eps
    that reproduces the observed surface temperature.

(3) Two-band single layer. The longwave spectrum is split at the edges of the
    atmospheric window (default 8-12 micron). Band fractions are obtained by
    numerically integrating the Planck function, so they depend on temperature.
    The window has optical depth tau_w and the rest of the spectrum tau_b, with
    emissivity 1 - exp(-tau) in each. The surface temperature is found by fixed
    point iteration because the band fractions move as the surface moves.

(4) Monte Carlo. Published uncertainties on the solar constant, the planetary
    albedo and the observed surface temperature are propagated to the inferred
    emissivity and the inferred fractional layer count. A second Monte Carlo
    propagates uncertainty in the window edges and the window optical depth to
    the grey-versus-two-band response ratio. Standard errors are computed from
    the runs themselves and convergence is recorded as trials accumulate.

ASSUMPTIONS, STATED PLAINLY
---------------------------
 * Global annual mean. One number for the whole planet, one number for the whole
   year. No latitude, no season, no day and night.
 * Pure radiative equilibrium. No convection, no latent heat, no sensible heat.
   The real troposphere moves roughly 100 W/m^2 upward by non-radiative means.
 * Layers are isothermal and infinitesimally thin, and there is no atmosphere
   between them. The real atmosphere is a continuum with a lapse rate.
 * Grey longwave (model 1 and 2). The real absorption spectrum of water vapour
   and CO2 is a forest of lines, not a constant.
 * The atmosphere is transparent to sunlight. It is not: roughly 23% of incoming
   solar is absorbed in the air, not at the ground.
 * Surface emissivity is exactly 1.
 * Albedo is held fixed while the greenhouse strength is varied, which is
   physically incoherent for large changes and is flagged where it matters.
 * sigma * mean(T)^4 is used where the honest quantity is mean(sigma * T^4).
   Section 8 quantifies the size of that error.

LIMITATIONS THE ARTICLE DISCUSSES AT LENGTH
-------------------------------------------
 * No integer number of perfectly absorbing layers reproduces 288 K.
 * The two-band single layer cannot simultaneously reproduce the observed
   surface temperature and the observed 22 W/m^2 of surface emission that
   escapes directly to space through the window. Section 12 measures the size
   of the contradiction.
 * The grey assumption badly overstates the temperature response to added
   absorber, because a grey atmosphere cannot saturate a band.

SEED
----
Master seed 20260321, hard coded below, fed to numpy's SeedSequence and spawned
into independent streams for each Monte Carlo. Generator is PCG64. The whole
output is deterministic.

RUN
---
    python layered-greenhouse.py > layered-greenhouse-output.txt
"""

import sys
import time

import numpy as np

T_START = time.time()

SEED = 20260321

# ---------------------------------------------------------------------------
# Physical constants (CODATA 2018 / SI defining constants)
# ---------------------------------------------------------------------------
SIGMA = 5.670374419e-8      # Stefan-Boltzmann, W m^-2 K^-4  (exact from SI defs)
H_PLANCK = 6.62607015e-34   # J s   (exact)
C_LIGHT = 2.99792458e8      # m s^-1 (exact)
K_B = 1.380649e-23          # J K^-1 (exact)

# ---------------------------------------------------------------------------
# Published inputs.  Sources are cited in the article; nothing here was measured
# by the club.
# ---------------------------------------------------------------------------
S0 = 1361.0        # total solar irradiance, W/m^2   (Kopp & Lean 2011: 1360.8 +/- 0.5)
S0_SD = 0.5
ALBEDO = 0.293     # CERES EBAF planetary albedo     (Loeb et al. 2018)
ALBEDO_SD = 0.005
ALBEDO_TEXTBOOK = 0.30      # the round number most textbooks use
T_OBS = 288.0      # global mean near-surface air temperature, K
T_OBS_SD = 0.5
SURFACE_LW_OBS = 398.2      # observed surface upward longwave, W/m^2 (Wild et al. 2013)
OLR_OBS = 239.0             # observed outgoing longwave radiation, W/m^2
WINDOW_ESCAPE_OBS = 22.0    # OLR that is directly transmitted surface emission,
                            # W/m^2 (Costa & Shine 2012)
ALBEDO_SURFACE_ONLY = 0.15  # an airless Earth keeps only its surface albedo
# Lunar values, all from the NASA NSSDC Moon fact sheet
MOON_ALBEDO = 0.11
MOON_TE_NASA = 270.4        # blackbody temperature listed on that sheet, K
MOON_EQ_MIN = 95.0          # equatorial diurnal minimum, K
MOON_EQ_MAX = 390.0         # equatorial diurnal maximum, K

WIN_LO_UM = 8.0    # atmospheric window, microns
WIN_HI_UM = 12.0
TAU_W_PRESENT = 0.30        # present-day clear-sky window optical depth, the
                            # club's working value; scanned in section 13


def rule(title):
    print()
    print("=" * 78)
    print(title)
    print("=" * 78)


def eff_temp(s, albedo):
    """Effective emitting temperature from a solar constant and an albedo."""
    return (s * (1.0 - albedo) / 4.0 / SIGMA) ** 0.25


# ===========================================================================
print("=" * 78)
print("A GREENHOUSE MODEL YOU CAN CHECK BY HAND, AND WHERE IT STOPS BEING HONEST")
print("Science Journaling Club, Volume 2 Issue 3, Spring 2026")
print("=" * 78)
print("python      : %s" % sys.version.split()[0])
print("numpy       : %s" % np.__version__)
print("master seed : %d" % SEED)
print()
print("This file computes. It does not observe. Every physical input below is a")
print("published number; every other number is arithmetic done here.")

rule("1. INPUTS AS USED")
print("  solar constant S          = %10.4f W/m^2   (sd %.2f)" % (S0, S0_SD))
print("  planetary albedo alpha    = %10.4f        (sd %.4f)" % (ALBEDO, ALBEDO_SD))
print("  textbook albedo           = %10.4f" % ALBEDO_TEXTBOOK)
print("  absorbed solar F=S(1-a)/4 = %10.4f W/m^2" % (S0 * (1 - ALBEDO) / 4.0))
print("  observed T_surface        = %10.4f K       (sd %.2f)" % (T_OBS, T_OBS_SD))
print("  observed surface LW up    = %10.4f W/m^2" % SURFACE_LW_OBS)
print("  observed OLR              = %10.4f W/m^2" % OLR_OBS)
print("  observed window escape    = %10.4f W/m^2" % WINDOW_ESCAPE_OBS)
print("  Stefan-Boltzmann sigma    = %.11e W m^-2 K^-4" % SIGMA)

F_ABS = S0 * (1.0 - ALBEDO) / 4.0
T_E = eff_temp(S0, ALBEDO)
F_ABS_TB = S0 * (1.0 - ALBEDO_TEXTBOOK) / 4.0
T_E_TB = eff_temp(S0, ALBEDO_TEXTBOOK)


# ===========================================================================
rule("2. VALIDATION A. THE ZERO-LAYER CASE AGAINST THE ACCEPTED 255 K")
print("With no atmosphere in the longwave, the surface alone balances the absorbed")
print("sunlight: sigma T^4 = S(1-alpha)/4. Every textbook quotes about 255 K.")
print()
print("  %-42s %12s %12s %10s" % ("case", "club", "accepted", "diff (K)"))
ACCEPTED_TE = 255.0
for name, alb, val in (("albedo 0.293 (CERES, Loeb 2018)", ALBEDO, eff_temp(S0, ALBEDO)),
                       ("albedo 0.30  (round textbook value)", ALBEDO_TEXTBOOK,
                        eff_temp(S0, ALBEDO_TEXTBOOK))):
    print("  %-42s %12.4f %12.1f %+10.4f" % (name, val, ACCEPTED_TE, val - ACCEPTED_TE))
print()
print("  Both land within %.2f K of the quoted 255 K. The textbook figure is a" %
      max(abs(eff_temp(S0, ALBEDO) - 255.0), abs(eff_temp(S0, ALBEDO_TEXTBOOK) - 255.0)))
print("  rounding of a number that depends on which albedo you adopt; the spread")
print("  between the two albedos above is %.3f K, which is larger than the" %
      abs(eff_temp(S0, ALBEDO) - eff_temp(S0, ALBEDO_TEXTBOOK)))
print("  rounding error in either. VALIDATION A PASSES.")
print()
print("  Working value used from here on: T_e = %.4f K, F = %.4f W/m^2" % (T_E, F_ABS))


# ===========================================================================
def solve_layers(n, f_abs):
    """Solve the N-layer grey radiative-equilibrium system numerically.

    Unknowns x[0..n] are the blackbody fluxes sigma*T^4 at the surface (0) and
    at each layer (1..n, bottom to top). Returns (x, residual_max, cond)."""
    m = n + 1
    a = np.zeros((m, m), dtype=float)
    b = np.zeros(m, dtype=float)
    # surface: x0 - x1 = F      (with no layers, x0 = F)
    a[0, 0] = 1.0
    if n >= 1:
        a[0, 1] = -1.0
    b[0] = f_abs
    # layer i (row i): x_{i-1} - 2 x_i + x_{i+1} = 0 , x_{n+1} = 0
    for i in range(1, m):
        a[i, i - 1] = 1.0
        a[i, i] = -2.0
        if i + 1 < m:
            a[i, i + 1] = 1.0
        b[i] = 0.0
    x = np.linalg.solve(a, b)
    resid = np.max(np.abs(a @ x - b))
    cond = np.linalg.cond(a)
    return x, resid, cond


rule("3. VALIDATION B. NUMERIC AGAINST THE EXACT ANALYTIC SOLUTION")
print("The system has a closed-form answer: T_s = T_e (N+1)^(1/4), and the i-th")
print("layer counting up from the ground sits at T_i = T_e (N+1-i)^(1/4).")
print("Below, the numerically solved surface temperature is printed beside it.")
print()
print("  %3s %16s %16s %14s %14s %12s" %
      ("N", "T_s numeric (K)", "T_s analytic (K)", "difference", "rel. diff", "solve resid"))
max_abs = 0.0
max_rel = 0.0
rows_validation = []
for n in range(0, 11):
    x, resid, cond = solve_layers(n, F_ABS)
    ts_num = (x[0] / SIGMA) ** 0.25
    ts_ana = T_E * (n + 1.0) ** 0.25
    d = ts_num - ts_ana
    rel = abs(d) / ts_ana
    max_abs = max(max_abs, abs(d))
    max_rel = max(max_rel, rel)
    rows_validation.append((n, ts_num, ts_ana, d, rel, resid))
    print("  %3d %16.10f %16.10f %+14.3e %14.3e %12.2e" %
          (n, ts_num, ts_ana, d, rel, resid))
print()
print("  max |numeric - analytic| over N=0..10 : %.3e K" % max_abs)
print("  max relative difference               : %.3e" % max_rel)
print("  double-precision epsilon              : %.3e" % np.finfo(float).eps)
print("  ratio of max rel. diff to eps         : %.2f" % (max_rel / np.finfo(float).eps))
print()
if max_rel < 20 * np.finfo(float).eps:
    print("  Agreement is at machine precision. VALIDATION B PASSES.")
else:
    print("  *** DISAGREEMENT BEYOND MACHINE PRECISION. INVESTIGATE. ***")


# ===========================================================================
rule("4. VALIDATION B2. EVERY LEVEL, NOT JUST THE SURFACE, AND LARGE N")
print("For N = 5 the whole column is printed. Level 0 is the ground.")
print()
x5, r5, c5 = solve_layers(5, F_ABS)
print("  %5s %16s %16s %14s %14s" %
      ("level", "T numeric (K)", "T analytic (K)", "difference", "sigma T^4"))
for i in range(6):
    tn = (x5[i] / SIGMA) ** 0.25
    ta = T_E * (5 + 1 - i) ** 0.25
    print("  %5d %16.10f %16.10f %+14.3e %14.4f" % (i, tn, ta, tn - ta, x5[i]))
print()
print("  top-of-atmosphere check: flux leaving the top layer = %.6f W/m^2" % x5[5])
print("  absorbed solar                                      = %.6f W/m^2" % F_ABS)
print("  imbalance                                           = %+.3e W/m^2" % (x5[5] - F_ABS))
print()
print("Pushing N further, to see where the dense linear solve starts to hurt.")
print()
print("  %5s %18s %18s %14s %14s" % ("N", "T_s numeric", "T_s analytic", "rel. diff", "cond(A)"))
big_rows = []
for n in (20, 50, 100, 200, 500):
    x, resid, cond = solve_layers(n, F_ABS)
    tn = (x[0] / SIGMA) ** 0.25
    ta = T_E * (n + 1.0) ** 0.25
    big_rows.append((n, tn, ta, abs(tn - ta) / ta, cond))
    print("  %5d %18.9f %18.9f %14.3e %14.3e" % (n, tn, ta, abs(tn - ta) / ta, cond))
print()
print("  The condition number grows roughly as N^2, so the relative error grows")
print("  with it. Even at N = 500 it stays far below anything that could matter")
print("  physically; we print it because it is the honest place where 'exact'")
print("  turns into 'exact enough'.")


# ===========================================================================
rule("5. SURFACE TEMPERATURE AGAINST LAYER COUNT")
print("Observed global mean surface temperature: %.2f K" % T_OBS)
print("No-atmosphere effective temperature      : %.4f K" % T_E)
print()
print("  %3s %14s %12s %14s %16s" %
      ("N", "T_s (K)", "T_s (C)", "T_s - 288 (K)", "T_s - T_e (K)"))
table5 = []
for n in range(0, 11):
    ts = T_E * (n + 1.0) ** 0.25
    table5.append((n, ts))
    print("  %3d %14.4f %12.4f %+14.4f %+16.4f" %
          (n, ts, ts - 273.15, ts - T_OBS, ts - T_E))
print()
n_star = (T_OBS / T_E) ** 4 - 1.0
print("  The observed 288 K sits between N = 0 (%.2f K, %.2f K too cold) and" %
      (T_E, T_E - T_OBS))
print("  N = 1 (%.2f K, %.2f K too hot)." % (T_E * 2 ** 0.25, T_E * 2 ** 0.25 - T_OBS))
print()
print("  Solving T_e (N+1)^(1/4) = %.2f for a real-valued N:" % T_OBS)
print("      N* = (T_obs/T_e)^4 - 1 = %.6f" % n_star)
print()
print("  So the answer to the question in the title is: about two thirds of one")
print("  layer. No whole number of perfectly absorbing layers reproduces Earth.")
print("  One layer overshoots by %.2f K, which is %.0f%% of the entire 33 K" %
      (T_E * 2 ** 0.25 - T_OBS, 100 * (T_E * 2 ** 0.25 - T_OBS) / (T_OBS - T_E)))
print("  greenhouse effect the model is trying to explain.")
print()
print("  Greenhouse effect as this model defines it: %.4f K" % (T_OBS - T_E))
print("  Greenhouse effect as a flux, sigma T_obs^4 - F: %.3f W/m^2" %
      (SIGMA * T_OBS ** 4 - F_ABS))
print("  Same quantity from observations, LW_up - OLR:   %.3f W/m^2" %
      (SURFACE_LW_OBS - OLR_OBS))
print("  difference                                      %+.3f W/m^2 (%.1f%%)" %
      (SIGMA * T_OBS ** 4 - F_ABS - (SURFACE_LW_OBS - OLR_OBS),
       100 * (SIGMA * T_OBS ** 4 - F_ABS - (SURFACE_LW_OBS - OLR_OBS)) /
       (SURFACE_LW_OBS - OLR_OBS)))
print()
print("  That gap is not a bug. sigma*288^4 = %.2f W/m^2 but the measured" %
      (SIGMA * T_OBS ** 4))
print("  surface emission is %.1f W/m^2. Section 8 explains where the missing" % SURFACE_LW_OBS)
print("  %.1f W/m^2 comes from." % (SURFACE_LW_OBS - SIGMA * T_OBS ** 4))


# ===========================================================================
rule("6. THE SAME QUESTION FOR TWO OTHER PLANETS")
print("The model is not about Earth. It is about any transparent-to-sunlight,")
print("grey-in-the-infrared atmosphere. Published solar constants, albedos and")
print("mean surface temperatures, run through the same one-line formula.")
print()
# Every column here comes from the NASA NSSDC planetary fact sheets, which give
# a solar irradiance, a Bond albedo and their own blackbody temperature. That
# last column lets us check our formula against theirs.
bodies = [
    ("Mars",   586.2,  0.250, 214.0, 209.8),
    ("Earth",  1361.0, 0.294, 288.0, 254.0),
    ("Venus",  2601.3, 0.770, 737.0, 226.6),
]
print("  %-8s %10s %8s %12s %12s %10s %12s %12s" %
      ("body", "S (W/m2)", "albedo", "T_e club", "T_e NASA", "diff", "T mean", "N* layers"))
body_rows = []
for nm, s_irr, a, tobs, te_nasa in bodies:
    te = eff_temp(s_irr, a)
    nn = (tobs / te) ** 4 - 1.0
    body_rows.append((nm, s_irr, a, te, te_nasa, tobs, nn))
    print("  %-8s %10.1f %8.3f %12.3f %12.1f %+10.3f %12.1f %12.4f" %
          (nm, s_irr, a, te, te_nasa, te - te_nasa, tobs, nn))
print()
print("  Two of the three agree with NASA's own blackbody temperature to better")
print("  than %.2f K, which is a real check: their number and ours were computed" %
      max(abs(body_rows[0][3] - body_rows[0][4]), abs(body_rows[2][3] - body_rows[2][4])))
print("  independently from the same two inputs.")
print()
print("  Earth does not agree, and the discrepancy is %+.2f K. That is not a bug" %
      (body_rows[1][3] - body_rows[1][4]))
print("  in our code. Solving backwards, the albedo that reproduces NASA's listed")
print("  254.0 K from their listed 1361.0 W/m^2 is:")
alb_implied = 1.0 - 4.0 * SIGMA * 254.0 ** 4 / 1361.0
print("      alpha = %.4f" % alb_implied)
print("  against the %.3f printed on the same page. The older Earth Bond albedo" % 0.294)
print("  was 0.306. The fact sheet appears to carry an updated albedo beside a")
print("  blackbody temperature computed from the previous one. We mention it")
print("  because it is exactly the kind of thing a reader can check in one line,")
print("  and because our own numbers deserve the same treatment.")
print()
print("  Venus needs about %.0f layers of perfect absorber; Mars needs %.3f." %
      (body_rows[2][6], body_rows[0][6]))
print("  The model spans two orders of magnitude in solar flux and a factor of")
print("  three in surface temperature, which is the reason it earns a blackboard.")
print("  The Mars figure is the weakest row: the fact sheet gives an average")
print("  atmospheric temperature rather than a mean surface temperature, and we")
print("  have used it as though the two were the same thing. They are not.")


# ===========================================================================
rule("7. WHAT 'NO ATMOSPHERE' ACTUALLY MEANS")
print("The 33 K figure quietly assumes that removing the greenhouse effect leaves")
print("the albedo untouched. It would not. Three different airless Earths:")
print()
cases = [
    ("albedo 0.293 held fixed (the standard 33 K claim)", ALBEDO),
    ("albedo 0.15, surface only, clouds gone with the air", ALBEDO_SURFACE_ONLY),
    ("albedo 0.00, a perfect blackbody sphere", 0.0),
]
print("  %-52s %10s %12s" % ("case", "T (K)", "'greenhouse'"))
airless_rows = []
for nm, a in cases:
    t = eff_temp(S0, a)
    airless_rows.append((nm, a, t, T_OBS - t))
    print("  %-52s %10.3f %+12.3f" % (nm, t, T_OBS - t))
print()
print("  The famous 33 K becomes %.1f K if you take the clouds away with the" %
      airless_rows[1][3])
print("  atmosphere, which you would have to. That is a %.0f%% change in the" %
      (100 * abs(airless_rows[1][3] - airless_rows[0][3]) / airless_rows[0][3]))
print("  headline number produced entirely by a modelling choice.")


# ===========================================================================
rule("8. TWO PLACES THE GLOBAL MEAN LIES, MEASURED")
print("(a) sigma*mean(T)^4 is not mean(sigma*T^4). Jensen's inequality guarantees")
print("    the second is larger. Surface temperature varies across the planet with")
print("    a spatial-plus-seasonal standard deviation of roughly 15 to 20 K.")
print()
print("  %10s %18s %18s %14s" % ("sd (K)", "mean(sigma T^4)", "sigma mean(T)^4", "excess"))
for sd in (10.0, 15.0, 20.0, 25.0):
    # exact for a normal distribution: E[T^4] = m^4 + 6 m^2 s^2 + 3 s^4
    m = T_OBS
    e_t4 = m ** 4 + 6 * m ** 2 * sd ** 2 + 3 * sd ** 4
    print("  %10.1f %18.3f %18.3f %+14.3f" %
          (sd, SIGMA * e_t4, SIGMA * m ** 4, SIGMA * (e_t4 - m ** 4)))
lo, hi = 0.0, 80.0
for _ in range(200):
    mid = 0.5 * (lo + hi)
    val = SIGMA * (T_OBS ** 4 + 6 * T_OBS ** 2 * mid ** 2 + 3 * mid ** 4)
    if val < SURFACE_LW_OBS:
        lo = mid
    else:
        hi = mid
sd_needed = 0.5 * (lo + hi)
print()
print("  The measured surface emission of %.1f W/m^2 is recovered at sd = %.2f K," %
      (SURFACE_LW_OBS, sd_needed))
print("  which is squarely inside the plausible range. So the %.1f W/m^2 gap in" %
      (SURFACE_LW_OBS - SIGMA * T_OBS ** 4))
print("  section 5 is mostly Jensen's inequality, not a broken model. Expressed as")
print("  a temperature, the planet radiates like a %.2f K surface while its mean" %
      ((SURFACE_LW_OBS / SIGMA) ** 0.25))
print("  temperature is %.2f K, a gap of %.2f K." %
      (T_OBS, (SURFACE_LW_OBS / SIGMA) ** 0.25 - T_OBS))
print()
print("(b) The same inequality, run the other way, on an airless body. The Moon")
print("    receives the same sunlight as Earth and has no atmosphere at all, so")
print("    it is the closest thing to a control the Solar System offers.")
moon_te = eff_temp(S0, MOON_ALBEDO)
t_subsolar = (S0 * (1 - MOON_ALBEDO) / SIGMA) ** 0.25
# zero thermal inertia, non-rotating: T = T_ss * mu^(1/4) on the dayside, 0 at night
mu = np.linspace(0.0, 1.0, 2000001)
t_local = t_subsolar * mu ** 0.25
_trap = getattr(np, "trapezoid", None) or np.trapz
t_mean_zero_inertia = 0.5 * _trap(t_local, mu)
print()
print("  Moon T_e from S and albedo, computed here : %10.3f K" % moon_te)
print("  Moon T_e listed by NASA                   : %10.1f K" % MOON_TE_NASA)
print("  difference                                : %+10.3f K" % (moon_te - MOON_TE_NASA))
print("  Moon subsolar temperature                 : %10.3f K" % t_subsolar)
print("  zero-thermal-inertia area-mean (computed) : %10.3f K" % t_mean_zero_inertia)
print("  analytic check, 0.4 * T_subsolar          : %10.3f K" % (0.4 * t_subsolar))
print("  difference                                : %+10.3e K" %
      (t_mean_zero_inertia - 0.4 * t_subsolar))
print("  NASA equatorial diurnal range             : %6.0f to %.0f K" %
      (MOON_EQ_MIN, MOON_EQ_MAX))
print()
print("  Our T_e reproduces NASA's to %.2f K, so the inputs are being used the" %
      abs(moon_te - MOON_TE_NASA))
print("  same way. Now look at the third line. A body with no atmosphere and no")
print("  thermal inertia would average %.1f K, which is %.1f K BELOW the" %
      (t_mean_zero_inertia, moon_te - t_mean_zero_inertia))
print("  effective temperature that supposedly describes it. No greenhouse gas is")
print("  involved. The whole gap is Jensen's inequality: the Moon balances its")
print("  energy budget in sigma T^4, and the mean of T^4 has almost nothing to do")
print("  with the fourth power of the mean of T when T runs from %.0f to %.0f K." %
      (MOON_EQ_MIN, MOON_EQ_MAX))
print("  The real Moon has thermal inertia and so sits above %.0f K, but well" % t_mean_zero_inertia)
print("  below %.0f K." % moon_te)
print()
print("  Which means the 33 K we keep quoting for Earth is the gap to a fictitious")
print("  uniform-temperature airless Earth, not to any airless planet that could")
print("  actually exist.")


# ===========================================================================
rule("9. THE SINGLE LAYER WITH A TUNABLE EMISSIVITY")
print("Drop the perfect-absorber assumption. One layer, emissivity eps in the")
print("longwave, still transparent to sunlight.")
print("    layer:   eps sigma T_s^4 = 2 eps sigma T_a^4   ->  T_a^4 = T_s^4 / 2")
print("    surface: F + eps sigma T_a^4 = sigma T_s^4")
print("    =>       sigma T_s^4 = F / (1 - eps/2),   T_s = T_e (1 - eps/2)^(-1/4)")
print()
print("Note what the first line says: the layer temperature does not depend on")
print("eps at all. A thin veil and an opaque slab sit at the same temperature.")
print("They differ only in how much of the ground they can see.")
print()


def ts_from_eps(eps, te=T_E):
    return te * (1.0 - eps / 2.0) ** -0.25


def eps_from_ts(ts, te=T_E):
    return 2.0 * (1.0 - (te / ts) ** 4)


eps_star_guess = 2.0 * (1.0 - (T_E / T_OBS) ** 4)
print("  %8s %14s %14s %14s" % ("eps", "T_s (K)", "T_s - T_e (K)", "T_a (K)"))
eps_rows = []
for eps in (0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, eps_star_guess, 0.8, 0.9, 1.0):
    ts = ts_from_eps(eps)
    ta = ts / 2 ** 0.25
    eps_rows.append((eps, ts, ta))
    print("  %8.4f %14.4f %+14.4f %14.4f" % (eps, ts, ts - T_E, ta))
print()
eps_star = eps_from_ts(T_OBS)
# independent numerical confirmation by bisection, not by the inverted formula
lo, hi = 0.0, 1.0
for _ in range(200):
    mid = 0.5 * (lo + hi)
    if ts_from_eps(mid) < T_OBS:
        lo = mid
    else:
        hi = mid
eps_bisect = 0.5 * (lo + hi)
print("  Emissivity that reproduces T_obs = %.2f K" % T_OBS)
print("    by inverting the formula : %.10f" % eps_star)
print("    by bisection on T_s(eps) : %.10f" % eps_bisect)
print("    difference               : %+.3e" % (eps_star - eps_bisect))
print()
print("  eps* = %.4f. The layer temperature it implies is T_a = %.3f K," %
      (eps_star, T_OBS / 2 ** 0.25))
print("  which corresponds to an emission altitude of about %.2f km at a lapse" %
      ((T_OBS - T_OBS / 2 ** 0.25) / 6.5))
print("  rate of 6.5 K/km. The real effective emission level is around 5 km, so")
print("  this is the right order and about %.0f%% off." %
      (100 * abs((T_OBS - T_OBS / 2 ** 0.25) / 6.5 - 5.0) / 5.0))
print()
d_eps = (ts_from_eps(eps_star + 1e-6) - ts_from_eps(eps_star - 1e-6)) / 2e-6
print("  Sensitivity: d T_s / d eps at eps* = %.3f K per unit emissivity, so" % d_eps)
print("  one percentage point of emissivity is worth %.3f K here." % (0.01 * d_eps))


# ===========================================================================
rule("10. MONTE CARLO ON eps* AND N*")
print("Published uncertainties propagated to the two inferred quantities.")
print("  S      ~ Normal(%.1f, %.2f)  W/m^2   Kopp & Lean 2011" % (S0, S0_SD))
print("  albedo ~ Normal(%.3f, %.3f)          Loeb et al. 2018" % (ALBEDO, ALBEDO_SD))
print("  T_obs  ~ Normal(%.1f, %.2f)  K" % (T_OBS, T_OBS_SD))
print()
N_MC = 200000
ss = np.random.SeedSequence(SEED)
child_a, child_b = ss.spawn(2)
rng = np.random.default_rng(child_a)
s_draw = rng.normal(S0, S0_SD, N_MC)
a_draw = rng.normal(ALBEDO, ALBEDO_SD, N_MC)
t_draw = rng.normal(T_OBS, T_OBS_SD, N_MC)
te_draw = (s_draw * (1.0 - a_draw) / 4.0 / SIGMA) ** 0.25
eps_draw = 2.0 * (1.0 - (te_draw / t_draw) ** 4)
nstar_draw = (t_draw / te_draw) ** 4 - 1.0
ghe_draw = t_draw - te_draw

print("  trials: %d   stream: PCG64 spawned from SeedSequence(%d), child 0" % (N_MC, SEED))
print()
print("  %-22s %12s %12s %12s %16s" % ("quantity", "mean", "sd", "SE", "95% interval"))
for nm, arr in (("eps*", eps_draw), ("N* (layers)", nstar_draw),
                ("T_e (K)", te_draw), ("greenhouse (K)", ghe_draw)):
    se = arr.std(ddof=1) / np.sqrt(N_MC)
    lo95, hi95 = np.percentile(arr, [2.5, 97.5])
    print("  %-22s %12.6f %12.6f %12.3e %7.4f-%7.4f" %
          (nm, arr.mean(), arr.std(ddof=1), se, lo95, hi95))
print()
se_eps = eps_draw.std(ddof=1) / np.sqrt(N_MC)
print("  Point estimate from the central inputs: eps* = %.6f" % eps_star)
print("  Monte Carlo mean                      : eps* = %.6f" % eps_draw.mean())
print("  difference                            : %+.3e (%.2f standard errors)" %
      (eps_draw.mean() - eps_star, abs(eps_draw.mean() - eps_star) / se_eps))
print()
print("  The Monte Carlo mean sits a few standard errors off the point estimate")
print("  because eps* is a nonlinear function of the inputs, so the mean of the")
print("  function is not the function of the mean. The offset, %.2e, is orders" %
      abs(eps_draw.mean() - eps_star))
print("  of magnitude below the spread, so it changes nothing physical. We print")
print("  it rather than round it away.")
print()
print("  Which input dominates? One-at-a-time, each varied by its own sd:")
for nm, ds, da, dt in (("S", S0_SD, 0.0, 0.0),
                       ("albedo", 0.0, ALBEDO_SD, 0.0),
                       ("T_obs", 0.0, 0.0, T_OBS_SD)):
    te_hi = eff_temp(S0 + ds, ALBEDO + da)
    e_hi = 2.0 * (1.0 - (te_hi / (T_OBS + dt)) ** 4)
    print("    %-8s +1 sd moves eps* by %+.5f" % (nm, e_hi - eps_star))
print()
print("  Convergence of the running mean of eps*:")
print()
print("  %12s %14s %14s %14s" % ("trials", "running mean", "running SE", "dev from final"))
checkpoints = [100, 300, 1000, 3000, 10000, 30000, 60000, 100000, 150000, 200000]
csum = np.cumsum(eps_draw)
csum2 = np.cumsum(eps_draw ** 2)
conv_rows = []
final_mean = eps_draw.mean()
for k in checkpoints:
    m = csum[k - 1] / k
    v = max(csum2[k - 1] / k - m * m, 0.0) * k / (k - 1)
    se = np.sqrt(v / k)
    conv_rows.append((k, m, se, m - final_mean))
    print("  %12d %14.6f %14.6f %+14.3e" % (k, m, se, m - final_mean))
print()
print("  SE falls as 1/sqrt(n) as it must: the ratio of the SE at %d trials to" % checkpoints[0])
print("  the SE at %d is %.2f, against sqrt(%d/%d) = %.2f." %
      (checkpoints[-1], conv_rows[0][2] / conv_rows[-1][2], checkpoints[-1],
       checkpoints[0], np.sqrt(checkpoints[-1] / checkpoints[0])))


# ===========================================================================
rule("11. PLANCK BAND FRACTIONS, COMPUTED")
print("The two-band model needs to know what fraction of a blackbody's emission")
print("falls in the atmospheric window. That is an integral of the Planck")
print("function, done here by Simpson's rule on a fine wavelength grid.")
print()


def planck_lambda(lam_m, t):
    """Spectral radiance B_lambda, W m^-3 sr^-1. lam_m in metres."""
    lam_m = np.asarray(lam_m, dtype=float)
    t = np.asarray(t, dtype=float)
    a = 2.0 * H_PLANCK * C_LIGHT ** 2 / lam_m ** 5
    xx = H_PLANCK * C_LIGHT / (lam_m * K_B * t)
    return a / np.expm1(xx)


def simpson(y, dx, axis=-1):
    """Composite Simpson on a uniform grid with an even number of intervals."""
    y = np.moveaxis(np.asarray(y, dtype=float), axis, -1)
    n = y.shape[-1] - 1
    if n % 2 != 0:
        raise ValueError("Simpson needs an even number of intervals")
    w = np.ones(n + 1, dtype=float)
    w[1:-1:2] = 4.0
    w[2:-1:2] = 2.0
    return np.tensordot(y, w, axes=([-1], [0])) * dx / 3.0


def band_fraction(t, lo_um=WIN_LO_UM, hi_um=WIN_HI_UM, npts=401):
    """Fraction of blackbody emission at T between lo_um and hi_um microns."""
    t = np.atleast_1d(np.asarray(t, dtype=float))
    lo_m = np.atleast_1d(np.asarray(lo_um, dtype=float)) * 1e-6
    hi_m = np.atleast_1d(np.asarray(hi_um, dtype=float)) * 1e-6
    t, lo_m, hi_m = np.broadcast_arrays(t, lo_m, hi_m)
    u = np.linspace(0.0, 1.0, npts)
    lam = lo_m[..., None] + (hi_m - lo_m)[..., None] * u
    b = planck_lambda(lam, t[..., None])
    integral = simpson(b, 1.0 / (npts - 1), axis=-1) * (hi_m - lo_m)
    return np.pi * integral / (SIGMA * t ** 4)


C2 = H_PLANCK * C_LIGHT / K_B      # second radiation constant, m K


def frac_below_series(lam_um, t, terms=60):
    """Exact fraction of blackbody emission below wavelength lam, by the standard
    series F(0->x) = (15/pi^4) sum_n e^{-nx}(x^3/n + 3x^2/n^2 + 6x/n^3 + 6/n^4).
    Used only to validate the Simpson integrator; nothing else calls it."""
    x = C2 / (lam_um * 1e-6 * t)
    tot = 0.0
    for n in range(1, terms + 1):
        tot += np.exp(-n * x) * (x ** 3 / n + 3 * x ** 2 / n ** 2
                                 + 6 * x / n ** 3 + 6.0 / n ** 4)
    return 15.0 / np.pi ** 4 * tot


print("  VALIDATION C. The Simpson integrator against the exact series solution")
print("  for the blackbody fraction function. The series is closed form; no")
print("  quadrature is involved in it at all.")
print()
print("  %10s %8s %8s %18s %18s %14s" %
      ("T (K)", "lo (um)", "hi (um)", "Simpson", "exact series", "difference"))
val_c_max = 0.0
for tv, l1, l2 in ((288.0, 8.0, 12.0), (255.0, 8.0, 12.0), (242.0, 8.0, 12.0),
                   (737.0, 8.0, 12.0), (93.7, 8.0, 12.0)):
    sim = band_fraction(tv, l1, l2)[0]
    exact = frac_below_series(l2, tv) - frac_below_series(l1, tv)
    val_c_max = max(val_c_max, abs(sim - exact))
    print("  %10.1f %8.1f %8.1f %18.12f %18.12f %+14.3e" % (tv, l1, l2, sim, exact, sim - exact))
print()
print("  max |Simpson - exact series| over the bands the model uses : %.3e" % val_c_max)
if val_c_max < 1e-9:
    print("  Exact to better than one part in 10^9. VALIDATION C PASSES.")
else:
    print("  *** INTEGRATOR FAILS. ***")
print()
print("  The default grid is 401 points across the band, which is plenty for a")
print("  4 um interval. It is not plenty for a wide one, and pretending otherwise")
print("  would be the easy dishonesty here. The same integrator over 4 to 100 um")
print("  at 288 K, refined:")
print()
print("  %10s %20s %16s" % ("grid points", "Simpson", "error vs series"))
wide_exact = frac_below_series(100.0, 288.0) - frac_below_series(4.0, 288.0)
for npt in (401, 1601, 6401, 25601, 102401):
    sw = band_fraction(288.0, 4.0, 100.0, npts=npt)[0]
    print("  %10d %20.12f %+16.3e" % (npt, sw, sw - wide_exact))
print("  %10s %20.12f" % ("exact", wide_exact))
print()
print("  Fourth-order convergence, as Simpson's rule promises: each fourfold")
print("  refinement cuts the error by about 256. Nothing in this study integrates")
print("  a range wider than 5.5 um, so the 401-point grid stands.")
print()
print("  Window fraction (%.1f to %.1f um) against temperature:" % (WIN_LO_UM, WIN_HI_UM))
print()
print("  %10s %16s %16s" % ("T (K)", "window frac", "absorbing frac"))
for t in (210.0, 240.0, 255.0, 270.0, 288.0, 300.0, 320.0):
    f = band_fraction(t)[0]
    print("  %10.1f %16.6f %16.6f" % (t, f, 1 - f))
F_WIN_SURF = band_fraction(T_OBS)[0]
print()
print("  At the observed surface temperature the window carries %.4f of the" % F_WIN_SURF)
print("  emission, which is %.2f W/m^2 out of sigma*288^4 = %.2f W/m^2." %
      (F_WIN_SURF * SIGMA * T_OBS ** 4, SIGMA * T_OBS ** 4))
print("  Wien peak at 288 K: %.3f um, inside the window." % (2897.771 / T_OBS))
print("  That is the awkward fact the grey model cannot represent: the window")
print("  sits on top of the emission peak, not off in the tail.")


# ===========================================================================
rule("12. THE TWO-BAND MODEL, AND THE CONTRADICTION IT EXPOSES")
print("One layer, two spectral bands. The window (%.0f-%.0f um) has optical depth" %
      (WIN_LO_UM, WIN_HI_UM))
print("tau_w, the rest of the spectrum has tau_b. Emissivity in each band is")
print("1 - exp(-tau). Band fractions come from section 11 and move with T_s, so")
print("the surface temperature is found by fixed-point iteration.")
print()


def two_band_ts(tau_w, tau_b, te=T_E, lo_um=WIN_LO_UM, hi_um=WIN_HI_UM,
                tol=1e-12, itmax=200):
    """Return (T_s, f_window, A_eff, iterations). A_eff is the effective grey
    absorptivity the two-band atmosphere presents to the surface."""
    ew = 1.0 - np.exp(-tau_w)
    eb = 1.0 - np.exp(-tau_b)
    ts = te
    it = 0
    for it in range(1, itmax + 1):
        f = band_fraction(ts, lo_um, hi_um)[0]
        a_eff = ew * f + eb * (1.0 - f)
        ts_new = te * (1.0 - a_eff / 2.0) ** -0.25
        if abs(ts_new - ts) < tol:
            ts = ts_new
            break
        ts = ts_new
    f = band_fraction(ts, lo_um, hi_um)[0]
    a_eff = ew * f + eb * (1.0 - f)
    return ts, f, a_eff, it


def grey_ts(tau, te=T_E):
    return te * (1.0 - (1.0 - np.exp(-tau)) / 2.0) ** -0.25


tau_grey_star = -np.log(1.0 - eps_star)
print("  Grey calibration. eps* = %.6f corresponds to tau = -ln(1-eps*) = %.6f" %
      (eps_star, tau_grey_star))
print("  check: T_s(tau) = %.6f K against target %.2f K, diff %+.3e" %
      (grey_ts(tau_grey_star), T_OBS, grey_ts(tau_grey_star) - T_OBS))
print()
print("  Two-band calibration with the window held at tau_w = %.2f. Solve for" % TAU_W_PRESENT)
print("  tau_b such that T_s = %.2f K." % T_OBS)


def solve_tau_b_direct(tau_w, target=T_OBS, lo_um=WIN_LO_UM, hi_um=WIN_HI_UM, te=T_E):
    """Closed-form calibration. At the target surface temperature the band
    fractions are fixed, so the required absorptivity follows from one line of
    algebra and no root finding is needed. Returns (tau_b, feasible)."""
    f = band_fraction(target, lo_um, hi_um)[0]
    ew = 1.0 - np.exp(-tau_w)
    eps_target = 2.0 * (1.0 - (te / target) ** 4)
    eb = (eps_target - ew * f) / (1.0 - f)
    if eb >= 1.0:
        return float("inf"), False
    return -np.log(1.0 - eb), True


def solve_tau_b_bisect(tau_w, target=T_OBS, lo=1e-6, hi=60.0):
    """Independent confirmation of the line above, by bisection on T_s."""
    for _ in range(100):
        mid = 0.5 * (lo + hi)
        if two_band_ts(tau_w, mid)[0] < target:
            lo = mid
        else:
            hi = mid
    return 0.5 * (lo + hi)


f288 = band_fraction(T_OBS)[0]
ew_min = (eps_star - (1.0 - f288)) / f288
print("  First, a feasibility question the grey model never has to face. Even")
print("  with the absorbing band made completely opaque, the band only covers")
print("  %.4f of the emission. To reach 288 K the window must supply the rest:" % (1 - f288))
print("    required window emissivity at least : %.6f" % ew_min)
print("    i.e. window optical depth at least  : %.6f" % (-np.log(1 - ew_min)))
print("  A perfectly clear window makes 288 K unreachable with one layer,")
print("  whatever you do to the rest of the spectrum. The classroom picture of")
print("  'a window that lets radiation straight out' is quantitatively too")
print("  generous by itself.")
print()

tau_b_star, feasible = solve_tau_b_direct(TAU_W_PRESENT)
tau_b_bis = solve_tau_b_bisect(TAU_W_PRESENT)
print("  Calibration cross-check, closed form against bisection:")
print("    tau_b* closed form : %.10f" % tau_b_star)
print("    tau_b* bisection   : %.10f" % tau_b_bis)
print("    difference         : %+.3e" % (tau_b_star - tau_b_bis))
print()
ts_2b, f_2b, a_2b, its = two_band_ts(TAU_W_PRESENT, tau_b_star)
print("    tau_b*         = %.6f" % tau_b_star)
print("    band emissivity= %.6f" % (1 - np.exp(-tau_b_star)))
print("    window emissiv.= %.6f" % (1 - np.exp(-TAU_W_PRESENT)))
print("    window fraction= %.6f" % f_2b)
print("    A_eff          = %.6f   (compare grey eps* = %.6f)" % (a_2b, eps_star))
print("    T_s            = %.6f K (target %.2f, diff %+.3e, %d iterations)" %
      (ts_2b, T_OBS, ts_2b - T_OBS, its))
print()
print("  The two models agree exactly on the present-day state, by construction.")
print("  A_eff and eps* are the same number to %.1e. The equilibrium equations" %
      abs(a_2b - eps_star))
print("  are identical once you collapse the bands into one absorptivity.")
print()
print("  Now the first thing the grey model cannot tell you.")
print()
trans = (f_2b * np.exp(-TAU_W_PRESENT) + (1 - f_2b) * np.exp(-tau_b_star)) * SIGMA * ts_2b ** 4
trans_win = f_2b * np.exp(-TAU_W_PRESENT) * SIGMA * ts_2b ** 4
trans_grey = np.exp(-tau_grey_star) * SIGMA * T_OBS ** 4
print("  Surface emission reaching space without being absorbed:")
print("    grey model, tau = %.4f          : %8.3f W/m^2" % (tau_grey_star, trans_grey))
print("    two-band model, through window  : %8.3f W/m^2" % trans_win)
print("    two-band model, total           : %8.3f W/m^2" % trans)
print("    OBSERVED (Costa & Shine 2012)   : %8.3f W/m^2" % WINDOW_ESCAPE_OBS)
print()
print("    grey overshoots observation by  : %+8.3f W/m^2 (%.1f times too much)" %
      (trans_grey - WINDOW_ESCAPE_OBS, trans_grey / WINDOW_ESCAPE_OBS))
print("    two-band overshoots by          : %+8.3f W/m^2 (%.1f times too much)" %
      (trans - WINDOW_ESCAPE_OBS, trans / WINDOW_ESCAPE_OBS))
print()
print("  Both models leak far too much radiation straight to space. Turn it")
print("  around: what tau_w reproduces the measured %.0f W/m^2 of window escape," % WINDOW_ESCAPE_OBS)
print("  with the absorbing band taken as completely opaque?")
print()
tau_b_opaque = 30.0
lo, hi = 1e-6, 20.0
for _ in range(200):
    mid = 0.5 * (lo + hi)
    ts_t, f_t, a_t, _ = two_band_ts(mid, tau_b_opaque)
    esc = f_t * np.exp(-mid) * SIGMA * ts_t ** 4
    if esc > WINDOW_ESCAPE_OBS:
        lo = mid
    else:
        hi = mid
tau_w_fit = 0.5 * (lo + hi)
ts_fit, f_fit, a_fit, _ = two_band_ts(tau_w_fit, tau_b_opaque)
esc_fit = f_fit * np.exp(-tau_w_fit) * SIGMA * ts_fit ** 4
print("    tau_w needed   = %.6f" % tau_w_fit)
print("    window escape  = %.4f W/m^2 (target %.1f)" % (esc_fit, WINDOW_ESCAPE_OBS))
print("    resulting T_s  = %.4f K" % ts_fit)
print("    against 288 K  = %+.4f K" % (ts_fit - T_OBS))
print()
print("  THE DOUBLE BIND. A single-layer two-band model can match the observed")
print("  surface temperature or the observed window escape. Not both.")
print("    tuned to T_s = 288 K -> window escape %.1f W/m^2, %.1f times measured" %
      (trans, trans / WINDOW_ESCAPE_OBS))
print("    tuned to escape 22 W/m^2 -> T_s = %.1f K, %.1f K too hot" %
      (ts_fit, ts_fit - T_OBS))
print()
print("  One more thing, which is sharper than either number above. Rerun the")
print("  calibration at several window opacities. Each time, tau_b is re-solved so")
print("  that the surface still lands on 288 K.")
print()
print("  %10s %12s %16s %16s %12s" %
      ("tau_w", "tau_b*", "window escape", "total escape", "T_s"))
pin_rows = []
for tw in (0.30, 0.60, 1.00, 1.50, 2.00):
    tb_i, ok_i = solve_tau_b_direct(tw)
    if not ok_i:
        continue
    b_i = two_band_ts(tw, tb_i)
    s4 = SIGMA * b_i[0] ** 4
    w_i = b_i[1] * np.exp(-tw) * s4
    t_i = w_i + (1 - b_i[1]) * np.exp(-tb_i) * s4
    pin_rows.append((tw, tb_i, w_i, t_i, b_i[0]))
    print("  %10.2f %12.4f %16.3f %16.3f %12.3f" % (tw, tb_i, w_i, t_i, b_i[0]))
print()
print("  The window escape falls by a factor of %.1f across that range and the" %
      (pin_rows[0][2] / pin_rows[-1][2]))
print("  total escape does not move at all: %.3f W/m^2 in every row, to the last" % pin_rows[0][3])
print("  printed digit. It cannot move. Once the model is pinned to 288 K the")
print("  effective absorptivity is fixed at eps* = %.6f, so the untouched" % eps_star)
print("  fraction leaving the top is (1 - eps*) sigma T_s^4 = %.3f W/m^2 however" %
      ((1 - eps_star) * SIGMA * T_OBS ** 4))
print("  you divide the opacity between the bands. Closing the window forces the")
print("  absorbing band open by exactly the compensating amount. At tau_w = 1.50")
print("  the window escape alone reads %.3f W/m^2, which is the measured figure," % pin_rows[3][2])
print("  and the model is still leaking %.1f W/m^2 in total." % pin_rows[3][3])
print()
print("  The reason is physical and the model cannot fix it. One layer has one")
print("  temperature. The real atmosphere emits from many altitudes: the window")
print("  leaks from near the ground where it is warm, the band centres radiate")
print("  from high and cold. A single slab has to pick one, and whichever it")
print("  picks, the other quantity is off by a factor of four, or by twelve")
print("  kelvin. Adding layers helps; making them grey does not.")


# ===========================================================================
rule("13. HOW WRONG IS GREY? THE RESPONSE TO ADDED ABSORBER")
print("The present-day state is not where the grey assumption fails. Both models")
print("were forced through the same point. It fails in the derivative.")
print()
print("Experiment: multiply the optical depth by a factor and recompute T_s.")
print("Grey scales one tau; the two-band scales tau_w and tau_b together.")
print()
print("  %10s %14s %14s %14s %14s %10s" %
      ("multiplier", "grey T_s", "two-band T_s", "grey dT", "2-band dT", "ratio"))
mult_rows = []
for mult in (1.0, 1.25, 1.5, 2.0, 3.0, 4.0, 8.0):
    tg = grey_ts(tau_grey_star * mult)
    t2 = two_band_ts(TAU_W_PRESENT * mult, tau_b_star * mult)[0]
    dg = tg - T_OBS
    d2 = t2 - T_OBS
    r = dg / d2 if abs(d2) > 1e-9 else float("nan")
    mult_rows.append((mult, tg, t2, dg, d2, r))
    print("  %10.2f %14.4f %14.4f %+14.4f %+14.4f %10s" %
          (mult, tg, t2, dg, d2, ("%.3f" % r) if r == r else "--"))
print()
d_grey_2x = mult_rows[3][3]
d_2b_2x = mult_rows[3][4]
print("  Doubling the absorber:")
print("    grey model warms by     %+.3f K" % d_grey_2x)
print("    two-band model warms by %+.3f K" % d_2b_2x)
print("    grey overstates by a factor of %.2f" % (d_grey_2x / d_2b_2x))
print()
print("  Why: the grey model can always absorb more, because a single band with")
print("  tau = %.2f is only %.1f%% opaque. The two-band model's absorbing band is" %
      (tau_grey_star, 100 * (1 - np.exp(-tau_grey_star))))
print("  already %.2f%% opaque at present day, so doubling it buys almost" %
      (100 * (1 - np.exp(-tau_b_star))))
print("  nothing and the warming has to come from the window alone.")
print()
print("  A harder version, closer to what CO2 actually does. CO2 absorbs in the")
print("  band and barely touches the 8-12 um window, so scale tau_b only:")
print()
print("  %10s %16s %14s" % ("tau_b mult", "two-band T_s", "dT (K)"))
band_only = []
for mult in (1.0, 2.0, 4.0, 8.0, 100.0):
    t2 = two_band_ts(TAU_W_PRESENT, tau_b_star * mult)[0]
    band_only.append((mult, t2, t2 - T_OBS))
    print("  %10.1f %16.4f %+14.4f" % (mult, t2, t2 - T_OBS))
print()
print("  A hundredfold increase in band opacity is worth %.4f K. The band is" % band_only[-1][2])
print("  saturated. This is the single most important thing the grey model hides,")
print("  and it is part of why the real CO2 forcing grows as the logarithm of the")
print("  concentration rather than linearly.")
print()
print("  Saturation ceiling of the two-band model (tau_b -> infinity):")
ts_ceiling = two_band_ts(TAU_W_PRESENT, 500.0)[0]
print("    T_s = %.4f K, i.e. %+.4f K above present day, no matter how much" %
      (ts_ceiling, ts_ceiling - T_OBS))
print("    absorber you add to that band.")
print()
print("  Sensitivity to the window optical depth, which is the softest number in")
print("  the whole exercise:")
print()
print("  %10s %12s %14s %14s %12s" %
      ("tau_w", "tau_b*", "grey dT (2x)", "2-band dT (2x)", "ratio"))
scan_rows = []
for tw in (0.10, 0.20, 0.30, 0.40, 0.60, 0.90):
    tb, _ok = solve_tau_b_direct(tw)
    t2 = two_band_ts(tw * 2, tb * 2)[0]
    d2 = t2 - T_OBS
    scan_rows.append((tw, tb, d_grey_2x, d2, d_grey_2x / d2))
    print("  %10.2f %12.4f %+14.4f %+14.4f %12.3f" % (tw, tb, d_grey_2x, d2, d_grey_2x / d2))
print()
print("  The ratio stays above 1 everywhere in this range, so the direction of")
print("  the conclusion is safe. Its size is not: the ratio runs from %.2f to %.2f" %
      (min(r[4] for r in scan_rows), max(r[4] for r in scan_rows)))
print("  across these six values of tau_w. The safe reading is that grey")
print("  overstates the response; the factor of two is an estimate that carries")
print("  the window optical depth's uncertainty along with it.")


# ===========================================================================
rule("14. MONTE CARLO ON THE GREY-VERSUS-TWO-BAND RATIO")
print("The window edges are not sharp and its optical depth is not well known.")
print("Sampling them:")
print("  window low edge  ~ Uniform(7.5, 8.5) um")
print("  window high edge ~ Uniform(11.5, 13.0) um")
print("  tau_w            ~ Uniform(0.15, 0.60)")
print()
N_MC2 = 4000
rng2 = np.random.default_rng(child_b)
lo_draw = rng2.uniform(7.5, 8.5, N_MC2)
hi_draw = rng2.uniform(11.5, 13.0, N_MC2)
tw_draw = rng2.uniform(0.15, 0.60, N_MC2)
print("  trials: %d   stream: PCG64 spawned from SeedSequence(%d), child 1" % (N_MC2, SEED))
print("  Each trial recalibrates tau_b to hit 288 K, then doubles both depths.")
print()

ratio = np.empty(N_MC2)
d2b = np.empty(N_MC2)
tb_draw = np.empty(N_MC2)
n_infeasible = 0
for i in range(N_MC2):
    tw = tw_draw[i]
    lo_i, hi_i = lo_draw[i], hi_draw[i]
    tb, ok = solve_tau_b_direct(tw, lo_um=lo_i, hi_um=hi_i)
    if not ok:
        n_infeasible += 1
        d2b[i] = np.nan
        ratio[i] = np.nan
        tb_draw[i] = np.nan
        continue
    tb_draw[i] = tb
    t2 = two_band_ts(tw * 2, tb * 2, lo_um=lo_i, hi_um=hi_i, itmax=80)[0]
    d2b[i] = t2 - T_OBS
    ratio[i] = d_grey_2x / (t2 - T_OBS)

print("  trials where no tau_b can reach 288 K: %d of %d" % (n_infeasible, N_MC2))
if n_infeasible:
    keep = ~np.isnan(ratio)
    ratio = ratio[keep]
    d2b = d2b[keep]
    N_MC2 = int(keep.sum())
    print("  those are dropped; %d usable trials remain" % N_MC2)
print()
se_ratio = ratio.std(ddof=1) / np.sqrt(N_MC2)
se_d2b = d2b.std(ddof=1) / np.sqrt(N_MC2)
print("  %-30s %12s %12s %12s" % ("quantity", "mean", "sd", "SE"))
print("  %-30s %12.4f %12.4f %12.4f" % ("two-band dT for doubling (K)",
                                        d2b.mean(), d2b.std(ddof=1), se_d2b))
print("  %-30s %12.4f %12.4f %12.4f" % ("grey / two-band ratio",
                                        ratio.mean(), ratio.std(ddof=1), se_ratio))
print()
lo95, hi95 = np.percentile(ratio, [2.5, 97.5])
print("  95%% of trials give a ratio between %.3f and %.3f." % (lo95, hi95))
print("  Fraction of trials with ratio > 2 : %.4f" % float((ratio > 2).mean()))
print("  Fraction of trials with ratio > 1 : %.4f" % float((ratio > 1).mean()))
print()
print("  Convergence of the running mean of the ratio:")
print()
print("  %10s %14s %14s %14s" % ("trials", "running mean", "running SE", "dev from final"))
cs = np.cumsum(ratio)
cs2 = np.cumsum(ratio ** 2)
final_ratio = ratio.mean()
conv2_rows = []
for k in sorted(set([c for c in (50, 100, 250, 500, 1000, 2000, 3000) if c <= N_MC2] + [N_MC2])):
    m = cs[k - 1] / k
    v = max(cs2[k - 1] / k - m * m, 0.0) * k / (k - 1)
    se = np.sqrt(v / k)
    conv2_rows.append((k, m, se, m - final_ratio))
    print("  %10d %14.5f %14.5f %+14.5f" % (k, m, se, m - final_ratio))
print()
print("  The grey model overstates the response to doubled absorber by")
print("  %.2f +/- %.2f (SE). The 2.5th percentile is %.2f, so even the most" %
      (ratio.mean(), se_ratio, lo95))
print("  forgiving window we sampled still leaves grey overstating by a quarter.")
print()
print("  One caution about the %d dropped trials. They were dropped because a" % n_infeasible)
print("  wide, clear window makes 288 K unreachable, and those are exactly the")
print("  cases where the grey model would have looked worst. Dropping them pulls")
print("  the reported ratio down, so %.2f is a conservative figure rather than a" % ratio.mean())
print("  flattering one.")


# ===========================================================================
rule("15. WHERE A DIFFERENT MODELLING CHOICE WOULD CHANGE THE ANSWER")
print("Each row swaps one decision and reports the headline numbers.")
print()
print("  %-46s %10s %10s" % ("choice", "T_e (K)", "eps*"))
alts = [
    ("baseline: S=1361, albedo=0.293, T_obs=288", S0, ALBEDO, T_OBS),
    ("textbook albedo 0.30", S0, 0.30, T_OBS),
    ("albedo 0.29 (Stephens et al. 2012)", S0, 0.29, T_OBS),
    ("older S=1366 (pre-SORCE value)", 1366.0, ALBEDO, T_OBS),
    ("T_obs = 287.0 K (a cooler baseline period)", S0, ALBEDO, 287.0),
    ("T_obs = 289.0 K (a warmer baseline period)", S0, ALBEDO, 289.0),
    ("T_obs = 289.4 K (from measured LW, sect. 8)", S0, ALBEDO,
     (SURFACE_LW_OBS / SIGMA) ** 0.25),
]
alt_rows = []
for nm, s, a, tt in alts:
    te = eff_temp(s, a)
    e = 2.0 * (1.0 - (te / tt) ** 4)
    alt_rows.append((nm, te, e, (tt / te) ** 4 - 1.0))
    print("  %-46s %10.3f %10.4f" % (nm, te, e))
print()
eps_lo = min(r[2] for r in alt_rows)
eps_hi = max(r[2] for r in alt_rows)
print("  eps* spread across these choices: %.4f to %.4f, a range of %.4f," %
      (eps_lo, eps_hi, eps_hi - eps_lo))
print("  which is %.1f times the Monte Carlo standard deviation of %.4f." %
      ((eps_hi - eps_lo) / eps_draw.std(ddof=1), eps_draw.std(ddof=1)))
print("  The modelling choices matter more than the measurement uncertainty.")
print("  Anyone quoting eps* to three decimals without saying which albedo they")
print("  used is quoting noise.")
print()
print("  And the choice that matters most of all is not on this list: whether the")
print("  atmosphere is allowed to convect. Pure radiative equilibrium with a")
print("  realistic absorber gives a surface far hotter than 288 K and a lapse")
print("  rate steeper than any real air column can hold, which is what Manabe and")
print("  Strickler found in 1964; convection has to be added by hand before the")
print("  answer becomes Earth-like. None of that is in this file.")


# ===========================================================================
rule("16. FIGURE DATA")
print("Every number the article plots, printed so the figures can be checked.")
print()
print("[FIG1] surface temperature against layer count")
print("  N, T_s(K)")
for n in range(0, 11):
    print("  %d, %.4f" % (n, T_E * (n + 1.0) ** 0.25))
print("  T_e = %.4f ; T_obs = %.4f ; N* = %.6f" % (T_E, T_OBS, n_star))
print()
print("[FIG2] vertical temperature profile for several N")
for n in (1, 2, 3, 5, 10):
    xs, _, _ = solve_layers(n, F_ABS)
    ts_col = [(v / SIGMA) ** 0.25 for v in xs]
    print("  N=%d: %s" % (n, ", ".join("%.3f" % v for v in ts_col)))
print()
print("[FIG3] surface temperature against single-layer emissivity")
print("  eps, T_s(K)")
for k in range(0, 21):
    e = k / 20.0
    print("  %.2f, %.4f" % (e, ts_from_eps(e)))
print("  eps* = %.6f -> T_s = %.4f" % (eps_star, T_OBS))
print()
print("[FIG4] response to scaled optical depth, grey against two-band")
print("  multiplier, grey T_s, two-band T_s, band-only T_s")
for mult in (0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 2.5, 3.0, 4.0, 6.0, 8.0):
    tg = grey_ts(tau_grey_star * mult)
    t2 = two_band_ts(TAU_W_PRESENT * mult, tau_b_star * mult)[0]
    tbo = two_band_ts(TAU_W_PRESENT, tau_b_star * mult)[0]
    print("  %.2f, %.4f, %.4f, %.4f" % (mult, tg, t2, tbo))
print()
print("[FIG5] Monte Carlo convergence of eps*")
print("  trials, running mean, running SE")
for k, m, se, dev in conv_rows:
    print("  %d, %.6f, %.6f" % (k, m, se))
print()
print("[FIG5b] Monte Carlo convergence of the grey/two-band ratio")
print("  trials, running mean, running SE")
for k, m, se, dev in conv2_rows:
    print("  %d, %.5f, %.5f" % (k, m, se))
print()
print("[TABLE] the results table in the article")
print("  N, T_s analytic, T_s numeric, difference, T_s - T_obs, equivalent eps")
for n in range(0, 11):
    xs, _, _ = solve_layers(n, F_ABS)
    tn = (xs[0] / SIGMA) ** 0.25
    ta = T_E * (n + 1.0) ** 0.25
    e_equiv = 2.0 * (1.0 - (T_E / ta) ** 4)
    print("  %d, %.6f, %.6f, %+.3e, %+.4f, %.4f" %
          (n, ta, tn, tn - ta, ta - T_OBS, e_equiv))


# ===========================================================================
rule("17. HEADLINE NUMBERS")
print("  effective temperature T_e                    %10.3f K" % T_E)
print("  observed surface temperature                 %10.3f K" % T_OBS)
print("  greenhouse effect (this model's definition)  %10.3f K" % (T_OBS - T_E))
print("  layers needed, N*                            %10.4f" % n_star)
print("  T_s at N = 1                                 %10.3f K (%+.2f K)" %
      (T_E * 2 ** 0.25, T_E * 2 ** 0.25 - T_OBS))
print("  single-layer emissivity eps*                 %10.4f +/- %.4f (sd)" %
      (eps_star, eps_draw.std(ddof=1)))
print("  grey optical depth tau*                      %10.4f" % tau_grey_star)
print("  two-band tau_b* at tau_w = %.2f              %10.4f" % (TAU_W_PRESENT, tau_b_star))
print("  window fraction at 288 K                     %10.4f" % F_WIN_SURF)
print("  window escape, model                         %10.2f W/m^2" % trans)
print("  window escape, observed                      %10.2f W/m^2" % WINDOW_ESCAPE_OBS)
print("  T_s if escape is forced to 22 W/m^2          %10.2f K (%+.2f K)" %
      (ts_fit, ts_fit - T_OBS))
print("  doubling response, grey                      %10.3f K" % d_grey_2x)
print("  doubling response, two-band                  %10.3f K" % d_2b_2x)
print("  grey overstatement factor                    %10.2f +/- %.2f (SE)" %
      (ratio.mean(), se_ratio))
print("  band-only 100x response                      %10.4f K" % band_only[-1][2])
print()
print("  Runtime: %.2f s" % (time.time() - T_START))
print("  Seed: %d. Rerunning this file reproduces every digit above." % SEED)