sidereon_core/broadcast/
mod.rs

1//! Broadcast-ephemeris orbit and clock evaluation (GPS LNAV, Galileo I/NAV,
2//! BeiDou D1/D2).
3//!
4//! Evaluates a broadcast navigation message into an ECEF satellite position and
5//! a satellite clock offset, by the standard Keplerian construction of
6//! IS-GPS-200 (Section 20.3.3.4.3.1, Table 20-IV) and the equivalent Galileo OS
7//! and BeiDou SIS ICD sections. The constellations share the algorithm and
8//! differ in the gravitational constant, the Earth-rotation rate, and the
9//! relativistic clock constant ([`ConstellationConstants`]). BeiDou
10//! geostationary satellites additionally take a custom-frame-to-ECEF rotation
11//! (the `is_geo` path of [`satellite_position_ecef`]); GPS, Galileo, and BeiDou
12//! MEO/IGSO satellites use the direct rotation.
13//!
14//! This is a 0-ULP parity target: the operation order reproduces the canonical
15//! reference recipe (`parity/generator/broadcast_eval.py`) bit-for-bit. The
16//! transcendentals are the libm scalar `sin`/`cos`/`sqrt`/`atan2` (matching the
17//! recipe's CPython `math` calls on the pinned Apple-libm target); separate
18//! `.sin()` and `.cos()` calls are used deliberately (never `.sin_cos()`, whose
19//! fused evaluation can differ in the last bit). Integer powers are explicit
20//! repeated multiplies and there is no fused multiply-add (Rust does not
21//! auto-contract `a * b + c`).
22
23use crate::astro::constants::models::broadcast::{
24    BEIDOU_OMEGA_E_RAD_S, GALILEO_BEIDOU_DTR_F, GALILEO_GM_M3_S2, GPS_DTR_F,
25    GPS_GALILEO_OMEGA_E_RAD_S, GPS_GM_M3_S2,
26};
27use crate::error::{Error, Result};
28use crate::frame::{FrameValueError, ItrfPositionM};
29
30/// Half a week, the fold threshold for a time difference against `toe`/`toc`.
31pub use crate::constants::HALF_WEEK_S;
32/// Seconds in one GPS/Galileo week.
33pub use crate::constants::SECONDS_PER_WEEK;
34
35/// Eccentric-anomaly fixed-point convergence threshold (radians).
36pub const KEPLER_TOL: f64 = 1.0e-12;
37/// Maximum eccentric-anomaly fixed-point iterations.
38pub const KEPLER_MAX_ITER: usize = 30;
39/// Satellite-clock time-argument refinement count (RTKLIB `eph2clk` convention).
40pub const CLOCK_MAX_ITER: usize = 2;
41
42/// Per-constellation physical constants used by the broadcast evaluation.
43///
44/// The literals match the values the broadcast reference recipe and the `rinex`
45/// crate use, so the Python and Rust sides share identical `f64` bit patterns.
46#[derive(Debug, Clone, Copy, PartialEq)]
47pub struct ConstellationConstants {
48    /// Gravitational constant GM (m^3 / s^2).
49    pub gm_m3_s2: f64,
50    /// Earth rotation rate used by the longitude-of-node (Sagnac) term (rad/s).
51    pub omega_e_rad_s: f64,
52    /// Relativistic clock constant `F = -2 * sqrt(GM) / c^2` (s / sqrt(m)).
53    pub dtr_f: f64,
54}
55
56impl ConstellationConstants {
57    /// GPS constants (IS-GPS-200).
58    pub const GPS: Self = Self {
59        gm_m3_s2: GPS_GM_M3_S2,
60        omega_e_rad_s: GPS_GALILEO_OMEGA_E_RAD_S,
61        dtr_f: GPS_DTR_F,
62    };
63    /// Galileo constants (OS SIS ICD); shares the GPS rotation rate.
64    pub const GALILEO: Self = Self {
65        gm_m3_s2: GALILEO_GM_M3_S2,
66        omega_e_rad_s: GPS_GALILEO_OMEGA_E_RAD_S,
67        dtr_f: GALILEO_BEIDOU_DTR_F,
68    };
69    /// BeiDou constants (BDS-SIS-ICD); its own Earth-rotation rate.
70    pub const BEIDOU: Self = Self {
71        gm_m3_s2: GALILEO_GM_M3_S2,
72        omega_e_rad_s: BEIDOU_OMEGA_E_RAD_S,
73        dtr_f: GALILEO_BEIDOU_DTR_F,
74    };
75}
76
77/// Broadcast Keplerian orbital elements (SI units; angles in radians; `toe_sow`
78/// in seconds of the constellation's week).
79#[derive(Debug, Clone, Copy, PartialEq)]
80pub struct KeplerianElements {
81    /// Square root of the semi-major axis (sqrt(m)).
82    pub sqrt_a: f64,
83    /// Eccentricity (dimensionless).
84    pub e: f64,
85    /// Mean anomaly at reference time (rad).
86    pub m0: f64,
87    /// Mean motion difference from computed value (rad/s).
88    pub delta_n: f64,
89    /// Longitude of ascending node at weekly epoch (rad).
90    pub omega0: f64,
91    /// Inclination at reference time (rad).
92    pub i0: f64,
93    /// Argument of perigee (rad).
94    pub omega: f64,
95    /// Rate of right ascension (rad/s).
96    pub omega_dot: f64,
97    /// Rate of inclination (rad/s).
98    pub idot: f64,
99    /// Latitude argument cosine correction (rad).
100    pub cuc: f64,
101    /// Latitude argument sine correction (rad).
102    pub cus: f64,
103    /// Orbit radius cosine correction (m).
104    pub crc: f64,
105    /// Orbit radius sine correction (m).
106    pub crs: f64,
107    /// Inclination cosine correction (rad).
108    pub cic: f64,
109    /// Inclination sine correction (rad).
110    pub cis: f64,
111    /// Ephemeris reference time, seconds of week.
112    pub toe_sow: f64,
113}
114
115/// Broadcast satellite-clock polynomial about `toc_sow`.
116#[derive(Debug, Clone, Copy, PartialEq)]
117pub struct ClockPolynomial {
118    /// Clock bias (s).
119    pub af0: f64,
120    /// Clock drift (s/s).
121    pub af1: f64,
122    /// Clock drift rate (s/s^2).
123    pub af2: f64,
124    /// Clock reference time, seconds of week.
125    pub toc_sow: f64,
126}
127
128/// A solved eccentric anomaly and the iteration count that produced it.
129#[derive(Debug, Clone, Copy, PartialEq)]
130pub struct EccentricAnomaly {
131    /// Eccentric anomaly (rad).
132    pub value: f64,
133    /// Fixed-point iterations performed.
134    pub iterations: usize,
135}
136
137/// The full intermediate substrate of a broadcast orbit evaluation.
138///
139/// Every field is exposed so a 0-ULP parity test can localize a mismatch to a
140/// single operation. [`OrbitState::position`] returns the ECEF position.
141#[derive(Debug, Clone, Copy, PartialEq)]
142pub struct OrbitState {
143    /// Semi-major axis (m).
144    pub a: f64,
145    /// Computed mean motion (rad/s).
146    pub n0: f64,
147    /// Corrected mean motion (rad/s).
148    pub n: f64,
149    /// Time from ephemeris reference epoch, half-week folded (s).
150    pub tk: f64,
151    /// Mean anomaly (rad).
152    pub mk: f64,
153    /// Eccentric anomaly (rad).
154    pub eccentric_anomaly: f64,
155    /// Number of Kepler iterations.
156    pub kepler_iterations: usize,
157    /// sin(E).
158    pub sin_e: f64,
159    /// cos(E).
160    pub cos_e: f64,
161    /// True anomaly (rad).
162    pub nu: f64,
163    /// Argument of latitude before correction (rad).
164    pub phi: f64,
165    /// sin(2*phi).
166    pub s2: f64,
167    /// cos(2*phi).
168    pub c2: f64,
169    /// Argument-of-latitude correction (rad).
170    pub du: f64,
171    /// Radius correction (m).
172    pub dr: f64,
173    /// Inclination correction (rad).
174    pub di: f64,
175    /// Corrected argument of latitude (rad).
176    pub u: f64,
177    /// Corrected radius (m).
178    pub r: f64,
179    /// Corrected inclination (rad).
180    pub i: f64,
181    /// Orbital-plane x (m).
182    pub xp: f64,
183    /// Orbital-plane y (m).
184    pub yp: f64,
185    /// Corrected longitude of ascending node (rad).
186    pub omega_k: f64,
187    /// ECEF x (m).
188    pub x_m: f64,
189    /// ECEF y (m).
190    pub y_m: f64,
191    /// ECEF z (m).
192    pub z_m: f64,
193}
194
195impl OrbitState {
196    /// The Earth-fixed (ITRF/ECEF) satellite position in meters.
197    pub const fn position(&self) -> core::result::Result<ItrfPositionM, FrameValueError> {
198        ItrfPositionM::new(self.x_m, self.y_m, self.z_m)
199    }
200}
201
202/// The satellite clock offset, split into its components.
203#[derive(Debug, Clone, Copy, PartialEq)]
204pub struct ClockOffset {
205    /// Polynomial term (s).
206    pub dt_clock_poly_s: f64,
207    /// Relativistic eccentricity term (s).
208    pub dt_rel_s: f64,
209    /// Group delay subtracted for the single-frequency user (s).
210    pub tgd_s: f64,
211    /// Total satellite clock offset (s).
212    pub dt_clock_total_s: f64,
213}
214
215/// Broadcast relativistic satellite-clock correction, seconds.
216///
217/// Evaluates the periodic eccentric-orbit term
218/// `F * e * sqrt(A) * sin(E)` from IS-GPS-200 and the equivalent Galileo/BeiDou
219/// broadcast-clock models. Use [`ConstellationConstants::dtr_f`] for `F` when
220/// evaluating a full broadcast record.
221pub fn relativistic_clock_correction_s(
222    dtr_f_s_sqrt_m: f64,
223    eccentricity: f64,
224    sqrt_a_m_sqrt: f64,
225    eccentric_anomaly_sin: f64,
226) -> Result<f64> {
227    validate_finite(dtr_f_s_sqrt_m, "dtr_f_s_sqrt_m")?;
228    validate_eccentricity(eccentricity)?;
229    validate_positive(sqrt_a_m_sqrt, "sqrt_a_m_sqrt")?;
230    validate_finite(eccentric_anomaly_sin, "eccentric_anomaly_sin")?;
231    let correction = relativistic_clock_correction_s_unchecked(
232        dtr_f_s_sqrt_m,
233        eccentricity,
234        sqrt_a_m_sqrt,
235        eccentric_anomaly_sin,
236    );
237    validate_finite(correction, "relativistic_clock_correction_s")?;
238    Ok(correction)
239}
240
241#[inline]
242pub(crate) fn relativistic_clock_correction_s_unchecked(
243    dtr_f_s_sqrt_m: f64,
244    eccentricity: f64,
245    sqrt_a_m_sqrt: f64,
246    eccentric_anomaly_sin: f64,
247) -> f64 {
248    dtr_f_s_sqrt_m * eccentricity * sqrt_a_m_sqrt * eccentric_anomaly_sin
249}
250
251/// Time difference `t - t_ref` (seconds), folded into the +/- half-week range.
252///
253/// Used for both `tk = t - toe` (orbit) and `t - toc` (clock); the fold handles
254/// a query that straddles the week rollover relative to the reference.
255pub fn time_from_reference_s(t_sow_s: f64, t_ref_sow_s: f64) -> f64 {
256    let mut dt = t_sow_s - t_ref_sow_s;
257    if dt > HALF_WEEK_S {
258        dt -= SECONDS_PER_WEEK;
259    }
260    if dt < -HALF_WEEK_S {
261        dt += SECONDS_PER_WEEK;
262    }
263    dt
264}
265
266/// Solve Kepler's equation by fixed-point iteration `E = M + e*sin(E)`, seeded
267/// at `E = M`; stops on `|dE| <= KEPLER_TOL` or after `KEPLER_MAX_ITER` steps.
268pub fn eccentric_anomaly(mean_anomaly_rad: f64, eccentricity: f64) -> Result<EccentricAnomaly> {
269    validate_finite(mean_anomaly_rad, "mean_anomaly_rad")?;
270    validate_eccentricity(eccentricity)?;
271
272    Ok(eccentric_anomaly_unchecked(mean_anomaly_rad, eccentricity))
273}
274
275pub(crate) fn eccentric_anomaly_unchecked(
276    mean_anomaly_rad: f64,
277    eccentricity: f64,
278) -> EccentricAnomaly {
279    let mut e_k = mean_anomaly_rad;
280    let mut iterations = 0usize;
281    while iterations < KEPLER_MAX_ITER {
282        let e_prev = e_k;
283        e_k = mean_anomaly_rad + eccentricity * e_prev.sin();
284        iterations += 1;
285        let delta = (e_k - e_prev).abs();
286        if delta <= KEPLER_TOL {
287            break;
288        }
289    }
290    EccentricAnomaly {
291        value: e_k,
292        iterations,
293    }
294}
295
296/// Evaluate the broadcast Keplerian orbit at `t_sow_s` (seconds of week).
297///
298/// `is_geo` selects the BeiDou geostationary path (the node omits the
299/// Earth-rotation-during-`tk` term and the position is rotated to ECEF by
300/// `Rz(omega_e*tk) . Rx(-5deg)`); GPS, Galileo, and BeiDou MEO/IGSO use `false`.
301/// The statement order reproduces `broadcast_eval.satellite_position_ecef`.
302pub fn satellite_position_ecef(
303    elements: &KeplerianElements,
304    consts: &ConstellationConstants,
305    t_sow_s: f64,
306    is_geo: bool,
307) -> Result<OrbitState> {
308    validate_elements(elements)?;
309    validate_constants(consts)?;
310    validate_finite(t_sow_s, "t_sow_s")?;
311
312    let state = satellite_position_ecef_unchecked(elements, consts, t_sow_s, is_geo);
313    validate_orbit_state(&state)?;
314    Ok(state)
315}
316
317pub(crate) fn satellite_position_ecef_unchecked(
318    elements: &KeplerianElements,
319    consts: &ConstellationConstants,
320    t_sow_s: f64,
321    is_geo: bool,
322) -> OrbitState {
323    let sqrt_a = elements.sqrt_a;
324    let e = elements.e;
325    let gm = consts.gm_m3_s2;
326    let omega_e = consts.omega_e_rad_s;
327
328    // 1. Semi-major axis and mean motion. a^3 as an explicit multiply chain.
329    let a = sqrt_a * sqrt_a;
330    let n0 = (gm / (a * a * a)).sqrt();
331    let n = n0 + elements.delta_n;
332
333    // 2. Time from ephemeris reference epoch (half-week folded).
334    let tk = time_from_reference_s(t_sow_s, elements.toe_sow);
335
336    // 3. Mean anomaly and eccentric anomaly.
337    let mk = elements.m0 + n * tk;
338    let kepler = eccentric_anomaly_unchecked(mk, e);
339    let ecc_anom = kepler.value;
340    let sin_e = ecc_anom.sin();
341    let cos_e = ecc_anom.cos();
342
343    // 4. True anomaly (atan2 form) and argument of latitude.
344    let e2 = e * e;
345    let nu = ((1.0 - e2).sqrt() * sin_e).atan2(cos_e - e);
346    let phi = nu + elements.omega;
347
348    // 5. Second-harmonic corrections (sine term first).
349    let two_phi = 2.0 * phi;
350    let s2 = two_phi.sin();
351    let c2 = two_phi.cos();
352    let du = elements.cus * s2 + elements.cuc * c2;
353    let dr = elements.crs * s2 + elements.crc * c2;
354    let di = elements.cis * s2 + elements.cic * c2;
355
356    // 6. Corrected argument of latitude, radius, inclination.
357    let u = phi + du;
358    let r = a * (1.0 - e * cos_e) + dr;
359    let i = elements.i0 + di + elements.idot * tk;
360
361    // 7. Position in the orbital plane.
362    let xp = r * u.cos();
363    let yp = r * u.sin();
364
365    // 8. Corrected longitude of ascending node. The BeiDou GEO node omits the
366    // Earth-rotation-during-tk term (applied by the final rotation instead).
367    let omega_k = if is_geo {
368        elements.omega0 + elements.omega_dot * tk - omega_e * elements.toe_sow
369    } else {
370        elements.omega0 + (elements.omega_dot - omega_e) * tk - omega_e * elements.toe_sow
371    };
372
373    // 9. Coordinates in the (custom) frame from the node rotation.
374    let sin_o = omega_k.sin();
375    let cos_o = omega_k.cos();
376    let sin_i = i.sin();
377    let cos_i = i.cos();
378    let xg = xp * cos_o - yp * cos_i * sin_o;
379    let yg = xp * sin_o + yp * cos_i * cos_o;
380    let zg = yp * sin_i;
381
382    // 10. Earth-fixed coordinates. The standard path is the identity; the BeiDou
383    // GEO path applies Rz(omega_e*tk) . Rx(-5deg) (BDS-SIS-ICD).
384    let (x, y, z) = if is_geo {
385        let deg5 = 5.0_f64.to_radians();
386        let cos_phi = deg5.cos();
387        let sin_phi = -deg5.sin();
388        let z_ang = omega_e * tk;
389        let cos_z = z_ang.cos();
390        let sin_z = z_ang.sin();
391        let yr = yg * cos_phi + zg * sin_phi;
392        let zr = -yg * sin_phi + zg * cos_phi;
393        (xg * cos_z + yr * sin_z, -xg * sin_z + yr * cos_z, zr)
394    } else {
395        (xg, yg, zg)
396    };
397
398    OrbitState {
399        a,
400        n0,
401        n,
402        tk,
403        mk,
404        eccentric_anomaly: ecc_anom,
405        kepler_iterations: kepler.iterations,
406        sin_e,
407        cos_e,
408        nu,
409        phi,
410        s2,
411        c2,
412        du,
413        dr,
414        di,
415        u,
416        r,
417        i,
418        xp,
419        yp,
420        omega_k,
421        x_m: x,
422        y_m: y,
423        z_m: z,
424    }
425}
426
427/// Evaluate the broadcast satellite clock offset (seconds).
428///
429/// `sin_e` is the eccentric-anomaly sine from the position evaluation at the
430/// same instant; `tgd_s` is the single-frequency group delay. The statement
431/// order reproduces `broadcast_eval.satellite_clock_offset_s`.
432pub fn satellite_clock_offset_s(
433    clock: &ClockPolynomial,
434    consts: &ConstellationConstants,
435    elements: &KeplerianElements,
436    sin_e: f64,
437    t_sow_s: f64,
438    tgd_s: f64,
439) -> Result<ClockOffset> {
440    validate_clock(clock)?;
441    validate_constants(consts)?;
442    validate_elements(elements)?;
443    validate_finite(sin_e, "sin_e")?;
444    validate_finite(t_sow_s, "t_sow_s")?;
445    validate_finite(tgd_s, "tgd_s")?;
446
447    let offset = satellite_clock_offset_s_unchecked(clock, consts, elements, sin_e, t_sow_s, tgd_s);
448    validate_clock_offset(&offset)?;
449    Ok(offset)
450}
451
452pub(crate) fn satellite_clock_offset_s_unchecked(
453    clock: &ClockPolynomial,
454    consts: &ConstellationConstants,
455    elements: &KeplerianElements,
456    sin_e: f64,
457    t_sow_s: f64,
458    tgd_s: f64,
459) -> ClockOffset {
460    let af0 = clock.af0;
461    let af1 = clock.af1;
462    let af2 = clock.af2;
463
464    // Time from clock reference, folded; then refine out the SV clock itself.
465    let dt0 = time_from_reference_s(t_sow_s, clock.toc_sow);
466    let mut dt = dt0;
467    let mut refine = 0usize;
468    while refine < CLOCK_MAX_ITER {
469        dt = dt0 - (af0 + af1 * dt + af2 * dt * dt);
470        refine += 1;
471    }
472    let dt_poly = af0 + af1 * dt + af2 * dt * dt;
473
474    // Relativistic eccentricity term (sqrt_a is the broadcast sqrt(A)).
475    let dt_rel =
476        relativistic_clock_correction_s_unchecked(consts.dtr_f, elements.e, elements.sqrt_a, sin_e);
477
478    let dt_total = dt_poly + dt_rel - tgd_s;
479
480    ClockOffset {
481        dt_clock_poly_s: dt_poly,
482        dt_rel_s: dt_rel,
483        tgd_s,
484        dt_clock_total_s: dt_total,
485    }
486}
487
488/// A satellite's broadcast orbit and clock evaluated together at one instant.
489#[derive(Debug, Clone, Copy, PartialEq)]
490pub struct SatelliteState {
491    /// The orbit evaluation (ECEF position and all intermediates).
492    pub orbit: OrbitState,
493    /// The clock offset evaluation.
494    pub clock: ClockOffset,
495}
496
497/// Evaluate the broadcast orbit and clock at the same instant.
498///
499/// This is the intended public entry point: it solves the orbit and feeds that
500/// solution's eccentric-anomaly sine into the clock's relativistic term, so a
501/// caller cannot accidentally pair a clock evaluation with `sin(E)` from a
502/// different epoch. [`satellite_position_ecef`] and [`satellite_clock_offset_s`]
503/// remain available for component-level parity testing.
504pub fn satellite_state(
505    elements: &KeplerianElements,
506    clock: &ClockPolynomial,
507    consts: &ConstellationConstants,
508    t_sow_s: f64,
509    tgd_s: f64,
510    is_geo: bool,
511) -> Result<SatelliteState> {
512    validate_elements(elements)?;
513    validate_clock(clock)?;
514    validate_constants(consts)?;
515    validate_finite(t_sow_s, "t_sow_s")?;
516    validate_finite(tgd_s, "tgd_s")?;
517
518    let state = satellite_state_unchecked(elements, clock, consts, t_sow_s, tgd_s, is_geo);
519    validate_orbit_state(&state.orbit)?;
520    validate_clock_offset(&state.clock)?;
521    Ok(state)
522}
523
524pub(crate) fn satellite_state_unchecked(
525    elements: &KeplerianElements,
526    clock: &ClockPolynomial,
527    consts: &ConstellationConstants,
528    t_sow_s: f64,
529    tgd_s: f64,
530    is_geo: bool,
531) -> SatelliteState {
532    let orbit = satellite_position_ecef_unchecked(elements, consts, t_sow_s, is_geo);
533    let clock =
534        satellite_clock_offset_s_unchecked(clock, consts, elements, orbit.sin_e, t_sow_s, tgd_s);
535    SatelliteState { orbit, clock }
536}
537
538fn validate_elements(elements: &KeplerianElements) -> Result<()> {
539    validate_positive(elements.sqrt_a, "elements.sqrt_a")?;
540    validate_eccentricity(elements.e)?;
541    validate_finite(elements.m0, "elements.m0")?;
542    validate_finite(elements.delta_n, "elements.delta_n")?;
543    validate_finite(elements.omega0, "elements.omega0")?;
544    validate_finite(elements.i0, "elements.i0")?;
545    validate_finite(elements.omega, "elements.omega")?;
546    validate_finite(elements.omega_dot, "elements.omega_dot")?;
547    validate_finite(elements.idot, "elements.idot")?;
548    validate_finite(elements.cuc, "elements.cuc")?;
549    validate_finite(elements.cus, "elements.cus")?;
550    validate_finite(elements.crc, "elements.crc")?;
551    validate_finite(elements.crs, "elements.crs")?;
552    validate_finite(elements.cic, "elements.cic")?;
553    validate_finite(elements.cis, "elements.cis")?;
554    validate_sow(elements.toe_sow, "elements.toe_sow")
555}
556
557fn validate_clock(clock: &ClockPolynomial) -> Result<()> {
558    validate_finite(clock.af0, "clock.af0")?;
559    validate_finite(clock.af1, "clock.af1")?;
560    validate_finite(clock.af2, "clock.af2")?;
561    validate_sow(clock.toc_sow, "clock.toc_sow")
562}
563
564fn validate_constants(consts: &ConstellationConstants) -> Result<()> {
565    validate_positive(consts.gm_m3_s2, "consts.gm_m3_s2")?;
566    validate_finite(consts.omega_e_rad_s, "consts.omega_e_rad_s")?;
567    validate_finite(consts.dtr_f, "consts.dtr_f")
568}
569
570fn validate_orbit_state(state: &OrbitState) -> Result<()> {
571    validate_finite(state.a, "orbit.a")?;
572    validate_finite(state.n0, "orbit.n0")?;
573    validate_finite(state.n, "orbit.n")?;
574    validate_finite(state.tk, "orbit.tk")?;
575    validate_finite(state.mk, "orbit.mk")?;
576    validate_finite(state.eccentric_anomaly, "orbit.eccentric_anomaly")?;
577    validate_finite(state.sin_e, "orbit.sin_e")?;
578    validate_finite(state.cos_e, "orbit.cos_e")?;
579    validate_finite(state.nu, "orbit.nu")?;
580    validate_finite(state.phi, "orbit.phi")?;
581    validate_finite(state.s2, "orbit.s2")?;
582    validate_finite(state.c2, "orbit.c2")?;
583    validate_finite(state.du, "orbit.du")?;
584    validate_finite(state.dr, "orbit.dr")?;
585    validate_finite(state.di, "orbit.di")?;
586    validate_finite(state.u, "orbit.u")?;
587    validate_finite(state.r, "orbit.r")?;
588    validate_finite(state.i, "orbit.i")?;
589    validate_finite(state.xp, "orbit.xp")?;
590    validate_finite(state.yp, "orbit.yp")?;
591    validate_finite(state.omega_k, "orbit.omega_k")?;
592    validate_finite(state.x_m, "orbit.x_m")?;
593    validate_finite(state.y_m, "orbit.y_m")?;
594    validate_finite(state.z_m, "orbit.z_m")
595}
596
597fn validate_clock_offset(clock: &ClockOffset) -> Result<()> {
598    validate_finite(clock.dt_clock_poly_s, "clock.dt_clock_poly_s")?;
599    validate_finite(clock.dt_rel_s, "clock.dt_rel_s")?;
600    validate_finite(clock.tgd_s, "clock.tgd_s")?;
601    validate_finite(clock.dt_clock_total_s, "clock.dt_clock_total_s")
602}
603
604fn validate_eccentricity(eccentricity: f64) -> Result<()> {
605    validate_finite(eccentricity, "eccentricity")?;
606    if (0.0..1.0).contains(&eccentricity) {
607        Ok(())
608    } else {
609        Err(invalid_input("eccentricity", "out of range"))
610    }
611}
612
613fn validate_sow(value: f64, field: &'static str) -> Result<()> {
614    validate_finite(value, field)?;
615    if (0.0..SECONDS_PER_WEEK).contains(&value) {
616        Ok(())
617    } else {
618        Err(invalid_input(field, "out of range"))
619    }
620}
621
622fn validate_positive(value: f64, field: &'static str) -> Result<()> {
623    validate_finite(value, field)?;
624    if value > 0.0 {
625        Ok(())
626    } else {
627        Err(invalid_input(field, "not positive"))
628    }
629}
630
631fn validate_finite(value: f64, field: &'static str) -> Result<()> {
632    if value.is_finite() {
633        Ok(())
634    } else {
635        Err(invalid_input(field, "not finite"))
636    }
637}
638
639fn invalid_input(field: &'static str, reason: &'static str) -> Error {
640    Error::InvalidInput(format!("{field} {reason}"))
641}
642
643#[cfg(test)]
644mod public_api_tests {
645    use super::*;
646
647    #[test]
648    fn relativistic_clock_correction_exposes_broadcast_formula() {
649        let dtr_f = ConstellationConstants::GPS.dtr_f;
650        let eccentricity = 0.013_456_789;
651        let sqrt_a = 5_153.795_477_5;
652        let sin_e = -0.625;
653        let got = relativistic_clock_correction_s(dtr_f, eccentricity, sqrt_a, sin_e)
654            .expect("valid relativistic correction");
655        let want = dtr_f * eccentricity * sqrt_a * sin_e;
656        assert_eq!(got.to_bits(), want.to_bits());
657    }
658
659    #[test]
660    fn relativistic_clock_correction_rejects_invalid_inputs() {
661        assert!(relativistic_clock_correction_s(f64::NAN, 0.01, 5_153.0, 0.5).is_err());
662        assert!(relativistic_clock_correction_s(
663            ConstellationConstants::GPS.dtr_f,
664            1.0,
665            5_153.0,
666            0.5
667        )
668        .is_err());
669        assert!(
670            relativistic_clock_correction_s(ConstellationConstants::GPS.dtr_f, 0.01, 0.0, 0.5)
671                .is_err()
672        );
673    }
674}
675
676#[cfg(all(test, sidereon_repo_tests))]
677mod tests;
sidereon_core/broadcast/mod.rs

sidereon_core/broadcast/
mod.rs
