"""Shared constants, paths, literature values and plot style for the impact_window study. Every literature number used anywhere in the analysis is defined here with its source, so a reader can check each one in a single place. Ages are in Ga (10^9 yr before present) unless stated; diameters in km. """ from __future__ import annotations import os import numpy as np HERE = os.path.dirname(os.path.abspath(__file__)) STUDY = os.path.abspath(os.path.join(HERE, "..")) DATA = os.path.join(STUDY, "data") FIGS = os.path.join(STUDY, "figures") RAW = os.environ.get("IMPACT_WINDOW_RAW", "fermi_raw/impact_window") for _d in (DATA, FIGS, RAW): os.makedirs(_d, exist_ok=True) SEED = 20261012 # --------------------------------------------------------------------------------------------- # Data sources (all re-downloaded by fetch_data.py) # --------------------------------------------------------------------------------------------- URL_NEUMANN = "https://www.ebi.ac.uk/europepmc/webservices/rest/PMC4646831/fullTextXML" URL_FASSETT = "https://ntrs.nasa.gov/api/citations/20130014881/downloads/20130014881.pdf" URL_ORGEL = "https://www.hou.usra.edu/meetings/lpsc2018/pdf/1395.pdf" URL_SBDB = "https://ssd-api.jpl.nasa.gov/sbdb_query.api" # --------------------------------------------------------------------------------------------- # Physical constants # --------------------------------------------------------------------------------------------- R_EARTH = 6371.0 # km R_MOON = 1737.4 # km VESC_EARTH = 11.186 # km/s VESC_MOON = 2.376 # km/s M_EARTH = 5.972e24 # kg RHO_IMP = 3000.0 # kg/m^3 (Marchi et al. 2014 use 3000) A_MOON_RE = (40.0, 60.3) # Earth-Moon distance in Earth radii, range used for 4.4-3.8 Ga # --------------------------------------------------------------------------------------------- # Lunar chronology functions. N1 = cumulative density of craters >= 1 km (km^-2), T in Gyr. # N1(T) = a (exp(bT) - 1) + cT # --------------------------------------------------------------------------------------------- CF_EXP = { # Neukum 1983 / Neukum & Ivanov 1994 / Neukum, Ivanov & Hartmann 2001; coefficients as quoted # by Marchi et al. 2009 (arXiv:0903.5137, eq. 11) and Morbidelli et al. 2012 (arXiv:1208.4624). "neukum1983": dict(a=5.44e-14, b=6.93, c=8.38e-4), # Marchi et al. 2009, AJ 137, 4936 (arXiv:0903.5137), eq. 11 and following text. "marchi2009": dict(a=1.23e-15, b=7.85, c=1.30e-3), # Yue et al. 2026, Sci. Adv. 12, eady9265 (PMC12871470), eq. 1: refit with the Chang'e-5 and # Chang'e-6 calibration points, SPA N(1) = 0.369 +/- 0.048 km^-2 at 4.247 +/- 0.005 Ga. "yue2026": dict(a=3.885e-15, b=7.595, c=7.377e-4), } # Nesvorny et al. 2023, Icarus 399, 115545 (arXiv:2303.17736), eq. 5: # N20(T) [per 10^6 km^2] = a exp(-(t/6 Myr)^0.45) + b exp(-(t/65 Myr)^0.6) + c T, # t = 4570 - T (Myr), a = 6.7e4, b = 40, c = 1.7e-3 Myr^-1. NESV = dict(a=6.7e4, tau1=6.0, al1=0.45, b=40.0, tau2=65.0, al2=0.6, c=1.7e-3, t0=4570.0) # Morbidelli et al. 2012, EPSL 355-356, 144 (arXiv:1208.4624), caption of their Fig. 3: # dN20/dt [km^-2 Gyr^-1] = 2.7e-16 exp(6.93 t) + 5.9e-7 for t < 4.1 Gyr # = 2.5e-2 exp(-((4.5 - t)/0.003)^0.34) for t > 4.1 Gyr SAW = dict(a=2.7e-16, b=6.93, c=5.9e-7, A=2.5e-2, t0=4.5, tau=0.003, be=0.34, tbreak=4.1) # Morbidelli et al. 2012 build the sawtooth with the Nectaris basin at 4.1 Ga (the onset of the # cataclysm), using N20(Nectaris) = 8.66e-5 km^-2 from Marchi et al. 2012. Orgel et al. 2018 # measure a higher density on Nectaris (4.17 Ga on Neukum's curve). The main-grid sawtooth is # therefore re-normalised so that Orgel's Nectaris density falls at SAW_NECTARIS_AGE; densities # up to Imbrium's (SAW_YOUNG_AGE on Neukum's curve) keep the published, sample-calibrated curve. SAW_NECTARIS_AGE = 4.10 SAW_NECTARIS_ORGEL = 4.17 # Orgel et al. 2018 model age of Nectaris on Neukum's curve SAW_YOUNG_AGE = 3.87 # Orgel et al. 2018 model age of Imbrium on Neukum's curve # N(1)/N(20) of the Neukum (1983) production function, as stated by Nesvorny et al. 2023 # ("N1/N20 ~ 1000" for the old production function; ~420 for the Neukum & Ivanov 1994 one). N1_OVER_N20 = 1000.0 # Radiometric anchors SPA_CE6 = (4.247, 0.005) # Su et al. 2025, Natl Sci Rev 12, nwaf103: 4247 +/- 5 Ma (Pb-Pb, norites) SPA_JOY = (4.33, 0.01) # Joy et al. 2025, Nat Astron 9, 55: ~4.33 Ga (NWA 2995); 0.01 Gyr assumed IMBRIUM_AGE = (3.922, 0.012) # Fra Mauro, as tabulated by Yue et al. 2026 (their Table 1) SPA_N1_YUE = (0.369, 0.048) # Yue et al. 2026 Table 1: N(1) of the SPA basin (km^-2) T_MOON = 4.51 # Moon-forming impact (calibration used by Moody et al. 2024) # --------------------------------------------------------------------------------------------- # Impactor-size thresholds on Earth (km). See report for the sources of each. # --------------------------------------------------------------------------------------------- THRESH = { "photic190": 190.0, # top ~200 m of ocean boiled (energy-scaled from Sleep et al. 1989; NOT checked at source) "ocean440": 440.0, # whole ocean vaporised: Sleep et al. 1989 (350-440 km per Zahnle et al. 2020) "ocean500": 500.0, # Marchi et al. 2014 definition of "global ocean vaporisation" "ocean740": 740.0, # Citron & Stewart 2022: >= 740 km from 3-D simulations (lower limit) "surface1000": 1000.0, # Marchi et al. 2014 "global sterilisation" threshold (citing Sleep et al. 1989) "melt2000": 2000.0, # Citron & Stewart 2022: 2000-2700 km to melt the entire surface } D_GRID = np.array([150., 190., 250., 300., 350., 440., 500., 600., 740., 1000., 1300., 2000., 2700.]) T_REPORT = np.array([4.4, 4.3, 4.2, 4.1, 4.0, 3.9, 3.8]) # Late-accretion (HSE) mass on Earth: (0.7-3.0)e22 kg, Walker 2009 as quoted by Marchi et al. 2014 HSE_MASS = (0.7e22, 3.0e22) # Treatments of that constraint (lower, upper bound in kg on the mass delivered by the modelled # impactor population since the Moon-forming impact). Every Monte Carlo realisation is kept or # discarded under each treatment, so all of them come from the same draws. # upper the veneer may have come from elsewhere (one or a few giant early bodies outside the # modelled population: Bottke et al. 2010; Brasser et al. 2016; Genda et al. 2017), so # the modelled population must only not exceed the budget [first-pass choice] # two the modelled population delivered the veneer: mass inside the HSE range # none no mass constraint # upper_x5 / two_x5 upper bound five times higher (metal of large differentiated impactors # partly lost to Earth's core, Marchi et al. 2018), without / with the lower bound # two_half lower bound halved HSE_CONDS = { "upper": (0.0, HSE_MASS[1]), "two": (HSE_MASS[0], HSE_MASS[1]), "none": (0.0, float("inf")), "upper_x5": (0.0, 5.0 * HSE_MASS[1]), "two_x5": (HSE_MASS[0], 5.0 * HSE_MASS[1]), "two_half": (0.5 * HSE_MASS[0], HSE_MASS[1]), } HSE_MAIN = ["upper", "two"] # the two treatments carried through every table HSE_LABEL = {"upper": "HSE mass as upper bound only", "two": "HSE mass, both bounds", "none": "no mass constraint", "upper_x5": "upper bound x 5 only", "two_x5": "lower bound and upper bound x 5", "two_half": "lower bound halved, upper bound"} # Archean spherule layers: Johnson et al. 2016, Icarus 271, 350, Table 1 # (name, age Ga, impactor diameter min, max km) SPHERULES = [ ("Dales Gorge & Kuruman", 2.49, 11, 39), ("Bee Gorge", 2.54, 13, 31), ("Reivilo & Paraburdoo", 2.55, 17, 29), ("Jeerinah, Carawine & Monteville", 2.63, 10, 67), ("S5", 3.23, 37, 79), ("S4", 3.24, 31, 49), ("S3", 3.24, 42, 67), ("S2", 3.26, 29, 88), ("S6", 3.28, 37, 79), ("S8", 3.30, 37, 79), ("S7", 3.42, 37, 79), ("S1 & Warrawoona", 3.47, 23, 39), ] SPHERULE_WINDOWS = [(2.44, 2.68), (3.2, 3.5)] # intervals Johnson et al. 2016 treat as searched # --------------------------------------------------------------------------------------------- # Biology-side dates, copied from the timing study (studies/timing/code/config.py, data/events.csv) # --------------------------------------------------------------------------------------------- T0_ZIRCON = 4.404 T0_LATE = 4.50 LUCA_PAPER = (4.20, 4.09, 4.33) # Moody et al. 2024 (centre, 2.5%, 97.5%) LIFE_DATES = { "LUCA 4.2 Ga (Moody+ 2024)": 4.20, "zircon graphite 4.1 Ga (contested)": 4.10, "Labrador 3.95 Ga (contested)": 3.95, "Isua 3.7 Ga (contested)": 3.70, "Dresser 3.48 Ga (accepted)": 3.48, } FUTURE_BASE = 0.9 # Gyr, Kipping's choice F_SHORT = (1.08, 0.14) F_LONG = (1.60, 1.87) # --------------------------------------------------------------------------------------------- # Plot style (same palette as the timing study) # --------------------------------------------------------------------------------------------- C = dict(blue="#2a78d6", orange="#eb6834", aqua="#1baf7a", yellow="#eda100", magenta="#e87ba4", green="#008300", violet="#6250d6", red="#e34948") SERIES = [C["blue"], C["orange"], C["aqua"], C["yellow"], C["magenta"], C["green"], C["violet"], C["red"]] SURFACE = "#fcfcfb" INK = "#0b0b0b" INK2 = "#52514e" GRID = "#e4e3df" def style(): import matplotlib as mpl mpl.use("Agg") mpl.rcParams.update({ "figure.facecolor": SURFACE, "axes.facecolor": SURFACE, "savefig.facecolor": SURFACE, "axes.edgecolor": "#b9b8b2", "axes.labelcolor": INK2, "text.color": INK, "xtick.color": INK2, "ytick.color": INK2, "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.6, "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, "figure.dpi": 130, "savefig.dpi": 160, "axes.axisbelow": True, }) def split_normal(rng, centre, q025, q975, size): """Two-piece normal with the given median and 2.5%/97.5% quantiles (as in the timing study).""" z = rng.standard_normal(size) s_lo = (centre - q025) / 1.959964 s_hi = (q975 - centre) / 1.959964 return centre + np.where(z < 0, s_lo, s_hi) * z