"""Detection probability of an inert object parked in near-Earth space, with the survey simulator calibrated in s02, plus tests against real objects. A. Synthetic parked populations (500 objects each): Sun-Earth L4/L5 (parked near the point; wide tadpoles), Earth horseshoes, quasi-satellites, Earth-like heliocentric 'drifters' (the orbits of spent upper stages), Earth-Moon L4/L5, high Earth orbits (60,000-300,000 km) and distant geocentric orbits (0.6-1.2 million km). -> P(discovered by 2026.8) versus absolute magnitude H for five survey-depth variants; published dedicated Earth-Trojan searches are added for the L4/L5 classes. Revision: (i) the elongation weight is the one measured in s02 for each depth variant (central, and its 95% lower / upper bounds), and the co-orbital classes are also run with the assumed ramps starting at 45, 60, 65 and 70 deg (sensitivity table); (ii) the geocentric classes are also run with a theoretical trailing-loss law, because the empirical law is fitted below 10 deg/day and these classes move at 2.5-214 deg/day. B. Tests on real objects: every catalogued asteroid on an Earth-like or co-orbital orbit (JPL Horizons daily ephemerides, s03). For each survey variant: the fraction of real discoveries the model calls impossible, the probability-integral-transform (PIT) of the discovery date, and - new, suggested by the referee - the PIT of the apparent magnitude on the night of first observation (a limit that is too deep passes the first two tests but fails this one). C. Human hardware in heliocentric orbit (GCAT hcocat): expected number of serendipitous survey detections 1998-2026 versus the number found. An order-of-magnitude check on the total only. D. Minimoon correction of the geocentric classes: the factor kappa on the nightly probability is SOLVED so that the simulator returns the two known minimoons (revision: the first version used 2/N, which is not self-consistent), with a range from the Poisson error on two objects combined with an assumed factor-4 uncertainty of the minimoon population. The lower end of the range is floored by human hardware: kappa may not be so small that the simulator contradicts (95%) the fraction of catalogued human objects in high Earth orbit that asteroid surveys have picked up (E). Projection for ten years at Rubin/LSST depth. E. (new) Check of the high-Earth-orbit class against catalogued human objects in orbit with apogee above 60,000 km (GCAT deepcat): how many are in B. Gray's list of objects whose orbits are maintained from asteroid-survey astrometry. Outputs: data/parked_pdet_curves.csv, data/parked_thresholds.csv, data/parked_elong_sensitivity.csv, data/minimoon_test.csv, data/geocentric_kappa.csv, data/parked_lsst_projection.csv, data/real_coorbital_test.csv, data/real_coorbital_objects.csv, data/elongation_check.csv, data/discovery_magnitude_test.csv, data/discovery_magnitude_hist.csv, data/human_hardware_prediction.csv, data/human_hardware_summary.json, data/heo_hardware_check.csv """ import json import os import re import numpy as np import pandas as pd import simlib as S from common import DATA, RAW, SEED, D_from_H, H_from_D, K_GAUSS, helio_xyz from s03_fetch_earthlike_ephem import designation_window N_CLASS = 500 H_GRID = np.arange(12.0, 38.01, 0.25) VARIANTS = ["base", "shallow", "deep", "vdeep", "hard22p5"] ALBEDOS = [0.2, 0.05, 0.01] def rayleigh_trunc(rng, sigma_deg, max_deg, n): x = rng.rayleigh(sigma_deg, n * 3) x = x[x < max_deg][:n] return np.radians(x) HELIO_CLASSES = ["Sun-Earth L4/L5, parked near the point (libration < 10 deg)", "Sun-Earth L4/L5, wide tadpole orbits", "Earth horseshoe co-orbitals", "Earth quasi-satellites", "Earth-like heliocentric orbits (0.9-1.1 au)"] # elongation-weight models: the measured table of each depth variant (central, 95% lower, 95% upper) and # assumed ramps (the same for every variant) RAMPS = {"ramp45": dict(e0=45.0, e1=75.0), "ramp60": dict(e0=60.0, e1=90.0), "ramp65": dict(e0=65.0, e1=95.0), "ramp70": dict(e0=70.0, e1=100.0)} LEVELS = {"measured": "c", "measured_lo": "lo", "measured_hi": "hi"} def weight_functions(ramps=True): wf = {(v, lv): S.make_wfun(level=lv, variant=v) for v in VARIANTS for lv in LEVELS.values()} if ramps: wf.update({name: S.make_wfun(**kw) for name, kw in RAMPS.items()}) return wf def helio_class_hists(rng): """Geometry histograms of the five heliocentric parked classes for every elongation-weight model. Returns {class: {(variant, level) or ramp name: hist}}, track statistics.""" n = N_CLASS two_pi = 2 * np.pi out, track_stats = {}, {} # --- Earth Trojans (L4; L5 is the mirror image and equivalent for visibility). phi_turn is the # turning point of the libration nearest Earth; L4 itself is at 60 deg. for name, lo, hi in [(HELIO_CLASSES[0], 50.0, 59.5), (HELIO_CLASSES[1], 30.0, 50.0)]: phi_turn = np.radians(rng.uniform(lo, hi, n)) PHI, X = S.coorbital_tracks(phi_turn, rng) e = rng.uniform(0, 0.1, n) inc = rayleigh_trunc(rng, 6.0, 30.0, n) Om, varpi = rng.uniform(0, two_pi, n), rng.uniform(0, two_pi, n) out[name] = S.hist_coorbital(PHI, X, e, inc, Om, varpi, wfun=weight_functions()) track_stats[name] = dict(phi_min_deg=float(np.degrees(PHI.min())), phi_max_deg=float(np.degrees(PHI.max()))) # --- horseshoes phi_turn = np.radians(rng.uniform(8.0, 22.0, n)) PHI, X = S.coorbital_tracks(phi_turn, rng) e = rng.uniform(0, 0.1, n) inc = rayleigh_trunc(rng, 6.0, 30.0, n) Om, varpi = rng.uniform(0, two_pi, n), rng.uniform(0, two_pi, n) out[HELIO_CLASSES[2]] = S.hist_coorbital(PHI, X, e, inc, Om, varpi, wfun=weight_functions()) dphi = np.abs(((PHI + np.pi) % two_pi) - np.pi) track_stats["horseshoe"] = dict(frac_time_within_30deg_of_earth=float((dphi < np.radians(30)).mean()), frac_time_within_60deg_of_earth=float((dphi < np.radians(60)).mean())) # --- quasi-satellites: guiding centre librates about Earth e = rng.uniform(0.05, 0.25, n) amp = rng.uniform(0, 0.5, n) * e per = rng.uniform(50, 300, n) * 365.25 psi = rng.uniform(0, two_pi, n) PHI = amp[:, None] * np.sin(two_pi * S.T_DAYS[None, :] / per[:, None] + psi[:, None]) dPHI = amp[:, None] * two_pi / per[:, None] * np.cos(two_pi * S.T_DAYS[None, :] / per[:, None] + psi[:, None]) X = -dPHI / (1.5 * K_GAUSS) inc = rayleigh_trunc(rng, 5.0, 25.0, n) Om, varpi = rng.uniform(0, two_pi, n), rng.uniform(0, two_pi, n) out[HELIO_CLASSES[3]] = S.hist_coorbital(PHI, X, e, inc, Om, varpi, wfun=weight_functions()) del PHI, X, dPHI # --- Earth-like drifters a = rng.uniform(0.9, 1.1, n * 2) a = a[np.abs(a - 1) > 0.006][:n] e = rng.uniform(0, 0.1, n) inc = rayleigh_trunc(rng, 2.0, 10.0, n) Om, w_, M0 = rng.uniform(0, two_pi, n), rng.uniform(0, two_pi, n), rng.uniform(0, two_pi, n) out[HELIO_CLASSES[4]] = S.hist_kepler(a, e, inc, Om, w_, M0, wfun=weight_functions(ramps=False)) return out, track_stats GEO_CLASSES = ["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)"] def geo_class_hists(rng): """Geometry histograms of the three geocentric classes for the two trailing-loss laws. extra_mag: allowance for moonlight / scattered light near the Moon and Earth. Returns {class: {trailing law: {variant: hist}}} (each variant with its own measured elongation weight) and the sky rates (deg/day) of each class.""" n = N_CLASS two_pi = 2 * np.pi spec = { GEO_CLASSES[0]: (10 ** rng.uniform(np.log10(6.0e4), np.log10(3.0e5), n), np.radians(rng.uniform(0, 30, n)), 0.5), GEO_CLASSES[1]: (np.full(n, 384400.0), np.radians(rng.uniform(0, 8, n)), 1.0), GEO_CLASSES[2]: (rng.uniform(0.6e6, 1.2e6, n), np.radians(rng.uniform(0, 30, n)), 0.0), } out, rates = {}, {} for name, (R, inc, extra) in spec.items(): ph = rng.uniform(0, two_pi, n) out[name] = {tr: S.hist_geocentric(R, inc, ph, extra_mag=extra, trail=tr, wfun={v: S.make_wfun(variant=v) for v in VARIANTS}) for tr in ("neomod2", "theory")} rate = 360.0 / (27.32166 * (R / 384400.0) ** 1.5) rates[name] = (float(rate.min()), float(rate.max())) return out, rates # Dedicated Earth-Trojan searches (published efficiencies; step functions of H). # L4: Lu et al. 2026 (WFST): 33.24% of the stable L4 cloud covered; cumulative 95% upper limit N(H < 19.1) <~ 19 # (abstract), differential limit 14.5 at H = 19.1 (their sect. 4.2) -> efficiency 3/19 ... 3/14.5 = # 0.16-0.21; 0.18 adopted for H <~ 19.4; # Cambioni et al. 2018 (OSIRIS-REx): upper limit 73 at H = 20.5 -> ~0.04; # Lifset et al. 2021 (DECam, 5.7 deg2): ~0.006 down to H ~ 21.8. # L5: Markwardt et al. 2020 (DECam): 24% of the L5 cloud covered; their own limit of 15-18 objects for # H <~ 19.7 implies an efficiency of 3/18 ... 3/15 = 0.17-0.20; 0.18 adopted (the first version of this # study used 0.12, which was too low); recovery 12.8% at H = 20.4 -> 0.24 x 0.128 = 0.03. def dedicated_trojan_eff(H): H = np.asarray(H, float) l4 = 1 - (1 - np.where(H <= 19.4, 0.18, 0.0)) * (1 - np.where(H <= 20.5, 0.04, 0.0)) * (1 - np.where(H <= 21.8, 0.006, 0.0)) l5 = np.where(H <= 19.7, 0.18, np.where(H <= 20.4, 0.031, 0.0)) return 0.5 * (l4 + l5) POISSON2_95 = (0.242, 7.225) # 95% interval for the mean of a Poisson variable observed as 2 MINIMOON_POP_FACTOR = 4.0 # assumed 95% uncertainty factor of the minimoon capture rate, either way def pcurve(h, f, k=1.0): return np.array([S.p_detect(h, H, f["vlim"], np.array(f["p"]) * k, f["rho"], vhard=f["vhard"]).mean() for H in H_GRID]) def Hat(P, level): ok = np.where(P >= level)[0] return float(H_GRID[ok.max()]) if len(ok) else np.nan def P_at_D(P, D_m, p_alb=0.2): return float(np.interp(H_from_D(D_m, p_alb), H_GRID, P)) def helio_curves(hh, fits): """Detection curves of the heliocentric classes. Columns: p_detect (measured central weight of the variant), p_wlo / p_whi (95% bounds of the measured weight); and the sensitivity table over weight models.""" rows, thr, sens = [], [], [] for cname, by_w in hh.items(): ded = dedicated_trojan_eff(H_GRID) if cname.startswith("Sun-Earth L4/L5") else np.zeros(len(H_GRID)) for v in VARIANTS: f = fits[v] P = {wm: pcurve(by_w[(v, lv)], f) for wm, lv in LEVELS.items()} P.update({wm: pcurve(by_w[wm], f) for wm in RAMPS if wm in by_w}) T = {wm: 1 - (1 - P[wm]) * (1 - ded) for wm in P} for i, H in enumerate(H_GRID): rows.append(dict(location=cname, variant=v, H=H, p_detect=T["measured"][i], p_detect_surveys_only=P["measured"][i], p_wlo=T["measured_lo"][i], p_whi=T["measured_hi"][i])) row = dict(location=cname, variant=v, kappa=1.0, P_max=float(T["measured"].max()), P_max_surveys_only=float(P["measured"].max()), H_90=Hat(T["measured"], 0.9), H_50=Hat(T["measured"], 0.5), H_10=Hat(T["measured"], 0.1), H_50_wlo=Hat(T["measured_lo"], 0.5), H_50_whi=Hat(T["measured_hi"], 0.5)) for tag in ["H_90", "H_50", "H_10", "H_50_wlo", "H_50_whi"]: for p_alb in ALBEDOS: row[f"D_m_{tag[2:]}_p{p_alb}"] = float(D_from_H(row[tag], p_alb)) if np.isfinite(row[tag]) else np.nan thr.append(row) for wm in P: sens.append(dict(location=cname, variant=v, weight_model=wm, P_max_surveys_only=float(P[wm].max()), P_max_with_dedicated=float(T[wm].max()), P_10km=P_at_D(T[wm], 1e4), P_1km=P_at_D(T[wm], 1e3), P_300m=P_at_D(T[wm], 300.0), P_100m=P_at_D(T[wm], 100.0), H_50=Hat(T[wm], 0.5), D50_m_p0p2=float(D_from_H(Hat(T[wm], 0.5), 0.2)) if np.isfinite(Hat(T[wm], 0.5)) else np.nan)) return pd.DataFrame(rows), pd.DataFrame(thr), pd.DataFrame(sens) # ---------------------------------------------------------------------------- # B. real Earth-like / co-orbital objects # ---------------------------------------------------------------------------- def real_object_test(fits): z = np.load(os.path.join(DATA, "earthlike_geom.npz")) g = z["g_centimag"].astype(float) / 100.0 gn = z["g_notrail_centimag"].astype(float) / 100.0 band = z["band"].astype(int) el = z["elong_decideg"].astype(float) / 10.0 dec_ok = z["dec_ok"].astype(float) obj = pd.read_csv(os.path.join(DATA, "earthlike_objects.csv")) nday = g.shape[1] t_year = 1998.0 + np.arange(nday) / 365.25 era = np.clip(np.searchsorted(S.ERA_EDGES, t_year, side="right") - 1, 0, S.N_ERA - 1) wins = [] for _, r in obj.iterrows(): des = r.prov if isinstance(r.prov, str) else r.pdes wins.append(designation_window(des)) obj["win0"] = [x[0] if x else np.nan for x in wins] obj["win1"] = [x[1] if x else np.nan for x in wins] use = obj.win0.notna() & (obj.win0 >= 1998.0) & (obj.win1 <= 2026.75) # night of first observation (the discovery night when it lies inside the designation half-month) fo = pd.to_datetime(obj.first_obs, errors="coerce") obj["k_first"] = (fo - pd.Timestamp("1998-01-01")).dt.days obj["t_first"] = 1998.0 + obj.k_first / 365.25 same = use & (obj.t_first >= obj.win0 - 4 / 365.25) & (obj.t_first <= obj.win1 + 4 / 365.25) & (obj.k_first >= 0) & (obj.k_first < nday) n_precov = int((use & (obj.t_first < obj.win0 - 1.0)).sum()) J = np.where(same.values)[0] kf = obj.k_first.values.astype(int) V_real = np.array([obj.H[j] + g[j, kf[j]] for j in J]) V_real_nt = np.array([obj.H[j] + gn[j, kf[j]] for j in J]) E_real = np.array([el[j, kf[j]] for j in J]) yr_real = obj.t_first.values[J] print(f"[parked] real Earth-like / co-orbital asteroids: {int(use.sum())} usable; first observation inside the designation " f"half-month for {len(J)}; first observation more than a year before designation for {n_precov} " f"({100 * n_precov / use.sum():.1f}%)") qs = [0.1, 0.5, 0.9, 0.99] print(f"[parked] REAL V on the night of first observation (with trailing term): quantiles 10/50/90/99% = " f"{np.round(np.quantile(V_real, qs), 2)}; fainter than 22.0: {100 * (V_real > 22).mean():.1f}%, than 22.5: " f"{100 * (V_real > 22.5).mean():.1f}%; without trailing term: median {np.median(V_real_nt):.2f}; " f"elongation below 90 deg: {100 * (E_real < 90).mean():.1f}%, below 75 deg: {100 * (E_real < 75).mean():.1f}%") vm_rows = [dict(variant="REAL (first-observation night)", n=len(J), V_median=float(np.median(V_real)), V_p90=float(np.quantile(V_real, 0.9)), frac_fainter_22p0=float((V_real > 22).mean()), frac_fainter_22p5=float((V_real > 22.5).mean()), pit_mean=np.nan, pit_se=np.nan, V_median_1998_2010=float(np.median(V_real[yr_real < 2011])), V_median_2011_2020=float(np.median(V_real[(yr_real >= 2011) & (yr_real < 2021)])), V_median_2021_2026=float(np.median(V_real[yr_real >= 2021])))] summ, per_obj = [], {} v_bins = np.arange(16.0, 24.51, 0.5) vh_rows = [dict(variant="REAL", V_lo=a_, V_hi=b_, n=float(x_)) for a_, b_, x_ in zip(v_bins[:-1], v_bins[1:], np.histogram(np.clip(V_real, 16.01, 24.49), bins=v_bins)[0])] el_bins = np.array([0, 60, 75, 90, 105, 120, 135, 150, 165, 181]) el_rows = [] h_real_true = np.histogram(E_real, bins=el_bins)[0].astype(float) for v, f in fits.items(): vlim, p, rho, vhard = map(np.array, (f["vlim"], f["p"], f["rho"], f["vhard"])) w = S.elong_weight(np.radians(el), variant=v) * dec_ok # measured weight of this variant (base table if none) u_list, imp_mag, imp_geom, P_now = [], 0, 0, {} h_mod = np.zeros(len(el_bins) - 1) pit_v, vals, wts = [], [], [] for j in np.where(use.values)[0]: V = obj.H[j] + g[j] Sg = np.where(V <= vhard[era], 1.0 / (1.0 + np.exp(np.clip((V - vlim[era]) / S.LOGISTIC_W, -60, 60))), 0.0) pn = p[era] * np.where(band[j] == 1, rho[era], 1.0) * Sg haz = -w[j] * np.log1p(-pn) Lam = np.concatenate([[0.0], np.cumsum(haz)]) k0 = int(np.clip((obj.win0[j] - 1998.0) * 365.25 - 3, 0, nday)) # +-3 d tolerance k1 = int(np.clip((obj.win1[j] - 1998.0) * 365.25 + 3, 0, nday)) lam_w = Lam[k1] - Lam[k0] tot = Lam[-1] P_now[j] = 1 - np.exp(-tot) if lam_w <= 0: # impossible under this variant: because of geometry cuts or because too faint? if (w[j, k0:k1] > 0).any(): imp_mag += 1 else: imp_geom += 1 continue km = (k0 + k1) // 2 u_list.append((1 - np.exp(-Lam[km])) / (1 - np.exp(-tot))) if same.values[j]: # model-predicted distribution of the discovery night for this object (over the whole survey) pdf = np.exp(-Lam[:-1]) * (1 - np.exp(-haz)) pdf = pdf / pdf.sum() h_mod += np.histogram(el[j], bins=el_bins, weights=pdf)[0] o = np.argsort(V) pit_v.append(float(np.interp(obj.H[j] + g[j, kf[j]], V[o], np.cumsum(pdf[o])))) vals.append(V) wts.append(pdf) u = np.array(u_list) n_use = int(use.sum()) us = np.sort(u) ks = float(np.max(np.abs(us - (np.arange(len(us)) + 0.5) / len(us)))) if len(us) else np.nan summ.append(dict(variant=v, n_objects=n_use, n_impossible_too_faint=imp_mag, n_impossible_geometry=imp_geom, frac_impossible_too_faint=imp_mag / n_use, frac_impossible_geometry=imp_geom / n_use, mean_u=float(u.mean()), median_u=float(np.median(u)), frac_u_gt_0p9=float((u > 0.9).mean()), ks_distance=ks, ks_critical_5pct=1.358 / np.sqrt(len(u)), mean_P_discovered_by_2026=float(np.mean(list(P_now.values()))))) per_obj[v] = P_now pit_v = np.array(pit_v) vals, wts = np.concatenate(vals), np.concatenate(wts) o = np.argsort(vals) cw = np.cumsum(wts[o]) / wts.sum() vm_rows.append(dict(variant=v, n=len(pit_v), V_median=float(vals[o][np.searchsorted(cw, 0.5)]), V_p90=float(vals[o][np.searchsorted(cw, 0.9)]), frac_fainter_22p0=float(wts[vals > 22.0].sum() / wts.sum()), frac_fainter_22p5=float(wts[vals > 22.5].sum() / wts.sum()), pit_mean=float(pit_v.mean()), pit_se=float(pit_v.std() / np.sqrt(len(pit_v))), vlim_last_era=float(vlim[-1]), nightly_p_last_era=float(p[-1]), rms_logit_history_fit=f["rms_logit"])) hv = np.histogram(np.clip(vals, 16.01, 24.49), bins=v_bins, weights=wts)[0] vh_rows += [dict(variant=v, V_lo=a_, V_hi=b_, n=float(x_)) for a_, b_, x_ in zip(v_bins[:-1], v_bins[1:], hv)] for a_, b_, x_, y_ in zip(el_bins[:-1], el_bins[1:], h_real_true, h_mod): el_rows.append(dict(variant=v, elong_lo=a_, elong_hi=b_, real_first_observations=x_, model_discoveries=y_)) pd.DataFrame(el_rows).to_csv(os.path.join(DATA, "elongation_check.csv"), index=False) vm = pd.DataFrame(vm_rows) vm.to_csv(os.path.join(DATA, "discovery_magnitude_test.csv"), index=False) pd.DataFrame(vh_rows).to_csv(os.path.join(DATA, "discovery_magnitude_hist.csv"), index=False) o2 = obj[use].copy() o2["first_obs_in_designation_halfmonth"] = same[use].values o2["V_first_obs_with_trailing"] = [obj.H[j] + g[j, kf[j]] if same.values[j] else np.nan for j in np.where(use.values)[0]] o2["elong_first_obs"] = [el[j, kf[j]] if same.values[j] else np.nan for j in np.where(use.values)[0]] for v in fits: o2["P_disc_by_2026_" + v] = [per_obj[v][j] for j in np.where(use.values)[0]] return pd.DataFrame(summ), o2, vm, n_precov # ---------------------------------------------------------------------------- # C. human hardware in heliocentric orbit # ---------------------------------------------------------------------------- MONTHS = {m: k + 1 for k, m in enumerate(["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"])} def gcat_date(s): m = re.match(r"\s*(\d{4})\s+([A-Za-z]{3})\s*(\d{1,2})?", str(s)) if not m: return np.nan d = int(m.group(3)) if m.group(3) else 15 return int(m.group(1)) + (MONTHS.get(m.group(2), 6) - 1) / 12.0 + (d - 1) / 365.25 def load_hcocat(): lines = open(os.path.join(RAW, "hcocat.tsv"), encoding="utf-8", errors="replace").read().split("\n") hdr = lines[0].lstrip("#").split("\t") rows = [l.split("\t") for l in lines[1:] if l and not l.startswith("#")] df = pd.DataFrame(rows, columns=hdr[:len(rows[0])]).apply(lambda c: c.str.strip()) for c in ["Length", "Diameter", "Span", "Perigee", "Apogee", "Inc", "DryMass"]: df[c] = pd.to_numeric(df[c], errors="coerce") df["t_start"] = [gcat_date(x) for x in df.SDate] df["t_launch"] = [gcat_date(x) for x in df.LDate] return df def human_hardware(fits, rng, n_mc=40): df = load_hcocat() n_rows = len(df) cur = df[(df.Status == "O") & (df.OpOrbit == "HCO")].copy() # phases that are current, in heliocentric orbit cur = cur.sort_values("t_start").drop_duplicates("JCAT", keep="last") cur = cur[cur.Type.str[0].isin(["P", "R"])] # payloads and rocket stages (not small components) cur = cur[cur.Perigee.notna() & cur.Apogee.notna() & (cur.Perigee > 0.3) & (cur.Apogee < 3.5) & (cur.Apogee >= cur.Perigee)] L = cur.Length.fillna(0).clip(lower=0) Dm = cur.Diameter.fillna(0).clip(lower=0) # mean projected area of a randomly oriented cylinder = (surface area)/4 area = (np.pi * Dm * L + np.pi * Dm ** 2 / 2) / 4 area = np.where(area > 0, area, np.pi * (np.maximum(L, Dm) / 2) ** 2) cur["D_eq_m"] = np.sqrt(4 * area / np.pi) cur = cur[cur.D_eq_m >= 1.0].copy() cur["a"] = (cur.Perigee + cur.Apogee) / 2 cur["e"] = (cur.Apogee - cur.Perigee) / (cur.Apogee + cur.Perigee) cur["earthlike"] = (cur.a > 0.9) & (cur.a < 1.1) & (cur.e < 0.12) cur["anchored"] = (cur.Parent == "Earth") | ((cur.t_start - cur.t_launch).abs() < 0.3) n = len(cur) a = cur.a.values e = np.clip(cur.e.values, 1e-4, 0.95) inc = np.radians(cur.Inc.fillna(1.0).clip(0, 60).values) Hmag = H_from_D(cur.D_eq_m.values, 0.2) ts = np.where(np.isfinite(cur.t_start.values), cur.t_start.values, cur.t_launch.values) anchored = cur.anchored.values acc = {v: np.zeros(n) for v in VARIANTS} for _ in range(n_mc): aj = a * (1 + rng.normal(0, 0.002, n)) # catalogue orbits are approximate # true anomaly at which r = 1 au (departure from Earth); perihelion if q>1, aphelion if Q<1 cosnu = np.clip((aj * (1 - e ** 2) - 1.0) / e, -1, 1) nu0 = np.arccos(cosnu) * rng.choice([-1.0, 1.0], n) E0 = 2 * np.arctan(np.sqrt((1 - e) / (1 + e)) * np.tan(nu0 / 2)) M_dep = E0 - e * np.sin(E0) t_dep_days = (ts - S.T0_YEAR) * 365.25 lamE_dep = S.LAM_E0 + K_GAUSS * t_dep_days varpi = lamE_dep - nu0 Om = lamE_dep.copy() nmot = K_GAUSS / aj ** 1.5 M0 = M_dep - nmot * t_dep_days # mean anomaly at T0 rnd = ~anchored | ~np.isfinite(ts) M0 = np.where(rnd, rng.uniform(0, 2 * np.pi, n), M0) varpi = np.where(rnd, rng.uniform(0, 2 * np.pi, n), varpi) Om = np.where(rnd, rng.uniform(0, 2 * np.pi, n), Om) # the object counts only from one year after departure (before that it was a tracked, known launch) h = S.hist_kepler(aj, e, inc, Om, varpi - Om, M0, t_start_days=np.nan_to_num(t_dep_days, nan=-1e9) + 365.25, wfun={v: S.make_wfun(variant=v) for v in VARIANTS}) for v in VARIANTS: f = fits[v] acc[v] += S.p_detect(h[v], Hmag, f["vlim"], f["p"], f["rho"], vhard=f["vhard"]) out = cur[["JCAT", "Type", "Name", "LDate", "Parent", "Length", "Diameter", "D_eq_m", "Perigee", "Apogee", "Inc", "a", "e", "earthlike", "anchored"]].copy() out["H_p0.2"] = Hmag for v in VARIANTS: out["P_detect_" + v] = acc[v] / n_mc summ = dict(n_catalogue_rows=int(n_rows), n_objects_current_heliocentric_ge1m=int(n), n_earthlike=int(cur.earthlike.sum()), n_ge3m=int((cur.D_eq_m >= 3).sum()), n_earthlike_ge3m=int((cur.earthlike & (cur.D_eq_m >= 3)).sum())) for v in VARIANTS: summ[f"expected_detections_{v}"] = float(out["P_detect_" + v].sum()) summ[f"expected_detections_ge3m_{v}"] = float(out.loc[out.D_eq_m >= 3, "P_detect_" + v].sum()) summ[f"expected_detections_earthlike_{v}"] = float(out.loc[out.earthlike, "P_detect_" + v].sum()) return out, summ # ---------------------------------------------------------------------------- # D. minimoon correction of the geocentric classes, and an LSST projection # ---------------------------------------------------------------------------- ERA_LEN = np.diff(np.minimum(S.ERA_EDGES, S.T1_YEAR)) # years in each era KGRID = 10 ** np.linspace(-4.0, 0.0, 81) def minimoon_rates(): """Temporarily captured natural satellites: Bolin et al. (2014), quoting Granvik et al. (2012): about one 1-m object (H ~ 33) is in orbit at any time, the mean capture lasts ~290 days, and a 5-10 m object (H 29.5-28) is captured about once per decade. A power law through these two statements gives the capture rate R(nk", hk, ln1m) # ln(1 - P) over the whole era out += ERA_LEN[k] * dr * (1 - np.exp(L * tau / ERA_LEN[k])).mean(axis=0) return out def solve_kappa(N_of_k, target, kgrid=KGRID): """kappa at which the expected number of discoveries equals target (1 if never reached).""" if N_of_k[-1] <= target: return 1.0 if N_of_k[0] >= target: return float(kgrid[0]) return float(10 ** np.interp(np.log10(target), np.log10(N_of_k), np.log10(kgrid))) def kappa_targets(): lo = 2.0 / np.exp(np.hypot(np.log(2.0 / POISSON2_95[0]), np.log(MINIMOON_POP_FACTOR))) hi = 2.0 * np.exp(np.hypot(np.log(POISSON2_95[1] / 2.0), np.log(MINIMOON_POP_FACTOR))) return 2.0, lo, hi def geo_curves(gh, fits, hw, hw_stats): """Detection curves of the geocentric classes with the minimoon correction. Columns: p_detect (empirical trailing law, central kappa), p_klo / p_khi (kappa range), p_theory, p_theory_klo, p_theory_khi (theoretical trailing law), p_detect_uncalibrated (kappa = 1).""" t_c, t_lo, t_hi = kappa_targets() mm_classes = [GEO_CLASSES[2], GEO_CLASSES[1]] krows, mrows, rows, thr = [], [], [], [] kap = {} for tr in ("neomod2", "theory"): for v in VARIANTS: f = fits[v] Ns = [] for c in mm_classes: N = minimoon_expected(gh[c][tr][v], f) Ns.append(N) mrows.append(dict(orbit_class=c, trailing_law=tr, variant=v, expected_discoveries_2005_2026_kappa1=float(N[-1]), known_minimoons=2)) N = np.mean(Ns, axis=0) k_c, k_mm, k_hi = solve_kappa(N, t_c), solve_kappa(N, t_lo), solve_kappa(N, t_hi) # floor from human hardware: the lowest kappa must not contradict (at 95%) the fraction of catalogued # human objects in high Earth orbit that asteroid surveys have picked up k_hw, p_hw1 = kappa_hardware_floor(gh[GEO_CLASSES[0]][tr][v], f, hw, hw_stats["fraction_lower95"]) k_lo = max(k_mm, k_hw) k_c = max(k_c, k_lo) k_old = min(1.0, 2.0 / N[-1]) # rule of the first version N_old = float(10 ** np.interp(np.log10(k_old), np.log10(KGRID), np.log10(N))) kap[(tr, v)] = (k_c, k_lo, k_hi) krows.append(dict(trailing_law=tr, variant=v, expected_at_kappa1=float(N[-1]), kappa=k_c, kappa_lo=k_lo, kappa_hi=k_hi, kappa_lo_minimoons_only=k_mm, kappa_floor_hardware=k_hw, hardware_mean_P_at_kappa1=p_hw1, target_lo=t_lo, target_hi=t_hi, kappa_first_version_rule=k_old, expected_with_first_version_rule=N_old)) for cname in GEO_CLASSES: for v in VARIANTS: f = fits[v] k_c, k_lo, k_hi = kap[("neomod2", v)] q_c, q_lo, q_hi = kap[("theory", v)] h, ht = gh[cname]["neomod2"][v], gh[cname]["theory"][v] P, Plo, Phi, Praw = pcurve(h, f, k_c), pcurve(h, f, k_lo), pcurve(h, f, k_hi), pcurve(h, f, 1.0) Q, Qlo, Qhi = pcurve(ht, f, q_c), pcurve(ht, f, q_lo), pcurve(ht, f, q_hi) for i, H in enumerate(H_GRID): rows.append(dict(location=cname, variant=v, H=H, p_detect=P[i], p_detect_surveys_only=P[i], p_detect_uncalibrated=Praw[i], p_klo=Plo[i], p_khi=Phi[i], p_theory=Q[i], p_theory_klo=Qlo[i], p_theory_khi=Qhi[i])) row = dict(location=cname, variant=v, kappa=k_c, kappa_lo=k_lo, kappa_hi=k_hi, P_max=float(P.max()), P_max_surveys_only=float(P.max()), H_90=Hat(P, 0.9), H_50=Hat(P, 0.5), H_10=Hat(P, 0.1), H_50_uncalibrated=Hat(Praw, 0.5), H_50_klo=Hat(Plo, 0.5), H_50_khi=Hat(Phi, 0.5), H_50_theory=Hat(Q, 0.5), H_50_theory_klo=Hat(Qlo, 0.5), H_50_theory_khi=Hat(Qhi, 0.5), P_1m=P_at_D(P, 1.0), P_1m_klo=P_at_D(Plo, 1.0), P_1m_khi=P_at_D(Phi, 1.0), P_1m_theory=P_at_D(Q, 1.0), P_3m=P_at_D(P, 3.0), P_3m_klo=P_at_D(Plo, 3.0), P_10m_klo=P_at_D(Plo, 10.0), P_10m_theory_klo=P_at_D(Qlo, 10.0)) for tag in ["H_90", "H_50", "H_10", "H_50_uncalibrated", "H_50_klo", "H_50_khi", "H_50_theory", "H_50_theory_klo", "H_50_theory_khi"]: for p_alb in ALBEDOS: row[f"D_m_{tag[2:]}_p{p_alb}"] = float(D_from_H(row[tag], p_alb)) if np.isfinite(row[tag]) else np.nan thr.append(row) return pd.DataFrame(rows), pd.DataFrame(thr), pd.DataFrame(krows), pd.DataFrame(mrows), kap def lsst_projection(hists, fits, kappa_geo=1.0, years=10.0, vlim=24.0, vhard=24.7): """Ten years of a survey with the geometry and nightly success probability of the 2021-2026 era but limiting magnitude `vlim` (Rubin single-visit depth ~24.5, minus ~0.5 mag for trailing and the need for repeated detections). A scaling exercise, not a simulation of the LSST cadence.""" f = fits["base"] rows = [] vl = np.array(f["vlim"], float) vh = np.array(f["vhard"], float) vl2, vh2 = vl.copy(), vh.copy() vl2[-1], vh2[-1] = vlim, vhard for cname, h in hists.items(): pk = np.array(f["p"]) * (kappa_geo if cname in GEO_CLASSES else 1.0) for H in H_GRID: L_now = S.survival_log(h, H, vl, pk, f["rho"], vh).sum(axis=1) L_new = S.survival_log(h, H, vl2, pk, f["rho"], vh2)[:, -1] * years / ERA_LEN[-1] rows.append(dict(location=cname, H=H, p_detect_now=float((1 - np.exp(L_now)).mean()), p_detect_after_lsst=float((1 - np.exp(L_now + L_new)).mean()))) return pd.DataFrame(rows) # ---------------------------------------------------------------------------- # E. catalogued human objects in high Earth orbit versus survey astrometry # ---------------------------------------------------------------------------- def load_heo_hardware(): """GCAT deepcat objects still in orbit (status O) with apogee above 60,000 km and an equivalent diameter of at least 1 m, cross-matched by international designation with B. Gray's list of high-orbit artificial objects whose orbits he maintains from astrometry reported by asteroid surveys ('objects over which asteroid hunters routinely stumble'). Being on the list is a lower bound on 'picked up by surveys': operating spacecraft with operator ephemerides (e.g. MMS) need no such orbit.""" lines = open(os.path.join(RAW, "deepcat.tsv"), encoding="utf-8", errors="replace").read().split("\n") hdr = lines[0].lstrip("#").split("\t") rows = [l.split("\t") for l in lines[1:] if l and not l.startswith("#")] df = pd.DataFrame([r[:len(hdr)] for r in rows], columns=hdr).apply(lambda c: c.str.strip()) for c in ["Length", "Diameter", "Perigee", "Apogee", "Inc"]: df[c] = pd.to_numeric(df[c], errors="coerce") cur = df[(df.Status == "O") & (df.Primary == "Earth") & df.Type.str[0].isin(["P", "R"]) & (df.Apogee > 60000) & (df.Apogee < 1.5e6)].copy() L, Dm = cur.Length.fillna(0).clip(lower=0), cur.Diameter.fillna(0).clip(lower=0) area = (np.pi * Dm * L + np.pi * Dm ** 2 / 2) / 4 area = np.where(area > 0, area, np.pi * (np.maximum(L, Dm) / 2) ** 2) cur["D_eq_m"] = np.sqrt(4 * area / np.pi) cur = cur[cur.D_eq_m >= 1.0].copy() txt = open(os.path.join(RAW, "gray_tle_list.txt"), encoding="utf-8", errors="replace").read() ids = set(c for _, c in re.findall(r"# ID:\s+(\S+)\s+(\S+)", txt)) n_unid = len([c for c in ids if c.startswith("00000")]) def key(piece): m = re.match(r"(\d{4})-(\d{3})([A-Z]+)$", piece) return (m.group(1)[2:] + m.group(2) + m.group(3)) if m else None cur["cospar_key"] = [key(x) for x in cur.Piece] cur["in_gray_list"] = cur.cospar_key.isin(ids) n, kk = len(cur), int(cur.in_gray_list.sum()) from scipy.stats import beta stats = dict(n_objects=n, n_in_gray_list=kk, fraction=kk / n if n else np.nan, fraction_lower95=float(beta.ppf(0.025, kk, n - kk + 1)) if kk else 0.0, fraction_upper95=float(beta.ppf(0.975, kk + 1, n - kk)) if kk < n else 1.0, n_ids_in_gray_list=len(ids), n_unidentified_in_gray_list=n_unid, deepcat_rows=len(df)) return cur, stats KGRID_HW = 10 ** np.linspace(-4.0, 0.0, 33) def kappa_hardware_floor(h, f, hw, target): """Smallest factor kappa on the nightly probability for which the high-Earth-orbit class of the simulator returns, averaged over the catalogued human objects in such orbits (sizes hw.D_eq_m, albedo 0.2), a detection probability of at least `target` (the 95% lower bound of the fraction of those objects that are in Gray's survey-astrometry list). Also returns the mean probability at kappa = 1.""" Hs, cnt = np.unique(np.round(H_from_D(hw.D_eq_m.values, 0.2), 3), return_counts=True) mean_p = np.array([sum(c * S.p_detect(h, H, f["vlim"], np.array(f["p"]) * k, f["rho"], vhard=f["vhard"]).mean() for H, c in zip(Hs, cnt)) / cnt.sum() for k in KGRID_HW]) if mean_p[-1] <= target: return 1.0, float(mean_p[-1]) if mean_p[0] >= target: return float(KGRID_HW[0]), float(mean_p[-1]) return float(10 ** np.interp(target, mean_p, np.log10(KGRID_HW))), float(mean_p[-1]) def heo_hardware_check(gh, fits, hw, stats): """Simulator probability (base variant, central kappa and lowest kappa) for each catalogued human object in high Earth orbit, next to whether it is in Gray's survey-astrometry list.""" f = fits["base"] kt = pd.read_csv(os.path.join(DATA, "geocentric_kappa.csv")).query("trailing_law == 'neomod2' and variant == 'base'").iloc[0] h = gh[GEO_CLASSES[0]]["neomod2"]["base"] out = hw[["JCAT", "Piece", "Type", "Name", "LDate", "Length", "Diameter", "D_eq_m", "Perigee", "Apogee", "Inc", "in_gray_list"]].copy() Hh = H_from_D(out.D_eq_m.values, 0.2) out["P_simulator_high_orbit_class"] = np.interp(Hh, H_GRID, pcurve(h, f, float(kt.kappa))) out["P_simulator_lowest_kappa"] = np.interp(Hh, H_GRID, pcurve(h, f, float(kt.kappa_lo))) out.to_csv(os.path.join(DATA, "heo_hardware_check.csv"), index=False) return out, dict(stats, mean_P_simulator=float(out.P_simulator_high_orbit_class.mean()), mean_P_simulator_lowest_kappa=float(out.P_simulator_lowest_kappa.mean())) def main(): rng = np.random.default_rng(SEED + 5) fits = json.load(open(os.path.join(DATA, "survey_fit.json"))) if not S.load_elong_table(): raise SystemExit("data/elong_weight_measured.csv missing: run s02 first") pd.set_option("display.width", 250) pd.set_option("display.max_colwidth", 60) # ---------------- A. heliocentric parked classes ---------------- hh, track_stats = helio_class_hists(rng) hc, hthr, sens = helio_curves(hh, fits) sens.to_csv(os.path.join(DATA, "parked_elong_sensitivity.csv"), index=False) print("[parked] co-orbital track statistics:", json.dumps(track_stats)) cols = ["location", "variant", "P_max", "P_max_surveys_only", "H_50", "D_m_90_p0.2", "D_m_50_p0.2", "D_m_50_wlo_p0.2", "D_m_50_whi_p0.2", "D_m_50_p0.05", "D_m_50_p0.01"] print(hthr[cols].round(2).to_string(index=False)) print("[parked] sensitivity of the co-orbital classes to the elongation weight (base variant; ramp60 = first version):") print(sens[sens.variant == "base"].drop(columns=["variant"]).round(3).to_string(index=False)) l4 = sens[sens.location == HELIO_CLASSES[0]] print("[parked] Sun-Earth L4/L5 parked: P(10 km) / P(1 km) / P(300 m), min-max over the five depth variants:") print(l4.groupby("weight_model")[["P_max_surveys_only", "P_10km", "P_1km", "P_300m", "D50_m_p0p2"]].agg(["min", "max"]).round(3).to_string()) # ---------------- D. geocentric classes with the minimoon correction ---------------- gh, rates = geo_class_hists(rng) print("[parked] sky rates of the geocentric classes (deg/day):", json.dumps(rates)) hw_heo, hw_stats = load_heo_hardware() gc, gthr, ktab, mm, kap = geo_curves(gh, fits, hw_heo, hw_stats) mm.to_csv(os.path.join(DATA, "minimoon_test.csv"), index=False) ktab.to_csv(os.path.join(DATA, "geocentric_kappa.csv"), index=False) print("[parked] minimoon correction factor kappa on the nightly probability of geocentric classes (targets: 2 known; range " f"{kappa_targets()[1]:.2f}-{kappa_targets()[2]:.1f} expected, from Poisson and a factor-{MINIMOON_POP_FACTOR:.0f} population uncertainty; " f"floor: {hw_stats['n_in_gray_list']} of {hw_stats['n_objects']} catalogued human objects in high Earth orbit are in Gray's " f"survey-astrometry list, 95% lower bound {hw_stats['fraction_lower95']:.2f}):") print(ktab.round(4).to_string(index=False)) gcols = ["location", "variant", "kappa", "kappa_lo", "D_m_50_p0.2", "D_m_50_klo_p0.2", "D_m_50_khi_p0.2", "D_m_50_theory_p0.2", "D_m_50_theory_klo_p0.2", "D_m_50_uncalibrated_p0.2", "P_1m", "P_1m_klo", "P_1m_khi", "P_1m_theory", "P_3m_klo", "P_10m_klo"] print(gthr[gcols].round(2).to_string(index=False)) curves = pd.concat([hc, gc], ignore_index=True) curves.to_csv(os.path.join(DATA, "parked_pdet_curves.csv"), index=False) thr = pd.concat([hthr, gthr], ignore_index=True) thr.to_csv(os.path.join(DATA, "parked_thresholds.csv"), index=False) # ---------------- LSST projection ---------------- central = {c: hh[c][("base", "c")] for c in HELIO_CLASSES} central.update({c: gh[c]["neomod2"]["base"] for c in GEO_CLASSES}) ls = lsst_projection(central, fits, kap[("neomod2", "base")][0]) ls.to_csv(os.path.join(DATA, "parked_lsst_projection.csv"), index=False) # ---------------- B. real objects ---------------- summ, o2, vm, n_precov = real_object_test({v: fits[v] for v in fits}) summ.to_csv(os.path.join(DATA, "real_coorbital_test.csv"), index=False) ec = pd.read_csv(os.path.join(DATA, "elongation_check.csv")) print(ec[ec.variant == "base"].round(1).to_string(index=False)) o2.to_csv(os.path.join(DATA, "real_coorbital_objects.csv"), index=False) print(summ.round(3).to_string(index=False)) print("[parked] discovery-magnitude test (real first-observation magnitudes of Earth-like / co-orbital asteroids versus each variant):") print(vm.round(3).to_string(index=False)) # ---------------- E. catalogued human objects in high Earth orbit ---------------- heo, hsum = heo_hardware_check(gh, fits, hw_heo, hw_stats) print("[parked] catalogued human objects in orbit with apogee > 60,000 km and D >= 1 m versus Gray's survey-astrometry list:") print(heo[["Piece", "Name", "D_eq_m", "Perigee", "Apogee", "in_gray_list", "P_simulator_high_orbit_class", "P_simulator_lowest_kappa"]].round(2).to_string(index=False)) print("[parked]", json.dumps(hsum)) # ---------------- C. human hardware in heliocentric orbit ---------------- hw, hs = human_hardware(fits, rng) hw.to_csv(os.path.join(DATA, "human_hardware_prediction.csv"), index=False) # hardware in heliocentric orbit that surveys later picked up as 'asteroids': J002E3 (Apollo 12 S-IVB), 2020 SO # (Surveyor 2 Centaur), 2024 BY15 (possibly the Helios B Centaur), 2018 CN41 (Falcon Heavy upper stage / Tesla) named = hw[hw.Name.isin(["Saturn S-IVB-507", "Centaur AC-7", "Centaur TC-5", "Falcon Heavy-001 Stage 2"])] hs["model_probabilities_of_recovered_objects"] = {str(r.Name): float(r.P_detect_base) for _, r in named.iterrows()} top = hw.sort_values("P_detect_base", ascending=False).head(5) hs["top5_contributors_base"] = {str(r.Name): float(r.P_detect_base) for _, r in top.iterrows()} json.dump(dict(hardware=hs, tracks=track_stats, heo_check=hsum, n_earthlike_first_obs_over_1yr_before_designation=n_precov), open(os.path.join(DATA, "human_hardware_summary.json"), "w"), indent=1) print("[parked] human hardware:", json.dumps(hs, indent=1)) show = ["Name", "LDate", "D_eq_m", "Perigee", "Apogee", "H_p0.2", "P_detect_base", "P_detect_shallow", "P_detect_deep"] print(hw.sort_values("P_detect_base", ascending=False).head(15)[show].round(3).to_string(index=False)) print("[parked] hardware later recovered by surveys (model probabilities):") print(named[show].round(3).to_string(index=False)) if __name__ == "__main__": main()