"""Step 3: (a) covariance forecast of sigma(A1) for the catalogued NEO population (and sigma(A2) for validation), using the real epochs/stations of the MPC astrometry and JPL radar astrometry; (b) [revision] a direct, linearised *measurement* of A1 from the real astrometric residuals about JPL's trajectory for every object treated with the numerical model. Outputs: data/fisher_results.csv (one row per analysed object; columns m_* hold the measurement), data/validation_kepler_vs_numeric.csv, data/vetting_numeric.csv, data/validation_seligman2024.csv, data/validation_oumuamua.csv Usage: python s03_fisher.py [--limit=N] [--no-fetch] """ import os import sys import time import warnings import numpy as np import pandas as pd import orbfisher as of import s01_download as dl from common import DATA, RAW warnings.filterwarnings("ignore") PLANETS_ON = True GAIA_CATS = set("UVWXYZ") # MPC catalogue codes of Gaia DR1/DR2/EDR3/DR3... LCO = {"K92", "K93", "Q63", "Q64", "V37", "W84", "W85", "W86", "W87", "K91", "E10", "F65"} FIX = {"704": 1.0, "G96": 0.5, "F51": 0.2, "F52": 0.2, "G45": 0.6, "699": 0.8, "D29": 0.75, "C51": 1.0, "E12": 0.75, "608": 0.6, "J75": 1.0, "645": 0.3, "673": 0.3, "689": 0.5, "950": 0.5, "H01": 0.3, "J04": 0.4, "Y28": 0.3} def assign_sigma(mjd, stn, note2, cat, lon): """Per-observation 1-sigma [arcsec] following the weighting scheme of Veres et al. (2017, Icarus 296, 139; Tables 2-5), including their sqrt(N/4) de-weighting of >4 observations per station per night.""" n = len(mjd) sig = np.full(n, 1.0) for k in range(n): s = stn[k]; t = mjd[k]; c = cat[k]; n2 = note2[k] if n2 in "PAN": v = 10.0 if t < 11368 else (5.0 if t < 33282 else 2.5) # 1890-01-01, 1950-01-01 elif n2 == "E" or n2 == "H": v = 0.2 elif n2 == "T": v = 0.5 elif n2 == "e": v = 0.75 elif n2 == "M": v = 2.0 elif s in FIX: v = FIX[s] elif n2 == "S": v = 1.5 elif s == "703": v = 1.0 if t < 56658 else 0.8 elif s == "691": v = 0.6 if t < 52640 else 0.5 elif s == "644": v = 0.6 if t < 52883 else 0.4 elif s == "568": v = 0.1 if c in GAIA_CATS else (0.2 if c == "t" else 0.5) elif s in ("T09", "T12", "T14"): v = 0.1 if c in GAIA_CATS else 1.0 elif s in ("309", "G83"): v = 0.2 if c in GAIA_CATS else 0.3 elif s in LCO: v = 0.5 else: v = 1.0 sig[k] = v night = np.floor(mjd + lon / 360.0 - 0.5).astype(np.int64) key = np.array(["%s%d" % (a, b) for a, b in zip(stn, night)]) _, inv, cnt = np.unique(key, return_inverse=True, return_counts=True) nn = cnt[inv] sig = sig * np.where(nn > 4, np.sqrt(nn / 4.0), 1.0) n_nights = len(np.unique(night)) return sig, n_nights, inv class ObsStore: def __init__(self): t0 = time.time() df = pd.read_csv(os.path.join(RAW, dl.OBS_FILE), dtype={"pdes": str, "note2": str, "cat": str, "stn": str}, keep_default_na=False, na_values={"ox": [""], "oy": [""], "oz": [""]}) df = df.sort_values(["pdes", "mjd_utc"], kind="stable").reset_index(drop=True) bounds = np.flatnonzero(np.r_[True, df.pdes.values[1:] != df.pdes.values[:-1], True]) names = df.pdes.values[bounds[:-1]] self.idx = {n: (bounds[i], bounds[i + 1]) for i, n in enumerate(names)} self.mjd = df.mjd_utc.values; self.ra = df.ra.values; self.dec = df.dec.values self.note2 = df.note2.values; self.cat = df.cat.values; self.stn = df.stn.values self.xyz = df[["ox", "oy", "oz"]].values.astype(float) oc = pd.read_csv(os.path.join(DATA, "mpc_obscodes.csv"), dtype={"code": str}) self.oc = {r.code: (r.lon_deg, r.rho_cos, r.rho_sin) for r in oc.itertuples()} rad = pd.read_csv(os.path.join(DATA, "jpl_radar_astrometry.csv"), dtype={"des": str}) rad["mjd"] = (pd.to_datetime(rad.epoch) - pd.Timestamp("1858-11-17")) / pd.Timedelta(days=1) self.radar = {k: g for k, g in rad.groupby("des")} print(" observation store: %d optical obs of %d objects (%d space-based), radar for %d (%.0f s)" % (len(df), len(self.idx), int(np.isfinite(self.xyz[:, 0]).sum()), len(self.radar), time.time() - t0)) def get(self, pdes, numbered=False, api=False, t_max=None, hst_sigma=None): """Astrometry of one object. Ground stations are placed with the MPC parallax constants; space-based observers with the position vector of their 's' record. Records from stations with neither are dropped, as are exact duplicates (same station, epochs within 2 s).""" if (numbered or api) or pdes not in self.idx: if not (numbered or api): return None rows = dl.mpc_api_obs(pdes) if not rows: return None mjd = np.array([r[0] for r in rows]); ra = np.array([r[1] for r in rows]) dec = np.array([r[2] for r in rows]); note2 = np.array([r[3] for r in rows]) cat = np.array([r[4] for r in rows]); stn = np.array([r[5] for r in rows]) xyz = np.array([r[6:9] for r in rows], float) o = np.argsort(mjd, kind="stable") mjd, ra, dec, note2, cat, stn, xyz = mjd[o], ra[o], dec[o], note2[o], cat[o], stn[o], xyz[o] else: a, b = self.idx[pdes] mjd, ra, dec = self.mjd[a:b], self.ra[a:b], self.dec[a:b] note2, cat, stn, xyz = self.note2[a:b], self.cat[a:b], self.stn[a:b], self.xyz[a:b] n_raw = len(mjd) is_sat = np.isfinite(xyz[:, 0]) known = np.array([s_ in self.oc for s_ in stn], bool) keep = is_sat | known # duplicates: same station, epochs within 2e-5 d o = np.lexsort((mjd, stn)) dup = np.zeros(n_raw, bool) dup[o[1:]] = (stn[o[1:]] == stn[o[:-1]]) & (np.abs(mjd[o[1:]] - mjd[o[:-1]]) < 2e-5) keep &= ~dup if t_max is not None: keep &= mjd <= t_max mjd, ra, dec, note2, cat, stn, xyz, is_sat = (v[keep] for v in (mjd, ra, dec, note2, cat, stn, xyz, is_sat)) if len(mjd) == 0: return None par = np.array([self.oc.get(s_, (0.0, 0.0, 0.0)) for s_ in stn]) sig, n_nights, grp = assign_sigma(mjd, stn, note2, cat, par[:, 0]) if hst_sigma is not None: sig = np.where(stn == "250", hst_sigma, sig) off = of.observer_offset(mjd, par[:, 0], par[:, 1], par[:, 2]) off[is_sat] = xyz[is_sat] return dict(mjd_utc=mjd, ra=ra, dec=dec, sig=sig, off=off, n_nights=n_nights, stn=stn, grp=grp, n_space=int(is_sat.sum()), n_raw=n_raw) def get_radar(self, pdes): g = self.radar.get(pdes) if g is None: return None return dict(mjd_utc=g.mjd.values.astype(float), kind=np.where(g.units.values == "us", "d", "D"), sig=g.sigma.values.astype(float), freq=g.freq.values.astype(float)) def hz_state(hz, pdes, jd): """Stored JPL Horizons state of `pdes` nearest in time to `jd` -> (state, jd_ref) or (None, None). Any epoch defines the same JPL trajectory, so the nearest stored one is used.""" lst = hz.get(pdes) if not lst: return None, None j = int(np.argmin([abs(x[0] - jd) for x in lst])) return lst[j][1], lst[j][0] def analyse(row, store, hz, force_numeric=False, use_api=False, t_max=None, measure=True, ng_override=None): """row: Series from the plan (+ A1,A2,A3 columns). Returns dict of results.""" pdes = row.pdes obs = store.get(pdes, numbered=(row.kind == "an"), api=use_api or bool(row.post_freeze), t_max=t_max) res = dict(pdes=pdes, model="", n_obs=0) if obs is None or len(obs["mjd_utc"]) < 3: res["model"] = "no_obs" return res radar = store.get_radar(pdes) if radar is not None and t_max is not None: m = radar["mjd_utc"] <= t_max radar = {k: v[m] for k, v in radar.items()} if m.any() else None numeric = force_numeric or row.ocat in ("short_inarc", "short_std") or row.sbin in ("1-10 m", "10-100 m") ng = [0.0 if not np.isfinite(x) else float(x) for x in (row.A1, row.A2, row.A3)] if ng_override is not None: ng = list(ng_override) t_first, t_last = obs["mjd_utc"][0], obs["mjd_utc"][-1] try: if numeric and (row.ocat != "short_inarc"): state0, jd0 = hz_state(hz, pdes, float(row.jd_hz)) if state0 is None: res["model"] = "no_state" return res t_ref = jd0 - 2400000.5 else: state0 = of.elements_to_state(row.a, row.e, row.i, row.om, row.w, row.ma) t_ref = float(row.epoch_mjd) if not np.all(np.isfinite(state0)): res["model"] = "no_state" return res if numeric: out = of.numeric_fisher(state0, t_ref, obs, radar, ng=ng, planets=PLANETS_ON, rtol=1e-10) res["model"] = "numeric" else: out = of.kepler_fisher(state0, t_ref, obs, radar) res["model"] = "kepler" except Exception as exc: # noqa res["model"] = "failed:" + type(exc).__name__ return res N, No = out["N"], out["N_opt"] base = [0, 1, 2, 3, 4, 5] res.update( n_obs=len(obs["mjd_utc"]), n_nights=obs["n_nights"], arc_obs=t_last - t_first, n_space=obs["n_space"], n_raw=obs["n_raw"], n_radar=0 if radar is None else len(radar["mjd_utc"]), sig_A1=of.marginal_sigma(N, base + [6], 6), sig_A1_opt=of.marginal_sigma(No, base + [6], 6), sig_A2=of.marginal_sigma(N, base + [7], 7), sig_A1_with_A2=of.marginal_sigma(N, base + [6, 7], 6), dist_min=float(out["dist"].min()), dist_med=float(np.median(out["dist"])), res_med=float(np.median(np.hypot(out["res_ra"], out["res_dec"]))), res_p90=float(np.percentile(np.hypot(out["res_ra"], out["res_dec"]), 90)), nfev=out["nfev"], int_ok=out["ok"]) if res["model"] == "numeric": rr = np.hypot(out["res_ra"], out["res_dec"]) sat = np.zeros(len(rr), bool) if obs["n_space"] == 0 else np.linalg.norm(obs["off"], axis=1) > 1.02 * of.R_EARTH_AU res["res_med_space"] = float(np.median(rr[sat])) if sat.any() else np.nan res["sig_A1_with_A2A3"] = of.marginal_sigma(N, base + [6, 7, 8], 6) # JPL's own parameter set (for the validation against catalogued sigmas) sel = [6 + j for j, x in enumerate((row.A1, row.A2, row.A3)) if np.isfinite(x)] if res["model"] == "kepler": sel = [s for s in sel if s < 8] if sel: for j, nm in ((6, "A1"), (7, "A2"), (8, "A3")): if j in sel: res["sig_%s_jplset" % nm] = of.marginal_sigma(N, base + sel, j) # ---- direct measurement: 6 state corrections + A1 fitted to the real optical residuals if measure and res["model"] == "numeric" and len(obs["mjd_utc"]) >= 5: try: m = of.fit_a1(out["Ha"], out["Hd"], out["res_ra"], out["res_dec"], obs["sig"], groups=obs["grp"], a1_ref=ng[0]) for k in ("m_A1", "m_A1_all"): if k in m: m[k] += ng[0] # total A1 = reference value + fitted correction res.update(m) res["m_ref_A1"] = ng[0] # significant linear estimates on single-apparition arcs are re-estimated with a full # (iterated, non-linear) fit: the linear step is biased for large signals if not any(ng) and (t_last - t_first) <= 200.0: res.update(refine_nonlinear(state0, t_ref, obs, m)) except Exception as exc: # noqa res["m_fail"] = type(exc).__name__ return res def load_plan(): plan = pd.read_csv(os.path.join(DATA, "analysis_plan.csv"), dtype={"pdes": str}, low_memory=False) plan["jd_hz"] = 0.5 * (plan.mjd_first + plan.mjd_last) + 2400000.5 return plan def load_hz(): path = os.path.join(DATA, "horizons_states.csv") hz = {} if os.path.exists(path): h = pd.read_csv(path, dtype={"pdes": str}) h = h[np.isfinite(h.x)] for r in h.itertuples(): hz.setdefault(r.pdes, []).append((float(r.jd_ref), np.array([r.x, r.y, r.z, r.vx, r.vy, r.vz]))) return hz def _median_epoch_jd(store, row): obs = store.get(row.pdes, numbered=(row.kind == "an"), api=bool(row.post_freeze)) if obs is None or len(obs["mjd_utc"]) < 3: return np.nan return round(float(np.median(obs["mjd_utc"])) + 2400000.5, 3) def _fetch_missing(need, fetch): """Fetch JPL Horizons states only for objects that have none stored yet.""" hz = load_hz() todo = need[~need.pdes.isin(set(hz))].dropna() if len(todo) and fetch: dl.fetch_horizons(todo) hz = load_hz() return hz # Seligman et al. 2024 (PNAS), Table 2: formal sigma(A1) [au/d^2] of the five multi-apparition # ("outer") dark comets, all fitted with A1, A2 and A3. Their data ended in 2024. SELIGMAN24_SIGA1 = {"139359": ("2001 ME1", 1.37e-10), "452639": ("2005 UY6", 4.28e-10), "152667": ("1998 FR11", 2.96e-10), "2012 UR158": ("2012 UR158", 6.66e-10), "523599": ("2003 RM", 1.10e-10)} SELIGMAN_TMAX = 60491.0 # 2024-07-01, approximate end of the data used in that paper def main(limit=None, fetch=True): print("[3] Fisher forecast + direct measurement of A1") plan = load_plan() store = ObsStore() sel = plan[plan.selected].copy() # The numeric model is also used for long-arc objects smaller than 100 m (observed during # close Earth approaches, where two-body partials fail) and for every object with a fitted # A1; their reference state is taken near the median observation epoch. ex = (sel.sbin.isin(["1-10 m", "10-100 m"]) | sel.A1.notna()) & ~sel.ocat.isin(["short_inarc", "short_std"]) for idx in sel.index[ex]: sel.loc[idx, "jd_hz"] = _median_epoch_jd(store, sel.loc[idx]) # numeric-vs-Kepler cross-check sample: 40 random long-arc unnumbered objects (>100 m) chk = sel[(sel.ocat == "long_unnumbered") & ~ex].sort_values("u").head(40).copy() for idx in chk.index: chk.loc[idx, "jd_hz"] = _median_epoch_jd(store, chk.loc[idx]) need = pd.concat([plan[plan.need_horizons][["pdes", "jd_ref"]], sel[ex][["pdes", "jd_hz"]].rename(columns={"jd_hz": "jd_ref"}), chk[["pdes", "jd_hz"]].rename(columns={"jd_hz": "jd_ref"})]).dropna() hz = _fetch_missing(need, fetch) if limit: sel = sel.sample(limit, random_state=1) rows = [] t0 = time.time() for k, row in enumerate(sel.itertuples(index=False)): r = analyse(row, store, hz, force_numeric=bool(np.isfinite(row.A1))) # positive control on real data: objects for which JPL fitted A1 are re-measured from a # reference trajectory with A1 switched off (JPL state, A2 and A3 kept) if np.isfinite(row.A1) and r.get("model") == "numeric": c = analyse(row, store, hz, force_numeric=True, ng_override=[0.0] + [0.0 if not np.isfinite(x) else float(x) for x in (row.A2, row.A3)]) for kk in ("m_A1", "m_sig", "m_rms6", "m_rms7"): r["ctrl_" + kk[2:]] = c.get(kk, np.nan) rows.append(r) if (k + 1) % 2000 == 0: print(" %d/%d (%.0f s)" % (k + 1, len(sel), time.time() - t0), flush=True) res = pd.DataFrame(rows) res.to_csv(os.path.join(DATA, "fisher_results.csv"), index=False) print(" wrote fisher_results.csv", len(res), "rows;", res.model.value_counts().to_dict(), " (%.0f s)" % (time.time() - t0)) if limit: return res # ---- cross-check: numeric model on long arcs that the main run treats with Kepler rows = [] for row in chk.itertuples(index=False): r = analyse(row, store, hz, force_numeric=True) r["sig_A1_kepler"] = float(res.loc[res.pdes == row.pdes, "sig_A1"].iloc[0]) if (res.pdes == row.pdes).any() else np.nan rows.append(r) pd.DataFrame(rows).to_csv(os.path.join(DATA, "validation_kepler_vs_numeric.csv"), index=False) # ---- vetting: two-body-model objects whose JPL fit is poor (rms >= 1.4) or for which JPL used # under 80% of the MPC astrometry are re-analysed with the numerical model so that A1 can # be measured from their real residuals (see s04_limits: fit-quality criterion) m = sel.merge(res[["pdes", "model", "n_obs"]], on="pdes") vet = m[(m.model == "kepler") & ((m.rms >= 1.4) | (m.n_obs_used < 0.8 * m.n_obs))].copy() for idx in vet.index: vet.loc[idx, "jd_hz"] = _median_epoch_jd(store, vet.loc[idx]) hz = _fetch_missing(vet[["pdes", "jd_hz"]].rename(columns={"jd_hz": "jd_ref"}), fetch) rows = [] for row in vet.drop(columns=["model", "n_obs"]).itertuples(index=False): r = analyse(row, store, hz, force_numeric=True) r.update(rms_jpl=row.rms, n_obs_used=row.n_obs_used, sbin=row.sbin, H=row.H) rows.append(r) pd.DataFrame(rows).to_csv(os.path.join(DATA, "vetting_numeric.csv"), index=False) print(" vetting run (numeric model on %d two-body objects): %.0f s" % (len(rows), time.time() - t0)) # ---- long-arc forecast against published sigma(A1) (Seligman et al. 2024, Table 2) sg = plan[plan.pdes.isin(SELIGMAN24_SIGA1)].copy() for idx in sg.index: sg.loc[idx, "jd_hz"] = _median_epoch_jd(store, sg.loc[idx]) hz = _fetch_missing(sg[["pdes", "jd_hz"]].rename(columns={"jd_hz": "jd_ref"}), fetch) rows = [] for row in sg.itertuples(index=False): name, pub = SELIGMAN24_SIGA1[row.pdes] r = dict(pdes=row.pdes, object=name, sig_A1_published=pub) for tag, tmax in (("now", None), ("to2024", SELIGMAN_TMAX)): k = analyse(row, store, hz, force_numeric=False, t_max=tmax, measure=False) nmr = analyse(row, store, hz, force_numeric=True, t_max=tmax, measure=False) r.update({"kepler_A1_%s" % tag: k.get("sig_A1"), "kepler_A1A2_%s" % tag: k.get("sig_A1_with_A2"), "numeric_A1A2A3_%s" % tag: nmr.get("sig_A1_with_A2A3"), "n_obs_%s" % tag: k.get("n_obs"), "model_%s" % tag: k.get("model")}) rows.append(r) pd.DataFrame(rows).to_csv(os.path.join(DATA, "validation_seligman2024.csv"), index=False) # ---- 1I/'Oumuamua (not a NEO-catalogue member): same machinery, JPL solution as reference try: oum = oumuamua(store) pd.DataFrame([oum]).to_csv(os.path.join(DATA, "validation_oumuamua.csv"), index=False) print(" 'Oumuamua: forecast sigma(A1) with A1,A2,A3 solved = %.2e (generic weights), %.2e (HST at 0.05 arcsec); " "JPL formal = %.2e au/d^2" % (oum["sig_A1_jplset"], oum["sig_A1_jplset_hst0p05"], oum["A1_sigma_jpl"])) print(" 'Oumuamua positive control: A1 measured from a gravity-only start = %.3e +- %.3e au/d^2 " "(Micheli et al. 2018: %.3e +- %.3e)" % (oum["A1_measured"], oum["A1_measured_sigma"], oum["A1_micheli2018"], oum["A1_sigma_micheli2018"])) except Exception as exc: # noqa print(" 'Oumuamua validation failed:", exc) return res NL_Z = 4.0 # linear estimates with |A1/sigma| >= NL_Z are refined by a full non-linear fit def refine_nonlinear(state0, t_ref, obs, m, n_iter=12): """Full Gauss-Newton fit (6 state components + A1, 3-sigma clipping on the residuals of the model WITH A1) started from the gravity-only reference orbit, for an object whose linearised estimate m (output of orbfisher.fit_a1) is significant. Returns m_nl_* entries, or {}.""" if "m_A1" not in m or not m.get("m_sig", 0) > 0: return {} z = m["m_A1"] / (m["m_sig"] * max(1.0, m["m_rms7"])) if not np.isfinite(z) or abs(z) < NL_Z: return {} try: g = gauss_newton(state0, t_ref, obs, a1=0.0, solve_a1=True, n_iter=n_iter, clip=3.0) except Exception: # noqa return dict(m_nl_ok=False) rms = float(np.sqrt(g["chi2"] / max(2 * g["n"] - 7, 1))) ok = bool(np.isfinite(g["sigma_a1"]) and np.isfinite(g["a1"]) and g["iterations"] < n_iter - 1 and rms <= max(1.5 * m["m_rms7"], 1.5) and g["n"] >= 5) return dict(m_nl_A1=float(g["a1"]), m_nl_sig=float(g["sigma_a1"]), m_nl_rms7=rms, m_nl_n=int(g["n"]), m_nl_iter=int(g["iterations"]), m_nl_ok=ok) def gauss_newton(state0, t_ref, obs, a1=0.0, solve_a1=False, n_iter=8, clip=None): """Iterated weighted least-squares orbit fit (6 state components, optionally A1). Returns dict(state, a1, sigma_a1, chi2, n, iterations).""" x = np.array(state0, float) npar = 7 if solve_a1 else 6 keep = np.ones(len(obs["sig"]), bool) chi_prev = None sig_a1 = np.nan for it in range(n_iter): out = of.numeric_fisher(x, t_ref, obs, None, ng=(a1, 0.0, 0.0)) sr = obs["sig"] * of.ARCSEC ra_, rd_ = out["res_ra"] / obs["sig"], out["res_dec"] / obs["sig"] if clip is not None and it >= 2: keep = (np.abs(ra_) < clip) & (np.abs(rd_) < clip) k2 = np.concatenate([keep, keep]) A = np.vstack([out["Ha"][:, :npar] / sr[:, None], out["Hd"][:, :npar] / sr[:, None]])[k2] b = np.concatenate([ra_, rd_])[k2] chi = float(b @ b) dx, cov = of._solve_normal(A.T @ A, A.T @ b) if dx is None: break if solve_a1: sig_a1 = float(np.sqrt(cov[6])) if chi_prev is not None and abs(chi_prev - chi) < 1e-6 * max(chi, 1e-12): break chi_prev = chi x = x + dx[:6] if solve_a1: a1 = a1 + dx[6] return dict(state=x, a1=a1, sigma_a1=sig_a1, chi2=chi, n=int(keep.sum()), iterations=it) def oumuamua(store=None): """1I/'Oumuamua with JPL's solution as reference. (i) Forecast with two weightings: the generic scheme (HST = satellite class, 1.5 arcsec) and HST astrometry at 0.05 arcsec (the HST positions dominate the published solutions). (ii) [revision] Positive control of the measurement on real data: starting from JPL's state with *no* non-gravitational term, a 6-parameter + A1 orbit is fitted to the real astrometry (HST at 0.05 arcsec, its position taken from the MPC satellite records) and compared with Micheli et al. (2018).""" ng = pd.read_csv(os.path.join(DATA, "sbdb_nongrav_all.csv"), dtype={"pdes": str}, low_memory=False) o = ng[ng.pdes == "2017 U1"].iloc[0] if store is None: store = ObsStore() jd = None res = dict(pdes="1I", A1_jpl=o.A1, A1_sigma_jpl=o.A1_sigma, A1_micheli2018=4.92e-6 / 20.0408, A1_sigma_micheli2018=0.16e-6 / 20.0408) b = [0, 1, 2, 3, 4, 5] for tag, hs in (("", None), ("_hst0p05", 0.05)): obs = store.get("1I", api=True, hst_sigma=hs) if jd is None: st, jd = hz_state(load_hz(), "A/2017 U1", float(np.median(obs["mjd_utc"])) + 2400000.5) if st is None: jd = round(float(np.median(obs["mjd_utc"])) + 2400000.5, 3) dl.fetch_horizons(pd.DataFrame({"pdes": ["A/2017 U1"], "jd_ref": [jd]})) st, jd = hz_state(load_hz(), "A/2017 U1", jd) res.update(n_obs=len(obs["mjd_utc"]), n_hst=int((obs["stn"] == "250").sum()), arc_obs=obs["mjd_utc"][-1] - obs["mjd_utc"][0]) out = of.numeric_fisher(st, jd - 2400000.5, obs, None, ng=[o.A1, o.A2, o.A3]) res["sig_A1" + tag] = of.marginal_sigma(out["N"], b + [6], 6) res["sig_A1_jplset" + tag] = of.marginal_sigma(out["N"], b + [6, 7, 8], 6) rr = np.hypot(out["res_ra"], out["res_dec"]) hst = obs["stn"] == "250" res["res_med_ground"] = float(np.median(rr[~hst])) res["res_med_hst"] = float(np.median(rr[hst])) res["dist_min"] = float(out["dist"].min()) # positive control: gravity-only start, fit state + A1 to the real data g0 = gauss_newton(st, jd - 2400000.5, obs, a1=0.0, solve_a1=False, clip=3.0) g1 = gauss_newton(st, jd - 2400000.5, obs, a1=0.0, solve_a1=True, clip=3.0) res.update(A1_measured=g1["a1"], A1_measured_sigma=g1["sigma_a1"], chi2_gravity_only=g0["chi2"], chi2_with_A1=g1["chi2"], n_fit_gravity_only=g0["n"], n_fit_with_A1=g1["n"]) # the same control without HST (ground-based astrometry only, generic weights) obs_g = store.get("1I", api=True) mg = obs_g["stn"] != "250" obs_g = {k: (v[mg] if isinstance(v, np.ndarray) and len(v) == len(mg) else v) for k, v in obs_g.items()} g2 = gauss_newton(st, jd - 2400000.5, obs_g, a1=0.0, solve_a1=True, clip=3.0) res.update(A1_measured_ground=g2["a1"], A1_measured_ground_sigma=g2["sigma_a1"]) return res if __name__ == "__main__": if "--only=oumuamua" in sys.argv: r_ = oumuamua() pd.DataFrame([r_]).to_csv(os.path.join(DATA, "validation_oumuamua.csv"), index=False) print(r_) sys.exit(0) lim = [int(a.split("=")[1]) for a in sys.argv if a.startswith("--limit=")] main(limit=lim[0] if lim else None, fetch="--no-fetch" not in sys.argv)