"""Where on the sky were near-Earth asteroids when they were first observed? Added in the revision (referee point 1): the survey simulator used to ASSUME that survey coverage starts at exactly 60 deg solar elongation. This step measures, from real data, the solar elongation and apparent magnitude of near-Earth asteroids on the night of their first observation, so that s02 can derive the coverage-versus-elongation function instead. Eligible objects: catalogued near-Earth asteroids with a provisional designation from 1998 onward whose first observation lies within 4 days of the half-month encoded in that designation (so that the first observation is the discovery apparition, not a precovery). Geometry on the first-observation date (08:00 UT): sample A orbits whose osculating epoch is within 60 days of the first observation (JPL gives such epochs to short-arc objects): two-body propagation from the JPL elements, Earth from the astropy built-in ephemeris. Every sample-A object with a two-body elongation below 90 deg, and 60 random ones, are re-computed with JPL Horizons (the validation is printed). sample B orbits given at the standard epoch (2026): a two-body propagation over years is NOT accurate enough, so it is used only to stratify. B1: two-body elongation < 105 deg -> all queried with JPL Horizons B2: two-body elongation >= 105 deg -> N_B2_SAMPLE random objects queried with Horizons and weighted by (size of stratum) / (number queried) Horizons results are cached in RAW/horizons_firstobs.csv (about 4,000 single-date queries, 45 minutes; re-queried only with --refresh). Output: data/discovery_geometry.csv (one row per object used, with its weight) """ import json import os import sys import time import urllib.parse import numpy as np import pandas as pd from common import DATA, RAW, SEED, K_GAUSS, fetch_text, helio_xyz, observe, phase_HG from s03_fetch_earthlike_ephem import designation_window, prov_from_fullname HOR = "https://ssd.jpl.nasa.gov/api/horizons.api" N_B2_SAMPLE = 1000 N_A_VALID = 60 SCREEN_DEG = 105.0 OBLIQ = np.radians(23.4392911) def earth_helio_ecl(jd): """Heliocentric ecliptic position of the Earth (au) from the astropy built-in ephemeris.""" import astropy.units as u from astropy.coordinates import get_body_barycentric, solar_system_ephemeris from astropy.time import Time t = Time(jd, format="jd", scale="tdb") with solar_system_ephemeris.set("builtin"): e = get_body_barycentric("earth", t) s = get_body_barycentric("sun", t) x = (e.x - s.x).to(u.au).value y = (e.y - s.y).to(u.au).value z = (e.z - s.z).to(u.au).value return x, y * np.cos(OBLIQ) + z * np.sin(OBLIQ), -y * np.sin(OBLIQ) + z * np.cos(OBLIQ) def two_body(df): jd = df.jd_first.values + 8.0 / 24.0 M = np.radians(df.ma.values) + K_GAUSS / df.a.values ** 1.5 * (jd - df.epoch.values) xyz = helio_xyz(df.a.values, df.e.values, np.radians(df.i.values), np.radians(df.om.values), np.radians(df.w.values), M) r, d, alpha, elong, unit = observe(xyz, earth_helio_ecl(jd)) sin_dec = unit[2] * np.cos(OBLIQ) + unit[1] * np.sin(OBLIQ) return pd.DataFrame(dict(r=r, delta=d, phase=np.degrees(alpha), elong=np.degrees(elong), eclat=np.degrees(np.arcsin(unit[2])), dec=np.degrees(np.arcsin(sin_dec))), index=df.index) def horizons_one(pdes, date): """r, delta, phase, elongation, ecliptic latitude, declination at `date` 08:00 UT (geocentric).""" cmd = f"'{pdes};'" if pdes.strip().isdigit() else f"'DES={pdes};'" d1 = (pd.Timestamp(date) + pd.Timedelta(days=1)).strftime("%Y-%m-%d") q = {"format": "text", "COMMAND": cmd, "OBJ_DATA": "NO", "MAKE_EPHEM": "YES", "EPHEM_TYPE": "OBSERVER", "CENTER": "'500@399'", "START_TIME": f"'{date} 08:00'", "STOP_TIME": f"'{d1} 08:00'", "STEP_SIZE": "'1 d'", "QUANTITIES": "'1,3,19,20,23,24,31'", "CSV_FORMAT": "YES", "ANG_FORMAT": "DEG"} txt = fetch_text(HOR + "?" + urllib.parse.urlencode(q), timeout=90, retries=3) if "$$SOE" not in txt: return None f = [x.strip() for x in txt.split("$$SOE")[1].split("$$EOE")[0].strip().split("\n")[0].split(",")] # 0 date, 1,2 flags, 3 RA, 4 DEC, 5 dRA*cosD, 6 dDEC, 7 r, 8 rdot, 9 delta, 10 deldot, 11 S-O-T, 12 /r, 13 S-T-O, 14 EcLon, 15 EcLat return dict(r=float(f[7]), delta=float(f[9]), phase=float(f[13]), elong=float(f[11]), eclat=float(f[15]), dec=float(f[4])) def horizons_many(todo, cache_file): """todo: DataFrame with pdes, first_obs. Appends to the cache as it goes (resumable).""" have = pd.read_csv(cache_file, dtype={"pdes": str}) if os.path.exists(cache_file) else pd.DataFrame(columns=["pdes"]) done = set(have.pdes) todo = todo[~todo.pdes.isin(done)] if len(todo): print(f"[discgeom] querying JPL Horizons for {len(todo)} objects (about {len(todo) * 0.7 / 60:.0f} min)", flush=True) new_file = not os.path.exists(cache_file) t0 = time.time() with open(cache_file, "a") as f: if new_file: f.write("pdes,date,r,delta,phase,elong,eclat,dec\n") for n, (_, r) in enumerate(todo.iterrows()): try: h = horizons_one(r.pdes, r.first_obs) except Exception as ex: # noqa print("[discgeom] failed", r.pdes, ex) h = None if h is None: f.write(f"{r.pdes},{r.first_obs},,,,,,\n") else: f.write(f"{r.pdes},{r.first_obs},{h['r']},{h['delta']},{h['phase']},{h['elong']},{h['eclat']},{h['dec']}\n") f.flush() if (n + 1) % 250 == 0: print(f"[discgeom] {n + 1}/{len(todo)} {time.time() - t0:.0f}s", flush=True) time.sleep(0.05) return pd.read_csv(cache_file, dtype={"pdes": str}).drop_duplicates("pdes", keep="last").set_index("pdes") def load_eligible(): d = json.load(open(os.path.join(RAW, "sbdb_neo_epoch.json"))) df = pd.DataFrame(d["data"], columns=d["fields"]) for c in ["H", "epoch", "a", "e", "i", "om", "w", "ma", "data_arc"]: df[c] = pd.to_numeric(df[c], errors="coerce") df["pdes"] = df.pdes.astype(str).str.strip() fo = pd.to_datetime(df.first_obs, errors="coerce") df["jd_first"] = [t.to_julian_date() if pd.notna(t) else np.nan for t in fo] df["t_first"] = fo.dt.year + (fo.dt.dayofyear - 1) / 365.25 wins = [designation_window(prov_from_fullname(n) or p) for n, p in zip(df.full_name, df.pdes)] df["win0"] = [w[0] if w else np.nan for w in wins] df["win1"] = [w[1] if w else np.nan for w in wins] n_all = len(df) n_1998 = int((df.win0 >= 1998.0).sum()) ok = (df.t_first >= df.win0 - 4 / 365.25) & (df.t_first <= df.win1 + 4 / 365.25) & (df.win0 >= 1998.0) & df.H.notna() & (df.e < 1) return df[ok].copy(), n_all, n_1998 def main(refresh=False): out_csv = os.path.join(DATA, "discovery_geometry.csv") cache = os.path.join(RAW, "horizons_firstobs.csv") if refresh and os.path.exists(cache): os.remove(cache) df, n_all, n_1998 = load_eligible() modal = df.epoch.mode()[0] print(f"[discgeom] {n_all} near-Earth asteroids; {n_1998} designated 1998 or later; {len(df)} with the first observation " f"inside the designation half-month (+-4 d)") rng = np.random.default_rng(SEED + 12) A = df[(df.epoch != modal) & ((df.epoch - df.jd_first).abs() <= 60)].copy() B = df[df.epoch == modal].copy() n_other = len(df) - len(A) - len(B) A = A.join(two_body(A)) B["elong_screen"] = two_body(B).elong B1 = B[B.elong_screen < SCREEN_DEG] B2pool = B[B.elong_screen >= SCREEN_DEG] B2 = B2pool.loc[rng.choice(B2pool.index.values, size=min(N_B2_SAMPLE, len(B2pool)), replace=False)] A_low = A[(A.elong < 90) | (A.delta < 0.003)] A_val = A.loc[rng.choice(A.index.values, size=N_A_VALID, replace=False)] print(f"[discgeom] sample A (epoch near first observation, two-body): {len(A)}; sample B (standard epoch): {len(B)} " f"-> B1 (screen < {SCREEN_DEG:.0f} deg, all queried): {len(B1)}, B2 pool: {len(B2pool)}, queried: {len(B2)}; " f"neither: {n_other}") todo = pd.concat([A_low, A_val, B1, B2])[["pdes", "first_obs"]].drop_duplicates("pdes") hz = horizons_many(todo, cache) # ---- validation of the two-body geometry of sample A ---- v = A_val.join(hz[["elong"]].rename(columns={"elong": "elong_hz"}), on="pdes").dropna(subset=["elong_hz"]) dv = (v.elong - v.elong_hz).abs() print(f"[discgeom] two-body versus Horizons for {len(v)} random sample-A objects: median |d elong| {dv.median():.2f} deg, " f"90th percentile {dv.quantile(0.9):.2f}, max {dv.max():.1f}") # ---- validation of the screen of sample B: how many B2 objects are really below 90 deg? ---- cols = ["r", "delta", "phase", "elong", "eclat", "dec"] parts = [] a = A.copy() a["source"] = "two-body" hit = a.pdes.isin(hz.index) & a.pdes.map(lambda p: p in hz.index and np.isfinite(hz.loc[p, "elong"])) for c in cols: a.loc[hit, c] = a.loc[hit, "pdes"].map(hz[c]) a.loc[hit, "source"] = "horizons" a["sample"], a["weight"] = "A", 1.0 parts.append(a) for name, sub, wgt_pool in [("B1", B1, None), ("B2", B2, len(B2pool))]: s = sub.join(hz[cols], on="pdes").dropna(subset=["elong"]) s["source"], s["sample"] = "horizons", name s["weight"] = 1.0 if wgt_pool is None else wgt_pool / len(s) if wgt_pool is None and len(s) < len(sub): s["weight"] = len(sub) / len(s) parts.append(s) out = pd.concat(parts) out["V"] = out.H + 5 * np.log10(out.r * out.delta) + phase_HG(np.radians(out.phase.values)) out["year"] = out.t_first b2 = out[out["sample"] == "B2"] print(f"[discgeom] B2 (screened as >= {SCREEN_DEG:.0f} deg): {len(b2)} queried, of which truly < 90 deg: {(b2.elong < 90).sum()}, " f"< 80 deg: {(b2.elong < 80).sum()}, < 70 deg: {(b2.elong < 70).sum()} (each stands for {b2.weight.iloc[0]:.1f} objects)") out[["pdes", "H", "year", "first_obs", "data_arc", "sample", "source", "weight", "elong", "V", "r", "delta", "phase", "eclat", "dec"]].to_csv( out_csv, index=False, float_format="%.5g") bins = [0, 45, 60, 70, 80, 90, 105, 120, 135, 150, 165, 180.01] out["bin"] = pd.cut(out.elong, bins) tab = out.groupby("bin", observed=False).agg(n_objects=("weight", "size"), n_weighted=("weight", "sum"), V_median=("V", "median")) tab["fraction"] = tab.n_weighted / tab.n_weighted.sum() pd.set_option("display.width", 200) print("[discgeom] solar elongation at first observation (all eligible near-Earth asteroids, weighted):") print(tab.round(4).to_string()) low = out[out.elong < 75].sort_values("elong") print(f"[discgeom] {len(low)} objects first observed below 75 deg elongation (weighted {low.weight.sum():.0f}); the 25 lowest:") print(low[["pdes", "H", "first_obs", "sample", "source", "elong", "delta", "V"]].head(25).round(2).to_string(index=False)) if __name__ == "__main__": main("--refresh" in sys.argv)