sidereon_core/
sidereal.rs

1//! Repeating-geometry residual filtering.
2//!
3//! This module is a pure post-processing kernel for residual arrays. It models
4//! local-environment repeat signatures as a phase-stacked template at a supplied
5//! repeat period, then subtracts only corrections backed by the requested prior
6//! coverage.
7
8use std::f64::consts::PI;
9
10use crate::astro::time::Duration;
11use crate::constants::SECONDS_PER_DAY;
12use crate::ephemeris::BroadcastEphemeris;
13use crate::estimation::{ewma_update, mad_spread, PrimitiveError};
14use crate::{GnssSatelliteId, GnssSystem};
15
16/// Length of the mean sidereal day used for constellation-default GPS phasing.
17///
18/// The value is 23 h 56 min 4.0905 s, stored as an integer nanosecond duration
19/// so [`repeat_period`] does not depend on floating-point rounding.
20pub const SIDEREAL_DAY_NANOS: i128 = 86_164_090_500_000;
21
22/// Length of the mean sidereal day in SI seconds.
23pub const SIDEREAL_DAY_SECONDS: f64 = 86_164.090_5;
24
25/// Errors returned by sidereal filtering and repeat-period helpers.
26#[derive(Debug, Clone, PartialEq, thiserror::Error)]
27pub enum SiderealFilterError {
28    /// A public input was not usable.
29    #[error("invalid sidereal filter {field}: {reason}")]
30    InvalidInput {
31        /// Field name for the failing input.
32        field: &'static str,
33        /// Human-readable reason for the failure.
34        reason: &'static str,
35    },
36    /// No broadcast record covered the requested satellite and epoch.
37    #[error("no broadcast record for {sat} at requested epoch")]
38    NoBroadcastRecord {
39        /// Satellite without a selected broadcast record.
40        sat: GnssSatelliteId,
41    },
42    /// The requested constellation has no Keplerian broadcast-repeat path.
43    #[error("unsupported orbit repeat constellation {system}")]
44    UnsupportedConstellation {
45        /// Unsupported constellation.
46        system: GnssSystem,
47    },
48}
49
50/// Template estimator used by [`sidereal_filter`].
51#[derive(Debug, Clone, Copy, PartialEq, Default)]
52pub enum SiderealTemplateMethod {
53    /// Arithmetic mean of the prior values in the same phase bin.
54    #[default]
55    Mean,
56    /// MAD-gated mean of the prior values in the same phase bin.
57    RobustMad,
58    /// Exponentially weighted mean of prior values in the same phase bin.
59    Ewma {
60        /// EWMA gain in `[0, 1]`.
61        alpha: f64,
62    },
63}
64
65/// Options controlling residual template stacking.
66#[derive(Debug, Clone, Copy, PartialEq)]
67pub struct SiderealFilterOptions {
68    /// Sampling interval of the residual array.
69    pub sample_interval: Duration,
70    /// Maximum number of prior repeats retained per phase bin.
71    pub prior_periods: usize,
72    /// Minimum prior samples required before a correction is applied.
73    pub min_coverage: usize,
74    /// Template estimator applied to prior samples.
75    pub template_method: SiderealTemplateMethod,
76}
77
78impl Default for SiderealFilterOptions {
79    fn default() -> Self {
80        Self {
81            sample_interval: Duration::from_nanos(1_000_000_000),
82            prior_periods: 1,
83            min_coverage: 1,
84            template_method: SiderealTemplateMethod::Mean,
85        }
86    }
87}
88
89/// Output of [`sidereal_filter`].
90#[derive(Debug, Clone, PartialEq)]
91pub struct SiderealFilterOutput {
92    /// Residuals after subtracting covered repeat-template corrections.
93    pub filtered: Vec<f64>,
94    /// Last correction template applied for each phase bin.
95    ///
96    /// Bins that never reached [`SiderealFilterOptions::min_coverage`] are
97    /// stored as `NaN` and marked in [`Self::under_covered`].
98    pub template: Vec<f64>,
99    /// Prior-sample count backing the last template decision for each phase bin.
100    pub coverage: Vec<usize>,
101    /// `true` for phase bins whose last decision had insufficient coverage.
102    pub under_covered: Vec<bool>,
103}
104
105/// Return the default ground-track repeat period for a GNSS constellation.
106///
107/// GPS, QZSS, NavIC, and SBAS default to one mean sidereal day. GLONASS,
108/// Galileo, and BeiDou use the design repeat cycles of 8, 10, and 7 sidereal
109/// days respectively. Use [`orbit_repeat_lag`] when a GPS-like per-satellite
110/// broadcast-orbit refinement is available.
111pub fn repeat_period(system: GnssSystem) -> Duration {
112    let sidereal_days = match system {
113        GnssSystem::Gps | GnssSystem::Qzss | GnssSystem::Navic | GnssSystem::Sbas => 1,
114        GnssSystem::Glonass => 8,
115        GnssSystem::Galileo => 10,
116        GnssSystem::BeiDou => 7,
117    };
118    Duration::from_nanos(SIDEREAL_DAY_NANOS * sidereal_days)
119}
120
121/// Compute a per-satellite orbital repeat lag from a broadcast ephemeris.
122///
123/// The selected broadcast record supplies the square root of semi-major axis
124/// and the mean-motion correction. The returned lag is `4*pi/n`, where
125/// `n = sqrt(GM / a^3) + delta_n`, matching the orbit-repeat construction used
126/// for modified sidereal filtering.
127pub fn orbit_repeat_lag(
128    ephemeris: &BroadcastEphemeris,
129    sat: GnssSatelliteId,
130    near_epoch_j2000_s: f64,
131) -> Result<Duration, SiderealFilterError> {
132    validate_finite(near_epoch_j2000_s, "near_epoch_j2000_s")?;
133    if !matches!(
134        sat.system,
135        GnssSystem::Gps | GnssSystem::Galileo | GnssSystem::BeiDou | GnssSystem::Qzss
136    ) {
137        return Err(SiderealFilterError::UnsupportedConstellation { system: sat.system });
138    }
139    let record = ephemeris
140        .select_record_at(sat, near_epoch_j2000_s)
141        .ok_or(SiderealFilterError::NoBroadcastRecord { sat })?;
142    let repeat_s = orbit_repeat_seconds(
143        record.elements.sqrt_a,
144        record.elements.delta_n,
145        record.constants().gm_m3_s2,
146    )?;
147    Duration::from_seconds(repeat_s).map_err(|_| invalid_input("orbit_repeat_lag", "out of range"))
148}
149
150/// Remove a repeating residual component by phase-stacking prior repeats.
151///
152/// Samples are assigned to phase bins by `(index * sample_interval) mod period`.
153/// A correction is subtracted only after a phase bin has at least
154/// [`SiderealFilterOptions::min_coverage`] prior samples. Under-covered samples
155/// are returned unchanged and the corresponding phase bins are flagged.
156pub fn sidereal_filter(
157    series: &[f64],
158    period: Duration,
159    options: SiderealFilterOptions,
160) -> Result<SiderealFilterOutput, SiderealFilterError> {
161    validate_series(series)?;
162    validate_options(options)?;
163    let period_s = validate_positive_duration(period, "period")?;
164    let sample_interval_s =
165        validate_positive_duration(options.sample_interval, "options.sample_interval")?;
166    let bins = phase_bin_count(period_s, sample_interval_s)?;
167
168    let mut histories = vec![Vec::<f64>::new(); bins];
169    let mut filtered = Vec::with_capacity(series.len());
170    let mut template = vec![f64::NAN; bins];
171    let mut coverage = vec![0usize; bins];
172    let mut under_covered = vec![true; bins];
173
174    for (sample_index, &sample) in series.iter().enumerate() {
175        let bin = phase_bin(sample_index, period_s, sample_interval_s, bins);
176        let history = &histories[bin];
177        let prior_count = history.len();
178        coverage[bin] = prior_count;
179        under_covered[bin] = prior_count < options.min_coverage;
180
181        if prior_count >= options.min_coverage {
182            let correction = estimate_template(history, options.template_method)?;
183            template[bin] = correction;
184            filtered.push(sample - correction);
185        } else {
186            filtered.push(sample);
187        }
188
189        let history = &mut histories[bin];
190        history.push(sample);
191        while history.len() > options.prior_periods {
192            history.remove(0);
193        }
194    }
195
196    Ok(SiderealFilterOutput {
197        filtered,
198        template,
199        coverage,
200        under_covered,
201    })
202}
203
204/// Score repeating components at candidate periods for 1 Hz samples.
205///
206/// The returned strength is a robust variance-reduction score in `[0, 1]`,
207/// where `1` means the period template removes all robust spread and `0` means
208/// no robust spread reduction was measured. Use
209/// [`periodicity_strength_with_sample_interval`] for other sample cadences.
210pub fn periodicity_strength(
211    series: &[f64],
212    candidate_periods: &[Duration],
213) -> Result<Vec<(Duration, f64)>, SiderealFilterError> {
214    periodicity_strength_with_sample_interval(
215        series,
216        candidate_periods,
217        Duration::from_nanos(1_000_000_000),
218    )
219}
220
221/// Score repeating components at candidate periods for an explicit cadence.
222///
223/// Each candidate period is converted to phase bins at `sample_interval`.
224/// Bins with fewer than two samples do not contribute a confident template.
225/// The score uses [`mad_spread`] for both the input and period-template
226/// residual spread.
227pub fn periodicity_strength_with_sample_interval(
228    series: &[f64],
229    candidate_periods: &[Duration],
230    sample_interval: Duration,
231) -> Result<Vec<(Duration, f64)>, SiderealFilterError> {
232    validate_series(series)?;
233    if series.len() < 2 {
234        return Err(invalid_input("series", "must contain at least two samples"));
235    }
236    let sample_interval_s = validate_positive_duration(sample_interval, "sample_interval")?;
237    let baseline = mad_spread(series, 0.0).map_err(map_primitive_error)?;
238    let mut scores = Vec::with_capacity(candidate_periods.len());
239
240    for &period in candidate_periods {
241        let period_s = validate_positive_duration(period, "candidate_periods")?;
242        let bins = phase_bin_count(period_s, sample_interval_s)?;
243        let strength = if baseline == 0.0 {
244            0.0
245        } else {
246            periodicity_strength_one(series, period_s, sample_interval_s, bins, baseline)?
247        };
248        scores.push((period, strength));
249    }
250
251    Ok(scores)
252}
253
254fn periodicity_strength_one(
255    series: &[f64],
256    period_s: f64,
257    sample_interval_s: f64,
258    bins: usize,
259    baseline: f64,
260) -> Result<f64, SiderealFilterError> {
261    let mut phases = vec![Vec::<f64>::new(); bins];
262    for (sample_index, &sample) in series.iter().enumerate() {
263        let bin = phase_bin(sample_index, period_s, sample_interval_s, bins);
264        phases[bin].push(sample);
265    }
266
267    let mut templates = vec![f64::NAN; bins];
268    let mut covered_bins = 0usize;
269    for (bin, values) in phases.iter().enumerate() {
270        if values.len() >= 2 {
271            templates[bin] = robust_mad_template(values)?;
272            covered_bins += 1;
273        }
274    }
275    if covered_bins == 0 {
276        return Ok(0.0);
277    }
278
279    let residuals = series
280        .iter()
281        .enumerate()
282        .map(|(sample_index, &sample)| {
283            let bin = phase_bin(sample_index, period_s, sample_interval_s, bins);
284            let correction = templates[bin];
285            if correction.is_finite() {
286                sample - correction
287            } else {
288                sample
289            }
290        })
291        .collect::<Vec<_>>();
292    let residual_spread = mad_spread(&residuals, 0.0).map_err(map_primitive_error)?;
293    let raw = 1.0 - (residual_spread / baseline).powi(2);
294    Ok(raw.clamp(0.0, 1.0))
295}
296
297fn orbit_repeat_seconds(
298    sqrt_a: f64,
299    delta_n: f64,
300    gm_m3_s2: f64,
301) -> Result<f64, SiderealFilterError> {
302    validate_positive(sqrt_a, "record.elements.sqrt_a")?;
303    validate_finite(delta_n, "record.elements.delta_n")?;
304    validate_positive(gm_m3_s2, "record.constants.gm_m3_s2")?;
305    let a = sqrt_a * sqrt_a;
306    let n0 = (gm_m3_s2 / (a * a * a)).sqrt();
307    let n = n0 + delta_n;
308    validate_positive(n, "record.mean_motion")?;
309    let repeat_s = 4.0 * PI / n;
310    validate_positive(repeat_s, "orbit_repeat_lag")
311}
312
313fn validate_options(options: SiderealFilterOptions) -> Result<(), SiderealFilterError> {
314    if options.prior_periods == 0 {
315        return Err(invalid_input("options.prior_periods", "must be positive"));
316    }
317    if options.min_coverage == 0 {
318        return Err(invalid_input("options.min_coverage", "must be positive"));
319    }
320    match options.template_method {
321        SiderealTemplateMethod::Mean | SiderealTemplateMethod::RobustMad => {}
322        SiderealTemplateMethod::Ewma { alpha } => {
323            validate_finite(alpha, "options.template_method.alpha")?;
324            if !(0.0..=1.0).contains(&alpha) {
325                return Err(invalid_input(
326                    "options.template_method.alpha",
327                    "must be in [0, 1]",
328                ));
329            }
330        }
331    }
332    Ok(())
333}
334
335fn validate_series(series: &[f64]) -> Result<(), SiderealFilterError> {
336    for &sample in series {
337        validate_finite(sample, "series")?;
338    }
339    Ok(())
340}
341
342fn validate_positive_duration(
343    duration: Duration,
344    field: &'static str,
345) -> Result<f64, SiderealFilterError> {
346    if duration.nanos <= 0 {
347        return Err(invalid_input(field, "must be positive"));
348    }
349    let seconds = duration.as_seconds();
350    validate_positive(seconds, field)
351}
352
353fn validate_finite(value: f64, field: &'static str) -> Result<f64, SiderealFilterError> {
354    if value.is_finite() {
355        Ok(value)
356    } else {
357        Err(invalid_input(field, "must be finite"))
358    }
359}
360
361fn validate_positive(value: f64, field: &'static str) -> Result<f64, SiderealFilterError> {
362    validate_finite(value, field)?;
363    if value > 0.0 {
364        Ok(value)
365    } else {
366        Err(invalid_input(field, "must be positive"))
367    }
368}
369
370fn phase_bin_count(period_s: f64, sample_interval_s: f64) -> Result<usize, SiderealFilterError> {
371    let bins = (period_s / sample_interval_s).ceil();
372    if !bins.is_finite() || bins < 1.0 || bins > usize::MAX as f64 {
373        return Err(invalid_input("period", "phase-bin count is out of range"));
374    }
375    Ok(bins as usize)
376}
377
378fn phase_bin(sample_index: usize, period_s: f64, sample_interval_s: f64, bins: usize) -> usize {
379    let phase_s = (sample_index as f64 * sample_interval_s).rem_euclid(period_s);
380    let bin = (phase_s / sample_interval_s).floor() as usize;
381    bin.min(bins - 1)
382}
383
384fn estimate_template(
385    values: &[f64],
386    method: SiderealTemplateMethod,
387) -> Result<f64, SiderealFilterError> {
388    match method {
389        SiderealTemplateMethod::Mean => mean(values),
390        SiderealTemplateMethod::RobustMad => robust_mad_template(values),
391        SiderealTemplateMethod::Ewma { alpha } => ewma_template(values, alpha),
392    }
393}
394
395fn mean(values: &[f64]) -> Result<f64, SiderealFilterError> {
396    if values.is_empty() {
397        return Err(invalid_input("values", "must not be empty"));
398    }
399    Ok(values.iter().sum::<f64>() / values.len() as f64)
400}
401
402fn ewma_template(values: &[f64], alpha: f64) -> Result<f64, SiderealFilterError> {
403    let Some((&first, rest)) = values.split_first() else {
404        return Err(invalid_input("values", "must not be empty"));
405    };
406    let mut template = first;
407    for &sample in rest {
408        template = ewma_update(template, sample, alpha).map_err(map_primitive_error)?;
409    }
410    Ok(template)
411}
412
413fn robust_mad_template(values: &[f64]) -> Result<f64, SiderealFilterError> {
414    let median = median(values)?;
415    let spread = mad_spread(values, 0.0).map_err(map_primitive_error)?;
416    if spread == 0.0 {
417        return Ok(median);
418    }
419    let gate = 3.0 * spread;
420    let mut sum = 0.0;
421    let mut count = 0usize;
422    for &value in values {
423        if (value - median).abs() <= gate {
424            sum += value;
425            count += 1;
426        }
427    }
428    if count == 0 {
429        Ok(median)
430    } else {
431        Ok(sum / count as f64)
432    }
433}
434
435fn median(values: &[f64]) -> Result<f64, SiderealFilterError> {
436    if values.is_empty() {
437        return Err(invalid_input("values", "must not be empty"));
438    }
439    let mut sorted = values.to_vec();
440    sorted.sort_by(|a, b| a.total_cmp(b));
441    let mid = sorted.len() / 2;
442    if sorted.len() % 2 == 1 {
443        Ok(sorted[mid])
444    } else {
445        Ok((sorted[mid - 1] + sorted[mid]) * 0.5)
446    }
447}
448
449fn map_primitive_error(error: PrimitiveError) -> SiderealFilterError {
450    match error {
451        PrimitiveError::InvalidInput { field, reason } => {
452            SiderealFilterError::InvalidInput { field, reason }
453        }
454    }
455}
456
457const fn invalid_input(field: &'static str, reason: &'static str) -> SiderealFilterError {
458    SiderealFilterError::InvalidInput { field, reason }
459}
460
461/// Return the 24-hour solar day as a duration for diagnostic callers.
462pub fn solar_day_period() -> Duration {
463    Duration::from_nanos((SECONDS_PER_DAY as i128) * 1_000_000_000)
464}
465
466#[cfg(test)]
467mod tests {
468    //! Validation provenance:
469    //!
470    //! - Sidereal day: the mean sidereal day constant is checked at the exact
471    //!   integer nanosecond value for 23 h 56 min 4.0905 s.
472    //! - Agnew and Larson, "Finding the Repeat Times of the GPS Constellation":
473    //!   the orbit-repeat formula `n = sqrt(GM / a^3) + delta_n`,
474    //!   `To = 4*pi/n` is checked against the Table 1 SV 1 orbit repeat time
475    //!   of 86152.95 s.
476    //! - Synthetic filtering tests use exact binary residual templates so the
477    //!   recovered template and current-cycle noise floor are bit-exact.
478
479    use super::*;
480    use crate::astro::time::GnssWeekTow;
481    use crate::broadcast::{ClockPolynomial, ConstellationConstants, KeplerianElements};
482    use crate::constants::{GPS_EPOCH_TO_J2000_S, SECONDS_PER_HOUR, SECONDS_PER_WEEK};
483    use crate::rinex_nav::{BroadcastGroupDelays, BroadcastIssue, BroadcastRecord, NavMessage};
484
485    #[test]
486    fn repeat_periods_use_documented_sidereal_day_multiples() {
487        assert_eq!(
488            repeat_period(GnssSystem::Gps),
489            Duration::from_nanos(SIDEREAL_DAY_NANOS)
490        );
491        assert_eq!(
492            repeat_period(GnssSystem::Glonass),
493            Duration::from_nanos(SIDEREAL_DAY_NANOS * 8)
494        );
495        assert_eq!(
496            repeat_period(GnssSystem::Galileo),
497            Duration::from_nanos(SIDEREAL_DAY_NANOS * 10)
498        );
499        assert_eq!(
500            repeat_period(GnssSystem::BeiDou),
501            Duration::from_nanos(SIDEREAL_DAY_NANOS * 7)
502        );
503    }
504
505    #[test]
506    fn orbit_repeat_lag_matches_agnew_larson_sv1_table_value() {
507        let sat = GnssSatelliteId::new(GnssSystem::Gps, 1).expect("valid satellite");
508        let target_repeat_s = 86_152.95;
509        let record = agnew_larson_table_record(sat, target_repeat_s);
510        let week = record.toe.week;
511        let toe = record.toe.tow_s;
512        let store = BroadcastEphemeris::new(vec![record]).expect("valid broadcast store");
513        let near_epoch_j2000_s = f64::from(week) * SECONDS_PER_WEEK + toe - GPS_EPOCH_TO_J2000_S;
514
515        let lag = orbit_repeat_lag(&store, sat, near_epoch_j2000_s)
516            .expect("orbit repeat lag from selected record");
517        let repeat_s = lag.as_seconds();
518
519        assert!(
520            (repeat_s - target_repeat_s).abs() <= 1.0e-9,
521            "repeat_s={repeat_s:.12}"
522        );
523    }
524
525    #[test]
526    fn sidereal_filter_recovers_exact_template_and_noise_floor() {
527        let injected = [1.0, 0.5, -0.25, -1.25, 2.0, -0.75, 0.125, 0.875];
528        let noise = [0.125, -0.25, 0.0625, 0.0, -0.125, 0.25, -0.0625, 0.1875];
529        let mut series = injected.to_vec();
530        series.extend(
531            injected
532                .iter()
533                .zip(noise)
534                .map(|(&signal, noise)| signal + noise),
535        );
536
537        let output = sidereal_filter(
538            &series,
539            Duration::from_nanos(8_000_000_000),
540            SiderealFilterOptions::default(),
541        )
542        .expect("sidereal filter");
543
544        assert_eq!(&output.template, &injected);
545        assert_eq!(&output.filtered[8..], &noise);
546        assert_eq!(variance(&output.filtered[8..]), variance(&noise));
547        assert_eq!(output.coverage, vec![1; injected.len()]);
548        assert_eq!(output.under_covered, vec![false; injected.len()]);
549    }
550
551    #[test]
552    fn short_or_undercovered_series_is_flagged_and_left_unchanged() {
553        let series = [0.25, -0.5, 0.75, -1.0, 1.25];
554        let output = sidereal_filter(
555            &series,
556            Duration::from_nanos(4_000_000_000),
557            SiderealFilterOptions {
558                min_coverage: 2,
559                ..SiderealFilterOptions::default()
560            },
561        )
562        .expect("sidereal filter");
563
564        assert_eq!(output.filtered, series);
565        assert_eq!(output.coverage, vec![1, 0, 0, 0]);
566        assert_eq!(output.under_covered, vec![true, true, true, true]);
567        assert!(output.template.iter().all(|value| value.is_nan()));
568    }
569
570    #[test]
571    fn periodicity_strength_separates_solar_and_sidereal_components() {
572        let solar = solar_day_period();
573        let sidereal = repeat_period(GnssSystem::Gps);
574        let cadence = Duration::from_nanos(60_000_000_000);
575        let samples = (5.0 * SECONDS_PER_DAY / cadence.as_seconds()) as usize;
576
577        let solar_only = patterned_period_series(samples, solar, cadence, 1.0);
578        let sidereal_only = patterned_period_series(samples, sidereal, cadence, 1.0);
579        let mixed = solar_only
580            .iter()
581            .zip(&sidereal_only)
582            .map(|(&solar, &sidereal)| solar + 0.5 * sidereal)
583            .collect::<Vec<_>>();
584        let candidates = [solar, sidereal];
585
586        let solar_scores =
587            periodicity_strength_with_sample_interval(&solar_only, &candidates, cadence)
588                .expect("solar scores");
589        let sidereal_scores =
590            periodicity_strength_with_sample_interval(&sidereal_only, &candidates, cadence)
591                .expect("sidereal scores");
592        let mixed_scores = periodicity_strength_with_sample_interval(&mixed, &candidates, cadence)
593            .expect("mixed scores");
594
595        assert_eq!(solar_scores[0].1.to_bits(), 1.0f64.to_bits());
596        assert_eq!(solar_scores[1].1.to_bits(), 0x3fd7_0a3d_70a3_d708);
597        assert_eq!(sidereal_scores[0].1.to_bits(), 0x3fcf_db97_530e_ca84);
598        assert_eq!(sidereal_scores[1].1.to_bits(), 1.0f64.to_bits());
599        assert_eq!(mixed_scores[0].1.to_bits(), 0x3fe9_fdb9_7530_eca9);
600        assert_eq!(mixed_scores[1].1.to_bits(), 0x3fd7_0a3d_70a3_d708);
601    }
602
603    fn agnew_larson_table_record(sat: GnssSatelliteId, repeat_s: f64) -> BroadcastRecord {
604        let sqrt_a = 5_153.795_477_5;
605        let a = sqrt_a * sqrt_a;
606        let n0 = (ConstellationConstants::GPS.gm_m3_s2 / (a * a * a)).sqrt();
607        let delta_n = 4.0 * PI / repeat_s - n0;
608        let toe_sow = 100_000.0;
609        let toe = GnssWeekTow::new(crate::astro::time::TimeScale::Gpst, 2_295, toe_sow)
610            .expect("valid toe");
611        BroadcastRecord {
612            satellite_id: sat,
613            message: NavMessage::GpsLnav,
614            issue_of_data: BroadcastIssue {
615                issue: 0,
616                message: NavMessage::GpsLnav,
617            },
618            week: toe.week,
619            toe,
620            toc: toe,
621            elements: KeplerianElements {
622                sqrt_a,
623                e: 0.01,
624                m0: 0.0,
625                delta_n,
626                omega0: 1.0,
627                i0: 0.94,
628                omega: 0.25,
629                omega_dot: -8.0e-9,
630                idot: 0.0,
631                cuc: 0.0,
632                cus: 0.0,
633                crc: 100.0,
634                crs: -50.0,
635                cic: 0.0,
636                cis: 0.0,
637                toe_sow,
638            },
639            clock: ClockPolynomial {
640                af0: 0.0,
641                af1: 0.0,
642                af2: 0.0,
643                toc_sow: toe_sow,
644            },
645            group_delays: BroadcastGroupDelays::default(),
646            cnav: None,
647            sv_health: 0.0,
648            sv_accuracy_m: 0.0,
649            fit_interval_s: Some(4.0 * SECONDS_PER_HOUR),
650        }
651    }
652
653    fn variance(values: &[f64]) -> f64 {
654        let mean = values.iter().sum::<f64>() / values.len() as f64;
655        values
656            .iter()
657            .map(|value| {
658                let residual = value - mean;
659                residual * residual
660            })
661            .sum::<f64>()
662            / values.len() as f64
663    }
664
665    fn patterned_period_series(
666        samples: usize,
667        period: Duration,
668        cadence: Duration,
669        amplitude: f64,
670    ) -> Vec<f64> {
671        let pattern = [1.0, -0.5, 0.25, -1.0, 0.75];
672        let period_s = period.as_seconds();
673        let cadence_s = cadence.as_seconds();
674        let bins = phase_bin_count(period_s, cadence_s).expect("phase bins");
675        (0..samples)
676            .map(|sample_index| {
677                let bin = phase_bin(sample_index, period_s, cadence_s, bins);
678                amplitude * pattern[bin % pattern.len()]
679            })
680            .collect()
681    }
682}
sidereon_core/sidereal.rs

sidereon_core/
sidereal.rs
