"""Figures (PNG, light surface). One y-axis per panel; categorical colours in fixed order."""
import os, json
import numpy as np, pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from common import DATA, FIGS, load_json

C = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100", "#e87ba4", "#008300", "#6250d6", "#e34948"]
SURF, INK, INK2, GRID = "#fcfcfb", "#0b0b0b", "#52514e", "#e4e3de"
plt.rcParams.update({"figure.facecolor": SURF, "axes.facecolor": SURF, "savefig.facecolor": SURF,
                     "axes.edgecolor": "#b9b8b2", "axes.labelcolor": INK2, "xtick.color": INK2, "ytick.color": INK2,
                     "text.color": INK, "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.7,
                     "axes.spines.top": False, "axes.spines.right": False, "font.size": 10,
                     "axes.titlesize": 11, "axes.titleweight": "bold", "axes.titlelocation": "left",
                     "legend.frameon": False, "lines.linewidth": 2.0})


def sci(v):
    """Two significant figures: plain below 1e-2 as 'a.b e-n', otherwise decimal."""
    if v >= 0.01:
        return f"{v:.2g}"
    m, e = f"{v:.1e}".split("e")
    return f"{m}e{int(e)}"


def grp(v):
    """Counts rounded to two significant figures with thousands separators."""
    r = float(f"{v:.2g}")
    return f"{r:,.0f}"


def fig1():
    cov = pd.read_csv(os.path.join(DATA, "macrostrat_coverage_1myr.csv"))
    b = pd.read_csv(os.path.join(DATA, "macrostrat_bins.csv"))
    fig, ax = plt.subplots(2, 1, figsize=(8.6, 8.2))
    a = ax[0]
    x = cov["age_ma"]
    def sm(y, k=5):
        return pd.Series(y).rolling(k, center=True, min_periods=1).mean()
    for col, lab, c in (("north_america_sed_farea", "North America, all sedimentary", C[0]),
                        ("north_america_nonmarine_farea", "North America, non-marine", C[1]),
                        ("north_america_marine_farea", "North America, marine", C[2]),
                        ("south_america_sed_farea", "South America, all sedimentary (58 columns)", C[3])):
        a.plot(x, sm(cov[col]).clip(lower=1e-4), color=c, label=lab)
    a.set_xscale("log"); a.set_yscale("log"); a.set_xlim(0.5, 3500); a.set_ylim(3e-4, 1.2)
    a.set_xlabel("Age (Ma, log scale)"); a.set_ylabel("Fraction of column area holding\nsedimentary rock of that age")
    a.set_title("(a) How much of a continent still holds sediment of a given age (Macrostrat)")
    a.legend(loc="lower left", fontsize=9)
    a = ax[1]
    bb = b.iloc[1:].reset_index(drop=True)          # skip the Holocene bin
    xs = np.arange(len(bb))
    for off, pre, lab, c in ((-0.14, "na_all", "North America", C[0]), (0.14, "ds_all", "Deep-sea drill sites", C[1])):
        ok = bb[f"{pre}_yr_per_cm_med"].notna()
        a.errorbar(xs[ok] + off, bb.loc[ok, f"{pre}_yr_per_cm_med"],
                   yerr=[bb.loc[ok, f"{pre}_yr_per_cm_med"] - bb.loc[ok, f"{pre}_yr_per_cm_p16"],
                         bb.loc[ok, f"{pre}_yr_per_cm_p84"] - bb.loc[ok, f"{pre}_yr_per_cm_med"]],
                   fmt="o", ms=6, color=c, ecolor=c, elinewidth=1.2, capsize=0, label=lab + ": time-weighted median, 16-84%")
        ok2 = bb[f"{pre}_yr_per_cm_med_unweighted"].notna()
        a.plot(xs[ok2] + off, bb.loc[ok2, f"{pre}_yr_per_cm_med_unweighted"], "o", ms=6, mfc=SURF, mec=c, mew=1.6,
               label=lab + ": plain median of units")
    a.axhline(300, color=INK2, lw=1, ls="--")
    a.text(1.01, 300, "a 300-yr\nepisode\n= 1 cm", transform=a.get_yaxis_transform(), ha="left", va="center", fontsize=9, color=INK2)
    a.set_yscale("log"); a.set_xticks(xs); a.set_xticklabels([s.replace(" (to 3.5 Ga)", "") for s in bb["bin"]], rotation=40, ha="right")
    a.set_ylabel("Years represented by one centimetre\n(unit duration / thickness)")
    a.set_title("(b) Long-term time resolution of the rock that survives")
    a.set_ylim(30, 4e6); a.legend(loc="upper left", fontsize=8.5, ncol=2, columnspacing=1.2, handletextpad=0.3)
    fig.tight_layout(); fig.savefig(os.path.join(FIGS, "fig1_rock_record.png"), dpi=150); plt.close(fig)


def fig2():
    dd = pd.read_csv(os.path.join(DATA, "pbdb_density_1myr.csv"))
    sl = pd.read_csv(os.path.join(DATA, "pbdb_psampled_slices.csv"))
    dm = pd.read_csv(os.path.join(DATA, "dmsed_bins.csv")).iloc[:12]
    fig, ax = plt.subplots(2, 1, figsize=(8.6, 8.0), sharex=True)
    a = ax[0]
    s = dd[dd["age_ma"] < 540]
    a.plot(s["age_ma"], s["colls_per_myr_all"].rolling(5, center=True, min_periods=1).mean(), color=C[0], label="Fossil collections, all environments (PBDB)")
    a.plot(s["age_ma"], s["colls_per_myr_terrestrial"].rolling(5, center=True, min_periods=1).mean().clip(lower=0.05), color=C[1], label="Fossil collections, terrestrial (PBDB)")
    xs = np.repeat(np.r_[dm["age_top_ma"].values, dm["age_base_ma"].values[-1]], 2)[1:-1]
    a.plot(xs, np.repeat(dm["samples_per_myr"].values, 2), color=C[2], label="Geochemical samples in one compilation (DM-SED; a lower bound)")
    cr = pd.read_csv(os.path.join(DATA, "cenogrid_resolution.csv"))
    a.plot(cr["age_ma"], cr["n_samples"].rolling(3, center=True, min_periods=1).mean(), color=C[3],
           label="Benthic isotope samples in one splice (CENOGRID)")
    a.set_yscale("log"); a.set_ylim(0.05, 3e5); a.set_ylabel("Records per million years\nof geological time")
    a.set_title("(a) Sampling effort per million years of Earth history"); a.legend(loc="upper right", fontsize=9)
    a = ax[1]
    s = sl[sl["age_ma"] >= 3]
    for col, lab, c in (("Pterr_A10000_tau0.01", "10,000 km$^2$", C[0]), ("Pterr_A100_tau0.01", "100 km$^2$", C[1]), ("Pterr_A1_tau0.01", "1 km$^2$", C[2])):
        y = s[col].rolling(5, center=True, min_periods=1).mean().clip(lower=1e-12)
        a.plot(s["age_ma"], y, color=c, label=lab)
        a.text(8, float(y.iloc[2]) * 1.6, lab, color=INK, fontsize=9, va="bottom")
    a.set_yscale("log"); a.set_ylim(1e-11, 1e-2); a.set_xlim(0, 540)
    a.set_xlabel("Age (Ma)"); a.set_ylabel("P(a terrestrial fossil collection lies inside\nthe footprint and the time window)")
    a.set_title("(b) Chance that palaeontologists have sampled a 10-kyr installation, by footprint")
    a.legend(loc="upper right", fontsize=9, title="Footprint")
    fig.tight_layout(); fig.savefig(os.path.join(FIGS, "fig2_sampling_effort.png"), dpi=150); plt.close(fig)


def fig3():
    d = pd.read_csv(os.path.join(DATA, "cenogrid_d13c_series.csv"))
    sc = pd.read_csv(os.path.join(DATA, "cenogrid_scan.csv"))
    r5 = pd.read_csv(os.path.join(DATA, "cenogrid_recovery.csv")); r3 = pd.read_csv(os.path.join(DATA, "cenogrid_recovery_k3.csv"))
    rs = pd.read_csv(os.path.join(DATA, "cenogrid_step_recovery.csv"))
    es = pd.read_csv(os.path.join(DATA, "cenogrid_step_events.csv"))
    e3 = pd.read_csv(os.path.join(DATA, "cenogrid_events_k3.csv"))
    n_step5 = int((es["criterion_sigma"] == 5).sum())
    fig = plt.figure(figsize=(8.6, 9.8)); gs = fig.add_gridspec(3, 1, height_ratios=[1, 1.05, 1.2])
    a = fig.add_subplot(gs[0])
    a.plot(d["age_ma"], d["d13c"], color="#8c8b85", lw=0.5)
    a.set_ylabel("Benthic $\\delta^{13}$C (permil)"); a.set_xlim(0, 67)
    a.set_title("(a) CENOGRID deep-sea carbon-isotope record, 23,939 samples (Westerhold et al. 2020)")
    a.tick_params(labelbottom=False)
    b = fig.add_subplot(gs[1], sharex=a)
    b.plot(sc["onset_age_ma"], sc["z"], color=C[0], lw=0.8)
    b.axhline(-5, color=INK2, lw=1, ls="--"); b.axhline(-3, color=INK2, lw=1, ls=":")
    b.text(0.6, -5.25, "5 sigma", ha="left", va="top", fontsize=9, color=INK2)
    b.text(0.6, -3.25, "3 sigma", ha="left", va="top", fontsize=9, color=INK2)
    # labels sit in the clear band below the 5-sigma line, joined to their dots by thin leaders
    place = {10.8: (10.8, "center"), 34.2: (34.2, "center"), 61.4: (60.4, "right"), 62.1: (62.6, "left")}
    for _, e in e3.iterrows():
        zz = max(e["z"], -9.2)
        b.plot([e["onset_age_ma"]], [zz], "o", ms=6, color=C[1], mec=SURF, mew=1.2, zorder=4)
        if e["z"] < -9:
            b.annotate(f"{e['onset_age_ma']:.1f} Ma (PETM, z = {e['z']:.0f}, off scale)", (e["onset_age_ma"], zz),
                       textcoords="offset points", xytext=(-8, 4), ha="right", va="center", fontsize=8.5, color=INK)
            continue
        xt, ha = place.get(round(e["onset_age_ma"], 1), (e["onset_age_ma"], "center"))
        b.annotate(f"{e['onset_age_ma']:.1f} Ma", (e["onset_age_ma"], zz), xytext=(xt, -6.6), ha=ha, va="center",
                   fontsize=8.5, color=INK, arrowprops=dict(arrowstyle="-", color=INK2, lw=0.7, shrinkA=1, shrinkB=3))
    b.set_ylim(-10, 6); b.set_xlabel("Age (Ma)"); b.set_ylabel("Matched-filter amplitude\n(local sigma units)")
    b.set_title("(b) Scan for an industrial-style pulse: abrupt drop, 100-kyr recovery")
    c = fig.add_subplot(gs[2])
    M = [200, 300, 500, 700, 1000, 1500, 2000, 3000, 5000, 10000]
    MS_ = [500, 700, 1000, 1500, 2000, 2500, 3000, 4000, 5000, 10000]
    row5 = r5[r5["interval"] == "Cenozoic (0.3-66 Ma)"].iloc[0]; row3 = r3[r3["interval"] == "Cenozoic (0.3-66 Ma)"].iloc[0]
    rowS = rs[(rs["interval"] == "Cenozoic (0.3-66 Ma)") & (rs["criterion_sigma"] == 5)].iloc[0]
    c.plot(M, [row5[f"P_flag_{m}GtC"] for m in M], "-o", ms=5, color=C[0],
           label="Matched filter, 5 sigma (natural events flagged: 1, the PETM)")
    c.plot(MS_, [rowS[f"P_flag_{m}GtC"] for m in MS_], "-s", ms=5, color=C[2],
           label=f"Step statistic, 5 sigma (natural events flagged: {n_step5})")
    c.plot(M, [row3[f"P_flag_{m}GtC"] for m in M], "-o", ms=5, color=C[1],
           label="Matched filter, 3 sigma (natural events flagged: 5)")
    for xm, lab, ha, dxf in ((500, "fossil carbon\nemitted to date", "right", 0.96), (5000, "all conventional\nfossil fuel", "left", 1.04)):
        c.axvline(xm, color=INK2, lw=1, ls="--"); c.text(xm * dxf, 0.42, lab, fontsize=9, color=INK2, va="center", ha=ha)
    c.set_xscale("log"); c.set_xlabel("Carbon released in a pulse (GtC)"); c.set_ylabel("Fraction of injected pulses flagged")
    c.set_title("(c) Injection-recovery with two detectors: what the best Cenozoic record would flag")
    c.legend(loc="upper left", fontsize=8.5, bbox_to_anchor=(0.0, -0.2), ncol=1); c.set_ylim(-0.03, 1.08)
    fig.tight_layout(); fig.savefig(os.path.join(FIGS, "fig3_cenogrid_injection.png"), dpi=150, bbox_inches="tight"); plt.close(fig)


def fig4():
    inst = pd.read_csv(os.path.join(DATA, "synthesis_installations.csv"))
    carb = pd.read_csv(os.path.join(DATA, "synthesis_carbon.csv"))
    excl = pd.read_csv(os.path.join(DATA, "synthesis_exclusion.csv"))
    eras = ["0-3 Ma", "3-66 Ma", "66-252 Ma", "252-539 Ma", "539-2400 Ma"]
    scen = ((1.0, 1e-2, "1 km$^2$\n10 kyr"), (100.0, 1e-2, "100 km$^2$\n10 kyr"), (1e4, 1e-2, "10,000 km$^2$\n10 kyr"),
            (1e6, 3e-4, "10$^6$ km$^2$ dispersed\n300 yr"))
    fig, ax = plt.subplots(2, 2, figsize=(12.4, 11.0))

    # (a) split of P(known), last 500 Myr: assumed surface part versus data-based deep-time part
    a = ax[0, 0]
    xs = np.arange(len(scen))
    for off, part, lab, c in ((-0.13, "surface", "Surface part: remains still at or near today's land surface (assumed)", C[1]),
                              (0.13, "deep", "Deep-time part: PBDB probe x E x recognition (data x assumptions)", C[0])):
        med = []; lo = []; hi = []
        for A, tau, _ in scen:
            r = inst[(inst["area_km2"] == A) & np.isclose(inst["occupation_myr"], tau)].iloc[0]
            med.append(r[f"P_{part}_part_last_500Myr_median"]); lo.append(r[f"P_{part}_part_last_500Myr_p05"]); hi.append(r[f"P_{part}_part_last_500Myr_p95"])
        med, lo, hi = map(np.array, (med, lo, hi))
        a.errorbar(xs + off, med, yerr=[med - lo, hi - med], fmt="o", ms=7, color=c, ecolor=c, elinewidth=1.4, capsize=0, label=lab)
        for x, m in zip(xs + off, med):
            a.text(x + (0.07 if off > 0 else -0.07), m, sci(m), ha="left" if off > 0 else "right", va="center", fontsize=8.5, color=INK)
    a.set_yscale("log"); a.set_ylim(3e-9, 2); a.set_xlim(-0.6, len(scen) - 0.4)
    a.set_xticks(xs); a.set_xticklabels([s_[2] for s_ in scen])
    a.set_ylabel("Contribution to P(known today),\nrandom time in the last 500 Myr")
    a.set_title("(a) Where the probability comes from\n(dot: median; bar: 5-95% of the assumed priors)")
    a.legend(loc="upper left", bbox_to_anchor=(0.0, -0.16), fontsize=8.5, ncol=1)

    # (b) by era, central assumptions; dashed = buried strata only
    a = ax[0, 1]
    xe = np.arange(len(eras))
    for (A, tau, lab), c in zip(((1e6, 3e-4, "10$^6$ km$^2$ dispersed built footprint, 300 yr"), (1e4, 1e-2, "10,000 km$^2$, 10 kyr"),
                                 (100.0, 1e-2, "100 km$^2$, 10 kyr"), (1.0, 1e-2, "1 km$^2$, 10 kyr")), C):
        r = inst[(inst["area_km2"] == A) & np.isclose(inst["occupation_myr"], tau)].iloc[0]
        y = np.clip([r[f"P_known_central_{e}"] for e in eras], 1e-11, 1)
        yd = np.clip([r[f"P_deep_only_central_{e}"] for e in eras], 1e-11, 1)
        a.plot(xe, y, "-o", ms=6, color=c, label=lab)
        a.plot(xe[:2], yd[:2], "--", lw=1.3, color=c)
        a.plot(xe[:1], yd[:1], "o", ms=6, mfc=SURF, mec=c, mew=1.5)
    a.plot([], [], "--o", color=INK2, lw=1.3, ms=6, mfc=SURF, mec=INK2, mew=1.5, label="hollow, dashed: buried strata only (no surface term)")
    a.set_yscale("log"); a.set_ylim(1e-11, 2); a.set_xticks(xe); a.set_xticklabels(eras, rotation=20)
    a.set_ylabel("Probability it would be known today\n(central assumptions)"); a.set_xlabel("When it happened")
    a.set_title("(b) Surface installations by age\n")
    a.legend(loc="upper left", bbox_to_anchor=(0.0, -0.24), fontsize=8.5, ncol=1)

    # (c) carbon pulse by era
    a = ax[1, 0]
    FLOOR = 1e-5
    cols = ["P_event_0-66Ma", "P_event_66-252Ma_median", "P_event_252-539Ma_median", "P_event_539-2400Ma_median"]
    xb = np.arange(len(cols))
    a.axvspan(0.5, 3.4, color="#efeee9", zorder=0)
    a.text(1.95, 9.0, "assumed coverage and noise factors\n(no injection test on real data)", ha="center", va="center", fontsize=8.5, color=INK2)
    a.text(0.0, 9.0, "measured\n(injection test)", ha="center", va="center", fontsize=8.5, color=INK2)
    for (M, lab), c in zip(((5000.0, "5,000 GtC"), (2000.0, "2,000 GtC"), (500.0, "500 GtC (fossil carbon emitted to date)")), C):
        for crit, ls, lw in (("5sigma_singular_outlier", "-", 2.0), ("3sigma_among_largest_excursions", "--", 1.3)):
            r = carb[(carb["total_GtC"] == M) & (carb["release_duration_yr"] == 300.0) & (carb["criterion"] == crit)].iloc[0]
            raw = np.array([r[k] for k in cols], float); y = np.clip(raw, FLOOR, 1)
            a.plot(xb, y, ls, lw=lw, color=c, label=lab if ls == "-" else None)
            a.scatter(xb, y, s=36, zorder=3, edgecolors=c, linewidths=1.5,
                      facecolors=[c if v >= FLOOR else SURF for v in raw])
    a.plot([], [], "-", color=INK2, lw=2.0, label="solid: singular outlier (5 sigma)")
    a.plot([], [], "--", color=INK2, lw=1.3, label="dashed: among the largest (3 sigma)")
    a.plot([], [], "o", ms=6, mfc=SURF, mec=INK2, mew=1.5, label="hollow: below 1e-5")
    a.set_yscale("log"); a.set_ylim(FLOOR * 0.6, 40); a.set_xlim(-0.4, 3.4)
    a.set_yticks([1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1])
    a.set_xticks(xb); a.set_xticklabels(["0-66 Ma", "66-252 Ma", "252-539 Ma", "539-2400 Ma"], rotation=20)
    a.set_xlabel("When it happened"); a.set_ylabel("Probability it is in today's catalogue of events")
    a.set_title("(c) Fossil-carbon pulse released in 300 yr: catalogued as an event\n(attribution to industry not included)")
    a.legend(loc="upper left", bbox_to_anchor=(0.0, -0.24), fontsize=8.5, ncol=2, columnspacing=1.2, handlelength=1.8)

    # (d) how many past episodes the null result allows (95%), last 2.4 Gyr
    a = ax[1, 1]
    e = excl[excl["window"] == "last_2400Myr"]
    nm = []; nv = []; lo = []; hi = []
    for A, tau, _ in scen:
        r = e[(e["area_km2"] == A) & np.isclose(e["occupation_myr"], tau)].iloc[0]
        nm.append(r["N95_marginalised"]); nv.append(r["N95_naive_3_over_median"]); lo.append(r["N95_at_p95"]); hi.append(r["N95_at_p05"])
    nm, nv, lo, hi = map(np.array, (nm, nv, lo, hi))
    a.axhspan(11000, 59000, color="#efeee9", zorder=0)
    a.axhline(21000, color=INK2, lw=1, ls="--",
              label="Star systems passing within 1 pc in 2.33 Gyr: 21,000\n(band: 11,000-59,000; sibling fly-by study)")
    a.vlines(xs, lo, hi, color=C[0], lw=1.4)
    a.plot(xs, nm, "o", ms=8, color=C[0], label="Uncertainty in P marginalised (this revision)")
    a.plot(xs, nv, "o", ms=7, mfc=SURF, mec=C[1], mew=1.6, label="3 / median P (first pass; too strong)")
    for x, v in zip(xs, nm):
        a.text(x + 0.1, v, grp(v), ha="left", va="center", fontsize=8.5, color=INK, bbox=dict(fc=SURF, ec="none", pad=1.0, alpha=0.8))
    a.set_yscale("log"); a.set_ylim(20, 3e6); a.set_xlim(-0.6, len(scen) - 0.4)
    a.set_xticks(xs); a.set_xticklabels([s_[2] for s_ in scen])
    a.set_ylabel("Largest number of past episodes compatible\nwith no known trace (95%), last 2.4 Gyr")
    a.set_title("(d) What the null result excludes\n(bar: 3/P at the 95th to 5th percentile of P)")
    a.legend(loc="upper left", bbox_to_anchor=(0.0, -0.16), fontsize=8.5, ncol=1)
    fig.tight_layout(); fig.savefig(os.path.join(FIGS, "fig4_p_known_by_age.png"), dpi=150, bbox_inches="tight"); plt.close(fig)


def fig5():
    rn = pd.read_csv(os.path.join(DATA, "signals_radionuclides.csv"))
    fu = pd.read_csv(os.path.join(DATA, "synthesis_funnel.csv"))
    fu = fu[(fu["area_km2"] == 1e4) & np.isclose(fu["occupation_myr"], 1e-2)].reset_index(drop=True)
    fig, ax = plt.subplots(1, 2, figsize=(11.0, 4.8))
    a = ax[0]
    ys = np.arange(len(rn))
    a.hlines(ys, rn["window_lo_yr"], rn["window_hi_yr"], color=C[0], lw=6, alpha=0.35)
    a.plot(rn["window_yr"], ys, "o", ms=8, color=C[0])
    for y, (nuc, wv) in enumerate(zip(rn["nuclide"], rn["window_yr"])):
        lab = f"{wv/1e6:.0f} Myr" if wv >= 1e6 else f"{wv/1e3:.0f} kyr"
        a.text(wv, y + 0.22, lab, ha="center", va="bottom", fontsize=9, color=INK, bbox=dict(fc=SURF, ec="none", pad=1.5))
    a.set_yticks(ys); a.set_yticklabels(["$^{240}$Pu", "$^{239}$Pu", "$^{129}$I", "$^{236}$U", "$^{244}$Pu"]); a.set_xscale("log"); a.set_xlim(3e4, 3e9)
    for xv, lab in ((6.6e7, "66 Myr"), (5e8, "500 Myr")):
        a.axvline(xv, color=INK2, lw=1, ls="--"); a.text(xv * 0.9, len(rn) - 0.1, lab, ha="right", va="center", fontsize=9, color=INK2)
    a.set_ylim(-0.6, len(rn) + 0.2)
    a.set_xlabel("Years for which a 20th-century-scale inventory stays measurable")
    a.set_title("(a) Nuclear fingerprints: measurability windows\n(dot: central floor; bar: floor x/÷ 10)")
    a = ax[1]
    xs = np.arange(len(fu)); wd = 0.26
    series = (("f_window_strata_survive", "strata of that 20-kyr window survive\n(fraction of footprint)", C[0]),
              ("f_horizon_is_mapped_bedrock_surface", "that horizon is today's mapped bedrock\n(fraction of footprint)", C[1]),
              ("P_pbdb_collection_inside_footprint_and_window", "a fossil collection samples it\n(probability)", C[2]))
    for k, (col, lab, c) in enumerate(series):
        a.bar(xs + (k - 1) * wd, np.clip(fu[col], 1e-9, 1), width=wd - 0.03, color=c, label=lab, bottom=0)
    a.set_yscale("log"); a.set_ylim(1e-8, 1); a.set_xticks(xs); a.set_xticklabels(fu["era"], rotation=20)
    a.set_xlabel("When it happened")
    a.set_title("(b) The funnel for a 10,000 km$^2$ installation occupied 10 kyr\n(first two bars are model outputs; the third is measured)")
    a.legend(loc="upper right", fontsize=8.5)
    fig.tight_layout(); fig.savefig(os.path.join(FIGS, "fig5_fingerprints_and_funnel.png"), dpi=150); plt.close(fig)


def main():
    fig1(); fig2(); fig3(); fig4(); fig5()
    print("figures written:", sorted(os.listdir(FIGS)))


if __name__ == "__main__":
    main()