From b598993a1addd4f8acf91accf9a02a14ebe15f59 Mon Sep 17 00:00:00 2001 From: matthewholman Date: Wed, 2 Sep 2026 18:17:00 -0400 Subject: [PATCH 1/4] Model the Shapiro delay in radar residuals A radar signal passes through the Sun's gravitational potential on each leg, and the resulting delay was not modelled. For a leg between A and B, dt = (2 GM / c^3) ln[(r_A + r_B + rho) / (r_A + r_B - rho)] with r_A and r_B heliocentric distances. At the geometries these observations are made it is one to two microseconds, which sounds negligible until you notice that radar delay uncertainties are 0.3 to 2 us: leaving it out biases every radar observation by several sigma. Validated against real Apophis radar astrometry, on the epochs near the reference orbit where propagation error is negligible. The delay residual falls from 2.16 sigma to 0.312 sigma, and the Doppler residual is unchanged at 0.28 sigma. The comparison is against an independent Python implementation of the same formula, so it checks the physics rather than this code; the harness is ISSUE_146_shapiro_validation.py. That confirms the conjecture recorded when the two-leg model was written, which attributed its residual to exactly this term. The ephemeris is passed into compute_radar_residuals rather than read from r->extras: the simulation used for residuals has no assist_extras attached, so that pointer is NULL there. The partials are left alone. The Jacobian is already a single-leg approximation by design, and a microsecond-scale, slowly varying term does not change the step direction. --- src/lib/orbit_fit/orbit_fit.cpp | 55 ++++++++++++++++++++++++++++++--- 1 file changed, 50 insertions(+), 5 deletions(-) diff --git a/src/lib/orbit_fit/orbit_fit.cpp b/src/lib/orbit_fit/orbit_fit.cpp index 4b4ea9f1..33d51e97 100644 --- a/src/lib/orbit_fit/orbit_fit.cpp +++ b/src/lib/orbit_fit/orbit_fit.cpp @@ -24,6 +24,7 @@ extern "C" { #include "rebound.h" +#include "assist.h" } namespace orbit_fit @@ -96,6 +97,10 @@ namespace orbit_fit // uncertainties from arcseconds to radians and scales residuals for display. static constexpr double ARCSEC_PER_RAD = 206265.0; + // Sun's gravitational parameter, au^3/day^2. Used only for the Shapiro delay, + // whose coefficient 2*GM/c^3 works out to 1.1402e-10 d = 9.851 us. + static constexpr double GM_SUN_AU3_D2 = 2.9591220828559115e-4; + // Geometry shared by all three observable residual paths, computed once by // compute_single_residuals after the light-time integration: the unit @@ -123,7 +128,8 @@ namespace orbit_fit // extrapolating the station state (pos+vel) to the transmit time t_obs - tau // using the observer acceleration supplied from Python. Shapiro (relativistic) // delay (~2 us here) is the next refinement. - void compute_radar_residuals(struct reb_simulation *r, const Observation &this_det, + void compute_radar_residuals(struct assist_ephem *ephem, struct reb_simulation *r, + const Observation &this_det, int var, int npar, const ResidualGeometry &g, residuals &resid, partials &parts) { @@ -156,12 +162,15 @@ namespace orbit_fit double tau_d = g.rho / SPEED_OF_LIGHT; double tau_u = tau_d; double rhu_x = g.rho_x, rhu_y = g.rho_y, rhu_z = g.rho_z, rho_u = g.rho; + // Kept outside the loop: the converged transmit-time station position is + // needed again below for the up leg's Shapiro term. + double rtx_x = xe, rtx_y = ye, rtx_z = ze; for (int it = 0; it < 3; it++) { double tau = tau_d + tau_u; - double rtx_x = xe - vox * tau - 0.5 * aox * tau * tau; - double rtx_y = ye - voy * tau - 0.5 * aoy * tau * tau; - double rtx_z = ze - voz * tau - 0.5 * aoz * tau * tau; + rtx_x = xe - vox * tau - 0.5 * aox * tau * tau; + rtx_y = ye - voy * tau - 0.5 * aoy * tau * tau; + rtx_z = ze - voz * tau - 0.5 * aoz * tau * tau; rhu_x = rbx - rtx_x; rhu_y = rby - rtx_y; rhu_z = rbz - rtx_z; @@ -177,6 +186,42 @@ namespace orbit_fit double vtx_z = voz - aoz * tau; double model_delay = tau_d + tau_u; // round-trip light time (days) + + // Shapiro (relativistic) delay. The signal traverses the Sun's potential + // on each leg; for a leg between A and B, + // dt = (2 GM / c^3) ln[(r_A + r_B + rho) / (r_A + r_B - rho)] + // with r_A and r_B heliocentric distances and rho the leg length. At the + // geometries these observations are made it is a few microseconds -- + // small, but radar uncertainties here are 0.3 to 2 us, so leaving it out + // biases the delay by several sigma per observation and inflates the + // chi-square enough to reject an otherwise good fit. + // + // The Sun is evaluated once, at the receive epoch. It moves of order + // 1e-5 au during a round trip, which shifts the logarithm's argument far + // below the microsecond level. + // The ephemeris is passed in rather than taken from r->extras: the + // simulation used for residuals has no assist_extras attached, so that + // pointer is NULL here. + if (ephem != NULL) + { + struct reb_particle sun = assist_get_particle(ephem, ASSIST_BODY_SUN, r->t); + auto heliocentric_distance = [&](double x, double y, double z) { + double dx = x - sun.x, dy = y - sun.y, dz = z - sun.z; + return sqrt(dx * dx + dy * dy + dz * dz); + }; + double r_bounce = heliocentric_distance(rbx, rby, rbz); + double r_receive = heliocentric_distance(xe, ye, ze); + double r_transmit = heliocentric_distance(rtx_x, rtx_y, rtx_z); + const double two_gm_over_c3 = + 2.0 * GM_SUN_AU3_D2 / (SPEED_OF_LIGHT * SPEED_OF_LIGHT * SPEED_OF_LIGHT); + double up = (r_transmit + r_bounce + rho_u) / (r_transmit + r_bounce - rho_u); + double down = (r_bounce + r_receive + g.rho) / (r_bounce + r_receive - g.rho); + // Guard the logarithms: the denominators vanish only for a signal + // grazing the Sun, which is not an observable geometry, but a + // non-positive argument must never reach log(). + 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 = @@ -429,7 +474,7 @@ namespace orbit_fit // for non-radar, and a streak adds its two rate rows on top. if (std::holds_alternative(this_det.observation_type)) { - compute_radar_residuals(r, this_det, var, npar, g, resid, parts); + compute_radar_residuals(ephem, r, this_det, var, npar, g, resid, parts); return; } compute_optical_residuals(r, this_det, var, npar, g, resid, parts); From 6f2dbcbd8510d9bf770ff13ae2c17f0ad1c836fa Mon Sep 17 00:00:00 2001 From: matthewholman Date: Wed, 2 Sep 2026 18:29:45 -0400 Subject: [PATCH 2/4] Carry the Shapiro delay in the radar test fixtures The synthetic radar truth was generated with delay = 2 rho/c and no Shapiro term. Once the fitter models it, that truth is inconsistent with the fitter by 19.2 us per observation -- on seven observations with a stated 1 us uncertainty, which walked the chi-square to 1.0e-2. Both generators now include the term: tests/data/radar_synthetic.json is regenerated (ISSUE_146_gen_radar_fixture.py, kept with the other harnesses), and test_radar_end_to_end builds its observables inline, so that copy is fixed too. The Sun is taken at the emission time, where the C++ light-time solution leaves the simulation. The end-to-end tests then pass unchanged. test_radar_validation needed its chi-square bound moved, because that bound changed meaning: the fixture comes from an independent implementation of the same model, and the two agree to about 0.3 ns per observation rather than exactly. That is 3e-4 of the fixture's own stated uncertainty, and a hundred times finer than the best real radar timing in the JPL database. The old 1e-6 bound was achievable only while fixture and fitter shared a formula. The state tolerances are tightened rather than relaxed -- 10x on position and 100x on velocity -- since those carry the physics the test exists to check, and recovery is 2.7e-9 au and 1.5e-11 au/day. --- tests/data/radar_synthetic.json | 16 +++++++-------- tests/layup/test_radar_end_to_end.py | 20 +++++++++++++++++-- tests/layup/test_radar_validation.py | 30 +++++++++++++++++++++++----- 3 files changed, 51 insertions(+), 15 deletions(-) diff --git a/tests/data/radar_synthetic.json b/tests/data/radar_synthetic.json index 08b8ea5b..322b1778 100644 --- a/tests/data/radar_synthetic.json +++ b/tests/data/radar_synthetic.json @@ -1,5 +1,5 @@ { - "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 (gen_radar_fixture.py).", + "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": 2451545.0, "epoch": 2459555.0, "true_state": [ @@ -25,7 +25,7 @@ 0.006821651639396558, 0.0029573124252989546 ], - "delay": 0.018636429721406632, + "delay": 0.018636429943384415, "doppler": -0.005527709239628141 }, { @@ -40,7 +40,7 @@ 0.005974675082844652, 0.002590684205279173 ], - "delay": 0.018547948073657488, + "delay": 0.0185479482944552, "doppler": -0.003644924931099627 }, { @@ -55,7 +55,7 @@ 0.005101795432686577, 0.002212605883097436 ], - "delay": 0.018496143174772404, + "delay": 0.01849614339489692, "doppler": -0.0017199773152476484 }, { @@ -70,7 +70,7 @@ 0.004207405088519609, 0.0018246951514722526 ], - "delay": 0.018481646660371547, + "delay": 0.018481646880337224, "doppler": 0.00022673915974917678 }, { @@ -85,7 +85,7 @@ 0.0032987925947447197, 0.0014302740766842501 ], - "delay": 0.018504627837237886, + "delay": 0.018504628057560944, "doppler": 0.0021694651973099667 }, { @@ -100,7 +100,7 @@ 0.002381224227046953, 0.0010320202556887116 ], - "delay": 0.01856481094800657, + "delay": 0.01856481116919987, "doppler": 0.004089017688450958 }, { @@ -115,7 +115,7 @@ 0.0014573943662972056, 0.0006313244406405992 ], - "delay": 0.018661620619189488, + "delay": 0.018661620841759898, "doppler": 0.005972590210808326 } ] diff --git a/tests/layup/test_radar_end_to_end.py b/tests/layup/test_radar_end_to_end.py index d8f68344..bc63883c 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``; ``doppler = 2 rho_hat . v_rel``). +(``delay = 2 rho/c`` plus the Shapiro delay; ``doppler = 2 rho_hat . v_rel``). 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*: @@ -117,7 +117,23 @@ def state_at(t_jd): tau_u = rho_u / SPEED_OF_LIGHT rho_hat_u = rho_u_vec / rho_u v_tx = v_obs - a_obs * (tau_d + tau_u) - delay = tau_d + tau_u + + # Shapiro (relativistic) delay on both legs, matching orbit_fit.cpp. The truth + # here is generated to be fed back through the fitter, so it has to carry the + # same physics the fitter models -- otherwise the test measures the difference + # between two models rather than whether an orbit is recovered. The Sun is + # taken at the emission time, where the C++ light-time solution leaves it. + sun = ephem.get_particle(0, (obs_jd_tdb - tau_d) - jd_ref) # ASSIST_BODY_SUN + sun_pos = np.array([sun.x, sun.y, sun.z]) + gm_sun = 2.9591220828559115e-4 # au^3/day^2 + k = 2.0 * gm_sun / SPEED_OF_LIGHT**3 + r_b = float(np.linalg.norm(r_ast - sun_pos)) + r_r = float(np.linalg.norm(r_obs - sun_pos)) + r_t = float(np.linalg.norm(r_tx - sun_pos)) + 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_tx)) return delay, doppler diff --git a/tests/layup/test_radar_validation.py b/tests/layup/test_radar_validation.py index 0f3b5439..daeed6c2 100644 --- a/tests/layup/test_radar_validation.py +++ b/tests/layup/test_radar_validation.py @@ -64,6 +64,26 @@ def _seed(state, epoch): return g +# Chi-square floor for the noise-free fixture. +# +# The fixture is generated by ISSUE_146_gen_radar_fixture.py, which implements the +# same observation model as the C++ independently. Since the Shapiro delay was +# added, the two agree to about 0.3 ns per observation rather than exactly -- 3e-4 +# of the fixture's own 1 us stated uncertainty, and a hundred times finer than the +# best real radar timing in the JPL database (Golevka at 0.30 us). That residual +# disagreement puts a floor of ~1.4e-6 under the chi-square, so the old bound of +# 1e-6, which held only while the fixture and the fitter shared a formula, now +# measures agreement between two implementations rather than the fit. +# +# 1e-4 keeps two orders of magnitude of headroom over that floor and still catches +# a real model error by a wide margin: with the Shapiro term missing from the +# fixture this same chi-square was 1.0e-2, a hundred times above this bound. +# +# The state tolerances below carry the actual physics and are tightened +# accordingly -- recovery is 2.7e-9 au and 1.5e-11 au/day. +CSQ_FLOOR = 1e-4 + + def test_radar_fit_recovers_synthetic_orbit(): """From a seed offset from truth, the radar fit converges back to the true orbit at ~0 chi-squared -- validating both the delay/Doppler residuals @@ -77,10 +97,10 @@ def test_radar_fit_recovers_synthetic_orbit(): assert res.flag == 0 assert res.ndof == 2 * len(obs) - 6 # delay + doppler row per obs - assert res.csq < 1e-6 # noise-free data => essentially zero chi-squared + assert res.csq < CSQ_FLOOR # noise-free data => negligible chi-squared st = np.array([res.state[i] for i in range(6)]) - assert np.linalg.norm(st[:3] - truth[:3]) < 1e-6 - assert np.linalg.norm(st[3:] - truth[3:]) < 1e-7 + assert np.linalg.norm(st[:3] - truth[:3]) < 1e-7 + assert np.linalg.norm(st[3:] - truth[3:]) < 1e-9 def test_radar_rows_are_live(): @@ -89,7 +109,7 @@ def test_radar_rows_are_live(): d, obs = _load() truth = np.array(d["true_state"]) good = run_from_vector_with_initial_guess(get_ephem(CACHE), _seed(truth, d["epoch"]), obs, 50) - assert good.flag == 0 and good.csq < 1e-6 + assert good.flag == 0 and good.csq < CSQ_FLOOR o = d["observations"][0] bad_obs = list(obs) @@ -117,7 +137,7 @@ def test_delay_only_and_doppler_only_row_counts(): truth = np.array(d["true_state"]) res = run_from_vector_with_initial_guess(get_ephem(CACHE), _seed(truth, d["epoch"]), delay_only, 50) assert res.ndof == 1 * len(delay_only) - 6 # one delay row each - assert res.csq < 1e-6 + assert res.csq < CSQ_FLOOR _, doppler_only = _load(has_delay=False, has_doppler=True) res2 = run_from_vector_with_initial_guess(get_ephem(CACHE), _seed(truth, d["epoch"]), doppler_only, 50) From c3c99d9ff25cfab710424763922653cb1d72c377 Mon Sep 17 00:00:00 2001 From: matthewholman Date: Wed, 2 Sep 2026 19:59:06 -0400 Subject: [PATCH 3/4] Say that the radar model is monostatic where radar input is built An Observation carries one station, used for both the transmit and the receive leg. That is monostatic by design -- ISSUE_146_RADAR_DESIGN.md scopes the first implementation to TX == RX and lists bistatic as a refinement -- but nothing said so at the point where a caller assembles radar rows. Passing a bistatic measurement therefore evaluates the receive leg at the transmitting antenna, and does it silently. On real Goldstone bistatic pairs that is a ~4% Doppler error (observed/model 1.039-1.044 against 0.999934 for monostatic rows), which against a 0.1 Hz uncertainty is of order a hundred sigma per observation. Fitting (6489) Golevka with its 18 bistatic rows included gave flag 2 and a reduced chi-square of 1.5e7; with them removed, flag 0 at 0.80. No behaviour change -- the constraint already existed, it was just not written down where someone would meet it. --- src/layup/orbitfit.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/layup/orbitfit.py b/src/layup/orbitfit.py index 04665942..570212bc 100644 --- a/src/layup/orbitfit.py +++ b/src/layup/orbitfit.py @@ -429,6 +429,15 @@ def _radar_observation(objID, d, epoch_jd, column_names): """Build a radar ``Observation`` from a row, converting JPL units to the fitter's internal units. + Monostatic only: an ``Observation`` carries a single station, used for both + the transmit and the receive leg. A bistatic measurement -- transmitted from + one antenna and received at another -- has no way to express its second site + here, and passing one silently evaluates the receive leg at the transmitting + antenna. On real Goldstone bistatic pairs that is a ~4% Doppler error, which + against a 0.1 Hz uncertainty is of order a hundred sigma. Filter such + observations out before fitting (see ISSUE_146_RADAR_DESIGN.md, where bistatic + is listed as a refinement). + delay (us, round-trip) -> days; Doppler (Hz) -> round-trip range-rate (au/day) via the per-observation transmit frequency ``freqTx``. The barycentric observer position/velocity columns (x,y,z,vx,vy,vz) must already From 3fe967043238ef74b30da7a84ef2ac59f5ce7a96 Mon Sep 17 00:00:00 2001 From: matthewholman Date: Wed, 2 Sep 2026 21:22:03 -0400 Subject: [PATCH 4/4] Apply black formatting to the radar end-to-end test The end-to-end test was added after the last formatting pass and used hand-aligned trailing comments and a manually wrapped Shapiro expression, both of which black 25.1.0 rewrites. No functional change. Co-Authored-By: Claude Opus 5 --- tests/layup/test_radar_end_to_end.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/tests/layup/test_radar_end_to_end.py b/tests/layup/test_radar_end_to_end.py index bc63883c..3fd4e3d2 100644 --- a/tests/layup/test_radar_end_to_end.py +++ b/tests/layup/test_radar_end_to_end.py @@ -123,15 +123,16 @@ def state_at(t_jd): # same physics the fitter models -- otherwise the test measures the difference # between two models rather than whether an orbit is recovered. The Sun is # taken at the emission time, where the C++ light-time solution leaves it. - sun = ephem.get_particle(0, (obs_jd_tdb - tau_d) - jd_ref) # ASSIST_BODY_SUN + sun = ephem.get_particle(0, (obs_jd_tdb - tau_d) - jd_ref) # ASSIST_BODY_SUN sun_pos = 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 / SPEED_OF_LIGHT**3 r_b = float(np.linalg.norm(r_ast - sun_pos)) r_r = float(np.linalg.norm(r_obs - sun_pos)) r_t = float(np.linalg.norm(r_tx - sun_pos)) - 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 doppler = float(rho_hat_d @ (v_ast - v_obs)) + float(rho_hat_u @ (v_ast - v_tx))