From 215fadb7445f41a7f5a05c5d53f6225ad99fbbd3 Mon Sep 17 00:00:00 2001 From: matthewholman Date: Thu, 3 Sep 2026 17:19:07 -0400 Subject: [PATCH 1/3] Move the radar fixture generator into the repository tests/data/radar_synthetic.json was produced by a script that lived outside the repository, so the fixture could not be regenerated by anyone else and quietly went stale whenever the radar model changed. Add it under tools/ as it stood. Co-Authored-By: Claude Opus 5 --- tools/gen_radar_fixture.py | 149 +++++++++++++++++++++++++++++++++++++ 1 file changed, 149 insertions(+) create mode 100644 tools/gen_radar_fixture.py diff --git a/tools/gen_radar_fixture.py b/tools/gen_radar_fixture.py new file mode 100644 index 00000000..2c6ab03b --- /dev/null +++ b/tools/gen_radar_fixture.py @@ -0,0 +1,149 @@ +"""Generate tests/data/radar_synthetic.json for the radar (delay/Doppler) fit. + +Reuses the streak fixture's true orbit, epoch, and per-observation observer +states, and computes the round-trip radar observables with the SAME light-time +convention as the C++ model (predict.cpp::integrate_light_time iterates to the +retarded time t_obs - rho/c): + + delay = 2 * rho / c + Shapiro round-trip light time, days + doppler = 2 * (rho_hat . v_rel) round-trip range rate, au/day + +rho/v are evaluated at the retarded emission time; v_rel = v_ast - v_obs. + +The Shapiro (relativistic) term is included, matching orbit_fit.cpp. It has to +be: this fixture is truth for a fit that models it, so omitting it would make the +test measure the difference between two models rather than the fitter's ability +to recover an orbit. It is 1-2 us here against a stated 1 us uncertainty, on +seven observations constraining six parameters, so the difference is not subtle. +""" +import json +from pathlib import Path + +import assist +import numpy as np +import pooch +import rebound + +AU_M = 149597870700.0 +C_AU_DAY = 2.99792458e8 * 86400.0 / AU_M # matches predict.cpp SPEED_OF_LIGHT + +CACHE = pooch.os_cache("layup") +STREAK = Path("tests/data/streak_synthetic.json") +OUT = Path("tests/data/radar_synthetic.json") + +ephem = assist.Ephem( + str(CACHE / "linux_p1550p2650.440"), + str(CACHE / "sb441-n16.bsp"), +) +JD_REF = ephem.jd_ref + + +def state_at(true_state, epoch, t_target_jd): + """Asteroid barycentric (r, v) at t_target_jd, integrating true_state@epoch.""" + sim = rebound.Simulation() + sim.t = epoch - JD_REF + sim.add( + x=true_state[0], y=true_state[1], z=true_state[2], + vx=true_state[3], vy=true_state[4], vz=true_state[5], + ) + ax = assist.Extras(sim, ephem) + sim.integrate(t_target_jd - JD_REF) + p = sim.particles[0] + r = np.array([p.x, p.y, p.z]) + v = np.array([p.vx, p.vy, p.vz]) + ax.detach(sim) + return r, v + + +def radar_observables(true_state, epoch, obs_epoch, r_obs, v_obs): + """Two-leg round-trip delay (days) and Doppler (au/day), observer accel = 0. + + Mirrors the C++ orbit_fit.cpp radar model with observer_acceleration left at + its default of zero (this fixture is fed to run_from_vector directly without + acceleration). Down leg: station at the receive epoch. Up leg: station linearly + extrapolated to the transmit time t - tau by v_obs (no accel term). + """ + r_obs = np.asarray(r_obs) + v_obs = np.asarray(v_obs) + # Down leg: retarded bounce time using the station at receive. + tau_d = 0.0 + for _ in range(4): + r_ast, v_ast = state_at(true_state, epoch, obs_epoch - tau_d) + rho_d_vec = r_ast - r_obs + rho_d = np.linalg.norm(rho_d_vec) + tau_d = rho_d / C_AU_DAY + rho_hat_d = rho_d_vec / rho_d + # Up leg: station at the transmit time t - (tau_d + tau_u), linear in v_obs. + tau_u = tau_d + for _ in range(5): + r_tx = r_obs - v_obs * (tau_d + tau_u) + rho_u_vec = r_ast - r_tx + rho_u = np.linalg.norm(rho_u_vec) + tau_u = rho_u / C_AU_DAY + rho_hat_u = rho_u_vec / rho_u + + # Shapiro delay on both legs -- the same formula and constant as + # orbit_fit.cpp::compute_radar_residuals. The Sun moves ~1e-5 au over a round + # trip, so a single evaluation at the receive epoch is ample. + # At the C++ residual, integrate_light_time has left the simulation at the + # emission (bounce) time, so the Sun is evaluated there and not at receive. + sun = ephem.get_particle(0, (obs_epoch - tau_d) - JD_REF) # ASSIST_BODY_SUN + S = np.array([sun.x, sun.y, sun.z]) + GM_SUN = 2.9591220828559115e-4 # au^3/day^2 + k = 2.0 * GM_SUN / C_AU_DAY**3 + r_b = np.linalg.norm(r_ast - S) + r_r = np.linalg.norm(r_obs - S) + r_t = np.linalg.norm(r_tx - S) + shapiro = k * (np.log((r_t + r_b + rho_u) / (r_t + r_b - rho_u)) + + np.log((r_b + r_r + rho_d) / (r_b + r_r - rho_d))) + + delay = tau_d + tau_u + shapiro + doppler = float(rho_hat_d @ (v_ast - v_obs)) + float(rho_hat_u @ (v_ast - v_obs)) + return delay, doppler + + +def main(): + d = json.loads(STREAK.read_text()) + true_state = d["true_state"] + epoch = d["epoch"] + + # 1-sigma uncertainties: realistic radar quality. + # JPL delay ~ a few us round-trip; Doppler ~ sub-Hz. Convert to internal units. + delay_unc_days = 1.0e-6 / 86400.0 # 1 us in days + doppler_unc_audy = 1.0e-9 # ~ mm/s-level range-rate, au/day + + out = { + "description": ( + "Synthetic noise-free radar (delay/Doppler) arc from the same MBA " + "orbit as streak_synthetic.json. delay = round-trip light time (days), " + "doppler = round-trip range rate (au/day); generated with ASSIST using " + "the C++ light-time convention including the Shapiro delay (gen_radar_fixture.py)." + ), + "jd_ref": JD_REF, + "epoch": epoch, + "true_state": true_state, + "delay_unc_days": delay_unc_days, + "doppler_unc_audy": doppler_unc_audy, + "observations": [], + } + for o in d["observations"]: + delay, doppler = radar_observables( + true_state, epoch, o["epoch"], o["observer_position"], o["observer_velocity"] + ) + out["observations"].append( + { + "epoch": o["epoch"], + "observer_position": o["observer_position"], + "observer_velocity": o["observer_velocity"], + "delay": delay, + "doppler": doppler, + } + ) + print(f" t={o['epoch']:.1f} delay={delay:.10e} d doppler={doppler:.10e} au/d") + + OUT.write_text(json.dumps(out, indent=2)) + print("wrote", OUT) + + +if __name__ == "__main__": + main() From aab669ab9f8c1956280d4aed6e8cca165c995a1e Mon Sep 17 00:00:00 2001 From: matthewholman Date: Thu, 3 Sep 2026 17:19:07 -0400 Subject: [PATCH 2/3] Carry the retardation denominators in the radar Doppler model The Doppler observable is c d(tau)/d(t_receive). Differentiating the two implicit light-time equations leaves a denominator on each leg: dt_bounce /dt_receive = (c + rho_down.v_receiver)/(c + rho_down.v_asteroid) dt_transmit/dt_receive = dt_bounce/dt_receive * (c - rho_up.v_asteroid)/(c - rho_up.v_transmitter) range rate = c (1 - dt_transmit/dt_receive) To first order that is the instantaneous sum of the two one-way range rates, which is what the model computed. The omitted term is a fractional error of order rho.v/c. That is negligible for most targets and not at all negligible for a fast one: on (6489) Golevka's 1995 apparition the line-of-sight rate reaches 6 km/s, making it 5 to 7 Hz against stated uncertainties of 0.09 to 0.40 Hz -- a coherent, one-signed bias on every Doppler row. The same denominator was already being applied to the range partials, as ltdenom, so the Jacobian and the residual disagreed. Both fixture generators encoded the old model and are updated with it, and tests/data/radar_synthetic.json is regenerated: only the Doppler values move, by 6e-7 to 1.7e-5 fractionally, and the delays are untouched. Co-Authored-By: Claude Opus 5 --- src/lib/orbit_fit/orbit_fit.cpp | 30 +++++++++++++++++++++++----- tests/data/radar_synthetic.json | 14 ++++++------- tests/layup/test_radar_end_to_end.py | 10 ++++++++-- tools/gen_radar_fixture.py | 13 ++++++++++-- 4 files changed, 51 insertions(+), 16 deletions(-) diff --git a/src/lib/orbit_fit/orbit_fit.cpp b/src/lib/orbit_fit/orbit_fit.cpp index 7a6ca171..216986b4 100644 --- a/src/lib/orbit_fit/orbit_fit.cpp +++ b/src/lib/orbit_fit/orbit_fit.cpp @@ -252,11 +252,31 @@ namespace orbit_fit if (up > 0.0 && down > 0.0) model_delay += two_gm_over_c3 * (log(up) + log(down)); } - // Round-trip range rate: down leg uses the station velocity at receive, - // up leg at transmit. - double model_doppler = - (g.rho_x * (vax - vox) + g.rho_y * (vay - voy) + g.rho_z * (vaz - voz)) + - (rhu_x * (vax - vtx_x) + rhu_y * (vay - vtx_y) + rhu_z * (vaz - vtx_z)); + // Round-trip range rate. The observable is c * d(tau)/d(t_receive), and + // differentiating the two implicit light-time equations leaves a + // retardation denominator on each leg: + // + // dt_bounce / dt_receive = (c + rho_down . v_receiver) + // / (c + rho_down . v_asteroid) + // dt_transmit / dt_receive = dt_bounce/dt_receive + // * (c - rho_up . v_asteroid) + // / (c - rho_up . v_transmitter) + // range rate = c * (1 - dt_transmit/dt_receive) + // + // To first order this is the instantaneous sum of the two one-way range + // rates, which is what this model used to compute. The omitted term is a + // fractional error of order rho.v/c: small, but on (6489) Golevka's 1995 + // apparition the line-of-sight rate reaches 6 km/s, making it 5 to 7 Hz + // against stated uncertainties of 0.09 to 0.40 Hz. The same denominator + // is already applied to the range partials below, as ltdenom. + double rd_v_ast = g.rho_x * vax + g.rho_y * vay + g.rho_z * vaz; + double rd_v_rcv = g.rho_x * vox + g.rho_y * voy + g.rho_z * voz; + double ru_v_ast = rhu_x * vax + rhu_y * vay + rhu_z * vaz; + double ru_v_tx = rhu_x * vtx_x + rhu_y * vtx_y + rhu_z * vtx_z; + double dt_bounce = (SPEED_OF_LIGHT + rd_v_rcv) / (SPEED_OF_LIGHT + rd_v_ast); + double dt_transmit = + dt_bounce * (SPEED_OF_LIGHT - ru_v_ast) / (SPEED_OF_LIGHT - ru_v_tx); + double model_doppler = SPEED_OF_LIGHT * (1.0 - dt_transmit); resid.delay_resid = rd.delay - model_delay; resid.doppler_resid = rd.doppler - model_doppler; diff --git a/tests/data/radar_synthetic.json b/tests/data/radar_synthetic.json index 322b1778..705d9e8e 100644 --- a/tests/data/radar_synthetic.json +++ b/tests/data/radar_synthetic.json @@ -26,7 +26,7 @@ 0.0029573124252989546 ], "delay": 0.018636429943384415, - "doppler": -0.005527709239628141 + "doppler": -0.00552779753855833 }, { "epoch": 2459548.3333333335, @@ -41,7 +41,7 @@ 0.002590684205279173 ], "delay": 0.0185479482944552, - "doppler": -0.003644924931099627 + "doppler": -0.0036449633375001795 }, { "epoch": 2459551.6666666665, @@ -56,7 +56,7 @@ 0.002212605883097436 ], "delay": 0.01849614339489692, - "doppler": -0.0017199773152476484 + "doppler": -0.0017199858776861766 }, { "epoch": 2459555.0, @@ -71,7 +71,7 @@ 0.0018246951514722526 ], "delay": 0.018481646880337224, - "doppler": 0.00022673915974917678 + "doppler": 0.00022673901388991867 }, { "epoch": 2459558.3333333335, @@ -86,7 +86,7 @@ 0.0014302740766842501 ], "delay": 0.018504628057560944, - "doppler": 0.0021694651973099667 + "doppler": 0.002169451630493268 }, { "epoch": 2459561.6666666665, @@ -101,7 +101,7 @@ 0.0010320202556887116 ], "delay": 0.01856481116919987, - "doppler": 0.004089017688450958 + "doppler": 0.004088969451188339 }, { "epoch": 2459565.0, @@ -116,7 +116,7 @@ 0.0006313244406405992 ], "delay": 0.018661620841759898, - "doppler": 0.005972590210808326 + "doppler": 0.005972487266112379 } ] } \ No newline at end of file diff --git a/tests/layup/test_radar_end_to_end.py b/tests/layup/test_radar_end_to_end.py index 3fd4e3d2..e83f77c7 100644 --- a/tests/layup/test_radar_end_to_end.py +++ b/tests/layup/test_radar_end_to_end.py @@ -17,7 +17,7 @@ were validated separately. The orbit is the same ~2.6 AU main-belt object near opposition as the streak fixture (``tests/data/streak_synthetic.json``); the truth observables are an independent ASSIST propagation at the C++ light-time convention -(``delay = 2 rho/c`` plus the Shapiro delay; ``doppler = 2 rho_hat . v_rel``). +(``delay = 2 rho/c`` plus the Shapiro delay; ``doppler = c d(tau)/d(t_receive)``). Radar over a short single-station arc weakly constrains the plane-of-sky position, so -- as in real radar astrometry -- the fit *refines a prior orbit*: @@ -135,7 +135,13 @@ def state_at(t_jd): ) delay = tau_d + tau_u + shapiro - doppler = float(rho_hat_d @ (v_ast - v_obs)) + float(rho_hat_u @ (v_ast - v_tx)) + # Round-trip range rate, carrying the same retardation denominators the fitter + # applies: the observable is c * d(tau)/d(t_receive), not the instantaneous sum + # of the two one-way range rates. + c = SPEED_OF_LIGHT + dt_bounce = (c + float(rho_hat_d @ v_obs)) / (c + float(rho_hat_d @ v_ast)) + dt_transmit = dt_bounce * (c - float(rho_hat_u @ v_ast)) / (c - float(rho_hat_u @ v_tx)) + doppler = c * (1.0 - dt_transmit) return delay, doppler diff --git a/tools/gen_radar_fixture.py b/tools/gen_radar_fixture.py index 2c6ab03b..e62963a5 100644 --- a/tools/gen_radar_fixture.py +++ b/tools/gen_radar_fixture.py @@ -1,12 +1,14 @@ """Generate tests/data/radar_synthetic.json for the radar (delay/Doppler) fit. +Run from the repository root: python tools/gen_radar_fixture.py + Reuses the streak fixture's true orbit, epoch, and per-observation observer states, and computes the round-trip radar observables with the SAME light-time convention as the C++ model (predict.cpp::integrate_light_time iterates to the retarded time t_obs - rho/c): delay = 2 * rho / c + Shapiro round-trip light time, days - doppler = 2 * (rho_hat . v_rel) round-trip range rate, au/day + doppler = c * d(tau)/d(t_receive) round-trip range rate, au/day rho/v are evaluated at the retarded emission time; v_rel = v_ast - v_obs. @@ -98,7 +100,14 @@ def radar_observables(true_state, epoch, obs_epoch, r_obs, v_obs): + np.log((r_b + r_r + rho_d) / (r_b + r_r - rho_d))) delay = tau_d + tau_u + shapiro - doppler = float(rho_hat_d @ (v_ast - v_obs)) + float(rho_hat_u @ (v_ast - v_obs)) + # Round-trip range rate with the same retardation denominators the fitter + # applies: the observable is c * d(tau)/d(t_receive), not the instantaneous + # sum of the two one-way range rates. Monostatic here, so the transmitting + # station is the receiving one, Taylor-extrapolated back to transmit. + v_tx = v_obs + dt_bounce = (C_AU_DAY + float(rho_hat_d @ v_obs)) / (C_AU_DAY + float(rho_hat_d @ v_ast)) + dt_transmit = dt_bounce * (C_AU_DAY - float(rho_hat_u @ v_ast)) / (C_AU_DAY - float(rho_hat_u @ v_tx)) + doppler = C_AU_DAY * (1.0 - dt_transmit) return delay, doppler From 9bcd0c517dca0482911fe2b564d48319549efa5c Mon Sep 17 00:00:00 2001 From: matthewholman Date: Thu, 3 Sep 2026 17:24:43 -0400 Subject: [PATCH 3/3] Apply black to the imported fixture generator It was written outside the repository and so never saw the formatter. Co-Authored-By: Claude Opus 5 --- tools/gen_radar_fixture.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/tools/gen_radar_fixture.py b/tools/gen_radar_fixture.py index e62963a5..11a4fde1 100644 --- a/tools/gen_radar_fixture.py +++ b/tools/gen_radar_fixture.py @@ -18,6 +18,7 @@ to recover an orbit. It is 1-2 us here against a stated 1 us uncertainty, on seven observations constraining six parameters, so the difference is not subtle. """ + import json from pathlib import Path @@ -45,8 +46,12 @@ def state_at(true_state, epoch, t_target_jd): sim = rebound.Simulation() sim.t = epoch - JD_REF sim.add( - x=true_state[0], y=true_state[1], z=true_state[2], - vx=true_state[3], vy=true_state[4], vz=true_state[5], + x=true_state[0], + y=true_state[1], + z=true_state[2], + vx=true_state[3], + vy=true_state[4], + vz=true_state[5], ) ax = assist.Extras(sim, ephem) sim.integrate(t_target_jd - JD_REF) @@ -89,15 +94,16 @@ def radar_observables(true_state, epoch, obs_epoch, r_obs, v_obs): # trip, so a single evaluation at the receive epoch is ample. # At the C++ residual, integrate_light_time has left the simulation at the # emission (bounce) time, so the Sun is evaluated there and not at receive. - sun = ephem.get_particle(0, (obs_epoch - tau_d) - JD_REF) # ASSIST_BODY_SUN + sun = ephem.get_particle(0, (obs_epoch - tau_d) - JD_REF) # ASSIST_BODY_SUN S = np.array([sun.x, sun.y, sun.z]) - GM_SUN = 2.9591220828559115e-4 # au^3/day^2 + GM_SUN = 2.9591220828559115e-4 # au^3/day^2 k = 2.0 * GM_SUN / C_AU_DAY**3 r_b = np.linalg.norm(r_ast - S) r_r = np.linalg.norm(r_obs - S) r_t = np.linalg.norm(r_tx - S) - shapiro = k * (np.log((r_t + r_b + rho_u) / (r_t + r_b - rho_u)) - + np.log((r_b + r_r + rho_d) / (r_b + r_r - rho_d))) + shapiro = k * ( + np.log((r_t + r_b + rho_u) / (r_t + r_b - rho_u)) + np.log((r_b + r_r + rho_d) / (r_b + r_r - rho_d)) + ) delay = tau_d + tau_u + shapiro # Round-trip range rate with the same retardation denominators the fitter @@ -118,8 +124,8 @@ def main(): # 1-sigma uncertainties: realistic radar quality. # JPL delay ~ a few us round-trip; Doppler ~ sub-Hz. Convert to internal units. - delay_unc_days = 1.0e-6 / 86400.0 # 1 us in days - doppler_unc_audy = 1.0e-9 # ~ mm/s-level range-rate, au/day + delay_unc_days = 1.0e-6 / 86400.0 # 1 us in days + doppler_unc_audy = 1.0e-9 # ~ mm/s-level range-rate, au/day out = { "description": (