diff --git a/src/lib/orbit_fit/orbit_fit.cpp b/src/lib/orbit_fit/orbit_fit.cpp index e3df663a..bdf7c67b 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 new file mode 100644 index 00000000..11a4fde1 --- /dev/null +++ b/tools/gen_radar_fixture.py @@ -0,0 +1,164 @@ +"""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 = 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. + +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 + # 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 + + +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()