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 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); 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..3fd4e3d2 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,24 @@ 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)