"""Lunar chronology models: crater density <-> absolute age. All models are expressed in one common unit, "Neukum-equivalent N(1)" (km^-2), so that a basin's measured crater density can be re-dated under each model: neukum1983, marchi2009, yue2026 : N1(T) = a (exp(bT) - 1) + cT (published in N1) nesvorny2023 : N20(T) of their eq. 5 times N1/N20 = 1000 (the conversion Nesvorny et al. use themselves to compare with Neukum's N1) sawtooth2012 : integral of the Morbidelli et al. (2012) rate (for T < 4.1 Ga it is exactly Neukum's curve divided by 1400, so N1eq = 1400 N20), re-normalised to the paper's own anchor, Nectaris at 4.1 Ga: N1eq = 1400 N20(T) s(T) with s = 1 up to the Imbrium density (3.87 Ga), rising log-linearly to s_N = 1.61 at 4.1 Ga and constant before, where s_N is the ratio of Orgel et al.'s Nectaris density to the curve's value at 4.1 Ga sawtooth_aspub : the same curve with s = 1 (first-pass implementation: published curve applied to Orgel's densities; Nectaris lands at 4.32 Ga) sawtooth_uniform : the same curve with s = s_N at all ages (referee's variant) anchor_* : no functional form. Ages are interpolated (linear in ln N1) between radiometric anchors: Imbrium 3.922 Ga, Nectaris at a sampled age, SPA at 4.247 Ga (Chang'e-6) or 4.33 Ga (Joy et al.); before SPA the flux is extrapolated as exp((T - T_SPA)/tau). Each model object gives n1(T) (vectorised, T in Ga, valid up to 4.51 Ga) and age(N1). """ from __future__ import annotations import numpy as np import config as cfg T_MAX = cfg.T_MOON # 4.51 Ga: nothing is extrapolated beyond the Moon-forming impact _TG = np.linspace(0.0, T_MAX, 9021) # 0.5 Myr grid used for numerical inversion class ExpCF: """N1 = a (exp(bT) - 1) + cT.""" def __init__(self, name, a, b, c): self.name, self.a, self.b, self.c = name, a, b, c self._n = self.n1(_TG) def n1(self, T): T = np.asarray(T, float) return self.a * np.expm1(self.b * T) + self.c * T def age(self, n1): return np.interp(np.asarray(n1, float), self._n, _TG) class TabCF: """Chronology given by a tabulated monotonic N1eq(T).""" def __init__(self, name, n_on_grid): self.name = name self._n = np.asarray(n_on_grid, float) def n1(self, T): return np.interp(np.asarray(T, float), _TG, self._n) def age(self, n1): return np.interp(np.asarray(n1, float), self._n, _TG) # clips at T_MAX def nesvorny_n20(T, scale_a=1.0): """Nesvorny et al. 2023 eq. 5; T in Ga; returns N20 per 10^6 km^2.""" p = cfg.NESV T = np.asarray(T, float) t = np.maximum(p["t0"] - 1000.0 * T, 1e-9) # Myr after first solids return (scale_a * p["a"] * np.exp(-(t / p["tau1"]) ** p["al1"]) + p["b"] * np.exp(-(t / p["tau2"]) ** p["al2"]) + p["c"] * 1000.0 * T) def sawtooth_rate(T): """Morbidelli et al. 2012 dN20/dt in km^-2 Gyr^-1 (solid curve of their Fig. 3).""" p = cfg.SAW T = np.asarray(T, float) young = p["a"] * np.exp(p["b"] * T) + p["c"] old = p["A"] * np.exp(-(np.maximum(p["t0"] - T, 0.0) / p["tau"]) ** p["be"]) return np.where(T < p["tbreak"], young, old) def _build(): models = {k: ExpCF(k, **v) for k, v in cfg.CF_EXP.items()} models["nesvorny2023"] = TabCF("nesvorny2023", nesvorny_n20(_TG) * 1e-6 * cfg.N1_OVER_N20) r = sawtooth_rate(_TG) n20 = np.concatenate([[0.0], np.cumsum(0.5 * (r[1:] + r[:-1]) * np.diff(_TG))]) aspub = n20 * 1400.0 nk = models["neukum1983"] s_n = float(nk.n1(cfg.SAW_NECTARIS_ORGEL) / np.interp(cfg.SAW_NECTARIS_AGE, _TG, aspub)) ramp = np.clip((_TG - cfg.SAW_YOUNG_AGE) / (cfg.SAW_NECTARIS_AGE - cfg.SAW_YOUNG_AGE), 0.0, 1.0) models["sawtooth2012"] = TabCF("sawtooth2012", aspub * np.exp(np.log(s_n) * ramp)) models["sawtooth_aspub"] = TabCF("sawtooth_aspub", aspub) models["sawtooth_uniform"] = TabCF("sawtooth_uniform", aspub * s_n) models["sawtooth2012"].s_nectaris = s_n return models MODELS = _build() FUNCTIONAL = ["neukum1983", "marchi2009", "yue2026", "nesvorny2023", "sawtooth2012"] ANCHORED = ["anchor_ce6", "anchor_joy433"] ALL = FUNCTIONAL + ANCHORED LABEL = { "neukum1983": "Neukum 1983 (classic)", "marchi2009": "Marchi+ 2009", "yue2026": "Yue+ 2026 (Chang'e-6 refit)", "nesvorny2023": "Nesvorny+ 2023 (accretion tail)", "sawtooth2012": "Morbidelli+ 2012 sawtooth, Nectaris pinned at 4.1 Ga", "sawtooth_aspub": "sawtooth curve as published, Orgel densities (Nectaris 4.32 Ga)", "sawtooth_uniform": "sawtooth curve x 1.61 at all ages (Nectaris 4.1 Ga)", "anchor_ce6": "anchors, SPA 4.247 Ga (CE-6)", "anchor_joy433": "anchors, SPA 4.33 Ga (Joy+)", } def anchored_ages(n1, n1_imb, n1_nec, n1_spa, t_imb, t_nec, t_spa): """Ages for densities n1 (array) by interpolation through three anchors (scalars for one draw). Younger than Imbrium: Neukum's curve shifted so that it passes through the Imbrium anchor. Imbrium..Nectaris..SPA: linear in ln N1. Denser than SPA: clipped to the SPA age (SPA is stratigraphically the oldest basin).""" nk = MODELS["neukum1983"] n1 = np.asarray(n1, float) out = np.empty_like(n1) young = n1 <= n1_imb out[young] = nk.age(n1[young]) + (t_imb - nk.age(n1_imb)) xs = np.log([n1_imb, n1_nec, n1_spa]) out[~young] = np.interp(np.log(n1[~young]), xs, [t_imb, t_nec, t_spa]) return out def anchored_n1_grid(Tg, n1_imb, n1_nec, n1_spa, t_imb, t_nec, t_spa, tau_pre): """N1eq on the age grid Tg for one draw of an anchored chronology (inverse of anchored_ages, plus the exponential extrapolation before SPA with e-folding time tau_pre).""" nk = MODELS["neukum1983"] Tg = np.asarray(Tg, float) out = np.empty_like(Tg) young = Tg <= t_imb out[young] = nk.n1(np.clip(Tg[young] - (t_imb - nk.age(n1_imb)), 0, None)) mid = (~young) & (Tg <= t_spa) out[mid] = np.exp(np.interp(Tg[mid], [t_imb, t_nec, t_spa], np.log([n1_imb, n1_nec, n1_spa]))) old = Tg > t_spa out[old] = n1_spa * np.exp((Tg[old] - t_spa) / tau_pre) return out def summary_table(): """Model ages of fixed densities, for the report (densities in Neukum-equivalent N1, km^-2).""" import pandas as pd nk = MODELS["neukum1983"] rows = [] marks = {"Orientale (Orgel 3.81)": nk.n1(3.81), "Imbrium (Orgel 3.87)": nk.n1(3.87), "Nectaris (Orgel 4.17)": nk.n1(4.17), "SPA (Orgel 4.31)": nk.n1(4.31), "SPA (Yue+ 2026 N1 = 0.369)": cfg.SPA_N1_YUE[0]} for k in FUNCTIONAL + ["sawtooth_aspub", "sawtooth_uniform"]: m = MODELS[k] row = dict(model=k) for lab, n in marks.items(): row[lab] = float(m.age(n)) for T in (4.0, 4.247, 4.4, 4.5): row[f"N1eq({T})"] = float(m.n1(T)) row["N1eq(4.5)/N1eq(4.247)"] = float(m.n1(4.5) / m.n1(4.247)) rows.append(row) return pd.DataFrame(rows) if __name__ == "__main__": import pandas as pd pd.set_option("display.width", 250) print(summary_table().round(4).T)