"""Step 4: validation of the covariance forecast, the direct search (measured A1), the recast fit-quality criterion, upper limits, vetting tables and figures. Inputs : data/analysis_plan.csv, data/fisher_results.csv (forecast + own-residual measurement), data/neocc_measurement.csv and data/neocc_nongrav.csv (s06_neocc), JPL SBDB snapshots, data/validation_*.csv written by s03 Outputs: data/limits.csv (measured search and fit-quality criterion; raw and weighted N), data/measured_flagged.csv, data/fit_quality_flagged.csv (every object that needed vetting), data/measured_5sigma_values.csv, data/followup_measured_candidates.csv, data/control_jpl_A1.csv, data/validation_A2_from_residuals.csv, data/validation_sigma.csv, data/completeness_by_bin.csv and completeness_grid.csv (forecast, not limits), data/robustness_limits.csv, data/poor_fits_catalogue.csv, data/sensitivity_per_object.csv, data/sigma_scaling_fit.csv, data/followup_list.csv, data/dark_comet_radial_bounds.csv, data/special_objects.csv, data/census_nongrav.csv, data/census_summary.csv, data/key_numbers.json, figures/fig1..fig5 """ import json import os import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import pandas as pd from common import (DATA, FIGS, A1_PER_AMR, AUD2_TO_MS2, MU_SUN, P_SRP, PV_NOM, PV_LO, PV_HI, H_of_diam, diam_m, amr_sphere, a2_yark_max, amr_of_A1) import s02_census as cen from orbfisher import upper_bound, LIN_MARGIN, RESOLVE_Z THRESH = [0.01, 0.1, 1.0] # m2/kg BINS = ["1-10 m", "10-100 m", ">100 m"] # palette (validated categorical order of the dataviz reference palette) C_BLUE, C_ORANGE, C_AQUA, C_YELLOW, C_MAGENTA, C_GREEN, C_VIOLET, C_RED = ( "#2a78d6", "#eb6834", "#1baf7a", "#eda100", "#e87ba4", "#008300", "#6250d6", "#e34948") INK, INK2, GRID, SURF = "#0b0b0b", "#52514e", "#e4e3df", "#fcfcfb" plt.rcParams.update({ "figure.facecolor": SURF, "axes.facecolor": SURF, "savefig.facecolor": SURF, "axes.edgecolor": "#b9b8b2", "axes.labelcolor": INK, "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", "legend.frameon": False, "figure.dpi": 110, "savefig.dpi": 160, }) def ul95(k, n): """Clopper-Pearson one-sided 95% upper limit on a binomial fraction.""" from scipy.stats import beta if n <= 0: return np.nan return float(beta.ppf(0.95, k + 1, n - k)) if k < n else 1.0 def size_bins(H, pv=PV_NOM): h10, h100 = H_of_diam(10.0, pv), H_of_diam(100.0, pv) return np.where(H > h10, "1-10 m", np.where(H > h100, "10-100 m", ">100 m")) # --------------------------------------------------------------- validation def validation(df): v1 = df[df.A1_sigma.notna() & df.sig_A1_jplset.notna()].copy() v1["ratio"] = v1.sig_A1_jplset / v1.A1_sigma v1["par"] = "A1"; v1["jpl"] = v1.A1_sigma; v1["fisher"] = v1.sig_A1_jplset v2 = df[df.A2_sigma.notna() & df.sig_A2_jplset.notna() & np.isfinite(df.sig_A2_jplset)].copy() v2["ratio"] = v2.sig_A2_jplset / v2.A2_sigma v2["par"] = "A2"; v2["jpl"] = v2.A2_sigma; v2["fisher"] = v2.sig_A2_jplset v = pd.concat([v1, v2])[["pdes", "par", "model", "sbin", "H", "arc_obs", "n_obs", "n_obs_used", "n_radar", "jpl", "fisher", "ratio", "res_med", "post_freeze"]] v.to_csv(os.path.join(DATA, "validation_sigma.csv"), index=False) out = {} for (par, model), g in v.groupby(["par", "model"]): lr = np.log10(g.ratio) out["%s_%s" % (par, model)] = dict(n=int(len(g)), median_ratio=float(10 ** lr.median()), p16=float(10 ** lr.quantile(0.16)), p84=float(10 ** lr.quantile(0.84)), frac_within_2x=float((lr.abs() < np.log10(2)).mean())) return v, out def kepler_check(): path = os.path.join(DATA, "validation_kepler_vs_numeric.csv") if not os.path.exists(path): return None, {} k = pd.read_csv(path, dtype={"pdes": str}) if "sig_A1" not in k.columns: return None, {} k = k[(k.model == "numeric") & np.isfinite(k.sig_A1) & np.isfinite(k.sig_A1_kepler)].copy() k["ratio"] = k.sig_A1_kepler / k.sig_A1 lr = np.log10(k.ratio) if len(k) < 5: return None, {} return k, dict(n=int(len(k)), median_ratio=float(10 ** lr.median()), p16=float(10 ** lr.quantile(0.16)), p84=float(10 ** lr.quantile(0.84)), max_ratio=float(k.ratio.max()), min_ratio=float(k.ratio.min())) # --------------------------------------------------------------- sensitivity R_LIST = [0.7, 1.0] # fixed rms ceilings of the fit-quality criterion (a sail would show rms >= 2R) SNR_T = 6.0 # measured sample: sigma(A/m) <= X / SNR_T K = A1_PER_AMR def add_sensitivity(df, cal_numeric=1.0, cal_kepler=1.0, scale=1.0, snr_t=SNR_T, clipped=True, own_multi=False, joint_multi=True): """Per-object quantities used by the three criteria. forecast amr_3sig : smallest A/m that would be a 3-sigma detection if A1 were fitted measured amr_T : smallest threshold X for which the object belongs to the tested sample T_X (formal sigma(A/m) of the optical data <= X / snr_t; this does not depend on the residuals) amr_hat, amr_msig : measured A/m and its uncertainty (scaled by the post-fit rms when that exceeds 1) amr_ul3 : 3-sigma upper bound of the primary fit (after outlier rejection); all bounds are max(estimate, 0) / 0.75 + 3 sigma: sunward estimates count as zero and 0.75 allows for the linearisation bias of a large signal found in the injection test amr_ul3_all : the same for the (linear) fit to all observations; infinite if the primary estimate is significant (>= 4 sigma), because the linear step is then biased and must not be used to clear the object amr_lono_ul : smallest 3-sigma upper bound when one station-night is left out; infinite if the estimate without that night is still >= 4 sigma amr_ul_best : the smallest of the three (used only to resolve flagged objects) Arcs <= 180 d: residuals about JPL's trajectory in this study's numerical model, 6 state corrections + A1. Arcs > 180 d: ESA NEOCC post-fit residuals (s06_neocc); with numerical partials 6 state corrections + A2 + A1 (joint_multi=True: the transverse Yarkovsky term is solved as a nuisance parameter), with two-body partials 6 + A1; unless own_multi=True. fit quality amr_fq_R : smallest X for which the forecast significance of X satisfies s >= 2 R sqrt(2 n); the sample is defined WITHOUT the observed rms fq_fail_rms_R, fq_fail_rej : members that fail (observed rms >= 2R; JPL used under 80% of the MPC astrometry) and therefore need vetting """ df = df.copy() cal = np.where(df.model == "kepler", cal_kepler, cal_numeric) * scale df["sigA1_cal"] = df.sig_A1 * cal bad = ~df.model.isin(["numeric", "kepler"]) | ~np.isfinite(df.sigA1_cal) df.loc[bad, "sigA1_cal"] = np.inf has = (df.A1.notna() & df.A1_sigma.notna()).values df["has_A1"] = has # ---- forecast df["amr_3sig"] = 3.0 * df.sigA1_cal / K df.loc[has, "amr_3sig"] = 3 * df.A1_sigma[has] / K # ---- measurement single = (df.arc_d <= 180).values sfx = "" if clipped else "_all" def pick(prefix): A, S, R = df[prefix + "A1" + sfx], df[prefix + "sig" + sfx], df[prefix + "rms7" + sfx] Aa, Sa, Ra = df[prefix + "A1_all"], df[prefix + "sig_all"], df[prefix + "rms7_all"] ok = np.isfinite(A) & np.isfinite(S) & (S > 0) hat = A / K msig = S * np.maximum(1.0, R.fillna(1.0)) * scale / K # upper bounds: max(estimate, 0) / LIN_MARGIN + 3 sigma (see orbfisher.upper_bound) ul = upper_bound(hat, msig) ul_all = upper_bound(Aa / K, Sa * np.maximum(1.0, Ra.fillna(1.0)) * scale / K) lono = upper_bound(df[prefix + "lono_A1"] / K, df[prefix + "lono_sig"] * scale / K) lono_z = df[prefix + "lono_A1"] / (df[prefix + "lono_sig"] * scale) return ok, hat, msig, ul, ul_all, lono, lono_z def usable(hat, msig, ul_all, lono, lono_z): """Bounds that may resolve a flagged object: the linear all-data fit only if the primary estimate is not significant, the leave-one-night-out fit only if what is left is not.""" z = hat / msig return (pd.Series(np.where(z < RESOLVE_Z, ul_all, np.inf), index=df.index), pd.Series(np.where(lono_z < RESOLVE_Z, lono, np.inf), index=df.index)) ok_o, hat_o, sg_o, ul_o, ula_o, lo_o, lz_o = pick("m_") ok_o &= (df.model == "numeric") if clipped: # significant values were re-estimated with a full non-linear fit in s03 (m_nl_*) nl = (df.m_nl_ok == True).values # noqa: E712 hat_o = pd.Series(np.where(nl, df.m_nl_A1 / K, hat_o), index=df.index) sg_o = pd.Series(np.where(nl, df.m_nl_sig * np.maximum(1.0, df.m_nl_rms7.fillna(1.0)) * scale / K, sg_o), index=df.index) ul_o = upper_bound(hat_o, sg_o) ula_o, lo_o = usable(hat_o, sg_o, ula_o, lo_o, lz_o) T_o = snr_t * df.sig_A1_opt * scale / K # optical-only forecast (the measurement uses no radar) # ESA residuals. Numerical partials (objects < 100 m): A2 is solved together with A1. Two-body # partials (objects > 100 m): A1 alone, because the two-body A2 column is only a secular # approximation (s05_injection.main_twobody: a joint fit is unreliable there) and the leak of # the largest plausible Yarkovsky drift of a 100 m body into A1 is 4e-5 m2/kg in the median and # 7e-3 m2/kg in the worst of 40 test objects (smaller in proportion to 1/D for larger bodies). jn = joint_multi & (df.model == "numeric").values p1, p2 = pick("nc_"), pick("nc2_") ok_n, hat_n, sg_n, ul_n, ula_n, lo_n, lz_n = [pd.Series(np.where(jn, b_, a_), index=df.index) for a_, b_ in zip(p1, p2)] ok_n = ok_n.astype(bool) ula_n, lo_n = usable(hat_n, sg_n, ula_n, lo_n, lz_n) okn_any = ok_n & (df.nc_status == "ok") ok_n = okn_any & df.nc_random T_n = snr_t * np.where(jn, df.nc2_sig_all, df.nc_sig_all) * scale / K # formal sigma, ESA's weights, all optical data use_own = single | own_multi ok = np.where(use_own, ok_o, ok_n) df["meas_source"] = np.where(ok, np.where(use_own, "own residuals (JPL trajectory)", np.where(jn, "ESA NEOCC residuals, A2 solved", "ESA NEOCC residuals")), "") df["amr_T"] = np.where(ok, np.where(use_own, T_o, T_n), np.inf) df["amr_hat"] = np.where(ok, np.where(use_own, hat_o, hat_n), np.nan) df["amr_msig"] = np.where(ok, np.where(use_own, sg_o, sg_n), np.nan) df["amr_ul3"] = np.where(ok, np.where(use_own, ul_o, ul_n), np.nan) df["amr_ul3_all"] = np.where(ok, np.where(use_own, ula_o, ula_n), np.nan) df["amr_lono_ul"] = np.where(ok, np.where(use_own, lo_o, lo_n), np.nan) df["m_rms6_used"] = np.where(use_own, df.m_rms6, np.where(jn, df.nc2_rms6, df.nc_rms6)) df["m_rms7_used"] = np.where(use_own, df.m_rms7, np.where(jn, df.nc2_rms7, df.nc_rms7)) df["lono_zmin"] = np.where(use_own, df.m_lono_zmin, np.where(jn, df.nc2_lono_zmin, df.nc_lono_zmin)) # objects with a catalogued A1 (JPL) or a fitted area-to-mass ratio in ESA's solution: # the catalogue estimate is the measurement esa = (np.isfinite(df.esa_amr) & np.isfinite(df.esa_amr_sigma) & ~has).values for m_, a_, s_, src in ((has, df.A1 / K, df.A1_sigma / K, "JPL catalogue A1"), (esa, df.esa_amr, df.esa_amr_sigma, "ESA NEOCC catalogue A/m")): df.loc[m_, "amr_T"] = snr_t * s_[m_] * scale df.loc[m_, "amr_hat"] = a_[m_] df.loc[m_, "amr_msig"] = s_[m_] * scale df.loc[m_, "amr_ul3"] = a_[m_] + 3 * scale * s_[m_] df.loc[m_, "amr_ul3_all"] = a_[m_] + 3 * scale * s_[m_] df.loc[m_, "amr_lono_ul"] = np.nan df.loc[m_, "meas_source"] = src df["amr_ul_best"] = df[["amr_ul3", "amr_ul3_all", "amr_lono_ul"]].min(axis=1) # ---- fit-quality criterion, recast so that membership does not depend on the observed rms n = df.n_obs_used df["fq_fail_rej"] = ~(df.n_obs_used >= 0.8 * df.n_obs).fillna(False) & ~has for R in R_LIST: col = "amr_fq_%g" % R df[col] = 2.0 * R * np.sqrt(2.0 * n) * df.sigA1_cal / K df.loc[has, col] = snr_t * df.A1_sigma[has] * scale / K df["fq_fail_rms_%g" % R] = (df.rms >= 2.0 * R) & ~has # measurement used to vet fit-quality failures (any object with a valid measurement, also the # flagged objects that are not part of the random measured sample) v_o = np.fmin(np.fmin(ul_o, ula_o), lo_o).where(ok_o) v_n = np.fmin(np.fmin(ul_n, ula_n), lo_n).where(okn_any) df["vet_ul3"] = np.where(single, v_o, v_n) df.loc[has | esa, "vet_ul3"] = df.amr_ul3[has | esa] return df def arc_mask(s, arcs): if arcs == "single-apparition (<=180 d)": return s.arc_d <= 180 if arcs == "multi-apparition (>180 d)": return s.arc_d > 180 return pd.Series(True, index=s.index) ARCS = ["single-apparition (<=180 d)", "multi-apparition (>180 d)", "all arcs"] def completeness(df, bins_col="sbin"): """Forecast completeness: weighted number of catalogued objects for which threshold X would be a 3-sigma detection if A1 were fitted. This is NOT a limit.""" rows = [] for b in BINS: g = df[df[bins_col] == b] n_cat = len(g) s = g[g.selected] for x in THRESH: hit = s.amr_3sig <= x n_eff = float((s.w * hit).sum()) var = float((s.w * (s.w - 1) * hit).sum()) # sampling variance of the weighted count rows.append(dict(size_bin=b, amr_threshold=x, n_catalogued=n_cat, n_analysed=int(len(s)), n_analysed_sensitive=int(hit.sum()), n_sensitive=n_eff, n_sensitive_err=np.sqrt(var), frac_of_catalogue=n_eff / n_cat if n_cat else np.nan)) return pd.DataFrame(rows) def _row(method, arcs, b, x, n_cat, T, flagged, unres, extra=None): nT = int(len(T)); nw = float(T.w.sum()) r = dict(method=method, arcs=arcs, size_bin=b, amr_threshold=x, n_catalogued=n_cat, n_tested=nT, n_tested_weighted=nw, frac_of_bin=nw / n_cat if n_cat else np.nan, k_flagged=int(len(flagged)), k_unresolved=int(len(unres)), ul95=ul95(len(unres), nT), ul95_if_flagged_count=ul95(len(flagged), nT), ul95_weightedN=ul95(len(unres), int(np.floor(nw))) if nw >= 1 else np.nan, flagged=";".join(flagged.pdes.astype(str)), unresolved=";".join(unres.pdes.astype(str))) if extra: r.update(extra) return r def measured_limits(d, bins_col="sbin", arcs_list=ARCS): """Direct search. T_X = measured objects whose formal sigma(A/m) <= X/6 (defined without looking at the residuals). A member is cleared when the 3-sigma upper bound of the fit after outlier rejection is below X; otherwise it is flagged. A flagged object is resolved only if the fit to all observations, or the data without one single station-night, exclude X at 3 sigma. Unresolved objects are counted as if they were detections. Limits use the raw number of tested objects (not the catalogue-weighted number).""" rows = [] for b in BINS: g = d[d[bins_col] == b] s0 = g[g.selected] for arcs in arcs_list: s = s0[arc_mask(s0, arcs)] for x in THRESH: T = s[s.amr_T <= x] flagged = T[~(T.amr_ul3 < x)] unres = flagged[~(flagged.amr_ul_best < x)] rows.append(_row("measured", arcs, b, x, len(g), T, flagged, unres, dict(n_tested_weighted=float(T.w_meas.sum()), frac_of_bin=float(T.w_meas.sum()) / len(g) if len(g) else np.nan, ul95_weightedN=ul95(len(unres), int(np.floor(T.w_meas.sum()))) if T.w_meas.sum() >= 1 else np.nan))) return pd.DataFrame(rows) def fq_limits(d, bins_col="sbin", arcs_list=ARCS, strict=False): """Fit-quality criterion (recast). Sample: forecast significance of X >= 2 R sqrt(2 n), R fixed. A member with A/m >= X would have a normalised rms >= 2R in JPL's gravity-only fit. Members that fail (rms >= 2R, or under 80% of the astrometry used) are flagged; a flagged object is resolved if a direct measurement of A1 excludes X at 3 sigma.""" rows = [] for b in BINS: g = d[d[bins_col] == b] s0 = g[g.selected] for arcs in arcs_list: s = s0[arc_mask(s0, arcs)] for R in R_LIST: for x in THRESH: T = s[s["amr_fq_%g" % R] <= x] if strict: # require that JPL used at least as many observations as the MPC lists T = T[(T.n_obs_used >= T.n_obs) | T.has_A1] fail = T[T["fq_fail_rms_%g" % R] | T.fq_fail_rej | (T.has_A1 & ~(T.amr_ul3 < x))] unres = fail[~(fail.vet_ul3 < x)] rows.append(_row("fit quality R=%g" % R, arcs, b, x, len(g), T, fail, unres, dict(k_fail_rms=int(T["fq_fail_rms_%g" % R].sum()), k_fail_rejections=int(T.fq_fail_rej.sum())))) return pd.DataFrame(rows) def cen_size_bin(H): return size_bins(np.asarray(H, float)) def flagged_tables(d): """Lists of every object that needed vetting (written to data/).""" s = d[d.selected] cols = ["pdes", "full_name", "sbin", "H", "arc_d", "n_obs", "n_obs_used", "n_space", "rms", "model", "meas_source", "amr_T", "amr_hat", "amr_msig", "amr_ul3", "amr_ul3_all", "amr_lono_ul", "amr_ul_best", "m_rms6_used", "m_rms7_used", "dist_min"] rows = [] for x in THRESH: T = s[s.amr_T <= x] f = T[~(T.amr_ul3 < x)].copy() f["amr_threshold"] = x f["resolved_by_all_data_fit"] = f.amr_ul3_all < x f["resolved_by_single_night_removal"] = f.amr_lono_ul < x f["unresolved"] = ~(f.amr_ul_best < x) rows.append(f[["amr_threshold"] + cols + ["resolved_by_all_data_fit", "resolved_by_single_night_removal", "unresolved"]]) meas = pd.concat(rows) rows = [] for R in R_LIST: for x in THRESH: T = s[s["amr_fq_%g" % R] <= x] f = T[T["fq_fail_rms_%g" % R] | T.fq_fail_rej | (T.has_A1 & ~(T.amr_ul3 < x))].copy() f["R"] = R; f["amr_threshold"] = x f["fail_rms"] = f["fq_fail_rms_%g" % R]; f["fail_rejections"] = f.fq_fail_rej f["resolved_by_measurement"] = f.vet_ul3 < x rows.append(f[["R", "amr_threshold"] + cols + ["nc_rms_header", "nc_n_optical", "nc_n_rejected_by_neocc", "fail_rms", "fail_rejections", "vet_ul3", "resolved_by_measurement"]]) fq = pd.concat(rows) return meas, fq # ------------------------------------------------------------------ figures def fig_census(census): fig, ax = plt.subplots(figsize=(9.2, 7.3)) D = np.logspace(-0.2, 4.2, 100) ax.plot(D, amr_sphere(D, 2000.0) * A1_PER_AMR, color=INK2, lw=1.4, ls="-", label="radiation pressure on a 2000 kg/m$^3$ sphere (A1)") ax.plot(D, a2_yark_max(D, 2000.0), color=INK2, lw=1.4, ls="--", label="maximum Yarkovsky for 2000 kg/m$^3$ (A2)") for amr, lab in ((1.0, "solar sail, 1 m$^2$/kg"), (1e-2, "spent rocket stage, 10$^{-2}$ m$^2$/kg")): ax.axhline(amr * A1_PER_AMR, color=C_RED, lw=1.0, ls=":") ax.text(1.2e4, amr * A1_PER_AMR * 1.25, lab, color=INK2, fontsize=8, ha="right", va="bottom") c = census[census.significant].copy() # one point per object: its largest significant component comp = [] for r in c.itertuples(): best = max(((abs(getattr(r, k)), k) for k in ("A1", "A2", "A3") if np.isfinite(getattr(r, k)) and getattr(r, "snr_" + k) >= 3), default=(np.nan, "")) comp.append(best) c["a_plot"] = [x[0] for x in comp]; c["comp"] = [x[1] for x in comp] groups = [ ("Yarkovsky-consistent", C_BLUE, "o", 10, 0.45, "A2 not excluded as Yarkovsky (%d)"), ("radiation pressure", C_AQUA, "s", 46, 1.0, "radiation pressure on natural body (%d)"), ("published dark comet", C_ORANGE, "D", 46, 1.0, "published dark comet / outgassing (%d)"), ("known active", C_VIOLET, "^", 60, 1.0, "known active or interstellar (%d)"), ("unpublished outgassing", C_RED, "X", 90, 1.0, "unpublished outgassing candidate, JPL-fitted A2 (%d)"), ("marginal", C_YELLOW, "P", 70, 1.0, "sunward A1, likely systematic (%d)"), ] for key, col, mk, sz, al, lab in groups: g = c[c.class_ng.str.startswith(key)] ax.scatter(g.D_nom, g.a_plot, s=sz, marker=mk, color=col, alpha=al, edgecolor=SURF if sz > 20 else "none", linewidth=0.8, label=lab % len(g), zorder=3 if sz > 20 else 2) art = cen.ARTIFICIAL ax.scatter([diam_m(a["H"]) for a in art], [a["A1"] for a in art], s=110, marker="*", color=INK, zorder=4, label="artificial (removed from catalogue) (%d)" % len(art)) lab = {"2017 U1": ("1I/'Oumuamua", 7, 4), "457175": ("362P (2008 GO98)", -8, 8), "883607": ("2016 TA56", 7, 4), "875163": ("P/1998 SH2", 7, 4), "523599": ("2003 RM", -8, -12), "2012 UR158": ("2012 UR158", 7, 4), "2016 NJ33": ("2016 NJ33", 8, -3), "2006 RH120": ("2006 RH120", -7, 8), "2020 CD3": ("2020 CD3", -2, -15), "4179": ("Toutatis", 7, 4), "2013 XY20": ("2013 XY20", 7, 4), "3200": ("Phaethon", -8, -13), "428209": ("2006 VC", -8, -11), "1998 KY26": ("1998 KY26", -8, -12), "2005 VL1": ("2005 VL1", -8, -2)} for r in c.itertuples(): if r.pdes in lab: t, dx, dy = lab[r.pdes] ax.annotate(t, (r.D_nom, r.a_plot), xytext=(dx, dy), textcoords="offset points", fontsize=8, color=INK, ha="left" if dx > 0 else "right") ax.annotate("2020 SO", (diam_m(art[0]["H"]), art[0]["A1"]), xytext=(8, 3), textcoords="offset points", fontsize=8) ax.annotate("J002E3", (diam_m(art[1]["H"]), art[1]["A1"]), xytext=(-8, 6), textcoords="offset points", fontsize=8, ha="right") ax.set_xscale("log"); ax.set_yscale("log") ax.set_xlim(0.6, 1.5e4); ax.set_ylim(2e-16, 1e-6) ax.set_xlabel("diameter from H (albedo 0.14 unless measured) [m]") ax.set_ylabel("largest significant non-gravitational parameter |A| at 1 au [au/d$^2$]") sec = ax.secondary_yaxis("right", functions=(lambda a: a / A1_PER_AMR, lambda m: m * A1_PER_AMR)) sec.set_ylabel("equivalent area-to-mass ratio [m$^2$/kg] (C$_R$ = 1)") ax.set_title("Census of non-gravitational accelerations on asteroid-designated objects (JPL SBDB, 2026-10-12)") fig.legend(*ax.get_legend_handles_labels(), loc="lower center", fontsize=8, ncol=3, bbox_to_anchor=(0.5, 0.005)) fig.tight_layout(rect=(0, 0.105, 1, 1)) fig.savefig(os.path.join(FIGS, "fig1_census.png")); plt.close(fig) def fig_validation(v, kchk, selt=None, pz=None, pnc=None): fig, axs = plt.subplots(2, 2, figsize=(11.6, 9.4)) ax = axs[0, 0] g = v[v.par == "A1"] ax.scatter(g.jpl, g.fisher, s=40, color=C_BLUE, edgecolor=SURF, linewidth=0.8, zorder=3) lim = [g.jpl.min() / 3, g.jpl.max() * 3] ax.plot(lim, lim, color=INK2, lw=1); ax.fill_between(lim, np.array(lim) / 2, np.array(lim) * 2, color=GRID, alpha=0.7, lw=0) ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlim(lim); ax.set_ylim(lim) ax.set_xlabel("JPL formal $\\sigma$(A1) [au/d$^2$]"); ax.set_ylabel("this work, forecast $\\sigma$(A1) [au/d$^2$]") ax.set_title("(a) Radial term: %d NEOs with a JPL-fitted A1" % len(g)) ax = axs[0, 1] g = v[v.par == "A2"] for model, col, mk in (("kepler", C_BLUE, "o"), ("numeric", C_ORANGE, "s")): h = g[g.model == model] ax.scatter(h.jpl, h.fisher, s=12 if model == "kepler" else 30, color=col, marker=mk, alpha=0.6 if model == "kepler" else 1, edgecolor="none" if model == "kepler" else SURF, linewidth=0.6, label="%s model (%d)" % ("two-body" if model == "kepler" else "numerical", len(h)), zorder=3) lim = [g.jpl.min() / 3, g.jpl.max() * 3] ax.plot(lim, lim, color=INK2, lw=1); ax.fill_between(lim, np.array(lim) / 2, np.array(lim) * 2, color=GRID, alpha=0.7, lw=0) ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlim(lim); ax.set_ylim(lim) ax.set_xlabel("JPL formal $\\sigma$(A2) [au/d$^2$]"); ax.set_ylabel("this work, forecast $\\sigma$(A2) [au/d$^2$]") ax.set_title("(b) Transverse term: %d objects with a JPL-fitted A2" % len(g)); ax.legend(fontsize=8, loc="upper left") ax = axs[1, 0] lo, hi = [], [] if kchk is not None and len(kchk): ax.scatter(kchk.sig_A1, kchk.sig_A1_kepler, s=30, color=C_BLUE, edgecolor=SURF, linewidth=0.6, zorder=3, label="two-body vs numerical forecast, same data (%d)" % len(kchk)) lo += [kchk.sig_A1.min(), kchk.sig_A1_kepler.min()]; hi += [kchk.sig_A1.max(), kchk.sig_A1_kepler.max()] if selt is not None and len(selt): ax.scatter(selt.sig_A1_published, selt.kepler_A1A2_to2024, s=70, marker="D", color=C_ORANGE, edgecolor=SURF, linewidth=0.8, zorder=4, label="two-body forecast vs published $\\sigma$(A1),\nSeligman et al. 2024 (%d dark comets)" % len(selt)) for r in selt.itertuples(): ax.annotate(r.object, (r.sig_A1_published, r.kepler_A1A2_to2024), xytext=(6, -3), textcoords="offset points", fontsize=7, color=INK2) lo += [selt.sig_A1_published.min(), selt.kepler_A1A2_to2024.min()]; hi += [selt.sig_A1_published.max(), selt.kepler_A1A2_to2024.max()] if lo: lim = [min(lo) / 2, max(hi) * 2] ax.plot(lim, lim, color=INK2, lw=1); ax.fill_between(lim, np.array(lim) / 2, np.array(lim) * 2, color=GRID, alpha=0.7, lw=0) ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlim(lim); ax.set_ylim(lim) ax.set_xlabel("reference $\\sigma$(A1) [au/d$^2$] (numerical model, or published)") ax.set_ylabel("two-body forecast $\\sigma$(A1) [au/d$^2$]") ax.set_title("(c) Multi-year arcs: two-body forecast"); ax.legend(fontsize=7.5, loc="upper left") ax = axs[1, 1] if pz is not None and len(pz): g = pz[pz.arc_d <= 180] edges = np.linspace(-5, 5, 51) ax.hist(np.clip(g.z.values, -4.99, 4.99), bins=edges, density=True, color=C_BLUE, alpha=0.85, edgecolor=SURF, linewidth=0.4, label="single apparition, residuals about JPL's\ntrajectory (%d NEOs)" % len(g)) if pnc is not None and len(pnc): ax.hist(np.clip(pnc.z.values, -4.99, 4.99), bins=edges, density=True, histtype="step", color=C_ORANGE, linewidth=2, label="multi-apparition, ESA NEOCC residuals\n(A2 solved below 100 m; %d NEOs)" % len(pnc)) xx = np.linspace(-5, 5, 400) ax.plot(xx, np.exp(-0.5 * xx ** 2) / np.sqrt(2 * np.pi), color=INK, lw=1.6, ls="--", label="unit Gaussian") ax.set_xlim(-5, 5) ax.legend(fontsize=7.5, loc="upper left") ax.set_xlabel("measured A1 / $\\sigma$(A1) (values beyond $\\pm$5 in the end bins)") ax.set_ylabel("probability density") ax.set_title("(d) Real residuals: pull of the measured A1") fig.suptitle("Validation of the forecast uncertainty (grey band: factor 2)", fontsize=11, fontweight="bold") fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig2_validation.png")); plt.close(fig) def fig_sensitivity(df): fig, ax = plt.subplots(figsize=(9.2, 5.8)) s = df[df.selected & np.isfinite(df.amr_3sig) & (df.arc_obs > 0)] for b, col, mk in zip(BINS[::-1], (C_AQUA, C_ORANGE, C_BLUE), ("^", "s", "o")): g = s[s.sbin == b] ax.scatter(g.arc_obs, g.amr_3sig, s=5 if b != "1-10 m" else 9, color=col, marker=mk, alpha=0.35 if b != "1-10 m" else 0.7, edgecolor="none", label="%s (%d analysed)" % (b, len(g)), rasterized=True) for x, lab in zip(THRESH, ("rocket stage 0.01", "0.1", "sail 1 m$^2$/kg")): ax.axhline(x, color=C_RED, lw=1.0, ls=":") ax.text(3.6e4, x * 1.25, lab, fontsize=8, color=INK, ha="right", bbox=dict(facecolor=SURF, edgecolor="none", pad=1.0, alpha=0.85)) ax.axhspan(1e-4, 1e-3, color=GRID, alpha=0.8, lw=0, zorder=0) ax.text(0.75, 3e-4, "metre-scale rocks", fontsize=8, color=INK2, va="center") ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlim(0.6, 4e4); ax.set_ylim(1e-7, 1e4) ax.set_xlabel("observed arc length [days]") ax.set_ylabel("forecast: smallest A/m detectable at 3$\\sigma$ if fitted [m$^2$/kg]") ax.set_title("Per-object sensitivity of existing astrometry to a radial 1/r$^2$ acceleration") leg = ax.legend(loc="upper right", fontsize=8, markerscale=3) for lh in leg.legend_handles: lh.set_alpha(1) fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig3_sensitivity_vs_arc.png")); plt.close(fig) def fig_completeness(df): fig, axs = plt.subplots(1, 2, figsize=(12.4, 4.8), sharey=False) edges = np.arange(16, 33.01, 1.0) mid = 0.5 * (edges[1:] + edges[:-1]) ax = axs[0] for x, col, mk in zip(THRESH, (C_BLUE, C_ORANGE, C_AQUA), ("o", "s", "^")): fr, er = [], [] for lo, hi in zip(edges[:-1], edges[1:]): g = df[(df.H > lo) & (df.H <= hi)] s = g[g.selected] n = len(g) hit = s.amr_3sig <= x f = (s.w * hit).sum() / n if n else np.nan fr.append(f); er.append(np.sqrt(max(f * (1 - f), 0) / max(len(s), 1)) if n else np.nan) ax.errorbar(mid, fr, yerr=er, color=col, marker=mk, ms=5, lw=2, capsize=0, label="%g m$^2$/kg" % x) ax.set_xlabel("absolute magnitude H"); ax.set_ylabel("fraction of catalogued NEOs") ax.set_title("Forecast: A/m detectable at 3$\\sigma$ if A1 were fitted, by size") ax.set_ylim(0, 1.02); ax.legend(title="area-to-mass threshold", fontsize=8, title_fontsize=8, loc="upper right") top = ax.secondary_xaxis("top", functions=(lambda h: h, lambda h: h)) hh = [17.75, 22.75, 27.75, 32.75] top.set_xticks(hh); top.set_xticklabels(["1 km", "100 m", "10 m", "1 m"]); top.set_xlabel("diameter (albedo 0.14)", fontsize=8) ax = axs[1] aedges = [1, 3, 10, 30, 100, 365, 3650, 40000] lab = ["1-3 d", "3-10 d", "10-30 d", "30-100 d", "100 d-1 yr", "1-10 yr", ">10 yr"] xx = np.arange(len(lab)) for k, (x, col, mk) in enumerate(zip(THRESH, (C_BLUE, C_ORANGE, C_AQUA), ("o", "s", "^"))): fr = [] for lo, hi in zip(aedges[:-1], aedges[1:]): s = df[df.selected & (df.arc_d >= lo) & (df.arc_d < hi)] fr.append((s.w * (s.amr_3sig <= x)).sum() / s.w.sum() if len(s) else np.nan) ax.plot(xx, fr, color=col, marker=mk, ms=6, lw=2, label="%g m$^2$/kg" % x) ax.set_xticks(xx); ax.set_xticklabels(lab, fontsize=8) ax.set_xlabel("data-arc length"); ax.set_ylabel("fraction of NEOs with that arc") ax.set_title("Forecast, by arc length"); ax.set_ylim(0, 1.02) fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig4_completeness.png")); plt.close(fig) def fig_limits(d, lim): from matplotlib.lines import Line2D from matplotlib.patches import Patch fig, axs = plt.subplots(1, 2, figsize=(13.4, 5.8), gridspec_kw=dict(width_ratios=[1.15, 1.0])) # ---- (a) the direct search ax = axs[0] s = d[d.selected & np.isfinite(d.amr_T) & np.isfinite(d.amr_ul3) & (d.amr_ul3 > 0)] cols = dict(zip(BINS, (C_BLUE, C_ORANGE, C_AQUA))) for b, mk in zip(BINS[::-1], ("^", "s", "o")): g = s[s.sbin == b] ax.scatter(g.amr_T / SNR_T, g.amr_ul3, s=5 if b != "1-10 m" else 9, color=cols[b], marker=mk, alpha=0.35 if b != "1-10 m" else 0.7, edgecolor="none", label="%s (%d measured)" % (b, len(g)), rasterized=True) for x, lab in zip(THRESH, ("rocket stage, 0.01", "0.1", "sail, 1 m$^2$/kg")): ax.plot([1e-7, x / SNR_T, x / SNR_T], [x, x, 1e-7], color=C_RED, lw=1.2, ls=":") ax.text(1.3e-6, x * 1.25, lab, fontsize=8, color=INK, bbox=dict(facecolor=SURF, edgecolor="none", pad=1.0, alpha=0.85)) fl = pd.concat([s[(s.amr_T <= x) & ~(s.amr_ul3 < x)] for x in THRESH]) if len(s) else s if len(fl): ax.scatter(fl.amr_T / SNR_T, fl.amr_ul3, s=70, facecolor="none", edgecolor=INK, linewidth=1.3, zorder=5, label="flagged: tested, not cleared by the primary fit (%d)" % fl.pdes.nunique()) opath = os.path.join(DATA, "validation_oumuamua.csv") if os.path.exists(opath): o = pd.read_csv(opath).iloc[0] if "A1_measured" in o: ax.scatter([o.A1_measured_sigma / K], [(o.A1_measured + 3 * o.A1_measured_sigma) / K], s=160, marker="*", color=INK, zorder=6, label="1I/'Oumuamua: same fit on its real astrometry") ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlim(1e-6, 1e3); ax.set_ylim(1e-6, 3e3) ax.set_xlabel("formal $\\sigma$(A/m) [m$^2$/kg] (an object is tested at a threshold if it lies left of that dotted corner)") ax.set_ylabel("measured 3$\\sigma$ upper bound on A/m [m$^2$/kg]") ax.set_title("(a) Direct search: A1 fitted to the real astrometric residuals") leg = ax.legend(loc="upper left", fontsize=7.5) for lh in leg.legend_handles[:3]: lh.set_alpha(1); lh.set_sizes([30]) # ---- (b) limits (upper limits: marker with a downward tick) ax = axs[1] series = [("measured", ARCS[0], "o", True, "measured, single apparition"), ("measured", ARCS[1], "s", True, "measured, multi-apparition (ESA residuals)"), ("fit quality R=1", "all arcs", "D", False, "fit quality (R = 1), all arcs")] tcol = dict(zip(THRESH, (C_BLUE, C_ORANGE, C_AQUA))) for j, x in enumerate(THRESH): for i, (meth, arcs, mk, filled, lab) in enumerate(series): g = lim[(lim.method == meth) & (lim.arcs == arcs) & (lim.amr_threshold == x)].set_index("size_bin").reindex(BINS) xs = np.arange(len(BINS)) + (j - 1) * 0.27 + (i - 1) * 0.075 ok = np.isfinite(g.ul95.values) & (g.n_tested.values >= 5) ax.scatter(xs[ok], g.ul95.values[ok], s=62, marker=mk, facecolor=tcol[x] if filled else SURF, edgecolor=tcol[x] if not filled else SURF, linewidth=1.6 if not filled else 0.8, zorder=3) for xi, yi in zip(xs[ok], g.ul95.values[ok]): ax.plot([xi, xi], [yi, yi / 1.9], color=tcol[x], lw=1.2, zorder=2) ax.set_yscale("log"); ax.set_xticks(np.arange(len(BINS))); ax.set_xticklabels(BINS) ax.set_ylabel("95% upper limit on the fraction of tested objects\nwith area-to-mass above threshold") ax.set_xlabel("size class") ax.set_title("(b) Upper limits among tested NEOs") ax.set_ylim(2e-4, 1.0) h = [Patch(facecolor=tcol[x], label="A/m $\\geq$ %g m$^2$/kg" % x) for x in THRESH] + \ [Line2D([], [], marker=mk, ls="none", markersize=7.5, markerfacecolor=INK2 if filled else SURF, markeredgecolor=INK2, label=lab) for (_, _, mk, filled, lab) in series] ax.legend(handles=h, fontsize=7.5, ncol=2, loc="upper center", bbox_to_anchor=(0.5, -0.14)) fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig5_limits.png")); plt.close(fig) # -------------------------------------------------------------------- main def scaling_fit(df): rows = [] for name, m in (("arc<=180d", df.arc_obs <= 180), ("arc>180d", df.arc_obs > 180)): s = df[df.selected & m & np.isfinite(df.sig_A1) & (df.sig_A1 > 0) & (df.n_nights > 0) & (df.arc_obs > 0)] if len(s) < 50: continue X = np.column_stack([np.ones(len(s)), np.log10(s.arc_obs), np.log10(s.n_obs), np.log10(s.dist_min), np.log10(s.q)]) y = np.log10(s.sig_A1) beta, *_ = np.linalg.lstsq(X, y, rcond=None) res = y - X @ beta rows.append(dict(sample=name, n=len(s), const=beta[0], slope_log_arc=beta[1], slope_log_nobs=beta[2], slope_log_distmin=beta[3], slope_log_q=beta[4], rms_dex=float(np.std(res)))) out = pd.DataFrame(rows) out.to_csv(os.path.join(DATA, "sigma_scaling_fit.csv"), index=False) return out def load(): plan = pd.read_csv(os.path.join(DATA, "analysis_plan.csv"), dtype={"pdes": str}, low_memory=False) fr = pd.read_csv(os.path.join(DATA, "fisher_results.csv"), dtype={"pdes": str}) df = plan.merge(fr, on="pdes", how="left") df["w"] = np.where(df.selected, 1.0 / df.frac, 0.0) # ESA NEOCC residual-based measurement (multi-apparition arcs) and ESA's fitted area-to-mass ratios nc = pd.read_csv(os.path.join(DATA, "neocc_measurement.csv"), dtype={"pdes": str}) df = df.merge(nc.drop(columns=[c for c in ("why",) if c in nc.columns]), on="pdes", how="left") u_max = float(nc.nc_u_max.iloc[0]) df["nc_status"] = df.nc_status.fillna("") # random part of the NEOCC sample: all analysed multi-apparition objects with the numerical model, # and the two-body objects with u <= u_max (uniform sampling fraction u_max of the catalogue) df["nc_random"] = (df.arc_d > 180) & ((df.model == "numeric") | ((df.model == "kepler") & (df.u <= u_max))) df["w_meas"] = np.where((df.arc_d > 180) & (df.model == "kepler"), 1.0 / u_max, df.w) ng = pd.read_csv(os.path.join(DATA, "neocc_nongrav.csv"), dtype={"name": str}) ng = ng[ng.amr_m2_per_kg.notna()][["name", "amr_m2_per_kg", "amr_sigma"]].rename( columns={"amr_m2_per_kg": "esa_amr", "amr_sigma": "esa_amr_sigma"}) df["name_nospace"] = df.pdes.str.replace(" ", "", regex=False) df = df.merge(ng, left_on="name_nospace", right_on="name", how="left").drop(columns=["name"]) for c in ("m_nl_ok", "m_nl_A1", "m_nl_sig", "m_nl_rms7", "m_nl_n", "m_nl_iter", "m_A1", "m_sig", "m_rms7", "m_rms6", "m_A1_all", "m_sig_all", "m_rms7_all", "m_lono_ul", "m_ref_A1", "n_space", "nc_A1", "nc_sig", "nc_rms7", "nc_rms6", "nc_A1_all", "nc_sig_all", "nc_rms7_all", "nc_lono_ul", "nc2_A1", "nc2_sig", "nc2_rms7", "nc2_rms6", "nc2_A1_all", "nc2_sig_all", "nc2_rms7_all", "nc2_lono_ul", "nc2_lono_zmin", "nc_lono_zmin", "m_lono_zmin", "m_lono_A1", "m_lono_sig", "nc_lono_A1", "nc_lono_sig", "nc2_lono_A1", "nc2_lono_sig", "nc2_n", "nc2_n_all", "nc_a2_val", "nc_a2_sig", "nc_a2_rms7", "nc_a2_dchi2", "nc_a3_val", "nc_a3_sig", "nc_a3_rms7", "nc_a3_dchi2", "nc2_dchi2", "nc_dchi2", "nc_night_z", "nc2_night_z", "nc_a3_night_z", "nc_night_val", "nc_night_sig", "nc2_night_val", "nc2_night_sig", "nc_n_nights_kept", "nc_rms_header", "nc_n_optical", "nc_n_rejected_by_neocc", "nc_nongrav_in_solution"): if c not in df.columns: df[c] = np.nan # primary own-residual estimate: the non-linear refit where it was made and converged nl = (df.m_nl_ok == True).values # noqa: E712 df["m_nonlinear"] = nl df["mp_A1"] = np.where(nl, df.m_nl_A1, df.m_A1) df["mp_sig"] = np.where(nl, df.m_nl_sig * np.maximum(1, df.m_nl_rms7), df.m_sig * np.maximum(1, df.m_rms7)) df["mp_rms7"] = np.where(nl, df.m_nl_rms7, df.m_rms7) return df def _pull_stats(z, extra=None): q = z.quantile([0.0228, 0.16, 0.5, 0.84, 0.9772]) out = dict(n=int(len(z)), z_median=float(q[0.5]), z_p16=float(q[0.16]), z_p84=float(q[0.84]), z_p2=float(q[0.0228]), z_p98=float(q[0.9772]), robust_width=float(0.5 * (q[0.84] - q[0.16])), frac_abs_z_gt3=float((z.abs() > 3).mean()), frac_abs_z_gt5=float((z.abs() > 5).mean()), n_z_ge3=int((z >= 3).sum()), n_z_le_m3=int((z <= -3).sum()), n_z_ge5=int((z >= 5).sum()), n_z_le_m5=int((z <= -5).sum())) if extra: out.update(extra) return out def pulls(d): """Empirical check of the measurement: distribution of measured A1 / sigma (sigma scaled by the post-fit rms when > 1), for (i) own residuals on single-apparition arcs, (ii) own residuals on multi-apparition arcs (dominated by this study's dynamical-model error; shown for comparison only), (iii) ESA NEOCC residuals on multi-apparition arcs with A1 alone and (iv) the same with A2 solved together with A1 (the primary multi-apparition measurement). Outward (z > 0) and sunward (z < 0, unphysical for radiation pressure) tails are counted separately.""" s = d[d.selected & ~d.has_A1].copy() own = s[(s.model == "numeric") & np.isfinite(s.m_A1) & (s.m_sig > 0)].copy() own["z"] = own.mp_A1 / own.mp_sig own["z_linear"] = own.m_A1 / (own.m_sig * np.maximum(1, own.m_rms7)) nc = s[(s.nc_status == "ok") & s.nc_random & np.isfinite(s.nc2_A1) & (s.nc2_sig > 0) & np.isfinite(s.nc_A1)].copy() nc["z_joint"] = nc.nc2_A1 / (nc.nc2_sig * np.maximum(1, nc.nc2_rms7)) nc["z_A1only"] = nc.nc_A1 / (nc.nc_sig * np.maximum(1, nc.nc_rms7)) nc["z"] = np.where(nc.model == "numeric", nc.z_joint, nc.z_A1only) # the primary measurement out = {} g = own[own.arc_d <= 180] out["own_single_apparition"] = _pull_stats(g.z, dict( median_rms6=float(g.m_rms6.median()), median_rms7=float(g.m_rms7.median()), median_rms_jpl=float(g.rms.median()), mean_frac_rejected=float(((g.m_n_all - g.m_n) / g.m_n_all).mean()), n_nonlinear_refit=int(g.m_nonlinear.sum()), n_refit_attempted=int(g.m_nl_ok.notna().sum()), n_z_ge5_robust_to_single_night=int(((g.z >= 5) & (g.m_lono_zmin >= 3)).sum()))) g2 = own[own.arc_d > 180] if len(g2) > 10: out["own_multi_apparition"] = _pull_stats(g2.z, dict(median_rms7=float(g2.m_rms7.median()))) for name, m in (("neocc_multi_numeric_partials", nc.model == "numeric"), ("neocc_multi_twobody_partials", nc.model == "kepler")): g3 = nc[m] if len(g3) > 10: out[name + "_A2_solved"] = _pull_stats(g3.z_joint, dict( median_rms7=float(g3.nc2_rms7.median()), median_rms_header=float(g3.nc_rms_header.median()), mean_frac_rejected=float(((g3.nc2_n_all - g3.nc2_n) / g3.nc2_n_all).mean()), median_sigma_inflation_from_A2=float((g3.nc2_sig / g3.nc_sig).median()), n_z_ge5_robust_to_single_night=int(((g3.z_joint >= 5) & (g3.nc2_lono_zmin >= 3)).sum()), night_compressed=_pull_stats(g3.nc2_night_z.dropna()))) out[name + "_A1_alone"] = _pull_stats(g3.z_A1only, dict( n_z_ge5_robust_to_single_night=int(((g3.z_A1only >= 5) & (g3.nc_lono_zmin >= 3)).sum()), night_compressed=_pull_stats(g3.nc_night_z.dropna()))) # cross-check on single-apparition arcs: own residuals against ESA residuals, same objects x = s[(s.arc_d <= 180) & (s.model == "numeric") & np.isfinite(s.m_A1) & (s.nc_status == "ok") & np.isfinite(s.nc_A1)].copy() if len(x) > 10: dz = (x.m_A1 - x.nc_A1) / np.hypot(x.m_sig, x.nc_sig) out["own_vs_neocc_single"] = dict(n=int(len(x)), median_diff_sigma=float(dz.median()), p16=float(dz.quantile(0.16)), p84=float(dz.quantile(0.84)), frac_abs_gt2=float((dz.abs() > 2).mean()), median_sigma_ratio=float((x.nc_sig / x.m_sig).median()), corr=float(np.corrcoef(x.m_A1 / x.m_sig, x.nc_A1 / x.nc_sig)[0, 1])) # every 5-sigma value, either sign cols = ["pdes", "full_name", "sbin", "H", "a", "e", "i", "arc_d", "n_obs", "dist_min", "rms", "A2"] a_ = own[(own.arc_d <= 180) & (own.z.abs() >= 5)].assign( source=lambda t: np.where(t.m_nonlinear, "own residuals, A1 (non-linear fit)", "own residuals, A1 (linear step)"), amr_hat=lambda t: t.mp_A1 / K, amr_sigma=lambda t: t.mp_sig / K, lono_ul=lambda t: t.m_lono_ul / K, lono_zmin=lambda t: t.m_lono_zmin, rms6=lambda t: t.m_rms6, rms7=lambda t: t.mp_rms7, z_night_compressed=np.nan, z_A3_night_compressed=np.nan, z_A1_alone=np.nan, A2_from_residuals=np.nan, A2_from_residuals_sigma=np.nan, A3_from_residuals=np.nan, A3_from_residuals_sigma=np.nan, dchi2_A1=lambda t: t.m_dchi2, dchi2_A2_alone=np.nan, dchi2_A3_alone=np.nan) jn = (nc.model == "numeric").values for c_ in ("A1", "sig", "rms7", "rms6", "lono_ul", "lono_zmin"): nc["p_" + c_] = np.where(jn, nc["nc2_" + c_], nc["nc_" + c_]) b_ = nc[(nc.z.abs() >= 5) | (nc.z_A1only.abs() >= 5)].assign( source=lambda t: np.where(t.model == "numeric", "ESA NEOCC residuals, A1 with A2 solved", "ESA NEOCC residuals, A1"), amr_hat=lambda t: t.p_A1 / K, amr_sigma=lambda t: t.p_sig * np.maximum(1, t.p_rms7) / K, lono_ul=lambda t: t.p_lono_ul / K, lono_zmin=lambda t: t.p_lono_zmin, rms6=lambda t: t.p_rms6, rms7=lambda t: t.p_rms7, z_night_compressed=lambda t: np.where(t.model == "numeric", t.nc2_night_z, t.nc_night_z), z_A3_night_compressed=lambda t: t.nc_a3_night_z, z_A1_alone=lambda t: t.z_A1only, A2_from_residuals=lambda t: t.nc_a2_val, A2_from_residuals_sigma=lambda t: t.nc_a2_sig, A3_from_residuals=lambda t: t.nc_a3_val, A3_from_residuals_sigma=lambda t: t.nc_a3_sig, dchi2_A1=lambda t: t.nc_dchi2, dchi2_A2_alone=lambda t: t.nc_a2_dchi2, dchi2_A3_alone=lambda t: t.nc_a3_dchi2) sigt = pd.concat([a_, b_])[cols + ["source", "amr_hat", "amr_sigma", "z", "lono_ul", "lono_zmin", "z_night_compressed", "z_A3_night_compressed", "rms6", "rms7", "z_A1_alone", "A2_from_residuals", "A2_from_residuals_sigma", "A3_from_residuals", "A3_from_residuals_sigma", "dchi2_A1", "dchi2_A2_alone", "dchi2_A3_alone"]] sigt = sigt.rename(columns={"A2": "A2_jpl"}).sort_values("z", ascending=False) sigt["robust_outward_candidate"] = (sigt.z >= 5) & (sigt.lono_zmin >= 3) sigt.to_csv(os.path.join(DATA, "measured_5sigma_values.csv"), index=False) # outward values of >= 5 sigma that survive the removal of any single station-night sigt[sigt.robust_outward_candidate].to_csv(os.path.join(DATA, "followup_measured_candidates.csv"), index=False) out["robust_outward_candidates"] = [str(x) for x in sigt[sigt.robust_outward_candidate].pdes] return own, nc, out def a2_check(d): """Validation of the residual-based method on real signals: the transverse parameter A2 estimated from ESA's post-fit residuals with this study's partials, against the A2 fitted by JPL, for analysed objects where JPL reports A2 at >= 3 sigma and ESA's own solution contains no non-gravitational parameter (otherwise the signal is absent from ESA's residuals).""" c = d[np.isfinite(d.A2) & np.isfinite(d.A2_sigma) & (d.A2_sigma > 0) & np.isfinite(d.nc_a2_val) & (d.nc_status == "ok")].copy() if len(c) == 0: return c, {} c["snr_jpl"] = (c.A2 / c.A2_sigma).abs() c["esa_solution_has_nongrav"] = c.nc_nongrav_in_solution.fillna(False).astype(bool) c["A2_res"] = c.nc_a2_val; c["A2_res_sigma"] = c.nc_a2_sig * np.maximum(1, c.nc_a2_rms7) c["ratio"] = c.A2_res / c.A2 c["diff_sigma"] = (c.A2_res - c.A2) / np.hypot(c.A2_res_sigma, c.A2_sigma) t = c[["pdes", "full_name", "sbin", "H", "arc_d", "model", "A2", "A2_sigma", "snr_jpl", "A2_res", "A2_res_sigma", "ratio", "diff_sigma", "esa_solution_has_nongrav"]] g = t[(t.snr_jpl >= 3) & ~t.esa_solution_has_nongrav] st = dict(n=int(len(g))) if len(g): st.update(median_ratio=float(g.ratio.median()), p16=float(g.ratio.quantile(0.16)), p84=float(g.ratio.quantile(0.84)), frac_within_30pct=float(((g.ratio > 0.7) & (g.ratio < 1.3)).mean()), frac_same_sign=float((g.ratio > 0).mean()), median_abs_diff_sigma=float(g.diff_sigma.abs().median()), median_sigma_ratio=float((g.A2_res_sigma / g.A2_sigma).median()), n_numeric=int((g.model == "numeric").sum()), n_twobody=int((g.model == "kepler").sum())) for mod, tag in (("numeric", "numeric"), ("kepler", "twobody")): h = g[g.model == mod] if len(h): st["median_ratio_" + tag] = float(h.ratio.median()) st["frac_within_30pct_" + tag] = float(((h.ratio > 0.7) & (h.ratio < 1.3)).mean()) g0 = t[(t.snr_jpl >= 3) & t.esa_solution_has_nongrav] st["n_esa_has_nongrav"] = int(len(g0)) if len(g0): st["median_ratio_when_esa_has_nongrav"] = float(g0.ratio.median()) return t, st def control_a1(d): """Positive control on real data: the NEOs with a JPL-fitted A1, re-measured here (i) from this study's own residuals about a reference trajectory with A1 switched off and (ii) from ESA NEOCC residuals (where ESA's solution does not itself include an area-to-mass ratio).""" if "ctrl_A1" not in d.columns: return pd.DataFrame(), {} c = d[d.has_A1].copy() c["snr_jpl"] = c.A1 / c.A1_sigma c["own_A1"] = c.ctrl_A1; c["own_sig"] = c.ctrl_sig * np.maximum(1, c.ctrl_rms7) # ESA-residual estimates are void where ESA's own solution already contains an area-to-mass # ratio, and where the linearisation fails (post-fit rms > 3: 2020 CD3, a temporarily captured # object whose heliocentric partials are chaotic) void = np.isfinite(c.esa_amr) | ~(c.nc_rms7 < 3) c["neocc_A1"] = np.where(void, np.nan, c.nc_A1) # A1 alone c["neocc_sig"] = np.where(void, np.nan, c.nc_sig * np.maximum(1, c.nc_rms7)) c["neocc2_A1"] = np.where(void, np.nan, c.nc2_A1) # A1 with A2 solved c["neocc2_sig"] = np.where(void, np.nan, c.nc2_sig * np.maximum(1, c.nc2_rms7)) c["neocc_linearisation_failed"] = ~(c.nc_rms7 < 3) & np.isfinite(c.nc_rms7) for k in ("A1", "A1_sigma", "own_A1", "own_sig", "neocc_A1", "neocc_sig", "neocc2_A1", "neocc2_sig"): c["amr_" + k] = c[k] / K c["jpl_solves"] = np.where(np.isfinite(c.A2), np.where(np.isfinite(c.A3), "A1,A2,A3", "A1,A2"), "A1") t = c[["pdes", "full_name", "H", "arc_d", "jpl_solves", "snr_jpl", "amr_A1", "amr_A1_sigma", "amr_own_A1", "amr_own_sig", "ctrl_rms7", "amr_neocc_A1", "amr_neocc_sig", "nc_rms7", "amr_neocc2_A1", "amr_neocc2_sig", "nc2_rms7", "nc_rms_header", "neocc_linearisation_failed", "esa_amr", "esa_amr_sigma"]].copy() sig = t[t.snr_jpl >= 3] st = dict(n=int(len(t)), n_jpl_3sigma_outward=int(len(sig))) for tag in ("own", "neocc", "neocc2"): g = sig[np.isfinite(sig["amr_%s_A1" % tag])] ratio = g["amr_%s_A1" % tag] / g.amr_A1 st[tag] = dict(n=int(len(g)), median_ratio=float(ratio.median()) if len(g) else np.nan, min_ratio=float(ratio.min()) if len(g) else np.nan, max_ratio=float(ratio.max()) if len(g) else np.nan, n_within_factor2=int(((ratio > 0.5) & (ratio < 2)).sum()), n_detected_3sigma=int((g["amr_%s_A1" % tag] / g["amr_%s_sig" % tag] >= 3).sum())) g = sig[np.isfinite(sig.esa_amr)] st["esa_catalogue"] = dict(n=int(len(g)), median_ratio=float((g.esa_amr / g.amr_A1).median()) if len(g) else np.nan, min_ratio=float((g.esa_amr / g.amr_A1).min()) if len(g) else np.nan, max_ratio=float((g.esa_amr / g.amr_A1).max()) if len(g) else np.nan) return t, st def neocc_summary(): t = pd.read_csv(os.path.join(DATA, "neocc_nongrav.csv"), dtype={"name": str}) a = t[t.amr_m2_per_kg.notna()] return dict(n_with_nongrav=int(len(t)), n_with_amr=int(len(a)), amr={r.name: [float(r.amr_m2_per_kg), float(r.amr_sigma)] for r in a.itertuples()}, amr_min=float(a.amr_m2_per_kg.min()), amr_max=float(a.amr_m2_per_kg.max())) def seligman_check(): path = os.path.join(DATA, "validation_seligman2024.csv") if not os.path.exists(path): return None, {} t = pd.read_csv(path, dtype={"pdes": str}) for c in ("kepler_A1A2_now", "kepler_A1A2_to2024", "numeric_A1A2A3_now", "numeric_A1A2A3_to2024", "kepler_A1_to2024"): t["ratio_" + c] = t[c] / t.sig_A1_published st = {c: dict(median=float(t["ratio_" + c].median()), min=float(t["ratio_" + c].min()), max=float(t["ratio_" + c].max())) for c in ("kepler_A1A2_now", "kepler_A1A2_to2024", "numeric_A1A2A3_now", "numeric_A1A2A3_to2024", "kepler_A1_to2024")} return t, st def main(): print("[4] validation, completeness, measured search, limits, figures") df = load() census = cen.classify(pd.read_csv(os.path.join(DATA, "sbdb_nongrav_all.csv"), dtype={"pdes": str}, low_memory=False) .query("kind.str.startswith('a')", engine="python")) census = refine_census(census) census.drop(columns=[c for c in census.columns if c in ("name", "spec_B", "spec_T")], errors="ignore") \ .to_csv(os.path.join(DATA, "census_nongrav.csv"), index=False) census.class_ng.value_counts().rename_axis("class").reset_index(name="n") \ .to_csv(os.path.join(DATA, "census_summary.csv"), index=False) v, vstats = validation(df) kchk, kstats = kepler_check() selt, selstats = seligman_check() cal_num = vstats.get("A1_numeric", {}).get("median_ratio", 1.0) cal_num = 1.0 / cal_num # JPL / forecast cal_kep = 1.0 / kstats["median_ratio"] if kstats else 1.0 # numeric / kepler cal_kep *= cal_num print(" validation:", json.dumps(vstats, indent=1)) print(" kepler vs numeric:", kstats) print(" forecast vs Seligman+2024 published sigma(A1):", json.dumps(selstats)) print(" calibration factors applied to forecast sigma(A1): numeric %.3f, kepler %.3f" % (cal_num, cal_kep)) d = add_sensitivity(df, cal_num, cal_kep) d[["pdes", "H", "sbin", "ocat", "arc_d", "n_obs_used", "rms", "selected", "w", "w_meas", "model", "n_obs", "n_space", "n_nights", "n_radar", "arc_obs", "dist_min", "res_med", "sig_A1", "sig_A1_opt", "sig_A1_with_A2", "sigA1_cal", "amr_3sig", "meas_source", "amr_T", "amr_hat", "amr_msig", "amr_ul3", "amr_ul3_all", "amr_lono_ul", "amr_ul_best", "m_rms6_used", "m_rms7_used", "nc_rms_header", "nc_n_optical", "nc_n_rejected_by_neocc", "amr_fq_0.7", "amr_fq_1", "fq_fail_rms_0.7", "fq_fail_rms_1", "fq_fail_rej", "vet_ul3"]] \ .to_csv(os.path.join(DATA, "sensitivity_per_object.csv"), index=False) pd.set_option("display.width", 250); pd.set_option("display.max_columns", 40); pd.set_option("display.max_colwidth", 60) # ---------------- forecast completeness (not a limit) comp = completeness(d) comp.to_csv(os.path.join(DATA, "completeness_by_bin.csv"), index=False) print(comp.round(3).to_string(index=False)) # ---------------- measured search and limits pz, pnc, pstats = pulls(d) print(" pulls of the measured A1:", json.dumps(pstats, indent=1)) a2t, a2stats = a2_check(d) if len(a2t): a2t.to_csv(os.path.join(DATA, "validation_A2_from_residuals.csv"), index=False) print(" A2 from ESA residuals vs JPL's fitted A2:", json.dumps(a2stats)) ctab, cstats = control_a1(d) if len(ctab): ctab.to_csv(os.path.join(DATA, "control_jpl_A1.csv"), index=False) print(ctab.round(6).to_string(index=False)) print(" positive control (JPL-fitted A1 re-measured):", json.dumps(cstats)) ml = measured_limits(d) fq = fq_limits(d) lim = pd.concat([ml, fq], ignore_index=True) lim.to_csv(os.path.join(DATA, "limits.csv"), index=False) show = ["method", "arcs", "size_bin", "amr_threshold", "n_catalogued", "n_tested", "n_tested_weighted", "frac_of_bin", "k_flagged", "k_unresolved", "ul95", "ul95_weightedN", "unresolved"] print(lim[show].round(5).to_string(index=False)) meas_f, fq_f = flagged_tables(d) meas_f.to_csv(os.path.join(DATA, "measured_flagged.csv"), index=False) fq_f.to_csv(os.path.join(DATA, "fit_quality_flagged.csv"), index=False) # every catalogued NEO with a poor gravity-only fit, whether or not it was analysed poor = d[d.rms >= 1.3][["pdes", "full_name", "sbin", "H", "arc_d", "n_obs_used", "rms", "ocat", "selected", "model", "n_obs", "amr_T", "amr_hat", "amr_msig", "amr_ul3", "m_rms6_used", "m_rms7_used"]].sort_values("rms", ascending=False) poor.to_csv(os.path.join(DATA, "poor_fits_catalogue.csv"), index=False) print(" catalogued NEOs with normalised rms >= 1.3:", len(poor)) print(poor.round(4).to_string(index=False)) # ---------------- robustness rob = [] def run(tag, dd, bins_col="sbin", strict=False, which=("forecast", "measured", "fq")): if "forecast" in which: c = completeness(dd, bins_col) c.insert(0, "method", "forecast (3 sigma if fitted; not a limit)"); c.insert(0, "variant", tag) c = c.rename(columns={"n_sensitive": "n_tested_weighted", "n_analysed_sensitive": "n_tested", "frac_of_catalogue": "frac_of_bin"}) rob.append(c) if "measured" in which: l = measured_limits(dd, bins_col); l.insert(0, "variant", tag); rob.append(l) if "fq" in which: l = fq_limits(dd, bins_col, strict=strict); l.insert(0, "variant", tag); rob.append(l) run("baseline", d) run("sigma x2 (pessimistic)", add_sensitivity(df, cal_num, cal_kep, scale=2.0)) run("sigma /2, or specular sail C_R = 2", add_sensitivity(df, cal_num, cal_kep, scale=0.5)) for pv, tag in ((PV_LO, "albedo 0.04 for size bins"), (PV_HI, "albedo 0.45 for size bins")): dd = d.copy(); dd["sb2"] = size_bins(dd.H.values, pv) run(tag, dd, "sb2") dd = df.copy(); dd["sig_A1"] = dd.sig_A1_opt run("optical only (no radar)", add_sensitivity(dd, cal_num, cal_kep), which=("forecast", "fq")) dd = df.copy(); dd["sig_A1"] = dd.sig_A1_with_A2 run("A1 and A2 solved together", add_sensitivity(dd, cal_num, cal_kep), which=("forecast", "fq")) dd = add_sensitivity(df, cal_num, cal_kep); dd["amr_3sig"] = dd.amr_3sig * 10.0 / 3.0 run("10-sigma instead of 3-sigma", dd, which=("forecast",)) run("JPL used all MPC astrometry (strict)", d, strict=True, which=("fq",)) run("measured sample at sigma <= X/3", add_sensitivity(df, cal_num, cal_kep, snr_t=3.0), which=("measured",)) run("measured sample at sigma <= X/10", add_sensitivity(df, cal_num, cal_kep, snr_t=10.0), which=("measured",)) run("measured without outlier rejection", add_sensitivity(df, cal_num, cal_kep, clipped=False), which=("measured",)) run("multi-apparition arcs measured from own residuals (dynamical-model error; comparison only)", add_sensitivity(df, cal_num, cal_kep, own_multi=True), which=("measured",)) run("multi-apparition arcs: A1 alone, A2 not solved (Yarkovsky leaks into A1; comparison only)", add_sensitivity(df, cal_num, cal_kep, joint_multi=False), which=("measured",)) dd = d.copy(); dd.loc[dd.n_space.fillna(0) > 0, ["amr_T", "amr_fq_0.7", "amr_fq_1"]] = np.inf run("objects with space-based astrometry excluded", dd, which=("measured", "fq")) rob = pd.concat(rob, ignore_index=True) rob.to_csv(os.path.join(DATA, "robustness_limits.csv"), index=False) # completeness grid (H x arc) hed = [15, 20, 22.75, 25, 27.75, 30, 35]; aed = [1, 3, 10, 30, 180, 3650, 40000] rows = [] for h0, h1 in zip(hed[:-1], hed[1:]): for a0, a1 in zip(aed[:-1], aed[1:]): g = d[(d.H > h0) & (d.H <= h1) & (d.arc_d >= a0) & (d.arc_d < a1)] s = g[g.selected] r = dict(H_lo=h0, H_hi=h1, arc_lo_d=a0, arc_hi_d=a1, n_catalogued=len(g), n_analysed=len(s)) for x in THRESH: r["frac_3sig_amr_%g" % x] = float((s.w * (s.amr_3sig <= x)).sum() / s.w.sum()) if len(s) else np.nan rows.append(r) pd.DataFrame(rows).to_csv(os.path.join(DATA, "completeness_grid.csv"), index=False) sc = scaling_fit(d) print(sc.to_string(index=False)) # follow-up list fu = followup(census, d) fu.to_csv(os.path.join(DATA, "followup_list.csv"), index=False) rb = radial_bounds(census) rb.to_csv(os.path.join(DATA, "dark_comet_radial_bounds.csv"), index=False) print(rb.round(6).to_string(index=False)) sp = special_objects(d, rb) sp.to_csv(os.path.join(DATA, "special_objects.csv"), index=False) print(sp.drop(columns=["note"]).to_string(index=False)) # figures fig_census(census); fig_validation(v, kchk, selt, pz, pnc); fig_sensitivity(d); fig_completeness(d); fig_limits(d, lim) # key numbers oum = census[census.pdes == "2017 U1"].iloc[0] s = d[d.selected] med = {} for b in BINS: g = d[d.sbin == b] an = g[g.selected & np.isfinite(g.amr_3sig)] # catalogue-weighted median with unanalysable objects (arc < 1 d etc.) counted as insensitive vals = np.concatenate([np.repeat(an.amr_3sig.values, np.round(an.w.values).astype(int)), np.full(max(len(g) - int(np.round(an.w.sum())), 0), np.inf)]) med[b] = dict(median_analysed=float(an.amr_3sig.median()), median_catalogue=float(np.median(vals)), median_arc_d=float(g.arc_d.median())) yk = census[census.class_ng == "Yarkovsky-consistent"] kn = dict( n_census=int(len(census)), census_classes=census.class_ng.value_counts().to_dict(), n_significant=int(census.significant.sum()), yarkovsky_over_nominal=int(yk.flag_A2_over_nominal.sum()), yarkovsky_over_nominal_3sig=[str(x) for x in yk[yk.flag_A2_over_nominal_3sig].full_name], yarkovsky_ratio_to_nominal_quantiles=[float(x) for x in yk.A2_over_nominal.quantile([0.05, 0.5, 0.95])], yarkovsky_frac_inward=float((yk.A2 < 0).mean()), n_neo_catalogued=int(len(df)), n_analysed=int(d.selected.sum()), model_counts=d[d.selected].model.value_counts().to_dict(), n_measured=int((s.amr_T < np.inf).sum()), n_measured_by_source=s[s.amr_T < np.inf].meas_source.value_counts().to_dict(), neocc=dict(n_catalogue=int(len(pd.read_csv(os.path.join(DATA, "neocc_measurement.csv")))), status=pd.read_csv(os.path.join(DATA, "neocc_measurement.csv")).nc_status.value_counts().to_dict(), nongrav=neocc_summary()), n_space_obs_objects=int((s.n_space.fillna(0) > 0).sum()), n_space_obs=int(s.n_space.fillna(0).sum()), n_obs_total=int(s.n_obs.fillna(0).sum()), validation=vstats, kepler_vs_numeric=kstats, seligman2024=selstats, cal_numeric=cal_num, cal_kepler=cal_kep, pulls=pstats, control_jpl_A1=cstats, a2_from_residuals_vs_jpl=a2stats, res_med_numeric_arcsec=float(d[d.model == "numeric"].res_med.median()), res_med_space_arcsec=float(d[d.model == "numeric"].res_med_space.median()) if "res_med_space" in d else None, A1_per_amr=A1_PER_AMR, oumuamua_A1_jpl=float(oum.A1), oumuamua_A1_sigma=float(oum.A1_sigma), oumuamua_amr_jpl=float(oum.A1 / A1_PER_AMR), oumuamua_amr_micheli=4.92e-6 / P_SRP, oumuamua_beta=float(oum.A1 / MU_SUN), dark_comet_amr_total_range=[float(x) for x in census[census.class_ng == "published dark comet"].amr_total_equiv.quantile([0, 0.5, 1])], artificial=[dict(name=a["name"], A1=a["A1"], amr=a["A1"] / A1_PER_AMR) for a in cen.ARTIFICIAL], max_amr_equiv_neo=float(census[(census.neo == "Y") & (census.snr_A1 >= 3)].amr_equiv.max()), frac_few_rejected=float((s.n_obs_used >= 0.8 * s.n_obs).mean()), median_sensitivity=med, overall_forecast_completeness={str(x): float((s.w * (s.amr_3sig <= x)).sum() / len(d)) for x in THRESH}, tested_fraction_measured={b: {str(x): float(ml[(ml.arcs == "all arcs") & (ml.size_bin == b) & (ml.amr_threshold == x)].frac_of_bin.iloc[0]) for x in THRESH} for b in BINS}, untested_fraction_fit_quality={b: {str(x): 1 - float(fq[(fq.method == "fit quality R=0.7") & (fq.arcs == "all arcs") & (fq.size_bin == b) & (fq.amr_threshold == x)].frac_of_bin.iloc[0]) for x in THRESH} for b in BINS}, n_poor_fits_rms_ge_1p3=int(len(poor)), linearisation_margin=LIN_MARGIN, ) if os.path.exists(os.path.join(DATA, "validation_oumuamua.csv")): o = pd.read_csv(os.path.join(DATA, "validation_oumuamua.csv")).iloc[0].to_dict() if "A1_measured" in o: o["amr_measured"] = o["A1_measured"] / A1_PER_AMR o["amr_measured_sigma"] = o["A1_measured_sigma"] / A1_PER_AMR o["amr_measured_ground"] = o["A1_measured_ground"] / A1_PER_AMR o["amr_measured_ground_sigma"] = o["A1_measured_ground_sigma"] / A1_PER_AMR kn["oumuamua_fisher"] = o with open(os.path.join(DATA, "key_numbers.json"), "w") as f: json.dump(kn, f, indent=1, default=float) return d, census, lim, rob, kn def refine_census(census): """Literature-based refinements of the automatic classification (see report).""" c = census.copy() # (875163) 1998 SH2: predicted from its orbit and then confirmed as a comet # (Farnocchia et al. 2026, Nat. Astron., doi:10.1038/s41550-026-02913-7) -> P/1998 SH2 m = c.pdes == "875163" c.loc[m, "class_ng"] = "published dark comet" c.loc[m, "lit_note"] = "P/1998 SH2: activity confirmed (Farnocchia et al. 2026, Nat. Astron.)" # (4179) Toutatis: 3.4-sigma *sunward* A1 of 1e-9 of solar gravity; a body of this size # cannot show SRP or outgassing at this level -> treated as model/systematic error. m = c.pdes == "4179" c.loc[m, "class_ng"] = "marginal sunward A1 (likely systematic)" # (883607) 2016 TA56: the transverse term was fitted by JPL (manual solution by D. Farnocchia, # 2026-06-17); comet-like orbit (T_J = 2.69). Not in the published dark-comet lists; no paper # found (web and arXiv search 2026-10-12). What is new here is only noticing that. m = c.pdes == "883607" c.loc[m, "class_ng"] = "unpublished outgassing candidate (JPL-fitted A2)" c.loc[m, "lit_note"] = ("A2 fitted by JPL (D. Farnocchia, 2026-06-17); Jupiter-family-like orbit; " "absent from Seligman+2023/2024 and Farnocchia+2023/2026") left = c.class_ng == "non-radiative, not in published lists" if left.any(): print(" WARNING: unreviewed non-radiative objects:", list(c[left].full_name)) return c # radial components of dark comets not fitted in SBDB: Seligman et al. 2024 (PNAS), Table 2 SELIGMAN24_A1 = {"2001 ME1": (-2.47e-13, 1.37e-10), "2005 UY6": (-2.92e-10, 4.28e-10), "1998 FR11": (-1.83e-10, 2.96e-10), "2012 UR158": (8.18e-11, 6.66e-10), "2003 RM": (-3.96e-11, 1.10e-10), "2010 VL65": (6.59e-10, 1.30e-9)} def radial_bounds(census): """Radial acceleration of each published dark comet (and of the unpublished candidate) expressed as an area-to-mass ratio, with the 3-sigma upper bound and the corresponding minimum mass per unit sunlit area.""" rows = [] c = census[census.class_ng.str.startswith(("published dark comet", "unpublished outgassing"))] seen = set() for r in c.itertuples(): prov = cen.PRIMARY_TO_PROV.get(r.pdes, r.pdes) prov = {"875163": "1998 SH2", "883607": "2016 TA56"}.get(prov, prov) seen.add(prov) if np.isfinite(r.A1): a1, s1, src = r.A1, r.A1_sigma, "JPL SBDB" elif prov in SELIGMAN24_A1: a1, s1 = SELIGMAN24_A1[prov]; src = "Seligman+2024 Table 2" else: a1, s1, src = np.nan, np.nan, "A1 not estimated" rows.append(dict(object=prov, H=r.H, D_nom_m=r.D_nom, A1=a1, A1_sigma=s1, source=src)) for prov in ("2010 VL65",): if prov not in seen: a1, s1 = SELIGMAN24_A1[prov] rows.append(dict(object=prov, H=29.22, D_nom_m=float(diam_m(29.22)), A1=a1, A1_sigma=s1, source="Seligman+2024 Table 2")) t = pd.DataFrame(rows) t["amr_equiv"] = t.A1 / A1_PER_AMR t["amr_sigma"] = t.A1_sigma / A1_PER_AMR t["amr_ul3"] = (t.A1 + 3 * t.A1_sigma).clip(lower=0) / A1_PER_AMR t["min_mass_per_area_kg_m2"] = 1.0 / t.amr_ul3 t["amr_rock_2000"] = amr_sphere(t.D_nom_m, 2000.0) t["excludes_rocket_stage_0p008"] = t.amr_ul3 < 0.0079 return t.sort_values("H", ascending=False) # Objects still carried as asteroids whose artificial origin has been discussed in the # literature, and published spacecraft hypotheses for two "dark comets". SUSPECTS = { "1991 VG": "Steel 1995 (Observatory 115, 78) discussed an artefact origin; Earth-like orbit", "2000 SG344": "Chodas & Chesley 2001 (BAAS): possible Apollo-era S-IVB stage; Earth-like orbit", "2019 EF3": "metre-class object with a 6-year linkage and a poor gravity-only fit (referee's flag)", "2007 AA2": "this study's only outward >= 5-sigma value that survives single-night removal (two apparitions, " "2007 and 2026; an out-of-plane A3 fits equally well); Earth-like orbit", } # object, hypothesis, variant, equivalent radial A/m [m2/kg] claimed or implied, reference. # Two variants each: (a) the full panel area permanently facing the Sun with C_R = 1 (the most the # hypothesis allows; NOT claimed by its authors), (b) the value the authors themselves adopt. SPACECRAFT_HYPOTHESES = [ ("2005 VL1", "Venera 2", "10 m2 Sun-facing, C_R = 1 (upper case, not the authors' claim)", 10.0 / 963.0, "Loeb & Cloete 2025, arXiv:2503.03552 (M = 0.963e6 g, S ~ 1e5 cm2)"), ("2005 VL1", "Venera 2", "authors' adopted |A| = 1.5e-9 au/d2 (their reflectance ~0.3), taken outward", 1.5e-9 / A1_PER_AMR, "Loeb & Cloete 2025, arXiv:2503.03552, their eq. 1 and text"), ("1998 KY26", "Phobos 1", "10 m2 Sun-facing, C_R = 1 (upper case, not the authors' claim)", 10.0 / 6220.0, "Hibberd et al. 2026, arXiv:2609.24283 (M = 6220 kg, panels ~10 m2)"), ("1998 KY26", "Phobos 1", "authors' effective area 4.5 m2 (from the 3.3e-9 m/s2 resultant of Seligman+2024)", 4.5 / 6220.0, "Hibberd et al. 2026, arXiv:2609.24283, their eq. 4"), ] def special_objects(d, rb): rows = [] for p, note in SUSPECTS.items(): g = d[d.pdes == p] if len(g) == 0: continue r = g.iloc[0] rows.append(dict(object=p, kind="suspect or flagged object still catalogued as asteroid", H=r.H, arc_d=r.arc_d, n_obs_used=r.n_obs_used, n_obs_mpc=r.n_obs, rms_norm=r.rms, sigma_A1_forecast=r.sigA1_cal, amr_3sigma_if_fitted=r.amr_3sig, amr_fit_quality_R0p7=r["amr_fq_0.7"], amr_measured=r.amr_hat, amr_measured_sigma=r.amr_msig, amr_measured_ul3=r.amr_ul3, measurement=r.meas_source, tested_down_to_amr=r.amr_T, amr_hypothesis=np.nan, tension_sigma=np.nan, note=note)) for obj, hyp, variant, amr_h, ref in SPACECRAFT_HYPOTHESES: g = rb[rb.object == obj] if len(g) == 0: continue r = g.iloc[0] rows.append(dict(object=obj, kind="published spacecraft hypothesis: %s; %s" % (hyp, variant), H=r.H, arc_d=np.nan, n_obs_used=np.nan, n_obs_mpc=np.nan, rms_norm=np.nan, sigma_A1_forecast=r.A1_sigma, amr_3sigma_if_fitted=3 * r.amr_sigma, amr_fit_quality_R0p7=np.nan, amr_measured=r.amr_equiv, amr_measured_sigma=r.amr_sigma, amr_measured_ul3=r.amr_ul3, measurement=r.source, tested_down_to_amr=np.nan, amr_hypothesis=amr_h, tension_sigma=(amr_h - r.amr_equiv) / r.amr_sigma, note=ref)) return pd.DataFrame(rows) def followup(census, d): c = census.copy() keep = c[c.class_ng.str.startswith(("unpublished", "published dark comet", "marginal")) | c.pdes.isin(["2020 GE", "2009 BD", "2020 CD3", "2011 CL50", "428209", "613995"])] cols = ["full_name", "H", "D_nom", "q", "e", "T_jup", "data_arc", "n_obs_used", "producer", "soln_date", "A1", "snr_A1", "A2", "snr_A2", "A3", "snr_A3", "amr_equiv", "amr_total_equiv", "rho_implied", "A2_yark_max", "A2_over_nominal", "class_ng", "lit_note"] return keep[cols].sort_values("amr_total_equiv", ascending=False) if __name__ == "__main__": main()