sidereon_core/
combinations.rs

1//! GNSS observable linear combinations.
2//!
3//! The ionosphere-free code and carrier-phase combinations are pure
4//! frequency-domain algebra. The operation order is intentionally simple and
5//! pinned to the original Sidereon/SciPy oracle recipe: square each carrier first,
6//! form `gamma = f1^2 / (f1^2 - f2^2)`, then evaluate
7//! `gamma * obs1 - (gamma - 1) * obs2` with normal Rust arithmetic.
8
9use std::collections::{BTreeMap, BTreeSet};
10
11use crate::constants::C_M_S;
12use crate::frequencies;
13pub use crate::frequencies::{CarrierBand, CarrierFrequency, CarrierPair};
14use crate::validate;
15use crate::GnssSystem;
16
17/// Error produced by ionosphere-free combination helpers.
18#[derive(Debug, Clone, PartialEq, Eq)]
19pub enum IonosphereFreeError {
20    /// The constellation has no standard ionosphere-free carrier pair.
21    UnknownSystem(char),
22    /// The requested band is not known for the constellation.
23    UnknownBand { system: char, band: String },
24    /// Equal carrier frequencies make the denominator vanish.
25    EqualFrequencies,
26    /// Cycle-to-meter conversion requires positive carrier frequencies.
27    InvalidFrequency,
28    /// Observation values must be finite, and the combined result must remain finite.
29    InvalidObservation,
30}
31
32impl core::fmt::Display for IonosphereFreeError {
33    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
34        match self {
35            Self::UnknownSystem(system) => {
36                write!(f, "unknown ionosphere-free constellation {system}")
37            }
38            Self::UnknownBand { system, band } => {
39                write!(f, "unknown carrier band {band} for constellation {system}")
40            }
41            Self::EqualFrequencies => write!(f, "equal carrier frequencies"),
42            Self::InvalidFrequency => write!(f, "carrier frequencies must be positive"),
43            Self::InvalidObservation => write!(f, "observations must be finite"),
44        }
45    }
46}
47
48impl std::error::Error for IonosphereFreeError {}
49
50/// Reason a satellite was not included in a paired pseudorange result.
51#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
52pub enum PseudorangeDropReason {
53    /// Present in band 2 only.
54    MissingBand1,
55    /// Present in band 1 only.
56    MissingBand2,
57    /// The satellite appeared more than once within at least one band.
58    DuplicateObservation,
59    /// The satellite's constellation or requested band pair is unsupported.
60    UnknownSystem,
61}
62
63/// A satellite-tokened pseudorange observation in meters.
64pub type PseudorangeObservation = (String, f64);
65
66/// Ionosphere-free pseudoranges paired by satellite token.
67pub type CombinedPseudoranges = Vec<PseudorangeObservation>;
68
69/// Per-satellite reasons for dropped pseudorange combinations.
70pub type DroppedPseudoranges = Vec<(String, PseudorangeDropReason)>;
71
72/// Result returned by [`ionosphere_free_pseudoranges`].
73pub type PseudorangeCombinationResult =
74    Result<(CombinedPseudoranges, DroppedPseudoranges), IonosphereFreeError>;
75
76/// The full carrier-frequency table for the supported standard pairs.
77pub const fn carrier_frequencies() -> [CarrierFrequency; 6] {
78    frequencies::iono_free_carrier_frequencies()
79}
80
81/// The standard carrier frequency in hertz for a constellation and band.
82pub fn carrier_frequency_hz(system: GnssSystem, band: CarrierBand) -> Option<f64> {
83    frequencies::frequency_hz(system, band)
84}
85
86/// Carrier frequency lookup by RINEX/IGS system letter and lower-case band name.
87pub fn frequency_hz(system: char, band: &str) -> Result<f64, IonosphereFreeError> {
88    let Some(system_id) = GnssSystem::from_letter(system) else {
89        return Err(IonosphereFreeError::UnknownSystem(system));
90    };
91    let Some(carrier_band) = CarrierBand::from_iono_free_name(band) else {
92        return Err(IonosphereFreeError::UnknownBand {
93            system,
94            band: band.to_owned(),
95        });
96    };
97    carrier_frequency_hz(system_id, carrier_band).ok_or_else(|| IonosphereFreeError::UnknownBand {
98        system,
99        band: band.to_owned(),
100    })
101}
102
103/// Standard dual-frequency ionosphere-free carrier pair for a constellation.
104pub fn default_pair(system: char) -> Result<CarrierPair, IonosphereFreeError> {
105    match GnssSystem::from_letter(system).and_then(frequencies::default_iono_free_pair) {
106        Some(pair) => Ok(pair),
107        None => Err(IonosphereFreeError::UnknownSystem(system)),
108    }
109}
110
111/// Ionosphere-free coefficient `gamma = f1^2 / (f1^2 - f2^2)`.
112pub fn gamma(f1_hz: f64, f2_hz: f64) -> Result<f64, IonosphereFreeError> {
113    let (f1_hz, f2_hz) = validate_distinct_frequencies(f1_hz, f2_hz)?;
114    let f1sq = finite_frequency_product(f1_hz * f1_hz)?;
115    let f2sq = finite_frequency_product(f2_hz * f2_hz)?;
116    let denominator = finite_frequency_product(f1sq - f2sq)?;
117    if denominator == 0.0 {
118        return Err(IonosphereFreeError::EqualFrequencies);
119    }
120    let gamma = f1sq / denominator;
121    validate::finite(gamma, "gamma").map_err(|_| IonosphereFreeError::InvalidFrequency)
122}
123
124/// Equal-variance noise amplification of the ionosphere-free combination.
125pub fn noise_amplification(f1_hz: f64, f2_hz: f64) -> Result<f64, IonosphereFreeError> {
126    let g = gamma(f1_hz, f2_hz)?;
127    let amplification = (g * g + (g - 1.0) * (g - 1.0)).sqrt();
128    validate::finite(amplification, "noise_amplification")
129        .map_err(|_| IonosphereFreeError::InvalidFrequency)
130}
131
132/// Ionosphere-free code or meter-valued phase combination.
133pub fn ionosphere_free(
134    obs1_m: f64,
135    obs2_m: f64,
136    f1_hz: f64,
137    f2_hz: f64,
138) -> Result<f64, IonosphereFreeError> {
139    let (f1_hz, f2_hz) = validate_distinct_frequencies(f1_hz, f2_hz)?;
140    let obs1_m = validate_observation(obs1_m, "obs1_m")?;
141    let obs2_m = validate_observation(obs2_m, "obs2_m")?;
142    validate_observation(
143        ionosphere_free_unchecked(obs1_m, obs2_m, f1_hz, f2_hz),
144        "ionosphere_free_m",
145    )
146}
147
148/// Ionosphere-free carrier-phase combination from meter-valued phase inputs.
149pub fn ionosphere_free_phase_m(
150    phase1_m: f64,
151    phase2_m: f64,
152    f1_hz: f64,
153    f2_hz: f64,
154) -> Result<f64, IonosphereFreeError> {
155    ionosphere_free(phase1_m, phase2_m, f1_hz, f2_hz)
156}
157
158/// Ionosphere-free carrier-phase combination from cycle-valued phase inputs.
159pub fn ionosphere_free_phase_cycles(
160    phi1_cycles: f64,
161    phi2_cycles: f64,
162    f1_hz: f64,
163    f2_hz: f64,
164) -> Result<f64, IonosphereFreeError> {
165    let (f1_hz, f2_hz) = validate_distinct_frequencies(f1_hz, f2_hz)?;
166    let phi1_cycles = validate_observation(phi1_cycles, "phi1_cycles")?;
167    let phi2_cycles = validate_observation(phi2_cycles, "phi2_cycles")?;
168    let phase1_m = C_M_S / f1_hz * phi1_cycles;
169    let phase2_m = C_M_S / f2_hz * phi2_cycles;
170    let phase1_m = validate_observation(phase1_m, "phase1_m")?;
171    let phase2_m = validate_observation(phase2_m, "phase2_m")?;
172    validate_observation(
173        ionosphere_free_unchecked(phase1_m, phase2_m, f1_hz, f2_hz),
174        "ionosphere_free_phase_m",
175    )
176}
177
178fn validate_distinct_frequencies(
179    f1_hz: f64,
180    f2_hz: f64,
181) -> Result<(f64, f64), IonosphereFreeError> {
182    let f1_hz = validate_frequency(f1_hz, "f1_hz")?;
183    let f2_hz = validate_frequency(f2_hz, "f2_hz")?;
184    if f1_hz == f2_hz {
185        Err(IonosphereFreeError::EqualFrequencies)
186    } else {
187        Ok((f1_hz, f2_hz))
188    }
189}
190
191fn validate_frequency(f_hz: f64, field: &'static str) -> Result<f64, IonosphereFreeError> {
192    validate::finite_positive(f_hz, field).map_err(|_| IonosphereFreeError::InvalidFrequency)
193}
194
195fn finite_frequency_product(value: f64) -> Result<f64, IonosphereFreeError> {
196    validate::finite(value, "frequency_product").map_err(|_| IonosphereFreeError::InvalidFrequency)
197}
198
199fn validate_observation(value: f64, field: &'static str) -> Result<f64, IonosphereFreeError> {
200    validate::finite(value, field).map_err(|_| IonosphereFreeError::InvalidObservation)
201}
202
203fn gamma_unchecked(f1_hz: f64, f2_hz: f64) -> f64 {
204    let f1sq = f1_hz * f1_hz;
205    let f2sq = f2_hz * f2_hz;
206    f1sq / (f1sq - f2sq)
207}
208
209fn ionosphere_free_unchecked(obs1_m: f64, obs2_m: f64, f1_hz: f64, f2_hz: f64) -> f64 {
210    let g = gamma_unchecked(f1_hz, f2_hz);
211    g * obs1_m - (g - 1.0) * obs2_m
212}
213
214/// Combine two satellite-keyed pseudorange bands into ionosphere-free ranges.
215///
216/// `overrides` is a list of `(system_letter, band1_name, band2_name)` entries.
217/// An invalid override band is treated the same way the original Sidereon wrapper
218/// treated a failed combination for a paired satellite: that satellite is
219/// reported as [`PseudorangeDropReason::UnknownSystem`].
220pub fn ionosphere_free_pseudoranges(
221    band1: &[PseudorangeObservation],
222    band2: &[PseudorangeObservation],
223    overrides: &[(char, String, String)],
224) -> PseudorangeCombinationResult {
225    let (m1, dups1) = map_and_duplicates(band1)?;
226    let (m2, dups2) = map_and_duplicates(band2)?;
227
228    let dups = dups1.union(&dups2).cloned().collect::<BTreeSet<_>>();
229    let ids1 = m1.keys().cloned().collect::<BTreeSet<_>>();
230    let ids2 = m2.keys().cloned().collect::<BTreeSet<_>>();
231
232    let mut combined = Vec::new();
233    let mut dropped = Vec::new();
234
235    for sat in ids1.intersection(&ids2) {
236        if dups.contains(sat) {
237            continue;
238        }
239        match combine_satellite(sat, m1[sat], m2[sat], overrides) {
240            Ok(range_m) => combined.push((sat.clone(), range_m)),
241            Err(IonosphereFreeError::UnknownSystem(_))
242            | Err(IonosphereFreeError::UnknownBand { .. }) => {
243                dropped.push((sat.clone(), PseudorangeDropReason::UnknownSystem))
244            }
245            Err(error) => return Err(error),
246        }
247    }
248
249    for sat in ids1.difference(&ids2) {
250        if !dups.contains(sat) {
251            dropped.push((sat.clone(), PseudorangeDropReason::MissingBand2));
252        }
253    }
254
255    for sat in ids2.difference(&ids1) {
256        if !dups.contains(sat) {
257            dropped.push((sat.clone(), PseudorangeDropReason::MissingBand1));
258        }
259    }
260
261    for sat in dups {
262        dropped.push((sat, PseudorangeDropReason::DuplicateObservation));
263    }
264
265    dropped.sort();
266    Ok((combined, dropped))
267}
268
269fn map_and_duplicates(
270    observations: &[(String, f64)],
271) -> Result<(BTreeMap<String, f64>, BTreeSet<String>), IonosphereFreeError> {
272    let mut counts = BTreeMap::<String, usize>::new();
273    let mut map = BTreeMap::<String, f64>::new();
274    for (sat, range_m) in observations {
275        let range_m = validate_observation(*range_m, "pseudorange_m")?;
276        *counts.entry(sat.clone()).or_insert(0) += 1;
277        map.insert(sat.clone(), range_m);
278    }
279    let dups = counts
280        .into_iter()
281        .filter_map(|(sat, count)| (count > 1).then_some(sat))
282        .collect();
283    Ok((map, dups))
284}
285
286fn combine_satellite(
287    sat: &str,
288    pr1_m: f64,
289    pr2_m: f64,
290    overrides: &[(char, String, String)],
291) -> Result<f64, IonosphereFreeError> {
292    let system = sat
293        .chars()
294        .next()
295        .ok_or(IonosphereFreeError::UnknownSystem('\0'))?;
296    let pair = pair_for(system, overrides)?;
297    let f1 = frequency_hz(system, pair.band1.name())?;
298    let f2 = frequency_hz(system, pair.band2.name())?;
299    ionosphere_free(pr1_m, pr2_m, f1, f2)
300}
301
302fn pair_for(
303    system: char,
304    overrides: &[(char, String, String)],
305) -> Result<CarrierPair, IonosphereFreeError> {
306    if let Some((_system, band1, band2)) = overrides.iter().find(|(s, _, _)| *s == system) {
307        let Some(band1) = CarrierBand::from_iono_free_name(band1) else {
308            return Err(IonosphereFreeError::UnknownBand {
309                system,
310                band: band1.clone(),
311            });
312        };
313        let Some(band2) = CarrierBand::from_iono_free_name(band2) else {
314            return Err(IonosphereFreeError::UnknownBand {
315                system,
316                band: band2.clone(),
317            });
318        };
319        Ok(CarrierPair::new(band1, band2))
320    } else {
321        default_pair(system)
322    }
323}
324
325#[cfg(test)]
326mod tests {
327    use super::*;
328    use crate::constants::{F_B1I_HZ, F_B3I_HZ, F_E1_HZ, F_E5A_HZ, F_L1_HZ, F_L2_HZ};
329
330    struct OracleCase {
331        f1_hz: u64,
332        f2_hz: u64,
333        pr1_m: u64,
334        pr2_m: u64,
335        gamma: u64,
336        noise: u64,
337        iono_free_m: u64,
338        phi1_cycles: u64,
339        phi2_cycles: u64,
340        phase1_m: u64,
341        phase2_m: u64,
342        phase_if_m: u64,
343        phase_if_cycles_m: u64,
344    }
345
346    fn f(bits: u64) -> f64 {
347        f64::from_bits(bits)
348    }
349
350    #[test]
351    fn frequency_table_matches_standard_carriers() {
352        assert_eq!(frequency_hz('G', "l1"), Ok(F_L1_HZ));
353        assert_eq!(frequency_hz('G', "l2"), Ok(F_L2_HZ));
354        assert_eq!(frequency_hz('E', "e1"), Ok(F_E1_HZ));
355        assert_eq!(frequency_hz('E', "e5a"), Ok(F_E5A_HZ));
356        assert_eq!(frequency_hz('C', "b1i"), Ok(F_B1I_HZ));
357        assert_eq!(frequency_hz('C', "b3i"), Ok(F_B3I_HZ));
358        assert_eq!(carrier_frequencies().len(), 6);
359    }
360
361    #[test]
362    fn frequency_lookup_classifies_system_and_band_errors() {
363        assert_eq!(
364            frequency_hz('X', "l1"),
365            Err(IonosphereFreeError::UnknownSystem('X'))
366        );
367        assert_eq!(
368            frequency_hz('G', "bad"),
369            Err(IonosphereFreeError::UnknownBand {
370                system: 'G',
371                band: "bad".to_owned(),
372            })
373        );
374        assert_eq!(frequency_hz('G', "l1"), Ok(F_L1_HZ));
375    }
376
377    #[test]
378    fn default_pairs_match_supported_systems() {
379        assert_eq!(
380            default_pair('G'),
381            Ok(CarrierPair::new(CarrierBand::L1, CarrierBand::L2))
382        );
383        assert_eq!(
384            default_pair('E'),
385            Ok(CarrierPair::new(CarrierBand::E1, CarrierBand::E5a))
386        );
387        assert_eq!(
388            default_pair('C'),
389            Ok(CarrierPair::new(CarrierBand::B1i, CarrierBand::B3i))
390        );
391        assert_eq!(
392            default_pair('X'),
393            Err(IonosphereFreeError::UnknownSystem('X'))
394        );
395    }
396
397    #[test]
398    fn scipy_oracle_cases_are_zero_ulp() {
399        let cases = [
400            OracleCase {
401                f1_hz: 0x41d779c018000000,
402                f2_hz: 0x41d24aec20000000,
403                pr1_m: 0x4175ef3c40772a36,
404                pr2_m: 0x4175ef3c6a2bcbb5,
405                gamma: 0x40045da686c28e3c,
406                noise: 0x4007d3777c503ebc,
407                iono_free_m: 0x4175ef3c00000000,
408                phi1_cycles: 0x419cd8990a6a993b,
409                phi2_cycles: 0x419682ad3bea73b9,
410                phase1_m: 0x4175f4f80ddd7ecd,
411                phase2_m: 0x4175fd37d057d184,
412                phase_if_m: 0x4175e837d93b3cba,
413                phase_if_cycles_m: 0x4175e837d93b3cba,
414            },
415            OracleCase {
416                f1_hz: 0x41d779c018000000,
417                f2_hz: 0x41d187ccf4000000,
418                pr1_m: 0x41775d7280ee546c,
419                pr2_m: 0x41775d72e7354588,
420                gamma: 0x400215b7b8bf1d8d,
421                noise: 0x4004b4e6a9e28198,
422                iono_free_m: 0x41775d71ffffffff,
423                phi1_cycles: 0x419eb9b5c28924dd,
424                phi2_cycles: 0x4196fa6edb5f553e,
425                phase1_m: 0x4177632debd8c260,
426                phase2_m: 0x41776c09259441ae,
427                phase_if_m: 0x41775803d8d388ae,
428                phase_if_cycles_m: 0x41775803d8d388ae,
429            },
430            OracleCase {
431                f1_hz: 0x41d7431dc4000000,
432                f2_hz: 0x41d2e70510000000,
433                pr1_m: 0x41753821627b0be3,
434                pr2_m: 0x417538219525d0a5,
435                gamma: 0x40078ca90724ddf1,
436                noise: 0x400c384adb005afd,
437                iono_free_m: 0x4175382100000000,
438                phi1_cycles: 0x419ba72a23131776,
439                phi2_cycles: 0x419680982c673023,
440                phase1_m: 0x41753deaa1cd1743,
441                phase2_m: 0x417545a9733bf939,
442                phase_if_m: 0x41752edcaa262834,
443                phase_if_cycles_m: 0x41752edcaa262834,
444            },
445        ];
446
447        for case in cases {
448            let f1 = f(case.f1_hz);
449            let f2 = f(case.f2_hz);
450            assert_eq!(gamma(f1, f2).unwrap().to_bits(), case.gamma);
451            assert_eq!(noise_amplification(f1, f2).unwrap().to_bits(), case.noise);
452            assert_eq!(
453                ionosphere_free(f(case.pr1_m), f(case.pr2_m), f1, f2)
454                    .unwrap()
455                    .to_bits(),
456                case.iono_free_m
457            );
458            assert_eq!(
459                ionosphere_free_phase_m(f(case.phase1_m), f(case.phase2_m), f1, f2)
460                    .unwrap()
461                    .to_bits(),
462                case.phase_if_m
463            );
464            assert_eq!(
465                ionosphere_free_phase_cycles(f(case.phi1_cycles), f(case.phi2_cycles), f1, f2)
466                    .unwrap()
467                    .to_bits(),
468                case.phase_if_cycles_m
469            );
470        }
471    }
472
473    #[test]
474    fn invalid_frequency_modes_are_tagged() {
475        assert_eq!(
476            gamma(F_L1_HZ, F_L1_HZ),
477            Err(IonosphereFreeError::EqualFrequencies)
478        );
479        assert_eq!(
480            gamma(f64::NAN, F_L2_HZ),
481            Err(IonosphereFreeError::InvalidFrequency)
482        );
483        assert_eq!(
484            gamma(f64::MAX, F_L2_HZ),
485            Err(IonosphereFreeError::InvalidFrequency)
486        );
487        assert_eq!(
488            noise_amplification(F_L1_HZ, f64::INFINITY),
489            Err(IonosphereFreeError::InvalidFrequency)
490        );
491        assert_eq!(
492            ionosphere_free(1.0, 2.0, F_L1_HZ, 0.0),
493            Err(IonosphereFreeError::InvalidFrequency)
494        );
495        assert_eq!(
496            ionosphere_free_phase_cycles(1.0, 2.0, F_L1_HZ, F_L1_HZ),
497            Err(IonosphereFreeError::EqualFrequencies)
498        );
499        assert_eq!(
500            ionosphere_free_phase_cycles(1.0, 2.0, 0.0, F_L2_HZ),
501            Err(IonosphereFreeError::InvalidFrequency)
502        );
503        assert_eq!(
504            ionosphere_free_phase_cycles(1.0, 2.0, -F_L1_HZ, F_L2_HZ),
505            Err(IonosphereFreeError::InvalidFrequency)
506        );
507    }
508
509    #[test]
510    fn invalid_observations_are_rejected() {
511        assert_eq!(
512            ionosphere_free(f64::NAN, 2.0, F_L1_HZ, F_L2_HZ),
513            Err(IonosphereFreeError::InvalidObservation)
514        );
515        assert_eq!(
516            ionosphere_free_phase_m(1.0, f64::INFINITY, F_L1_HZ, F_L2_HZ),
517            Err(IonosphereFreeError::InvalidObservation)
518        );
519        assert_eq!(
520            ionosphere_free_phase_cycles(f64::NAN, 2.0, F_L1_HZ, F_L2_HZ),
521            Err(IonosphereFreeError::InvalidObservation)
522        );
523        assert_eq!(
524            ionosphere_free(f64::MAX, -f64::MAX, F_L1_HZ, F_L2_HZ),
525            Err(IonosphereFreeError::InvalidObservation)
526        );
527    }
528
529    #[test]
530    fn pseudorange_pairing_reports_missing_unknown_and_duplicate_satellites() {
531        let band1 = vec![
532            ("G01".to_string(), 23_000_000.0),
533            ("G01".to_string(), 23_000_010.0),
534            ("G02".to_string(), 22_000_000.0),
535            ("G03".to_string(), 21_000_000.0),
536            ("X01".to_string(), 20_000_000.0),
537        ];
538        let band2 = vec![
539            ("G01".to_string(), 23_000_000.0),
540            ("G02".to_string(), 22_000_000.0),
541            ("G04".to_string(), 24_000_000.0),
542            ("X01".to_string(), 20_000_000.0),
543        ];
544
545        let (combined, dropped) =
546            ionosphere_free_pseudoranges(&band1, &band2, &[]).expect("valid pseudoranges");
547
548        assert_eq!(combined.len(), 1);
549        assert_eq!(combined[0].0, "G02");
550        assert_eq!(
551            dropped,
552            vec![
553                (
554                    "G01".to_string(),
555                    PseudorangeDropReason::DuplicateObservation
556                ),
557                ("G03".to_string(), PseudorangeDropReason::MissingBand2),
558                ("G04".to_string(), PseudorangeDropReason::MissingBand1),
559                ("X01".to_string(), PseudorangeDropReason::UnknownSystem),
560            ]
561        );
562    }
563
564    #[test]
565    fn pseudorange_pairing_rejects_non_finite_values() {
566        let band1 = vec![("G01".to_string(), f64::NAN)];
567        let band2 = vec![("G01".to_string(), 23_000_000.0)];
568
569        assert_eq!(
570            ionosphere_free_pseudoranges(&band1, &band2, &[]),
571            Err(IonosphereFreeError::InvalidObservation)
572        );
573    }
574}
sidereon_core/combinations.rs

sidereon_core/
combinations.rs
