"""Figures (static PNG).  Styling: shared plotstyle.py (validated reference palette, fixed categorical
order; one y-axis per panel; legends plus selective direct labels)."""
import os

import matplotlib.ticker
import numpy as np
import pandas as pd

import plotstyle as ps
from common import DATA_DIR, FIG_DIR, AGE_SUN
from abundances import ism_factor
import population as pop
import thermal

plt = ps.apply()

F_COL = {"F1_classical": ps.BLUE, "F2h_hybrid": ps.ORANGE, "F2_flat": ps.AQUA, "F3_stagnant": ps.YELLOW}
F_SHORT = {"F1_classical": "F1 classical, mobile lid", "F2h_hybrid": "F2h hybrid, mobile lid",
           "F2_flat": "F2 flat, mobile lid", "F3_stagnant": "F3 stagnant lid"}


def _line(ax, x, a, b, color, x0=AGE_SUN, lw=2.0, ls="-"):
    ax.plot(x, a + b * (x - x0), color=color, lw=lw, ls=ls)


def fig1_abundances(S, summary, emp):
    tw, dm, bb = S["solar_twins"], S["harps_gto"], S["battistini_bensby"]
    fig, axes = plt.subplots(1, 3, figsize=(15.5, 4.9))
    xs = np.linspace(0, 13, 50)

    ax = axes[0]
    m = tw["use_th"]
    ax.errorbar(tw.loc[m, "age"], tw.loc[m, "logf_th_direct"], yerr=tw.loc[m, "e_logf_th_direct"], fmt="o", ms=5,
                color=ps.BLUE, ecolor=ps.AXIS, elinewidth=0.8, mec=ps.SURFACE, mew=0.8,
                label=f"Th measured, corrected for decay in the star\n(Botelho+ 2019), N = {int(m.sum())}")
    m = tw["use_eu"]
    ax.errorbar(tw.loc[m, "age"], tw.loc[m, "logf_eu_rock"], yerr=tw.loc[m, "e_logf_eu_rock"], fmt="s", ms=4.5,
                color=ps.ORANGE, ecolor=ps.AXIS, elinewidth=0.8, mec=ps.SURFACE, mew=0.8,
                label=f"Eu measured (Spina+ 2018), N = {int(m.sum())}")
    _line(ax, xs[xs <= 9.2], summary["fit_th_age"]["a"], summary["fit_th_age"]["b"], ps.BLUE)
    _line(ax, xs[xs <= 9.2], summary["fits_eu_raw"]["solar_twins"]["a"], summary["fits_eu_raw"]["solar_twins"]["b"], ps.ORANGE)
    ax.plot([AGE_SUN], [0.0], marker="*", ms=15, color=ps.INK, mec=ps.SURFACE, mew=0.8, ls="none", label="Sun (definition)")
    ax.set_title("Solar twins (|[Fe/H]| < 0.15)")
    ax.set_xlabel("stellar age today (Gyr)")
    ax.set_ylabel("[Th/rock] at birth or [Eu/rock] (dex), rock = (Mg + Si)/2")
    ax.set_xlim(0, 13)
    ax.set_ylim(-0.3, 0.45)
    ax.legend(loc="upper right", fontsize=8.5)

    ax = axes[1]
    for t, col, mk, lab, key in ((dm, ps.AQUA, "o", "Eu measured (Delgado Mena+ 2017)", "harps_gto"),
                                 (bb, ps.MAGENTA, "s", "Eu measured (Battistini & Bensby 2016)", "battistini_bensby")):
        m = t["use_eu"]
        ax.plot(t.loc[m, "age"], t.loc[m, "logf_eu_rock"], mk, ms=4, color=col, mec=ps.SURFACE, mew=0.5, alpha=0.9,
                label=f"{lab}, N = {int(m.sum())}")
        _line(ax, xs, summary["fits_eu_raw"][key]["a"], summary["fits_eu_raw"][key]["b"], col)
    xc = np.clip(xs, 0, 12)
    raw = emp["a"] + emp["b"] * (xc - AGE_SUN)
    ax.plot(xs, raw, color=ps.INK, lw=2.2, ls="-", label="adopted trend, data only")
    ax.fill_between(xs, raw - emp["s"], raw + emp["s"], color=ps.NEUTRAL, alpha=0.18, lw=0, label="adopted intrinsic scatter (1 sigma)")
    for sc, ls, lab in ((pop.ISM_SCALE_BASE, "--", f"initial Th: data + decay model at strength {pop.ISM_SCALE_BASE:g} (baseline)"),
                        (1.0, ":", "initial Th: data + full one-zone decay model")):
        ax.plot(xs, raw + sc * ism_factor("Th232", xc), color=ps.INK, lw=1.6, ls=ls, label=lab)
    ax.plot([AGE_SUN], [0.0], marker="*", ms=15, color=ps.INK, mec=ps.SURFACE, mew=0.8, ls="none")
    ax.set_title("Broad samples (-1 < [Fe/H] < 0.55)")
    ax.set_xlabel("stellar age today (Gyr)")
    ax.set_ylabel("[Eu/rock] (dex); dashed and dotted lines add a model")
    ax.set_xlim(0, 13.5)
    ax.set_ylim(-0.3, 0.68)
    ax.legend(loc="upper left", fontsize=7.8)

    ax = axes[2]
    gk = pd.read_csv(os.path.join(DATA_DIR, "galah_k_binned_MgSi.csv"))
    for key, col, lab in (("galah_dr3", ps.VIOLET, "GALAH DR3"), ("galah_dr4", ps.GREEN, "GALAH DR4")):
        g = gk[gk["sample"] == key]
        xm = 0.5 * (g["age_lo"] + g["age_hi"])
        ax.fill_between(xm, g["p16"], g["p84"], color=col, alpha=0.15, lw=0)
        ax.plot(xm, g["median"], "o-", color=col, ms=5, mec=ps.SURFACE, mew=0.8,
                label=f"{lab}: median and 16-84% (N = {summary['fits_k'][key]['n_pass_cuts']} pass the cuts)")
    ax.axhline(0.0, color=ps.MUTED, lw=1.0)
    ax.set_title("Potassium (stable K), GALAH")
    ax.set_xlabel("stellar age today (Gyr)")
    ax.set_ylabel("[K/rock] (dex)")
    ax.set_xlim(0, 13.5)
    ax.set_ylim(-0.3, 0.45)
    ax.legend(loc="upper left", fontsize=8.5)
    fig.suptitle("Long-lived heat sources per unit rock (Mg, Si) in Sun-like stars, by formation epoch.  Points and solid lines are data; "
                 "dashed and dotted lines include the decay model", x=0.01, ha="left", fontsize=11.5, fontweight="bold")
    fig.tight_layout(rect=(0, 0, 1, 0.95))
    fig.savefig(os.path.join(FIG_DIR, "fig1_heat_sources_vs_age.png"))
    plt.close(fig)


def fig2_earth_histories():
    hist = pd.read_csv(os.path.join(DATA_DIR, "earth_twin_histories.csv"))
    cal = pd.read_csv(os.path.join(DATA_DIR, "earth_calibration.csv")).set_index("formulation")
    fig, axes = plt.subplots(1, 2, figsize=(12.5, 4.8))
    ax = axes[0]
    ax.axhspan(thermal.T_SOL0_LOW, thermal.T_SOL0_HIGH, color=ps.NEUTRAL, alpha=0.2, lw=0)
    ax.text(0.6, thermal.T_SOL0_LOW - 12, "band: dry solidus at the surface\n(Katz et al. 2003 to Foley & Smye 2018)",
            ha="left", va="top", fontsize=8.5, color=ps.INK2)
    for name, col in F_COL.items():
        h = hist[hist["formulation"] == name]
        ax.plot(h["t_gyr"], h["T_p_K"], color=col, label=F_SHORT[name])
    ax.plot([4.5], [thermal.T_REF], marker="*", ms=14, color=ps.INK, mec=ps.SURFACE, mew=0.8, ls="none", label="Earth today")
    ax.set_xlabel("planet age (Gyr)")
    ax.set_ylabel("mantle potential temperature (K)")
    ax.set_xlim(0, 30)
    ax.set_ylim(1150, 2150)
    ax.set_title("Thermal history of an Earth twin")
    ax.legend(loc="upper right", fontsize=9)

    ax = axes[1]
    ax.axhline(0.1, color=ps.INK2, lw=1.2, ls="--")
    ax.text(29.6, 0.115, "threshold used: 10% of modern Earth", ha="right", va="bottom", fontsize=8.5, color=ps.INK2)
    ax.axhline(1.0, color=ps.MUTED, lw=1.0, ls=":")
    ax.text(29.6, 1.15, "modern Earth", ha="right", va="bottom", fontsize=8.5, color=ps.INK2)
    for name, col in F_COL.items():
        h = hist[(hist["formulation"] == name) & (hist["degassing_rel"] > 0)]
        ax.plot(h["t_gyr"], h["degassing_rel"], color=col)
        life = cal.loc[name, "lifetime_Gyr"]
        ax.plot([life], [0.1], "o", ms=7, color=col, mec=ps.SURFACE, mew=1.0)
        off, ha = {"F1_classical": ((6, 6), "left"), "F2h_hybrid": ((7, -14), "left"), "F2_flat": ((-6, 7), "right"),
                   "F3_stagnant": ((-6, -14), "right")}[name]
        ax.annotate(f"{life:.1f} Gyr", (life, 0.1), xytext=off, textcoords="offset points", fontsize=8.5, color=ps.INK2, ha=ha)
    ax.set_yscale("log")
    ax.set_ylim(0.003, 3000)
    ax.set_xlim(0, 30)
    ax.set_xlabel("planet age (Gyr)")
    ax.set_ylabel("CO2 degassing per unit area / modern Earth")
    ax.set_title("Degassing of an Earth twin (ridge law for mobile lids)")
    fig.tight_layout()
    fig.savefig(os.path.join(FIG_DIR, "fig2_earth_twin_histories.png"))
    plt.close(fig)


def fig3_survival():
    cur = pd.read_csv(os.path.join(DATA_DIR, "survival_curves_baseline.csv"))
    hyp = pd.read_csv(os.path.join(DATA_DIR, "hyper_draws.csv"))
    hz = pd.read_csv(os.path.join(DATA_DIR, "stellar_hz_lifetimes.csv")).set_index("mass_msun")
    ages = [5, 8, 10, 15, 20, 30, 45]
    fig, axes = plt.subplots(1, 3, figsize=(17.5, 4.9), sharey=True)
    for ax, crit, title in ((axes[0], "degassing>=10%", "Degassing at least 10% of modern Earth, ridge law"),
                            (axes[1], "platespeed>=10%", "Degassing at least 10%, plate-speed law"),
                            (axes[2], "any_melting", "Any decompression melting (dry solidus)")):
        for name, col in F_COL.items():
            crit_n = "degassing>=10%" if (name == "F3_stagnant" and crit == "platespeed>=10%") else crit
            c = cur[(cur["formulation"] == name) & (cur["criterion"] == crit_n) & np.isclose(cur["stellar_age_today_gyr"], AGE_SUN)]
            ax.plot(c["t_gyr"], c["P_hab"], color=col, label=F_SHORT[name])
            g = hyp[(hyp["formulation"] == name) & (hyp["criterion"] == crit_n)]
            lo = [np.percentile(g[f"P_hab_{a}Gyr"], 16) for a in ages]
            hi = [np.percentile(g[f"P_hab_{a}Gyr"], 84) for a in ages]
            dx = {"F1_classical": 0.45, "F2h_hybrid": 0.15, "F2_flat": -0.15, "F3_stagnant": -0.45}[name]
            ax.vlines(np.array(ages) + dx, lo, hi, color=col, lw=3, alpha=0.45)
        for mass, ls in ((0.8, ":"), (0.7, "--")):
            x = hz.loc[mass, "hz_exit_median_gyr"]
            ax.axvline(x, color=ps.MUTED, lw=1.1, ls=ls)
            ax.text(x + 0.4, 0.99, f"{mass:g} Msun star\nleaves its\nhabitable zone\n(median orbit)", ha="left", va="top",
                    fontsize=8, color=ps.INK2)
        ax.set_xlim(0, 46)
        ax.set_ylim(0, 1.02)
        ax.set_xlabel("planet age (Gyr)")
        ax.set_title(title, fontsize=10.5)
    axes[0].set_ylabel("P(still geologically habitable)")
    axes[0].legend(loc="center right", fontsize=9)
    fig.suptitle("Probability that a rocky planet (0.5-2 Earth masses, formed with the Sun) is still geologically habitable, for three "
                 f"end-of-habitability rules.  Lines: baseline; bars: 16-84% over {hyp['draw'].nunique()} hyper-parameter draws",
                 x=0.01, ha="left", fontsize=11, fontweight="bold")
    fig.tight_layout(rect=(0, 0, 1, 0.94))
    fig.savefig(os.path.join(FIG_DIR, "fig3_survival_curves.png"))
    plt.close(fig)


def fig4_controls(emp):
    sens = pd.read_csv(os.path.join(DATA_DIR, "lifetime_vs_th_budget.csv"))
    var = pd.read_csv(os.path.join(DATA_DIR, "lifetime_variance_decomposition.csv"))
    pct = pd.read_csv(os.path.join(DATA_DIR, "lifetime_percentiles_baseline.csv"))
    fig, axes = plt.subplots(1, 3, figsize=(15, 4.8))
    ax = axes[0]
    lo, hi = 10 ** (emp["a"] - 2 * emp["s"]), 10 ** (emp["a"] + 2 * emp["s"])
    ax.axvspan(lo, hi, color=ps.NEUTRAL, alpha=0.2, lw=0)
    ax.text(np.sqrt(lo * hi), 0.6, "measured star-to-star\nrange (2 sigma)", ha="center", va="bottom", fontsize=8.5, color=ps.INK2)
    for name, col in F_COL.items():
        s = sens[sens["formulation"] == name]
        ax.plot(s["f_Th_U"], s["lifetime_degassing"], "o-", color=col, ms=5, mec=ps.SURFACE, mew=0.8, label=F_SHORT[name])
        ax.plot(s["f_Th_U"], s["lifetime_melting"], "--", color=col, lw=1.4)
    ax.plot([], [], "--", color=ps.INK2, lw=1.4, label="dashed: any-melting criterion")
    ax.set_xscale("log")
    ax.set_xticks([0.5, 0.7, 1, 1.4, 2, 3])
    ax.set_xticklabels(["0.5", "0.7", "1", "1.4", "2", "3"])
    ax.xaxis.set_minor_locator(matplotlib.ticker.NullLocator())
    ax.set_xlim(0.48, 3.1)
    ax.set_ylim(0, 28)
    ax.set_xlabel("initial Th and U per unit rock / Earth")
    ax.set_ylabel("habitable lifetime of an Earth-mass planet (Gyr)")
    ax.set_title("Lifetime versus Th and U budget")
    ax.legend(loc="upper left", fontsize=8.5)

    ax = axes[1]
    v = var[(var["criterion"] == "degassing>=10%") & (var["input_varied"] != "all together") &
            ~var["input_varied"].str.startswith("SHIFT")]
    labels = list(dict.fromkeys(v["input_varied"]))
    y0 = np.arange(len(labels))[::-1]
    for j, (name, col) in enumerate(F_COL.items()):
        g = v[v["formulation"] == name].set_index("input_varied").loc[labels]
        ax.barh(y0 + (1.5 - j) * 0.19, g["width_5_95_gyr"], height=0.17, color=col, label=F_SHORT[name])
    ax.set_yticks(y0)
    ax.set_yticklabels([l.replace(" (", "\n(") for l in labels], fontsize=8.5)
    ax.set_xlabel("5-95% width of lifetime when only this input varies (Gyr)")
    ax.set_title("What spreads the lifetimes (degassing criterion)")
    ax.grid(axis="y", visible=False)

    ax = axes[2]
    p = pct[pct["criterion"] == "degassing>=10%"]
    sc = pd.read_csv(os.path.join(DATA_DIR, "scenarios.csv"))
    sc = sc[sc["criterion"] == "degassing-type"]
    for name, col in F_COL.items():
        g = p[p["formulation"] == name].sort_values("stellar_age_today_gyr")
        band = []
        for key in ("decay-model strength 0", "decay-model strength 1"):
            r = sc[(sc["formulation"] == name) & sc["scenario"].str.startswith(key)].iloc[0]
            band.append(np.array([float(x) for x in r["median_lifetime_by_epoch"].split("/")]))
        ax.fill_between(pop.NODE_AGES, band[0], band[1], color=col, alpha=0.18, lw=0)
        ax.plot(pop.NODE_AGES, band[0], color=col, lw=1.2, ls=":")
        ax.plot(g["stellar_age_today_gyr"], g["median"], "o-", color=col, ms=5, mec=ps.SURFACE, mew=0.8, label=F_SHORT[name])
    xx = np.linspace(0, 12.5, 20)
    ax.plot(xx, xx, color=ps.MUTED, lw=1.2, ls="--")
    ax.text(0.4, 21.3, "points: baseline (decay-model strength 0.5)\nbands: strength 0 (dotted edge, data only) to 1\n"
            "grey dashed: lifetime = age today", fontsize=8.5, color=ps.INK2, ha="left", va="top")
    ax.set_xlim(0, 12.5)
    ax.set_ylim(0, 22)
    ax.set_xlabel("stellar age today (Gyr); older stars formed earlier")
    ax.set_ylabel("median habitable lifetime (Gyr)")
    ax.set_title("Lifetime versus formation epoch")
    fig.tight_layout()
    fig.savefig(os.path.join(FIG_DIR, "fig4_what_controls_lifetime.png"))
    plt.close(fig)


LADDER = [
    ("baseline (0.5-2 Earth masses", "baseline: 0.5-2 Earth masses, ridge law,\n10% threshold, decay strength 0.5"),
    ("planet mass 0.5-3", "planet mass 0.5-3 Earth masses"),
    ("planet mass 0.5-5", "planet mass 0.5-5 Earth masses"),
    ("5% of planets on the melting", "5% of planets on the melting criterion"),
    ("10% of planets on the melting", "10% of planets on the melting criterion"),
    ("threshold per planet log-uniform", "threshold per planet, log-uniform 3-100%"),
    ("plate-speed degassing law", "plate-speed degassing law"),
    ("decay-model strength 0", "heat budget from data only (strength 0)"),
    ("decay-model strength 1", "full one-zone decay model (strength 1)"),
    ("no thermostat", "no thermostat (first-pass behaviour)"),
    ("thermostat 2100 K", "thermostat 2100 K"),
    ("COMBINED highest", "combined, highest: 0.5-2, strength 0,\nthermostat 2100 K"),
    ("COMBINED: mass 0.5-5 + plate-speed", "combined: 0.5-5 + plate-speed law"),
    ("COMBINED low:", "combined, low: 0.5-5 + 10% on melting"),
    ("COMBINED lowest", "combined, lowest: 0.5-5 + 10% on melting\n+ strength 1"),
    ("MELT", "all planets on the melting criterion"),
]
F_MARK = {"F1_classical": "o", "F2h_hybrid": "s", "F2_flat": "D"}


def fig5_rank():
    rk = pd.read_csv(os.path.join(DATA_DIR, "rank_baseline.csv"))
    ref = pd.read_csv(os.path.join(DATA_DIR, "rank_reference_fixed_caps.csv")).set_index("case")
    hyp = pd.read_csv(os.path.join(DATA_DIR, "hyper_draws.csv"))
    sc = pd.read_csv(os.path.join(DATA_DIR, "scenarios.csv"))
    ns = [1, 2, 3, 4, 6, 8]
    fig = plt.figure(figsize=(14.5, 11.6))
    gs = fig.add_gridspec(2, 2, height_ratios=[1.0, 1.25], hspace=0.28, wspace=0.08)
    top = [fig.add_subplot(gs[0, 0]), None]
    top[1] = fig.add_subplot(gs[0, 1], sharey=top[0])
    for ax, crit, title in ((top[0], "degassing>=10%", "Baseline: degassing at least 10% of modern Earth (ridge law)"),
                            (top[1], "any_melting", "Baseline: any decompression melting (dry solidus)")):
        for case, ls, lab in (("fixed cap 10 Gyr, all stars", ":", "fixed 10 Gyr cap"),
                              ("fixed cap 20 Gyr, all stars", "-.", "fixed 20 Gyr cap"),
                              ("no cap, hosts > 0.6 Msun (K dwarfs and up)", "--", "no cap, K dwarfs and hotter")):
            ax.plot(ns, [ref.loc[case, f"rank_n{n}"] for n in ns], color=ps.INK2, lw=1.4, ls=ls, label=lab)
        for name, col in F_COL.items():
            g = rk[(rk["formulation"] == name) & (rk["criterion"] == crit) & (rk["hosts"] == "all stars") &
                   (rk["lid_mix"] == "all planets this regime")].iloc[0]
            ax.plot(ns, [g[f"rank_n{n}"] for n in ns], "o-", color=col, ms=6, mec=ps.SURFACE, mew=0.9, label=F_SHORT[name])
            h = hyp[(hyp["formulation"] == name) & (hyp["criterion"] == crit)]
            for n in (2, 4, 6):
                lo, hi = np.percentile(h[f"rank_n{n}"], [16, 84])
                dx = {"F1_classical": 0.09, "F2h_hybrid": 0.03, "F2_flat": -0.03, "F3_stagnant": -0.09}[name]
                ax.vlines(n + dx, lo, hi, color=col, lw=3, alpha=0.45)
        ax.axhline(0.1, color=ps.MUTED, lw=1.0)
        ax.axhline(0.01, color=ps.MUTED, lw=1.0)
        ax.set_yscale("log")
        ax.set_ylim(3e-5, 1.5)
        ax.set_xticks(ns)
        ax.set_xlabel("number of hard steps n")
        ax.set_title(title, fontsize=10.5)
    top[0].set_ylabel("humanity's birth rank\n(fraction of civilisations born earlier)")
    top[0].legend(loc="lower left", fontsize=8.5, ncol=2)
    plt.setp(top[1].get_yticklabels(), visible=False)

    # scenario ladder: the same rank when one assumption about the long-lived tail is changed
    y = np.arange(len(LADDER))[::-1]
    bot = [fig.add_subplot(gs[1, 0]), None]
    bot[1] = fig.add_subplot(gs[1, 1], sharey=bot[0])
    for ax, n in ((bot[0], 4), (bot[1], 6)):
        for j, (key, _lab) in enumerate(LADDER):
            for name in ("F1_classical", "F2h_hybrid", "F2_flat"):
                if key == "MELT":
                    r = sc[(sc["formulation"] == name) & (sc["criterion"] == "any_melting") & sc["scenario"].str.startswith("baseline")]
                else:
                    r = sc[(sc["formulation"] == name) & (sc["criterion"] == "degassing-type") & sc["scenario"].str.startswith(key)]
                dy = {"F1_classical": 0.2, "F2h_hybrid": 0.0, "F2_flat": -0.2}[name]
                ax.plot([r.iloc[0][f"rank_n{n}"]], [y[j] + dy], F_MARK[name], ms=8, color=F_COL[name], mec=ps.SURFACE, mew=1.0,
                        label=F_SHORT[name] if j == 0 else None)
        for x in (0.1, 0.01):
            ax.axvline(x, color=ps.MUTED, lw=1.0)
        nocap = ref.loc["no cap, hosts > 0.6 Msun (K dwarfs and up)", f"rank_n{n}"]
        ax.axvline(nocap, color=ps.INK2, lw=1.4, ls="--")
        ax.text(nocap * 1.12, y[0] + 0.45, "no planetary cap,\nK dwarfs and hotter", fontsize=8, color=ps.INK2, va="top", ha="left")
        for yy in y[:-1] - 0.5:
            ax.axhline(yy, color=ps.GRID, lw=0.8)
        ax.set_xscale("log")
        ax.set_xlim(1.5e-4, 1.5)
        ax.set_ylim(-0.6, len(LADDER) - 0.4)
        ax.grid(axis="y", visible=False)
        ax.set_xlabel(f"humanity's birth rank for n = {n} hard steps")
        ax.set_title(f"One assumption changed at a time, n = {n}", fontsize=10.5)
    bot[0].set_yticks(y)
    bot[0].set_yticklabels([lab for _k, lab in LADDER], fontsize=8.8)
    plt.setp(bot[1].get_yticklabels(), visible=False)
    bot[0].legend(loc="center left", fontsize=8.5)
    fig.suptitle("Humanity's birth rank with geophysical habitable-lifetime distributions (all stars, Madau-Dickinson star formation).\n"
                 f"Top: baseline by number of hard steps; bars: 16-84% over {hyp['draw'].nunique()} hyper-parameter draws.  "
                 "Bottom: dependence on assumptions about the long-lived tail",
                 x=0.01, ha="left", fontsize=11, fontweight="bold")
    fig.savefig(os.path.join(FIG_DIR, "fig5_birth_rank.png"))
    plt.close(fig)


def make_all(S, summary, emp):
    fig1_abundances(S, summary, emp)
    fig2_earth_histories()
    fig3_survival()
    fig4_controls(emp)
    fig5_rank()