"""Figures of the impact_window study (called by run_all.py with the in-memory results). Colour use (palette of the timing study): the three size-distribution families always take the first three categorical hues (blue, orange, aqua); chronologies take the seven hues in fixed order; the HSE treatment is carried by line style or marker fill (solid / filled = upper bound only, dashed / open = both bounds), never by colour; sterilisation thresholds are an ordered blue ramp. Text is always ink-coloured. """ from __future__ import annotations import os import numpy as np import chronology as chron import config as cfg import model cfg.style() import matplotlib.pyplot as plt # noqa: E402 from matplotlib.lines import Line2D # noqa: E402 CH_COL = dict(zip(chron.ALL, cfg.SERIES)) # seven categorical hues in fixed order CH_SHORT = {"neukum1983": "Neukum 1983", "marchi2009": "Marchi+ 2009", "yue2026": "Yue+ 2026 (CE-6 refit)", "nesvorny2023": "Nesvorny+ 2023", "sawtooth2012": "sawtooth, Nectaris at 4.1 Ga", "anchor_ce6": "anchors, SPA 4.247 Ga (CE-6)", "anchor_joy433": "anchors, SPA 4.33 Ga (Joy+)"} SFD_COL = {"mainbelt": cfg.C["blue"], "lunar_pl300": cfg.C["orange"], "lunar_pl600": cfg.C["aqua"]} SFD_SHORT = {"mainbelt": "main-belt shape", "lunar_pl300": "shallow power law (basins >= 300 km)", "lunar_pl600": "steep power law (basins >= 600 km)"} HSE_LS = {"upper": "-", "two": (0, (4, 2))} HSE_SHORT = {"upper": "HSE mass as upper bound only", "two": "HSE mass, both bounds"} THR_COL = {740.: "#86b6ef", 440.: "#2a78d6", 190.: "#104281"} # ordinal blue ramp, light -> dark EV = [("LUCA", 4.20), ("4.1", 4.10), ("3.95", 3.95), ("3.7", 3.70), ("3.48", 3.48)] def _dj(res, D): return int(np.argmin(np.abs(res["D"] - D))) def _out(fig, name): path = os.path.join(cfg.FIGS, name) fig.savefig(path, bbox_inches="tight") plt.close(fig) return path # --------------------------------------------------------------------------------------------- def fig_chronologies(res): """Fig. 1: the same crater densities under different chronologies.""" bt = res["basin_table"] bt = bt[bt.dated] nk = chron.MODELS["neukum1983"] T = np.linspace(3.7, 4.51, 400) fig, ax = plt.subplots(1, 2, figsize=(11.5, 4.7)) a = ax[0] for k in chron.FUNCTIONAL: a.plot(T, chron.MODELS[k].n1(T), color=CH_COL[k], label=CH_SHORT[k]) a.plot(T, chron.MODELS["sawtooth_aspub"].n1(T), color=CH_COL["sawtooth2012"], lw=1.2, ls=(0, (4, 2)), label="sawtooth as published (not re-normalised)") a.errorbar([cfg.SPA_CE6[0]], [cfg.SPA_N1_YUE[0]], yerr=[cfg.SPA_N1_YUE[1]], fmt="o", color=cfg.INK, ms=5, capsize=2, zorder=5) a.annotate("SPA, Chang'e-6\n4.247 Ga (Su+ 2025;\nN(1) from Yue+ 2026)", (cfg.SPA_CE6[0], cfg.SPA_N1_YUE[0]), xytext=(4.27, 0.045), fontsize=8, color=cfg.INK2, arrowprops=dict(arrowstyle="-", color=cfg.INK2, lw=0.7)) for lab, age in (("Imbrium", 3.87), ("Nectaris", 4.17), ("SPA", 4.31)): y = float(nk.n1(age)) a.axhline(y, color=cfg.INK2, lw=0.7, ls=":") if lab == "Imbrium": a.text(4.5, y * 1.08, f"{lab} crater density (Orgel+ 2018)", fontsize=7.5, color=cfg.INK2, ha="right") else: a.text(3.705, y * 1.08, f"{lab} crater density (Orgel+ 2018)", fontsize=7.5, color=cfg.INK2) a.set_yscale("log") a.set_ylim(0.012, 7) a.set_xlim(3.7, 4.51) a.set_xlabel("age (Ga)") a.set_ylabel("craters >= 1 km per km$^2$ (Neukum-equivalent)") a.set_title("a Five published chronology curves") a.legend(fontsize=7.3, loc="upper left", bbox_to_anchor=(0.0, 1.0)) b = ax[1] Tg = np.linspace(3.7, 4.51, 812) for k in chron.ALL: ages = bt[f"age_{k}"].to_numpy() n = (ages[None, :] <= Tg[:, None]).sum(axis=1) b.step(Tg, n, where="post", color=CH_COL[k], lw=1.6, label=CH_SHORT[k]) b.axvspan(cfg.LUCA_PAPER[1], cfg.LUCA_PAPER[2], color=cfg.GRID, alpha=0.8, lw=0) b.text(4.21, 0.4, "LUCA 95% (Moody+ 2024)", fontsize=7.5, color=cfg.INK2, ha="center") b.axvline(cfg.T0_ZIRCON, color=cfg.INK2, lw=0.9, ls="--") b.text(cfg.T0_ZIRCON + 0.006, 3.0, "oldest zircon\n4.404 Ga", fontsize=8, color=cfg.INK2) b.set_xlim(3.7, 4.51) b.set_ylim(0, 31) b.set_xlabel("age (Ga)") b.set_ylabel("dated lunar basins younger than this age (of 30)") b.set_title("b The same 30 basins, re-dated") b.legend(fontsize=7.5, loc="upper left") fig.suptitle("The lunar basin record is fixed; its calendar is not", x=0.01, ha="left", fontsize=13, fontweight="bold") fig.tight_layout() return _out(fig, "fig1_basin_chronologies.png") def fig_sfd(res): """Fig. 2: the size distribution above the largest lunar basin and its one terrestrial check.""" B = res["basins"] grid = res["grid"] a_sc, b_sc = model.scaling_fit() d = np.sort(B.d_imp_km.to_numpy())[::-1] n = np.arange(1, len(d) + 1) mb = model.MainBeltSFD() dd = np.logspace(np.log10(11), np.log10(3000), 300) d300, d600 = a_sc * 300 ** b_sc, a_sc * 600 ** b_sc n300, n600 = int((B.D_main_km >= 300).sum()), int((B.D_main_km >= 600).sum()) fig, axs = plt.subplots(1, 2, figsize=(12.6, 5.2), gridspec_kw=dict(width_ratios=[1.25, 1])) ax = axs[0] ax.step(np.concatenate([d, [d[-1]]]), np.concatenate([n, [n[-1]]]), where="pre", color=cfg.INK, lw=1.8, label=f"lunar basins (Neumann+ 2015, {len(B)}), sizes from basin diameters") y = n300 * mb.N(dd) / mb.N(d300) ax.plot(dd[dd <= mb.dmax], y[dd <= mb.dmax], color=SFD_COL["mainbelt"], label="main-belt shape, ends at Ceres (939 km)") for fam, d0, n0 in (("lunar_pl300", d300, n300), ("lunar_pl600", d600, n600)): q = grid[("yue2026", fam)]["q_lunar"] x = dd[dd >= d0] lo, md, hi = (n0 * (x / d0) ** (-np.quantile(q, p)) for p in (0.84, 0.5, 0.16)) ax.fill_between(x, lo, hi, color=SFD_COL[fam], alpha=0.18, lw=0) ax.plot(x, md, color=SFD_COL[fam], label=SFD_SHORT[fam] + f", slope {np.median(q):.2f}") for lab, D, y0, ha in (("190 km\nphotic zone", 190, 90, "left"), ("440 km\nocean (1-D)", 440, 90, "left"), ("740 km\nocean (3-D)", 740, 22, "right"), ("1000 km", 1000, 22, "left"), ("2000 km\nsurface melts", 2000, 6, "right")): ax.axvline(D, color=cfg.INK2, lw=0.7, ls=":") ax.text(D * (1.03 if ha == "left" else 0.97), y0, lab, fontsize=8, color=cfg.INK2, va="top", ha=ha) ax.annotate("SPA impactor: the largest\nin the lunar record", (d[0], 1.0), xytext=(40, 0.17), fontsize=8, color=cfg.INK2, arrowprops=dict(arrowstyle="-", color=cfg.INK2, lw=0.7)) mg = res["summary"]["megabasins_ge1000km"] ax.annotate(f"basins >= 1000 km: {mg['observed']} observed,\n{mg['mainbelt_predicted']:.1f} predicted by the main-belt shape", (a_sc * 1000 ** b_sc, mg["mainbelt_predicted"]), xytext=(12.5, 0.55), fontsize=8, color=cfg.INK2, arrowprops=dict(arrowstyle="-", color=cfg.INK2, lw=0.7)) ax.axhline(1.0 / res["ratio_stats"]["median"], color=cfg.INK2, lw=0.7, ls="--") ax.text(12, 1.0 / res["ratio_stats"]["median"] * 1.15, "below the dashed line: fewer than one such impactor expected\n" f"on Earth while the lunar record formed (Earth/Moon = {res['ratio_stats']['median']:.0f})", fontsize=7.5, color=cfg.INK2, va="bottom") ax.set_xscale("log") ax.set_yscale("log") ax.set_xlim(11, 3000) ax.set_ylim(0.004, 120) ax.set_xticks([20, 50, 100, 200, 500, 1000, 2000]) ax.set_xticklabels(["20", "50", "100", "200", "500", "1000", "2000"]) ax.set_yticks([0.01, 0.1, 1, 10, 100]) ax.set_yticklabels(["0.01", "0.1", "1", "10", "100"]) ax.set_xlabel("impactor diameter d (km)") ax.set_ylabel("number of lunar basins from impactors >= d") ax.set_title("a Every sterilisation threshold lies above the largest lunar basin") ax.legend(fontsize=7.6, loc="lower left", bbox_to_anchor=(0.0, 0.0)) # realised late-accreted mass by family against the HSE range b = axs[1] edges = np.logspace(20, 24.3, 44) ctr = np.sqrt(edges[1:] * edges[:-1]) b.axvspan(cfg.HSE_MASS[0], cfg.HSE_MASS[1], color=cfg.GRID, lw=0) b.text(np.sqrt(cfg.HSE_MASS[0] * cfg.HSE_MASS[1]), 0.985, "HSE late veneer\n(0.7-3.0) x 10$^{22}$ kg", fontsize=8, color=cfg.INK2, ha="center", va="top", transform=b.get_xaxis_transform()) hs = res["hse"].set_index("sfd") for fam in model.SFD_FAMILIES: m = np.concatenate([r["mass"] for (c, s), r in grid.items() if s == fam]) h, _ = np.histogram(m, bins=edges) b.step(ctr, h / h.sum() / np.diff(np.log10(edges)), where="mid", color=SFD_COL[fam], lw=1.8, label=f"{SFD_SHORT[fam].split(' (')[0]}: {100 * hs.loc[fam, 'frac_below_lower_0p7e22']:.0f} / " f"{100 * hs.loc[fam, 'frac_within_0p7_3e22']:.0f} / " f"{100 * hs.loc[fam, 'frac_above_upper_3e22']:.0f}") b.set_xscale("log") b.set_xlim(1e20, 2e24) b.set_ylim(0, 1.6) b.set_xlabel("mass delivered to Earth by the modelled impactors since 4.51 Ga (kg)") b.set_ylabel("probability per decade of mass") b.set_title("b Most realisations deliver less than the late veneer") b.legend(fontsize=7.6, loc="upper left", title="% of realisations below / inside /\nabove the HSE range", title_fontsize=7.6, alignment="left") fig.suptitle("The tail of the size distribution is assumed, and the late-veneer mass is its only terrestrial check", x=0.01, ha="left", fontsize=12.5, fontweight="bold") fig.tight_layout() return _out(fig, "fig2_size_distribution.png") def fig_survival(res): """Fig. 3: probability that no impactor above a size struck Earth after a given age.""" G = res["G"] Tg = model.T_GRID sel = Tg >= 3.3 fig, axs = plt.subplots(2, 2, figsize=(11.8, 8.6)) panels = [(190., "a >= 190 km: photic zone boiled"), (440., "b >= 440 km: ocean vaporised (Sleep+ 1989)"), (740., "c >= 740 km: ocean vaporised (Citron & Stewart 2022)")] def decorate(ax, title): for lab, t in EV: ax.axvline(t, color=cfg.INK2, lw=0.6, ls=":") ax.text(t, 1.03, lab, fontsize=7.5, color=cfg.INK2, ha="center") ax.axvline(cfg.T0_ZIRCON, color=cfg.INK2, lw=0.9, ls="--") ax.text(cfg.T0_ZIRCON, 1.03, "water", fontsize=7.5, color=cfg.INK2, ha="center") ax.set_xlim(4.51, 3.3) ax.set_ylim(0, 1.0) ax.set_xlabel("age T (Ga); time runs to the right") ax.set_ylabel("P(no such impact after T)") ax.set_title(title, pad=16) for ax, (D, title) in zip(axs.ravel()[:3], panels): for s in model.SFD_FAMILIES: for cond in cfg.HSE_MAIN: r = G[cond][f"sfd:{s}"] ax.plot(Tg[sel], r["surv"][_dj(r, D)][sel], color=SFD_COL[s], lw=1.3, ls=HSE_LS[cond]) for cond in cfg.HSE_MAIN: r = G[cond]["ALL"] ax.plot(Tg[sel], r["surv"][_dj(r, D)][sel], color=cfg.INK, lw=2.4, ls=HSE_LS[cond]) decorate(ax, title) handles = [Line2D([], [], color=cfg.INK, lw=2.4, label="all models")] handles += [Line2D([], [], color=SFD_COL[s], lw=1.6, label=SFD_SHORT[s]) for s in model.SFD_FAMILIES] handles += [Line2D([], [], color=cfg.INK2, lw=1.4, ls=HSE_LS[c], label=HSE_SHORT[c]) for c in cfg.HSE_MAIN] axs[0, 0].legend(handles=handles, fontsize=7.3, loc="lower right", bbox_to_anchor=(1.0, 0.02)) ax = axs[1, 1] for k in chron.ALL: r = G["upper"][f"chron:{k}"] ax.plot(Tg[sel], r["surv"][_dj(r, 440.)][sel], color=CH_COL[k], lw=1.4, label=CH_SHORT[k]) for cond in cfg.HSE_MAIN: r = G[cond]["ALL"] ax.plot(Tg[sel], r["surv"][_dj(r, 440.)][sel], color=cfg.INK, lw=2.4, ls=HSE_LS[cond], label="all models, " + ("upper bound only" if cond == "upper" else "both bounds")) decorate(ax, "d >= 440 km by chronology (HSE mass as upper bound only)") ax.legend(fontsize=7.3, loc="lower right", bbox_to_anchor=(1.0, 0.02)) fig.suptitle("When was the last large impact? The size-distribution tail and the late-veneer treatment " "matter more than the chronology", x=0.01, ha="left", fontsize=12.5, fontweight="bold") fig.tight_layout() return _out(fig, "fig3_last_sterilising_impact.png") def fig_bayes(res): """Fig. 4: what the bombardment does to the fast-versus-slow abiogenesis Bayes factor.""" bf = res["bf"] en = res["bf_enabled"] labs = list(cfg.LIFE_DATES) short = ["LUCA\n4.2 Ga", "zircon C\n4.1 Ga", "Labrador\n3.95 Ga", "Isua\n3.7 Ga", "Dresser\n3.48 Ga"] sfdg = [f"sfd:{s}" for s in model.SFD_FAMILIES] fig, ax = plt.subplots(1, 2, figsize=(12.4, 5.1), gridspec_kw=dict(width_ratios=[1.6, 1])) a = ax[0] series = [(740., "sterilised by >= 740 km"), (440., "sterilised by >= 440 km"), (190., "sterilised by >= 190 km")] for i, lab in enumerate(labs): b0 = bf[(bf.group == "ALL") & (bf.evidence == lab) & (bf.m == 1)].B_survive.iloc[0] a.plot([i - 0.40, i + 0.40], [b0, b0], color=cfg.INK, lw=2.2, zorder=2) for k, (D, name) in enumerate(series): col = THR_COL[D] x = i - 0.27 + 0.18 * k q = bf[(bf.D_km == D) & (bf.evidence == lab) & (bf.m == 1)] fam = q[q.group.isin(sfdg)].B_sterilise a.plot([x, x], [fam.min(), fam.max()], color=col, lw=4.5, alpha=0.45, solid_capstyle="butt", zorder=3) for cond, dx in (("upper", -0.035), ("two", 0.035)): v = q[(q.group == "ALL") & (q.hse == cond)].B_sterilise.iloc[0] a.plot([x + dx], [v], "o", ms=7, mew=1.6, mec=col, mfc=col if cond == "upper" else cfg.SURFACE, zorder=5) e = en[(en.group == "ALL") & (en.evidence == lab) & (en.m == 1)] e2 = e[e.W_Gyr <= 0.2].B_enabled e5 = e.B_enabled a.plot([i + 0.33, i + 0.33], [max(e5.min(), 0.82), e5.max()], color=cfg.C["orange"], lw=1.2, zorder=3) a.plot([i + 0.33, i + 0.33], [max(e2.min(), 0.82), e2.max()], color=cfg.C["orange"], lw=4.5, solid_capstyle="butt", zorder=4) a.axhline(1.0, color=cfg.INK2, lw=0.8) handles = [Line2D([], [], color=cfg.INK, lw=2.2, label="life survives every impact; also the value if only\n" "surface-melting impacts (>= 2000 km) sterilise")] handles += [Line2D([], [], color=THR_COL[D], lw=0, marker="o", ms=7, mec=THR_COL[D], mfc=THR_COL[D], label=name) for D, name in series[::-1]] handles += [Line2D([], [], color=cfg.INK2, lw=0, marker="o", ms=7, mew=1.6, mec=cfg.INK2, mfc=cfg.INK2, label="filled: HSE mass as upper bound only"), Line2D([], [], color=cfg.INK2, lw=0, marker="o", ms=7, mew=1.6, mec=cfg.INK2, mfc=cfg.SURFACE, label="open: HSE mass, both bounds"), Line2D([], [], color=cfg.INK2, lw=4.5, alpha=0.45, label="bar: range over size distributions, both treatments"), Line2D([], [], color=cfg.C["orange"], lw=4.5, label="origin only within 10-200 Myr of a >= 440-1000 km impact\n(thin line: 500 Myr window)")] a.set_yscale("log") a.set_ylim(0.8, 400) a.set_yticks([1, 2, 3, 5, 10, 20, 50, 100, 200]) a.set_yticklabels(["1", "2", "3", "5", "10", "20", "50", "100", "200"]) a.set_xticks(range(len(labs))) a.set_xticklabels(short) a.set_xlim(-0.6, len(labs) - 0.4) a.set_xlabel("earliest evidence of life adopted") a.set_ylabel("Bayes factor, fast : slow abiogenesis") a.set_title("a One later hard step", fontsize=10.5) a.legend(handles=handles, fontsize=7.3, loc="upper right", ncol=2, columnspacing=1.2) b = ax[1] lab = labs[0] ms = np.array(sorted(bf.m.unique())) q0 = bf[(bf.group == "ALL") & (bf.evidence == lab) & (bf.hse == "upper")] b.plot(ms, q0[q0.D_km == 440.].sort_values("m").B_survive, color=cfg.INK, lw=2.2, marker="o", ms=4, label="survives every impact") for D, name in series[::-1]: for cond in cfg.HSE_MAIN: q = bf[(bf.group == "ALL") & (bf.evidence == lab) & (bf.hse == cond) & (bf.D_km == D)].sort_values("m") b.plot(ms, q.B_sterilise, color=THR_COL[D], ls=HSE_LS[cond], lw=1.8, marker="o", ms=4, mfc=THR_COL[D] if cond == "upper" else cfg.SURFACE, label=name if cond == "upper" else None) e = en[(en.group == "ALL") & (en.hse == "upper") & (en.D_km == 740.) & (en.evidence == lab) & (en.W_Gyr == 0.1)].sort_values("m") b.plot(e.m, e.B_enabled, color=cfg.C["orange"], marker="s", ms=4, label="impact-enabled origin, 100 Myr window") b.plot([], [], color=cfg.INK2, ls=HSE_LS["two"], lw=1.6, label="dashed: HSE mass, both bounds") b.axhline(1.0, color=cfg.INK2, lw=0.8) b.set_yscale("log") b.set_ylim(0.8, 400) b.set_yticks([1, 2, 3, 5, 10, 20, 50, 100, 200]) b.set_yticklabels(["1", "2", "3", "5", "10", "20", "50", "100", "200"]) b.set_xlabel("hard steps after abiogenesis") b.set_ylabel("Bayes factor, fast : slow abiogenesis") b.set_title("b LUCA at 4.2 Ga, by number of later hard steps", fontsize=10.5) b.legend(fontsize=7.5, loc="upper right") fig.suptitle("If large impacts sterilise, the evidence for fast abiogenesis rises; " "if they are needed for the origin, it vanishes", x=0.01, ha="left", fontsize=12.5, fontweight="bold") fig.tight_layout() return _out(fig, "fig4_bayes_factor.png") def make(res): return [fig_chronologies(res), fig_sfd(res), fig_survival(res), fig_bayes(res)]