Summary
In the velocity computation of calc_satellite_position_and_velocity, the true-anomaly rate true_anomaly_dot is computed by dividing by sin(true_anomaly), which is 0 at perigee (ν = 0) and apogee (ν = π). At those points sin(eccentric_anomaly) is also 0, so the expression is 0/0 → NaN.
Location
src/sat_position.jl, in calc_satellite_position_and_velocity:
true_anomaly_dot =
sin(eccentric_anomaly) *
eccentric_anomaly_dot *
(1.0 + data.e * cos(true_anomaly)) /
(sin(true_anomaly) * (1.0 - data.e * cos(eccentric_anomaly)))
Why this is a problem
The expression is algebraically correct, but numerically singular. Using the
identity sin ν = √(1−e²)·sin E / (1 − e·cos E), the sin E / sin ν ratio is the
finite quantity (1 − e·cos E)/√(1−e²). But evaluating it as sin(E)/sin(ν):
- at
ν = 0 (perigee): E = 0 ⇒ sin(E) = sin(ν) = 0 ⇒ 0/0 = NaN,
- at
ν = π (apogee): E = π ⇒ sin(E) = sin(ν) = 0 ⇒ 0/0 = NaN.
A NaN here propagates into the entire satellite velocity vector (it feeds
corrected_argument_of_latitude_dot, corrected_radius_dot, and
corrected_inclination_dot), and from there into the receiver
velocity / clock-drift solve.
GNSS orbits have small eccentricity, but ν sweeps through 0 and π every
orbital period, so a satellite passes near perigee/apogee twice per orbit. An
exact hit (E landing on 0/π in floating point) gives a hard NaN; this is
easy to trigger with synthetic ephemerides and possible with real data.
Suggested fix
Use the singularity-free closed form (reuses the already-computed
eccentric_anomaly_dot and true_anomaly):
true_anomaly_dot =
eccentric_anomaly_dot * (1.0 + data.e * cos(true_anomaly)) / sqrt(1.0 - data.e^2)
This is exactly equal to the current expression away from the singular points
(Ė·(1+e·cos ν)/√(1−e²) = Ė·√(1−e²)/(1−e·cos E)) and is finite everywhere. An
equivalent E-only form is eccentric_anomaly_dot * sqrt(1 - data.e^2) / (1 - data.e * cos(eccentric_anomaly)).
Test
Add a regression test that propagates a satellite positioned at (or very near)
perigee/apogee and asserts the returned velocity is finite.
Summary
In the velocity computation of
calc_satellite_position_and_velocity, the true-anomaly ratetrue_anomaly_dotis computed by dividing bysin(true_anomaly), which is0at perigee (ν = 0) and apogee (ν = π). At those pointssin(eccentric_anomaly)is also0, so the expression is0/0 → NaN.Location
src/sat_position.jl, incalc_satellite_position_and_velocity:Why this is a problem
The expression is algebraically correct, but numerically singular. Using the
identity
sin ν = √(1−e²)·sin E / (1 − e·cos E), thesin E / sin νratio is thefinite quantity
(1 − e·cos E)/√(1−e²). But evaluating it assin(E)/sin(ν):ν = 0(perigee):E = 0⇒sin(E) = sin(ν) = 0⇒0/0 = NaN,ν = π(apogee):E = π⇒sin(E) = sin(ν) = 0⇒0/0 = NaN.A
NaNhere propagates into the entire satellite velocity vector (it feedscorrected_argument_of_latitude_dot,corrected_radius_dot, andcorrected_inclination_dot), and from there into the receivervelocity / clock-drift solve.
GNSS orbits have small eccentricity, but
νsweeps through0andπeveryorbital period, so a satellite passes near perigee/apogee twice per orbit. An
exact hit (
Elanding on0/πin floating point) gives a hardNaN; this iseasy to trigger with synthetic ephemerides and possible with real data.
Suggested fix
Use the singularity-free closed form (reuses the already-computed
eccentric_anomaly_dotandtrue_anomaly):This is exactly equal to the current expression away from the singular points
(
Ė·(1+e·cos ν)/√(1−e²) = Ė·√(1−e²)/(1−e·cos E)) and is finite everywhere. Anequivalent
E-only form iseccentric_anomaly_dot * sqrt(1 - data.e^2) / (1 - data.e * cos(eccentric_anomaly)).Test
Add a regression test that propagates a satellite positioned at (or very near)
perigee/apogee and asserts the returned velocity is finite.