Skip to content

Use the transmitting antenna for the radar up leg (bistatic) #534

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
matthewholman merged 1 commit into from
Sep 3, 2026
+127 −19
Merged

Use the transmitting antenna for the radar up leg (bistatic) #534

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 55 additions & 0 deletions src/layup/orbitfit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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,
)


Expand Down Expand Up @@ -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:
Expand Down
33 changes: 28 additions & 5 deletions src/lib/detection.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,12 +121,26 @@ namespace orbit_fit
ObservationType observation_type;
std::array<double, 3> observer_position;
std::array<double, 3> 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<double, 3> 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<double, 3> transmitter_position{{0.0, 0.0, 0.0}};
std::array<double, 3> transmitter_velocity{{0.0, 0.0, 0.0}};
bool has_transmitter{false};

// Computed unit direction vector
Eigen::Vector3d rho_hat;
// tangent plane vectors
Expand Down Expand Up @@ -301,12 +315,18 @@ namespace orbit_fit
const std::array<double, 3> &obs_velocity,
double delay_uncy = 1.0,
double doppler_uncy = 1.0,
const std::array<double, 3> &obs_acceleration = {{0.0, 0.0, 0.0}})
const std::array<double, 3> &obs_acceleration = {{0.0, 0.0, 0.0}},
const std::array<double, 3> &tx_position = {{0.0, 0.0, 0.0}},
const std::array<double, 3> &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;
}
};
Expand Down Expand Up @@ -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<double, 3>{{0.0, 0.0, 0.0}},
py::arg("transmitter_position") = std::array<double, 3>{{0.0, 0.0, 0.0}},
py::arg("transmitter_velocity") = std::array<double, 3>{{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")
Expand Down
58 changes: 44 additions & 14 deletions src/lib/orbit_fit/orbit_fit.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
Loading