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)