Model the Shapiro delay in radar residuals

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

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<RadarObservation>(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);

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
     }
   ]

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
 

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)