"""Synthesis: the completeness function C(location, size) for an inert artefact in three tiers, and the probability that a single artefact would already have been found. Tiers (revision: the first version mixed the first two under the word 'detected'): 1. in archived data surfaces: fraction of the surface imaged with pixel scale <= D/2 orbits: same as tier 2 (unlinked detections are not modelled, so this is a lower bound for orbits) 2. catalogued orbits: a designated catalogue object surfaces: imaged AND examined in a way that would log a discrete feature (tier 1 x examined fraction; no global catalogue of metre-scale lunar or martian features exists) 3. recognised catalogued AND recognised as a non-human artefact Provenance of the tier-1/2 completeness per location: measured Moon and Mars image coverage (s04, s07); other near-Earth orbits (catalogue / population model, s01); main belt and Jupiter Trojans (catalogue turnover, s06) simulated the three geocentric shells (minimoon-corrected) and the four co-orbital / Earth-like classes (survey simulator, s05, with the elongation weight measured in s02) literature Earth orbit below 60,000 km; Mercury / Venus / Ceres surfaces bound outer Solar System (photometric bound) Recognition is a separate factor R(location, size) built from measured quantities (fraction of objects with long arcs or any physical characterisation - for the belts now as a function of absolute magnitude -, fraction of imagery examined) and a small number of explicit judgement parameters, each with a low / central / high value (RECOG below). The low / high columns set every parameter low or high together: they are scenario envelopes, not confidence intervals. Outputs: data/completeness_matrix.csv, data/completeness_table.csv, data/detection_thresholds.csv, data/found_probability.csv, data/recognition_model.json, data/weights.csv, data/pan2026_comparison.csv """ import json import os import numpy as np import pandas as pd from common import DATA, RAW, D_from_H, H_from_D SIZES = np.round(10 ** np.arange(0.0, 4.0001, 0.1), 6) # 1 m ... 10 km TABLE_SIZES = [1, 3, 10, 30, 100, 300, 1000, 3000, 10000] ALBEDOS = [0.2, 0.05, 0.01] MOON_AREA_KM2 = 4 * np.pi * 1737.4 ** 2 # ---- judgement parameters of the recognition model: (low, central, high) ---- RECOG = { # P(attributed to a non-human origin | flagged as artificial). Human hardware beyond low Earth orbit is # at most ~13 m across (Saturn S-IVB; computed from GCAT below); smaller unidentified objects near Earth # are filed as human debris (95 unidentified geocentric objects in the Project Pluto index). "attrib_small": (0.0, 0.01, 0.1), # D <= largest human deep-space object "attrib_large": (0.5, 0.9, 1.0), # D >= 30 m # P(flagged | physically characterised) for an object with rock-like density (no radiation-pressure signal) "flag_given_characterised": (0.1, 0.3, 1.0), # weight of the 'hollow, spacecraft-like' hypothesis (radiation pressure betrays it once the arc is >= 45 d) "hollow_weight": (0.0, 0.5, 1.0), # fraction of metre-scale surface imagery (objects < 50 m) examined in a way that would notice an artefact; # low value = area of the only published lunar anomaly search (Lesnikowski et al. 2020), set in code "exam_small": (None, 0.01, 0.3), # same for objects >= 1 km (visible in global mosaics used for crater catalogues and geological maps) "exam_large": (0.5, 0.9, 1.0), # HiRISE frames are targeted and inspected by people "exam_hirise": (0.1, 0.5, 1.0), # P(flagged | albedo / lightcurve known) for a belt object "flag_belt": (0.001, 0.1, 1.0), "px_detect": (3.0, 2.0, 1.5), # pixels across the object needed to register as a feature "px_recognise": (13.0, 8.0, 6.0), # pixels needed to recognise an unnatural shape (Johnson criteria) } LV = {"lo": 0, "c": 1, "hi": 2} LOCATIONS = [ ("Lunar surface", "imaged everywhere at 0.5-2 m per pixel; almost none of it examined for artefacts"), ("Mars surface", "HiRISE covers 3.6%; elsewhere CTX at 5-6 m"), ("Other solid surfaces (Mercury, Venus, Ceres)", "global imaging at 250-500 m (literature)"), ("Earth orbit below 60,000 km", "tracked by radar/optical surveillance (literature); lost among human debris"), ("High Earth orbit (60,000-300,000 km)", "survey depth, trailing; presumed human debris"), ("Earth-Moon L4/L5 and lunar-distance orbits", "survey depth, trailing, moonlight; presumed human debris"), ("Distant geocentric orbits (0.6-1.2 million km)", "survey depth; presumed human debris"), ("Earth quasi-satellites", "survey depth; arcs too short to characterise"), ("Sun-Earth L4/L5 (parked)", "solar elongation 54-65 deg: measured survey coverage there is about 1% of that at opposition"), ("Earth horseshoe co-orbitals", "decades spent behind the Sun"), ("Earth-like heliocentric orbits (0.9-1.1 au)", "long synodic period; survey depth"), ("Other near-Earth orbits", "survey depth; brief apparitions"), ("Main asteroid belt", "survey depth; no characterisation"), ("Jupiter Trojans", "survey depth; no characterisation"), ("Outer Solar System (beyond 5.5 au)", "reflected light falls as r^-4"), ] PROVENANCE = { "Lunar surface": "measured", "Mars surface": "measured", "Other solid surfaces (Mercury, Venus, Ceres)": "literature", "Earth orbit below 60,000 km": "literature", "High Earth orbit (60,000-300,000 km)": "simulated", "Earth-Moon L4/L5 and lunar-distance orbits": "simulated", "Distant geocentric orbits (0.6-1.2 million km)": "simulated", "Earth quasi-satellites": "simulated", "Sun-Earth L4/L5 (parked)": "simulated", "Earth horseshoe co-orbitals": "simulated", "Earth-like heliocentric orbits (0.9-1.1 au)": "simulated", "Other near-Earth orbits": "measured", "Main asteroid belt": "measured", "Jupiter Trojans": "measured", "Outer Solar System (beyond 5.5 au)": "bound", } SURFACES = {"Lunar surface", "Mars surface", "Other solid surfaces (Mercury, Venus, Ceres)"} def provenance(loc, tier): """Provenance of a tier: the surface rows of tiers 2 and 3 rest on the examined fraction (a judgement parameter), and every tier-3 number rests on the recognition scenario.""" if tier == "in archived data": return PROVENANCE[loc] if tier == "catalogued": return PROVENANCE[loc] + " x examined fraction (judgement)" if loc in SURFACES else PROVENANCE[loc] return "scenario model (judgement parameters)" SIM_CLASS = { "High Earth orbit (60,000-300,000 km)": "High Earth orbit (60,000-300,000 km)", "Earth-Moon L4/L5 and lunar-distance orbits": "Earth-Moon L4/L5 and lunar-distance orbits", "Distant geocentric orbits (0.6-1.2 million km)": "Distant geocentric orbits (0.6-1.2 million km)", "Earth quasi-satellites": "Earth quasi-satellites", "Sun-Earth L4/L5 (parked)": "Sun-Earth L4/L5, parked near the point (libration < 10 deg)", "Earth horseshoe co-orbitals": "Earth horseshoe co-orbitals", "Earth-like heliocentric orbits (0.9-1.1 au)": "Earth-like heliocentric orbits (0.9-1.1 au)", } GEOCENTRIC = {"Earth orbit below 60,000 km", "High Earth orbit (60,000-300,000 km)", "Earth-Moon L4/L5 and lunar-distance orbits", "Distant geocentric orbits (0.6-1.2 million km)"} NEAR_EARTH_HELIO = {"Earth quasi-satellites", "Sun-Earth L4/L5 (parked)", "Earth horseshoe co-orbitals", "Earth-like heliocentric orbits (0.9-1.1 au)", "Other near-Earth orbits"} # ---- weighting schemes over locations (each sums to 1) ---- WEIGHTS = { # A is Earth-centric by construction: 10 of the 15 location bins are the Moon, geocentric shells or # near-Earth co-orbital classes. A volume- or surface-weighted prior would resemble C. "A: equal weight per location": {n: 1.0 for n, _ in LOCATIONS}, "B: Earth-watching (Moon, cislunar, co-orbital)": { "Lunar surface": 0.25, "Earth-Moon L4/L5 and lunar-distance orbits": 0.15, "Earth quasi-satellites": 0.15, "Sun-Earth L4/L5 (parked)": 0.10, "High Earth orbit (60,000-300,000 km)": 0.10, "Earth orbit below 60,000 km": 0.05, "Distant geocentric orbits (0.6-1.2 million km)": 0.05, "Earth horseshoe co-orbitals": 0.05, "Earth-like heliocentric orbits (0.9-1.1 au)": 0.05, "Other near-Earth orbits": 0.01, "Mars surface": 0.01, "Main asteroid belt": 0.01, "Jupiter Trojans": 0.01, "Other solid surfaces (Mercury, Venus, Ceres)": 0.005, "Outer Solar System (beyond 5.5 au)": 0.005}, "C: resource-rich and long-lived (belt, outer system)": { "Main asteroid belt": 0.40, "Outer Solar System (beyond 5.5 au)": 0.20, "Jupiter Trojans": 0.10, "Mars surface": 0.05, "Other solid surfaces (Mercury, Venus, Ceres)": 0.05, "Lunar surface": 0.05, "Other near-Earth orbits": 0.05, "Earth-like heliocentric orbits (0.9-1.1 au)": 0.02, "Earth quasi-satellites": 0.02, "Sun-Earth L4/L5 (parked)": 0.02, "Earth horseshoe co-orbitals": 0.02, "Earth-Moon L4/L5 and lunar-distance orbits": 0.005, "High Earth orbit (60,000-300,000 km)": 0.005, "Distant geocentric orbits (0.6-1.2 million km)": 0.005, "Earth orbit below 60,000 km": 0.005}, } # Residence weights of Pan et al. (2026, doi 10.21203/rs.3.rs-10268655/v1, Table 2) and their detection # probabilities (Table 3), for a like-for-like comparison. Mapping of their seven regions onto ours: NEO -> # other near-Earth orbits; Earth co-orbital -> quasi-satellites, Sun-Earth L4/L5, horseshoes and Earth-like # orbits in equal parts; Earth-Moon Lagrange -> Earth-Moon L4/L5 and lunar-distance orbits; lunar surface; # main belt; Kuiper belt + outer/misc -> outer Solar System. PAN_WEIGHTS = {"default (outer-dominated)": (0.020, 0.005, 0.002, 0.003, 0.200, 0.770), "diffuse (uniform)": (0.05, 0.01, 0.01, 0.01, 0.30, 0.62), "lurker (Benford)": (0.10, 0.40, 0.20, 0.10, 0.05, 0.15)} PAN_PDET = {"default (outer-dominated)": {10: 0.0019, 100: 0.0209, 1000: 0.184}, "diffuse (uniform)": {10: 0.0053, 100: 0.0461, 1000: 0.304}, "lurker (Benford)": {10: 0.0351, 100: 0.2401, 1000: 0.4343}} PAN_LAMBDA95 = {"default (outer-dominated)": {10: 14.3, 100: 3.16, 1000: 0.68}, "diffuse (uniform)": {10: 7.97, 100: 1.80, 1000: 0.46}, "lurker (Benford)": {10: 2.27, 100: 0.55, 1000: 0.34}} def pan_weights(name): neo, co, eml, moon, belt, outer = PAN_WEIGHTS[name] w = {n: 0.0 for n, _ in LOCATIONS} w["Other near-Earth orbits"] = neo for c in ["Earth quasi-satellites", "Sun-Earth L4/L5 (parked)", "Earth horseshoe co-orbitals", "Earth-like heliocentric orbits (0.9-1.1 au)"]: w[c] = co / 4.0 w["Earth-Moon L4/L5 and lunar-distance orbits"] = eml w["Lunar surface"] = moon w["Main asteroid belt"] = belt w["Outer Solar System (beyond 5.5 au)"] = outer return w def loginterp(D, pts): """Piecewise-linear in log10(D) through pts = [(D, value), ...]; constant outside.""" x = np.log10([p[0] for p in pts]) y = [p[1] for p in pts] return np.interp(np.log10(D), x, y) class Model: def __init__(self): d = DATA self.moon = pd.read_csv(os.path.join(d, "moon_nac_samples.csv")) self.ctx = pd.read_csv(os.path.join(d, "mars_ctx_samples.csv")) surf = pd.read_csv(os.path.join(d, "surface_coverage_summary.csv")) g = lambda m: float(surf[surf.metric == m].value.iloc[0]) # noqa: E731 self.hirise_any = g("union fraction of surface, any scale") self.hirise_03 = g("union fraction at <=0.3 m/pixel map scale") self.curves = pd.read_csv(os.path.join(d, "parked_pdet_curves.csv")) self.neo = pd.read_csv(os.path.join(d, "neo_completeness.csv")) self.belt = pd.read_csv(os.path.join(d, "belt_H_histograms.csv")) self.beltsum = pd.read_csv(os.path.join(d, "belt_completeness.csv")).set_index("region") rs = json.load(open(os.path.join(d, "recognition_stats.json"))) self.f_belt_phys = rs["all_asteroid_characterisation"]["MBA_frac_albedo"] # catalogue averages (first version; kept for the record) self.f_tjn_phys = rs["all_asteroid_characterisation"]["TJN_frac_albedo"] self.f_tno_phys = rs["all_asteroid_characterisation"]["TNO_frac_albedo"] # (revision) fraction with any physical characterisation as a function of absolute magnitude self.phys_by_H = {} for g_, rows_ in rs["characterisation_by_H"].items(): Hc_ = np.array([0.5 * (max(r["H_lo"], 4) + min(r["H_hi"], 22)) for r in rows_]) fr_ = np.array([r["n_any_physical"] / r["n"] if r["n"] >= 20 else np.nan for r in rows_]) ok_ = np.isfinite(fr_) self.phys_by_H[g_] = (Hc_[ok_], fr_[ok_]) hh = pd.read_csv(os.path.join(d, "human_hardware_prediction.csv")) self.d_human_max = float(hh.D_eq_m.max()) # measured recognition inputs: arcs and physical characterisation of real near-Earth asteroids el = pd.read_csv(os.path.join(d, "earthlike_objects.csv")) el["arc"] = pd.to_numeric(el.data_arc, errors="coerce") sb = json.load(open(os.path.join(RAW, "sbdb_neo.json"))) nea = pd.DataFrame(sb["data"], columns=sb["fields"]) nea = nea[nea.kind.str.startswith("a")].copy() for c in ["H", "data_arc"]: nea[c] = pd.to_numeric(nea[c], errors="coerce") ph = json.load(open(os.path.join(RAW, "sbdb_neo_phys.json"))) phys = pd.DataFrame(ph["data"], columns=ph["fields"]) phys["H"] = pd.to_numeric(phys.H, errors="coerce") phys["any"] = phys.rot_per.notna() | phys.spec_B.notna() | phys.spec_T.notna() | phys.albedo.notna() | phys.A1.notna() self.hb = np.array([-10.0, 19.0, 22.0, 25.0, 28.0, 99.0]) rows = [] for a, b in zip(self.hb[:-1], self.hb[1:]): e_ = el[(el.H >= a) & (el.H < b)] n_ = nea[(nea.H >= a) & (nea.H < b)] p_ = phys[(phys.H >= a) & (phys.H < b)] rows.append(dict(H_lo=a, H_hi=b, n_earthlike=len(e_), f_arc45_earthlike=float((e_.arc >= 45).mean()) if len(e_) else np.nan, n_nea=len(n_), f_arc45_nea=float((n_.data_arc >= 45).mean()), f_phys_nea=float(p_["any"].mean()))) self.arc = pd.DataFrame(rows) # too few Earth-like objects brighter than H = 22 to measure the fraction: use the all-NEA value there self.arc["f_arc45_earthlike"] = np.where(self.arc.n_earthlike >= 10, self.arc.f_arc45_earthlike, self.arc.f_arc45_nea) # Lesnikowski et al. (2020): 255,000 tiles covering ~64 km^2 around two Apollo sites self.lesnikowski_frac = 64.0 / MOON_AREA_KM2 self.params = {k: list(v) for k, v in RECOG.items()} self.params["exam_small"][0] = self.lesnikowski_frac def par(self, key, lv): return self.params[key][LV[lv]] # ------------------------------------------------------------------ detection def sim(self, cls, H, lv): """Survey-simulator completeness. Central: base depth variant, measured elongation weight, solved minimoon factor, empirical trailing law. Low / high: minimum / maximum over the five depth variants, the 95% bounds of the measured elongation weight (heliocentric classes), the range of the minimoon factor and the theoretical trailing law (geocentric classes).""" c = self.curves[self.curves.location == cls] alt = [k for k in ["p_wlo", "p_whi", "p_klo", "p_khi", "p_theory", "p_theory_klo", "p_theory_khi"] if k in c.columns] out, cen = [], None for v in sorted(c.variant.unique()): cv = c[c.variant == v].sort_values("H") for col in ["p_detect"] + alt: if cv[col].notna().all(): y = np.interp(H, cv.H.values, cv[col].values, left=cv[col].values[0], right=0.0) out.append(y) if v == "base" and col == "p_detect": cen = y out = np.array(out) return {"c": cen, "lo": out.min(axis=0), "hi": out.max(axis=0)}[lv] def f_phys(self, group, H): Hc_, fr_ = self.phys_by_H[group] return np.interp(H, Hc_, fr_) def neo_catalogue(self, H, lv): Hc = 0.5 * (self.neo.H1 + self.neo.H2).values res = [] for col in ["compl_diff_neomod2", "compl_diff_hc21"]: c = np.clip(self.neo[col].values, 1e-9, 1.0) # running maximum from the faint end so that the curve is monotonic (bright bins have small-number noise) c = np.maximum.accumulate(c[::-1])[::-1] slope = np.polyfit(Hc[-4:], np.log10(c[-4:]), 1)[0] lc = np.interp(H, Hc, np.log10(c), left=0.0) # brighter than H = 15.5: taken as complete lc = np.where(H > Hc[-1], np.log10(c[-1]) + slope * (H - Hc[-1]), lc) # extrapolated beyond H = 27.5 res.append(np.clip(10 ** lc, 0, 1)) res = np.array(res) return {"c": res[0], "lo": res.min(axis=0), "hi": res.max(axis=0)}[lv] def belt_curve(self, region, H): b = self.belt[self.belt.region == region].sort_values("H") Hmode = self.beltsum.loc[region, "H_mode"] r = np.clip(b.ratio_to_powerlaw.values, 0, 1) r = np.where(b.H.values < Hmode - 1.0, 1.0, r) r = np.minimum.accumulate(r) # monotonic decline return np.interp(H, b.H.values, r, left=1.0, right=0.0) def detect(self, loc, D, p, lv): D = np.asarray(D, float) H = H_from_D(D, p) npx = self.par("px_detect", lv) if loc == "Lunar surface": res = self.moon.best_res_goodlight.values if lv != "hi" else self.moon.best_res.values res = np.where(np.isfinite(res), res, np.inf) return np.array([(res <= d / npx).mean() for d in D]) if loc == "Mars surface": f_h = np.where(D >= npx * 0.5, self.hirise_any, np.where(D >= npx * 0.3, self.hirise_03, 0.0)) f_c = np.array([(self.ctx.best_res.values <= d / npx).mean() for d in D]) return 1 - (1 - f_h) * (1 - f_c) if loc.startswith("Other solid surfaces"): return np.mean([D >= npx * r for r in (250.0, 300.0, 500.0)], axis=0) # literature (Lazio 2026, Table 2) if loc == "Earth orbit below 60,000 km": # literature: objects > 10 cm tracked in low orbit by radar; ~0.3-1 m at geosynchronous altitude (optical). d_min = {0.2: 1.0, 0.05: 1.0, 0.01: 3.0}[p] * {"lo": 2.0, "c": 1.0, "hi": 0.5}[lv] return np.where(D >= d_min, 1.0, 0.5) # 0.5: the radar-tracked low-orbit half only if loc in SIM_CLASS: return self.sim(SIM_CLASS[loc], H, lv) if loc == "Other near-Earth orbits": return self.neo_catalogue(H, lv) if loc == "Main asteroid belt": cs = np.array([self.belt_curve(r, H) for r in ["Inner main belt (2.0-2.5 au)", "Middle main belt (2.5-2.82 au)", "Outer main belt (2.82-3.27 au)"]]) return {"c": cs.mean(axis=0), "lo": cs.min(axis=0), "hi": cs.max(axis=0)}[lv] if loc == "Jupiter Trojans": return self.belt_curve("Jupiter Trojans (4.8-5.5 au)", H + {"lo": 0.5, "c": 0.0, "hi": -0.5}[lv]) if loc.startswith("Outer Solar System"): # photometric bound: fraction of heliocentric distances 5.5-50 au (uniform in r, an assumption) at # which the object is brighter than V = 22 (21.5 / 22.5) at opposition vlim = {"lo": 21.5, "c": 22.0, "hi": 22.5}[lv] x = 10 ** ((vlim - H) / 5.0) # r (r - 1) <= x rmax = 0.5 * (1 + np.sqrt(1 + 4 * x)) return np.clip((rmax - 5.5) / (50.0 - 5.5), 0, 1) raise KeyError(loc) # ------------------------------------------------------------------ recognition def attribution(self, D, lv): return loginterp(D, [(self.d_human_max, self.par("attrib_small", lv)), (30.0, self.par("attrib_large", lv))]) def exam(self, D, lv): return 10 ** loginterp(D, [(50.0, np.log10(self.par("exam_small", lv))), (1000.0, np.log10(self.par("exam_large", lv)))]) def flag_near_earth(self, H, lv, earthlike): ib = np.clip(np.searchsorted(self.hb, H, side="right") - 1, 0, len(self.arc) - 1) f_arc = self.arc["f_arc45_earthlike" if earthlike else "f_arc45_nea"].values[ib] f_phys = self.arc.f_phys_nea.values[ib] w = self.par("hollow_weight", lv) return w * f_arc + (1 - w) * self.par("flag_given_characterised", lv) * f_phys def catalogued(self, loc, D, p, lv): """Tier 2. Orbits: a designated catalogue object (= detect). Surfaces: imaged at the detection pixel count AND examined in a way that would log a discrete feature.""" D = np.asarray(D, float) det = self.detect(loc, D, p, lv) if loc == "Lunar surface": return det * self.exam(D, lv) if loc == "Mars surface": npx = self.par("px_detect", lv) f_h = np.where(D >= npx * 0.5, self.hirise_any, np.where(D >= npx * 0.3, self.hirise_03, 0.0)) * self.par("exam_hirise", lv) f_c = np.array([(self.ctx.best_res.values <= d / npx).mean() for d in D]) * self.exam(D, lv) return 1 - (1 - f_h) * (1 - f_c) if loc.startswith("Other solid surfaces"): return det * self.par("exam_large", lv) return det def recognised(self, loc, D, p, lv): """Completeness for 'found AND recognised as a non-human artefact'.""" D = np.asarray(D, float) H = H_from_D(D, p) det = self.detect(loc, D, p, lv) npx = self.par("px_recognise", lv) if loc == "Lunar surface": res = self.moon.best_res_goodlight.values if lv != "hi" else self.moon.best_res.values res = np.where(np.isfinite(res), res, np.inf) return np.array([(res <= d / npx).mean() for d in D]) * self.exam(D, lv) if loc == "Mars surface": f_h = np.where(D >= npx * 0.5, self.hirise_any, np.where(D >= npx * 0.3, self.hirise_03, 0.0)) * self.par("exam_hirise", lv) f_c = np.array([(self.ctx.best_res.values <= d / npx).mean() for d in D]) * self.exam(D, lv) return 1 - (1 - f_h) * (1 - f_c) if loc.startswith("Other solid surfaces"): return np.mean([D >= npx * r for r in (250.0, 300.0, 500.0)], axis=0) * self.par("exam_large", lv) if loc in GEOCENTRIC: return det * self.attribution(D, lv) if loc in NEAR_EARTH_HELIO: return det * self.flag_near_earth(H, lv, earthlike=(loc != "Other near-Earth orbits")) * self.attribution(D, lv) grp = {"Main asteroid belt": "main_belt", "Jupiter Trojans": "trojans"}.get(loc, "tno_centaur") return det * self.f_phys(grp, H) * self.par("flag_belt", lv) * self.attribution(D, lv) TIERS = [("C_arch_", "in archived data"), ("C_cat_", "catalogued"), ("C_rec_", "recognised")] def main(): m = Model() names = [n for n, _ in LOCATIONS] lim = dict(LOCATIONS) rows = [] for loc in names: for p in ALBEDOS: r = dict(location=loc, albedo=p, provenance=PROVENANCE[loc], provenance_catalogued=provenance(loc, "catalogued")) cols = {} for lv in LV: cols["C_arch_" + lv] = m.detect(loc, SIZES, p, lv) cols["C_cat_" + lv] = m.catalogued(loc, SIZES, p, lv) cols["C_rec_" + lv] = m.recognised(loc, SIZES, p, lv) # make sure lo <= c <= hi after all the piecewise rules for k, _ in TIERS: st = np.sort(np.array([cols[k + "lo"], cols[k + "c"], cols[k + "hi"]]), axis=0) cols[k + "lo"], cols[k + "hi"] = st[0], st[2] for i, D in enumerate(SIZES): rr = dict(r, size_m=D, H=float(H_from_D(D, p))) rr.update({k: float(v[i]) for k, v in cols.items()}) rr["limitation"] = lim[loc] rows.append(rr) M = pd.DataFrame(rows) M.to_csv(os.path.join(DATA, "completeness_matrix.csv"), index=False) # ---- table: completeness versus size (albedo 0.2) in the three tiers ---- def at(df, D): return df.iloc[(np.abs(np.log10(df.size_m.values) - np.log10(D))).argmin()] trows, thr = [], [] for loc in names: for p in ALBEDOS: sub = M[(M.location == loc) & (M.albedo == p)] t = dict(location=loc, albedo=p, provenance=PROVENANCE[loc]) for key, tag in [("C_arch_", "arch"), ("C_cat_", "cat")]: for lv in ["c", "lo", "hi"]: for level in (0.5, 0.9): ok = sub[sub[key + lv] >= level] t[f"D{int(level * 100)}_{tag}_{lv}_m"] = float(ok.size_m.min()) if len(ok) else np.nan ok = sub[sub.C_rec_c >= 0.5] t["D50_recognised_c_m"] = float(ok.size_m.min()) if len(ok) else np.nan t["plateau_cat_c"] = float(sub.C_cat_c.max()) t["plateau_cat_lo"] = float(sub.C_cat_lo.max()) t["plateau_cat_hi"] = float(sub.C_cat_hi.max()) t["limitation"] = lim[loc] thr.append(t) sub = M[(M.location == loc) & (M.albedo == 0.2)] for key, name in TIERS: for lv in ["c", "lo", "hi"]: t = dict(location=loc, provenance=provenance(loc, name), quantity=name, level={"c": "central", "lo": "low", "hi": "high"}[lv]) for D in TABLE_SIZES: t[f"{D} m"] = float(at(sub, D)[key + lv]) trows.append(t) T = pd.DataFrame(trows) T.to_csv(os.path.join(DATA, "completeness_table.csv"), index=False) TH = pd.DataFrame(thr) TH.to_csv(os.path.join(DATA, "detection_thresholds.csv"), index=False) # ---- probability that a single artefact is in the archive / catalogued / recognised ---- def aggregate(w, scheme): out = [] tot = sum(w.values()) for p in ALBEDOS: for D in SIZES: r = dict(scheme=scheme, albedo=p, size_m=D) sub = M[(M.albedo == p) & (np.isclose(M.size_m, D))].set_index("location") for key, _ in TIERS: for lv in ["lo", "c", "hi"]: r[(key + lv).replace("C_", "P_")] = float(sum(w[loc] / tot * sub.loc[loc, key + lv] for loc in names)) # the same without the three surface rows (their weight is counted as not found) r["P_cat_c_orbits_only"] = float(sum(w[loc] / tot * sub.loc[loc, "C_cat_c"] for loc in names if loc not in SURFACES)) out.append(r) return out wrows, prow = [], [] for scheme, w in WEIGHTS.items(): tot = sum(w.values()) for loc in names: wrows.append(dict(scheme=scheme, location=loc, weight=w[loc] / tot)) prow += aggregate(w, scheme) P = pd.DataFrame(prow) # largest number of such artefacts compatible (95%) with none having been recognised: (1-P)^N = 0.05. # N95_cat is the same number if cataloguing implied recognition (the limit of perfect recognition): it is NOT # a statement that no artefact is in the catalogues, which nobody can make. for lv in ["lo", "c", "hi"]: for key in ["rec", "cat"]: pr = np.clip(P[f"P_{key}_{lv}"].values, 1e-300, 1 - 1e-12) P[f"N95_{key}_{lv}"] = np.log(0.05) / np.log1p(-pr) P.to_csv(os.path.join(DATA, "found_probability.csv"), index=False) pd.DataFrame(wrows).to_csv(os.path.join(DATA, "weights.csv"), index=False) # ---- comparison with Pan et al. (2026) under their residence weights ---- crow = [] for name in PAN_WEIGHTS: A = pd.DataFrame(aggregate(pan_weights(name), name)) A = A[A.albedo == 0.2] for D in (10, 100, 1000): a_ = A.iloc[(np.abs(np.log10(A.size_m.values) - np.log10(D))).argmin()] crow.append(dict(residence_model=name, size_m=D, pan_p_det=PAN_PDET[name][D], pan_lambda95=PAN_LAMBDA95[name][D], this_P_archive=a_.P_arch_c, this_P_catalogued=a_.P_cat_c, this_P_catalogued_lo=a_.P_cat_lo, this_P_catalogued_hi=a_.P_cat_hi, this_P_recognised=a_.P_rec_c, this_P_recognised_lo=a_.P_rec_lo, this_P_recognised_hi=a_.P_rec_hi, this_N95_catalogued=np.log(0.05) / np.log1p(-min(a_.P_cat_c, 1 - 1e-12)), this_N95_recognised=np.log(0.05) / np.log1p(-max(a_.P_rec_c, 1e-300)))) CMP = pd.DataFrame(crow) CMP.to_csv(os.path.join(DATA, "pan2026_comparison.csv"), index=False) json.dump(dict(parameters_low_central_high=m.params, largest_human_deep_space_object_m=m.d_human_max, lesnikowski_fraction_of_lunar_surface=m.lesnikowski_frac, measured_arc_and_characterisation=m.arc.to_dict(orient="records"), f_phys_main_belt_catalogue_average=m.f_belt_phys, f_phys_trojans_catalogue_average=m.f_tjn_phys, f_phys_tno_catalogue_average=m.f_tno_phys, f_phys_by_H={g: dict(H=list(map(float, v[0])), fraction=list(map(float, v[1]))) for g, v in m.phys_by_H.items()}, f_phys_main_belt_at_1km_10km=[float(m.f_phys("main_belt", H_from_D(1e3, 0.2))), float(m.f_phys("main_belt", H_from_D(1e4, 0.2)))], note="low/high set every parameter low or high together: scenario envelopes, not confidence intervals"), open(os.path.join(DATA, "recognition_model.json"), "w"), indent=1) # ---- print ---- pd.set_option("display.width", 260) pd.set_option("display.max_colwidth", 48) pd.set_option("display.max_rows", 400) fmt = lambda x: ("%.2f" % x if x >= 0.095 else ("%.0e" % x if x > 0 else "0")) if isinstance(x, float) else x # noqa: E731 print("[synth] measured recognition inputs:\n", m.arc.round(3).to_string(index=False)) print(f"[synth] largest human object in heliocentric orbit: {m.d_human_max:.1f} m equivalent diameter; " f"Lesnikowski et al. fraction of lunar surface: {m.lesnikowski_frac:.2e}; main-belt characterised fraction at 1 km / 10 km: " f"{m.f_phys('main_belt', H_from_D(1e3, 0.2)):.3f} / {m.f_phys('main_belt', H_from_D(1e4, 0.2)):.3f}") print("[synth] completeness versus size, albedo 0.2 (central, low, high):") print(T.to_string(index=False, formatters={c: fmt for c in T.columns if c.endswith(" m")})) print("[synth] smallest size with probability >= 50% / 90% (central; low-high):") print(TH[["location", "albedo", "D50_arch_c_m", "D50_cat_c_m", "D50_cat_lo_m", "D50_cat_hi_m", "D90_cat_c_m", "D50_recognised_c_m", "plateau_cat_c", "plateau_cat_lo", "plateau_cat_hi"]].to_string(index=False, float_format=lambda x: "%.3g" % x)) for alb in (0.2, 0.01): print(f"[synth] probability that a single artefact is in archived data / catalogued / recognised (albedo {alb}):") for scheme in WEIGHTS: sx = P[(P.scheme == scheme) & (P.albedo == alb)] sx = sx[sx.size_m.isin(TABLE_SIZES)] print(scheme) print(sx[["size_m", "P_arch_c", "P_cat_lo", "P_cat_c", "P_cat_hi", "P_cat_c_orbits_only", "P_rec_lo", "P_rec_c", "P_rec_hi", "N95_cat_c", "N95_rec_c", "N95_rec_hi"]].to_string(index=False, float_format=lambda x: "%.3g" % x)) print("[synth] comparison with Pan et al. (2026) under their residence weights:") print(CMP.to_string(index=False, float_format=lambda x: "%.3g" % x)) if __name__ == "__main__": main()