From b5573668ba2908afb06a8ba3204d765c5e9f6b5f Mon Sep 17 00:00:00 2001 From: matthewholman Date: Thu, 3 Sep 2026 13:50:21 -0400 Subject: [PATCH] Use the transmitting antenna for the radar up leg (#528) An Observation carried one station, used for both legs, so a bistatic measurement had its up leg evaluated at the antenna that received rather than the one that transmitted. On (6489) Golevka that is about 125 sigma per Doppler point -- observed/model 1.039-1.044 on bistatic rows against 0.999934 on monostatic ones -- and 20 of its 40 rows are bistatic. Observation now carries the transmitting antenna's state at the transmit epoch, and the up leg uses it directly. Where it is absent the receive station is extrapolated as before, so monostatic behaviour is unchanged. The transmit epoch does not have to be solved for: the round-trip delay is the observable, so it is t_receive - tau. Golevka's bistatic rows are all Doppler-only and carry no delay of their own, so tau is interpolated from the object's delay rows, which run smoothly from 46.5 to 47.2 s across the 1995 apparition. Evaluating the station rather than extrapolating to it also removes the Taylor truncation, measured at 2.2e-3 m/s in velocity at that round trip. Measured on Golevka, all 40 radar observations: bistatic, receive station extrapolated flag 2, reduced chi-square 1.487e7 bistatic, transmitter supplied flag 0, reduced chi-square 6.68 monostatic rows only (the old workaround) flag 0, reduced chi-square 0.80 So the fit converges on the full set for the first time, where the previous workaround required discarding 18 observations. Aten and Nyx are unchanged controls: neither has bistatic rows, and both move by less than 10%. A residual remains. All of it sits in the 20 bistatic Doppler rows, at about 5 Hz -- a round-trip range-rate error near 0.1 m/s. Ruled out so far: the bounce point (every Golevka delay row is already bp=C, Chesley et al. 2003's re-reduced values), the peak-power reference, station identity and coordinates (DSS-13 to DSS-14 checks out at 21.6 km), the velocity truncation (0.12 Hz, ~1 sigma), and uniform weight inflation (the two tight-sigma rows and the eighteen loose ones agree in Hz, so it is not a weights problem). Closes #528. --- src/layup/orbitfit.py | 55 +++++++++++++++++++++++++++++++ src/lib/detection.cpp | 33 ++++++++++++++++--- src/lib/orbit_fit/orbit_fit.cpp | 58 +++++++++++++++++++++++++-------- 3 files changed, 127 insertions(+), 19 deletions(-) diff --git a/src/layup/orbitfit.py b/src/layup/orbitfit.py index 3736a736..32d96b51 100644 --- a/src/layup/orbitfit.py +++ b/src/layup/orbitfit.py @@ -476,6 +476,13 @@ def _radar_observation(objID, d, epoch_jd, column_names): else: observer_acc = [0.0, 0.0, 0.0] + # Transmitting antenna at the transmit epoch (issue #528). When absent the + # C++ falls back to extrapolating the receive station, as it did before. + tx_cols = ("txx", "txy", "txz", "txvx", "txvy", "txvz") + has_tx = all(c in column_names for c in tx_cols) and not np.isnan(d["txx"]) + tx_pos = [float(d["txx"]), float(d["txy"]), float(d["txz"])] if has_tx else [0.0, 0.0, 0.0] + tx_vel = [float(d["txvx"]), float(d["txvy"]), float(d["txvz"])] if has_tx else [0.0, 0.0, 0.0] + return Observation.from_radar_with_id( str(objID), delay_days, @@ -488,6 +495,53 @@ def _radar_observation(objID, d, epoch_jd, column_names): delay_unc, doppler_unc, observer_acc, # Barycentric observer acceleration (au/day^2) + tx_pos, # Transmitting antenna state at the transmit epoch + tx_vel, + has_tx, + ) + + +def _append_transmitter_state(data, observatory): + """Append the transmitting antenna's state at the transmit epoch (issue #528). + + Columns ``txx, txy, txz, txvx, txvy, txvz``. The radar up leg uses these + directly instead of extrapolating the receive station back to the transmit + time. That fixes two things at once: the antenna is the one that actually + transmitted, which for a bistatic measurement is not the receiver, and the + Taylor truncation disappears. + + The transmit epoch is ``t_receive - tau``. Rows carrying a delay give tau + from the observable. Rows without one -- Doppler-only, which is what all of + (6489) Golevka's bistatic observations are -- take tau interpolated from + this object's delay rows, which vary smoothly (46.5 to 47.2 s across its + 1995 apparition). With no delay rows at all there is nothing to interpolate + from, so those rows are left without a transmitter state and the model falls + back to extrapolating. + + The transmitting antenna is named by ``trx``, the ADES field for it, so a + dataset carrying ADES radar columns needs no renaming; ``stnTx`` is accepted + as an alias. Without either, the receiving station is used, which is the + monostatic case and leaves behaviour unchanged. + """ + names = data.dtype.names + tx_col = "trx" if "trx" in names else ("stnTx" if "stnTx" in names else None) + if tx_col is None or "delay" not in names: + return data + has_delay = np.isfinite(data["delay"]) + if not has_delay.any(): + return data + order = np.argsort(data["et"][has_delay]) + tau_s = np.interp(data["et"], data["et"][has_delay][order], (data["delay"][has_delay] * 1e-6)[order]) + shifted = data.copy() + shifted["stn"] = np.array([str(t).strip() or str(v) for t, v in zip(data[tx_col], data["stn"])]) + shifted["et"] = data["et"] - tau_s + pv = np.atleast_1d(observatory.obscodes_to_barycentric(shifted)) + return rfn.append_fields( + data, + ["txx", "txy", "txz", "txvx", "txvy", "txvz"], + [pv["x"], pv["y"], pv["z"], pv["vx"], pv["vy"], pv["vz"]], + usemask=False, + asrecarray=True, ) @@ -1720,6 +1774,7 @@ def orbitfit( # columns are present so optical/streak fits are unaffected. if any(col in data.dtype.names for col in ("delay", "doppler")): data = _append_observer_acceleration(data, layup_observatory) + data = _append_transmitter_state(data, layup_observatory) bias_dict = None if debias: diff --git a/src/lib/detection.cpp b/src/lib/detection.cpp index 84668bdb..49e23b01 100644 --- a/src/lib/detection.cpp +++ b/src/lib/detection.cpp @@ -121,12 +121,26 @@ namespace orbit_fit ObservationType observation_type; std::array observer_position; std::array observer_velocity; - // Barycentric observer acceleration (au/day^2). Only used by the radar - // two-leg light-time model, which Taylor-extrapolates the station state - // back to the signal transmit time (~one round-trip earlier). Defaults to - // zero, so it has no effect on optical/streak observations. + // Barycentric observer acceleration (au/day^2). Retained for the radar + // two-leg light-time model's fallback path, which Taylor-extrapolates the + // station state back to the signal transmit time. Defaults to zero, so it + // has no effect on optical/streak observations. std::array observer_acceleration{{0.0, 0.0, 0.0}}; + // Barycentric state of the TRANSMITTING antenna AT THE TRANSMIT EPOCH + // (issue #528), so the up leg neither uses the wrong antenna nor + // extrapolates. For a bistatic measurement the transmitter is not the + // receiver: on (6489) Golevka that gives observed/model Doppler + // 1.039-1.044 against 0.999934 monostatic, ~125 sigma at a 0.11 Hz + // uncertainty. Evaluating rather than extrapolating also removes the + // truncation, measured at 2.2e-3 m/s in velocity over Golevka's ~47 s + // round trip. When has_transmitter is false the receive station is + // extrapolated instead, which is correct when TX == RX and is what the + // model did before. + std::array transmitter_position{{0.0, 0.0, 0.0}}; + std::array transmitter_velocity{{0.0, 0.0, 0.0}}; + bool has_transmitter{false}; + // Computed unit direction vector Eigen::Vector3d rho_hat; // tangent plane vectors @@ -301,12 +315,18 @@ namespace orbit_fit const std::array &obs_velocity, double delay_uncy = 1.0, double doppler_uncy = 1.0, - const std::array &obs_acceleration = {{0.0, 0.0, 0.0}}) + const std::array &obs_acceleration = {{0.0, 0.0, 0.0}}, + const std::array &tx_position = {{0.0, 0.0, 0.0}}, + const std::array &tx_velocity = {{0.0, 0.0, 0.0}}, + bool has_tx = false) { Observation obs = from_radar(delay, doppler, has_delay, has_doppler, epoch_val, obs_position, obs_velocity, delay_uncy, doppler_uncy, obs_acceleration); obs.objID = objID; + obs.transmitter_position = tx_position; + obs.transmitter_velocity = tx_velocity; + obs.has_transmitter = has_tx; return obs; } }; @@ -378,6 +398,9 @@ namespace orbit_fit py::arg("epoch"), py::arg("observer_position"), py::arg("observer_velocity"), py::arg("delay_unc") = 1.0, py::arg("doppler_unc") = 1.0, py::arg("observer_acceleration") = std::array{{0.0, 0.0, 0.0}}, + py::arg("transmitter_position") = std::array{{0.0, 0.0, 0.0}}, + py::arg("transmitter_velocity") = std::array{{0.0, 0.0, 0.0}}, + py::arg("has_transmitter") = false, "Construct a Radar observation (delay in days, doppler in au/day)") .def_readwrite("epoch", &Observation::epoch, "Observation epoch (as a double)") .def_readwrite("observation_type", &Observation::observation_type, "Variant holding the observation data") diff --git a/src/lib/orbit_fit/orbit_fit.cpp b/src/lib/orbit_fit/orbit_fit.cpp index 33d51e97..7a6ca171 100644 --- a/src/lib/orbit_fit/orbit_fit.cpp +++ b/src/lib/orbit_fit/orbit_fit.cpp @@ -156,34 +156,64 @@ namespace orbit_fit double aoy = this_det.observer_acceleration[1]; double aoz = this_det.observer_acceleration[2]; - // Up-leg light time: iterate tau_u with the station at the transmit - // time t_obs - (tau_d + tau_u), Taylor-extrapolated from the receive - // epoch. rho is the converged down-leg distance from the caller. + // Up-leg light time. rho is the converged down-leg distance from the + // caller. Two ways to place the transmitting station: + // + // supplied (issue #528): the caller evaluated the transmitting antenna + // at the transmit epoch, t_receive - tau, and passed its state. Exact, + // and correct for a bistatic measurement, where the transmitter is not + // the antenna that received. Feeding a bistatic row through the + // fallback below is a large error: on (6489) Golevka the modelled + // Doppler is wrong by ~125 sigma per point, and all 20 of its bistatic + // rows are Doppler-only. + // + // fallback: Taylor-extrapolate the receive station back by the round + // trip. Correct only when TX == RX, and truncated -- the neglected + // velocity term is 2.2e-3 m/s over Golevka's ~47 s round trip. + // + // Either way rtx_* holds the transmit-time station position, which the + // up leg's Shapiro term below needs. 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 vtx_x, vtx_y, vtx_z; + if (this_det.has_transmitter) { - double tau = tau_d + tau_u; - 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; + rtx_x = this_det.transmitter_position[0]; + rtx_y = this_det.transmitter_position[1]; + rtx_z = this_det.transmitter_position[2]; rhu_x = rbx - rtx_x; rhu_y = rby - rtx_y; rhu_z = rbz - rtx_z; rho_u = sqrt(rhu_x * rhu_x + rhu_y * rhu_y + rhu_z * rhu_z); tau_u = rho_u / SPEED_OF_LIGHT; + vtx_x = this_det.transmitter_velocity[0]; + vtx_y = this_det.transmitter_velocity[1]; + vtx_z = this_det.transmitter_velocity[2]; + } + else + { + for (int it = 0; it < 3; it++) + { + double tau = tau_d + tau_u; + 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; + rho_u = sqrt(rhu_x * rhu_x + rhu_y * rhu_y + rhu_z * rhu_z); + tau_u = rho_u / SPEED_OF_LIGHT; + } + const double tau = tau_d + tau_u; + vtx_x = vox - aox * tau; // station velocity at transmit + vtx_y = voy - aoy * tau; + vtx_z = voz - aoz * tau; } rhu_x /= rho_u; // up-leg unit vector (bounce -> transmit station) rhu_y /= rho_u; rhu_z /= rho_u; - double tau = tau_d + tau_u; - double vtx_x = vox - aox * tau; // station velocity at transmit - double vtx_y = voy - aoy * tau; - double vtx_z = voz - aoz * tau; double model_delay = tau_d + tau_u; // round-trip light time (days)