sidereon_core/
rinex_clock.rs

1//! RINEX clock (`.CLK`) satellite-clock parser and interpolation.
2//!
3//! The parser owns the product grammar for `AS` satellite clock-bias records.
4//! The strict parser reports malformed `AS` rows. Use
5//! [`RinexClock::parse_lossy`] only when best-effort input recovery is intended.
6
7use std::cmp::Ordering;
8use std::collections::BTreeMap;
9use std::fmt;
10
11use crate::astro::time::model::{Instant, InstantRepr, JulianDateSplit, TimeScale};
12use crate::constants::{GPS_EPOCH_TO_J2000_S, J2000_JD, SECONDS_PER_DAY};
13use crate::validate::{self, FieldError};
14
15const INSTANT_SCALE_ORDER_STRIDE_S: f64 = 1.0e15;
16
17/// One satellite clock-bias sample.
18#[derive(Debug, Clone, Copy, PartialEq)]
19pub struct ClockPoint {
20    /// Scale-tagged epoch from the RINEX clock file's declared time system.
21    pub epoch: Instant,
22    /// Satellite clock bias in seconds.
23    pub bias_s: f64,
24}
25
26impl ClockPoint {
27    /// This sample's epoch as GPS seconds, when the sample is actually GPST.
28    pub fn gps_seconds(&self) -> Option<f64> {
29        instant_to_gps_seconds(&self.epoch)
30    }
31}
32
33/// Civil epoch tag used by RINEX clock records, interpreted in the file's time scale.
34#[derive(Debug, Clone, Copy, PartialEq)]
35pub struct ClockEpoch {
36    /// Four-digit calendar year.
37    pub year: i32,
38    /// Calendar month, 1..=12.
39    pub month: u8,
40    /// Calendar day of month, 1..=31.
41    pub day: u8,
42    /// Hour of day, 0..=23.
43    pub hour: u8,
44    /// Minute of hour, 0..=59.
45    pub minute: u8,
46    /// Seconds of minute, including fractional seconds.
47    pub second: f64,
48}
49
50/// Parsed RINEX clock product.
51#[derive(Debug, Clone, PartialEq)]
52pub struct RinexClock {
53    /// Time scale declared by the RINEX clock header. Missing headers default to GPST.
54    pub time_scale: TimeScale,
55    /// Per-satellite, strictly time-ordered clock-bias series.
56    pub series: BTreeMap<String, Vec<ClockPoint>>,
57}
58
59/// RINEX clock parse error.
60#[derive(Debug, Clone, PartialEq, Eq)]
61pub enum RinexClockError {
62    /// An `AS` satellite clock row is too short to carry the required bias.
63    MalformedAsRecord {
64        /// One-based input line number.
65        line: usize,
66        /// Human-readable parse failure.
67        reason: &'static str,
68        /// The full record text.
69        record: String,
70    },
71    /// A required `AS` field could not be parsed or was out of range.
72    BadField {
73        /// One-based input line number.
74        line: usize,
75        /// Field name.
76        field: &'static str,
77        /// Source field value.
78        value: String,
79    },
80    /// Public manual input or query parameter was invalid.
81    InvalidInput {
82        /// Field name.
83        field: &'static str,
84        /// Human-readable validation failure.
85        reason: &'static str,
86    },
87}
88
89impl fmt::Display for RinexClockError {
90    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
91        match self {
92            RinexClockError::MalformedAsRecord {
93                line,
94                reason,
95                record,
96            } => write!(
97                f,
98                "malformed RINEX AS clock record at line {line}: {reason}: {record}"
99            ),
100            RinexClockError::BadField { line, field, value } => write!(
101                f,
102                "bad RINEX AS clock field at line {line}: {field}={value}"
103            ),
104            RinexClockError::InvalidInput { field, reason } => {
105                write!(f, "invalid RINEX clock input {field}: {reason}")
106            }
107        }
108    }
109}
110
111impl std::error::Error for RinexClockError {}
112
113impl RinexClock {
114    /// Parse a RINEX clock text into per-satellite `AS` records.
115    pub fn parse(text: &str) -> Result<Self, RinexClockError> {
116        let time_scale = parse_time_scale(text)?;
117        let lines = data_lines(text);
118        let mut by_sat = BTreeMap::<String, Vec<(ClockPoint, usize)>>::new();
119
120        for (line_number, line) in lines {
121            if let Some((sat, point)) = parse_record(line_number, line, time_scale)? {
122                by_sat.entry(sat).or_default().push((point, line_number));
123            }
124        }
125
126        Ok(Self {
127            time_scale,
128            series: build_series(by_sat),
129        })
130    }
131
132    /// Parse a RINEX clock text while skipping malformed and non-`AS` records.
133    pub fn parse_lossy(text: &str) -> Self {
134        let time_scale = parse_time_scale(text).unwrap_or(TimeScale::Gpst);
135        let lines = data_lines(text);
136        let mut by_sat = BTreeMap::<String, Vec<(ClockPoint, usize)>>::new();
137
138        for (line_number, line) in lines {
139            if let Ok(Some((sat, point))) = parse_record(line_number, line, time_scale) {
140                by_sat.entry(sat).or_default().push((point, line_number));
141            }
142        }
143
144        Self {
145            time_scale,
146            series: build_series(by_sat),
147        }
148    }
149
150    /// Rebuild a GPST product from the legacy public GPS-second row shape.
151    pub fn from_series_rows(rows: Vec<(String, Vec<(f64, f64)>)>) -> Result<Self, RinexClockError> {
152        let rows = rows
153            .into_iter()
154            .map(|(sat, points)| {
155                validate::require_strictly_increasing(
156                    points.iter().map(|&(gps_seconds, _)| gps_seconds),
157                    "gps_seconds",
158                )
159                .map_err(map_manual_order_error)?;
160                let points = points
161                    .into_iter()
162                    .map(|(gps_seconds, bias_s)| {
163                        validate_finite(bias_s, "bias_s")?;
164                        Ok((gps_seconds_to_instant(gps_seconds), bias_s))
165                    })
166                    .collect::<Result<Vec<_>, RinexClockError>>()?;
167                Ok((sat, points))
168            })
169            .collect::<Result<Vec<_>, RinexClockError>>()?;
170        Self::from_instant_series_rows(TimeScale::Gpst, rows)
171    }
172
173    /// Rebuild a parsed product from scale-tagged instant rows.
174    pub fn from_instant_series_rows(
175        time_scale: TimeScale,
176        rows: Vec<(String, Vec<(Instant, f64)>)>,
177    ) -> Result<Self, RinexClockError> {
178        let mut series = BTreeMap::new();
179        for (sat, points) in rows {
180            let mut indexed = points
181                .into_iter()
182                .enumerate()
183                .map(|(idx, (epoch, bias_s))| {
184                    let point = ClockPoint { epoch, bias_s };
185                    validate_clock_point(point)?;
186                    Ok((point, idx))
187                })
188                .collect::<Result<Vec<_>, RinexClockError>>()?;
189            validate_instant_series_order(&indexed)?;
190            indexed.sort_by(|(a, ai), (b, bi)| {
191                compare_instants(&a.epoch, &b.epoch).then_with(|| ai.cmp(bi))
192            });
193            series.insert(sat, dedup_by_time(indexed));
194        }
195        Ok(Self { time_scale, series })
196    }
197
198    /// Export GPST samples as `[(satellite, [(gps_seconds, bias_s), ...]), ...]`.
199    ///
200    /// Non-GPST samples are not coerced into GPS seconds and are omitted.
201    pub fn series_rows(&self) -> Vec<(String, Vec<(f64, f64)>)> {
202        self.series
203            .iter()
204            .map(|(sat, points)| {
205                (
206                    sat.clone(),
207                    points
208                        .iter()
209                        .filter_map(|point| Some((point.gps_seconds()?, point.bias_s)))
210                        .collect(),
211                )
212            })
213            .collect()
214    }
215
216    /// Export the product as scale-tagged instant rows.
217    pub fn instant_series_rows(&self) -> Vec<(String, Vec<(Instant, f64)>)> {
218        self.series
219            .iter()
220            .map(|(sat, points)| {
221                (
222                    sat.clone(),
223                    points
224                        .iter()
225                        .map(|point| (point.epoch, point.bias_s))
226                        .collect(),
227                )
228            })
229            .collect()
230    }
231
232    /// Interpolate one satellite clock bias at a civil epoch in this file's scale.
233    pub fn clock_s(
234        &self,
235        satellite_id: &str,
236        epoch: ClockEpoch,
237    ) -> Result<Option<f64>, RinexClockError> {
238        let epoch = civil_to_clock_instant(
239            self.time_scale,
240            epoch.year,
241            epoch.month,
242            epoch.day,
243            epoch.hour,
244            epoch.minute,
245            epoch.second,
246        )
247        .ok_or_else(|| invalid_input("epoch", "invalid civil clock epoch"))?;
248        self.clock_s_at_instant(satellite_id, epoch)
249    }
250
251    /// Interpolate one satellite clock bias at a scale-tagged instant.
252    pub fn clock_s_at_instant(
253        &self,
254        satellite_id: &str,
255        epoch: Instant,
256    ) -> Result<Option<f64>, RinexClockError> {
257        validate_instant(epoch, "epoch")?;
258        let Some(records) = self.series.get(satellite_id) else {
259            return Ok(None);
260        };
261        Ok(interpolate(records, epoch))
262    }
263
264    /// Interpolate one satellite clock bias at GPS seconds.
265    pub fn clock_s_at_gps_seconds(
266        &self,
267        satellite_id: &str,
268        gps_seconds: f64,
269    ) -> Result<Option<f64>, RinexClockError> {
270        validate_finite(gps_seconds, "gps_seconds")?;
271        self.clock_s_at_instant(satellite_id, gps_seconds_to_instant(gps_seconds))
272    }
273}
274
275/// Convert a civil clock tag in the given scale into a scale-tagged instant.
276pub fn civil_to_clock_instant(
277    scale: TimeScale,
278    year: i32,
279    month: u8,
280    day: u8,
281    hour: u8,
282    minute: u8,
283    second: f64,
284) -> Option<Instant> {
285    let civil = validate::civil_datetime_with_fractional_second_policy(
286        i64::from(year),
287        i64::from(month),
288        i64::from(day),
289        i64::from(hour),
290        i64::from(minute),
291        second,
292        civil_second_policy_for_time_scale(scale),
293    )
294    .ok()?;
295    civil_microsecond_to_instant(scale, civil).ok()
296}
297
298/// Convert a civil GPS-time tag into seconds since 1980-01-06 00:00:00.
299pub fn civil_to_gps_seconds(
300    year: i32,
301    month: u8,
302    day: u8,
303    hour: u8,
304    minute: u8,
305    second: f64,
306) -> Option<f64> {
307    let civil = validate::civil_datetime_with_fractional_second_policy(
308        i64::from(year),
309        i64::from(month),
310        i64::from(day),
311        i64::from(hour),
312        i64::from(minute),
313        second,
314        validate::CivilSecondPolicy::Continuous,
315    )
316    .ok()?;
317    gps_seconds_from_civil(civil)
318}
319
320fn parse_time_scale(text: &str) -> Result<TimeScale, RinexClockError> {
321    let mut time_scale = TimeScale::Gpst;
322    for (idx, line) in text.lines().enumerate() {
323        if line.contains("END OF HEADER") {
324            break;
325        }
326        if line.contains("TIME SYSTEM ID") {
327            let label = line
328                .split("TIME SYSTEM ID")
329                .next()
330                .unwrap_or(line)
331                .split_whitespace()
332                .next()
333                .unwrap_or("");
334            if label.is_empty() {
335                time_scale = TimeScale::Gpst;
336            } else {
337                time_scale = crate::parse::time_scale_label(label).ok_or_else(|| {
338                    RinexClockError::BadField {
339                        line: idx + 1,
340                        field: "time_system",
341                        value: label.to_string(),
342                    }
343                })?;
344            }
345        }
346    }
347    Ok(time_scale)
348}
349
350fn gps_seconds_to_instant(gps_seconds: f64) -> Instant {
351    let gps_epoch_jd = J2000_JD - GPS_EPOCH_TO_J2000_S / SECONDS_PER_DAY;
352    let days = (gps_seconds / SECONDS_PER_DAY).floor();
353    let seconds_of_day = gps_seconds - days * SECONDS_PER_DAY;
354    Instant::from_julian_date(
355        TimeScale::Gpst,
356        JulianDateSplit::new(gps_epoch_jd + days, seconds_of_day / SECONDS_PER_DAY)
357            .expect("valid split Julian date"),
358    )
359}
360
361fn validate_clock_point(point: ClockPoint) -> Result<(), RinexClockError> {
362    validate_instant(point.epoch, "epoch")?;
363    validate_finite(point.bias_s, "bias_s")
364}
365
366fn validate_instant(epoch: Instant, field: &'static str) -> Result<(), RinexClockError> {
367    match epoch.repr {
368        InstantRepr::JulianDate(split) => {
369            validate_finite(split.jd_whole, field)?;
370            validate_finite(split.fraction, field)?;
371            if !(-1.0..=1.0).contains(&split.fraction) {
372                return Err(invalid_input(field, "Julian-date fraction out of range"));
373            }
374            Ok(())
375        }
376        InstantRepr::Nanos(_) => Ok(()),
377    }
378}
379
380fn validate_finite(value: f64, field: &'static str) -> Result<(), RinexClockError> {
381    if value.is_finite() {
382        Ok(())
383    } else {
384        Err(invalid_input(field, "must be finite"))
385    }
386}
387
388fn invalid_input(field: &'static str, reason: &'static str) -> RinexClockError {
389    RinexClockError::InvalidInput { field, reason }
390}
391
392fn map_manual_order_error(error: FieldError) -> RinexClockError {
393    match error {
394        FieldError::NonFinite { field } => invalid_input(field, "must be finite"),
395        FieldError::OutOfRange { field, .. } => invalid_input(field, "must be strictly increasing"),
396        _ => invalid_input(error.field(), error.reason()),
397    }
398}
399
400fn validate_instant_series_order(points: &[(ClockPoint, usize)]) -> Result<(), RinexClockError> {
401    validate::require_strictly_increasing(
402        points
403            .iter()
404            .map(|(point, _)| instant_order_key(&point.epoch)),
405        "epoch",
406    )
407    .map_err(map_manual_order_error)
408}
409
410fn instant_order_key(epoch: &Instant) -> f64 {
411    let offset_s = time_scale_rank(epoch.scale) as f64 * INSTANT_SCALE_ORDER_STRIDE_S;
412    let instant_s = match epoch.repr {
413        InstantRepr::JulianDate(split) => {
414            split.jd_whole * SECONDS_PER_DAY + split.fraction * SECONDS_PER_DAY
415        }
416        InstantRepr::Nanos(nanos) => nanos as f64 / 1.0e9,
417    };
418    offset_s + instant_s
419}
420
421fn instant_to_gps_seconds(epoch: &Instant) -> Option<f64> {
422    if epoch.scale != TimeScale::Gpst {
423        return None;
424    }
425    instant_to_j2000_seconds(epoch).map(|seconds| seconds + GPS_EPOCH_TO_J2000_S)
426}
427
428fn instant_to_j2000_seconds(epoch: &Instant) -> Option<f64> {
429    match epoch.repr {
430        InstantRepr::JulianDate(split) => {
431            Some((split.jd_whole - J2000_JD) * SECONDS_PER_DAY + split.fraction * SECONDS_PER_DAY)
432        }
433        InstantRepr::Nanos(_) => None,
434    }
435}
436
437fn data_lines(text: &str) -> Vec<(usize, &str)> {
438    drop_header(
439        text.lines()
440            .enumerate()
441            .map(|(idx, line)| (idx + 1, line))
442            .collect(),
443    )
444}
445
446fn drop_header(lines: Vec<(usize, &str)>) -> Vec<(usize, &str)> {
447    match lines
448        .iter()
449        .position(|(_, line)| line.contains("END OF HEADER"))
450    {
451        Some(idx) => lines.into_iter().skip(idx + 1).collect(),
452        None => lines,
453    }
454}
455
456#[derive(Debug, Clone, Copy)]
457struct ClockEpochFields<'a> {
458    year: i32,
459    month: u8,
460    day: u8,
461    hour: u8,
462    minute: u8,
463    second: &'a str,
464}
465
466fn parse_record(
467    line_number: usize,
468    line: &str,
469    time_scale: TimeScale,
470) -> Result<Option<(String, ClockPoint)>, RinexClockError> {
471    let fields = line.split_whitespace().collect::<Vec<_>>();
472    if fields.first() != Some(&"AS") {
473        return Ok(None);
474    }
475    if fields.len() < 10 {
476        return Err(RinexClockError::MalformedAsRecord {
477            line: line_number,
478            reason: "expected at least 10 fields",
479            record: line.trim().to_string(),
480        });
481    }
482
483    let sat = validate::strict_gnss_satellite_id(fields[1], "satellite")
484        .map_err(|error| map_field_error(line_number, error, fields[1]))?
485        .to_string();
486    let year = parse_int_field::<i32>(line_number, "year", fields[2])?;
487    let month = parse_int_field::<u8>(line_number, "month", fields[3])?;
488    let day = parse_int_field::<u8>(line_number, "day", fields[4])?;
489    let hour = parse_int_field::<u8>(line_number, "hour", fields[5])?;
490    let minute = parse_int_field::<u8>(line_number, "minute", fields[6])?;
491    let epoch = ClockEpochFields {
492        year,
493        month,
494        day,
495        hour,
496        minute,
497        second: fields[7],
498    };
499    let bias_s = parse_f64_field(line_number, "bias", fields[9])?;
500    let epoch = civil_decimal_second_to_instant(time_scale, epoch)
501        .map_err(|error| map_epoch_error(line_number, error, epoch))?;
502
503    Ok(Some((sat, ClockPoint { epoch, bias_s })))
504}
505
506fn parse_int_field<T>(
507    line_number: usize,
508    field: &'static str,
509    value: &str,
510) -> Result<T, RinexClockError>
511where
512    T: std::str::FromStr,
513{
514    validate::strict_int(value, field).map_err(|error| map_field_error(line_number, error, value))
515}
516
517fn parse_f64_field(
518    line_number: usize,
519    field: &'static str,
520    value: &str,
521) -> Result<f64, RinexClockError> {
522    validate::strict_f64(value, field).map_err(|error| map_field_error(line_number, error, value))
523}
524
525fn civil_decimal_second_to_instant(
526    scale: TimeScale,
527    epoch: ClockEpochFields<'_>,
528) -> Result<Instant, FieldError> {
529    let civil = validate::civil_datetime_with_decimal_second_policy(
530        i64::from(epoch.year),
531        i64::from(epoch.month),
532        i64::from(epoch.day),
533        i64::from(epoch.hour),
534        i64::from(epoch.minute),
535        epoch.second,
536        civil_second_policy_for_time_scale(scale),
537    )?;
538    civil_microsecond_to_instant(scale, civil)
539}
540
541fn civil_microsecond_to_instant(
542    scale: TimeScale,
543    civil: validate::ValidCivilMicrosecond,
544) -> Result<Instant, FieldError> {
545    let split = civil_microsecond_to_julian_split(scale, civil)?;
546    Ok(Instant::from_julian_date(scale, split))
547}
548
549fn civil_microsecond_to_julian_split(
550    scale: TimeScale,
551    civil: validate::ValidCivilMicrosecond,
552) -> Result<JulianDateSplit, FieldError> {
553    if civil.year < 1 {
554        return Err(FieldError::InvalidCivilDate {
555            field: "civil datetime",
556            year: civil.year,
557            month: i64::from(civil.month),
558            day: i64::from(civil.day),
559        });
560    }
561
562    let month = i64::from(civil.month);
563    let day = i64::from(civil.day);
564    let a = (14 - month) / 12;
565    let y = civil.year + 4800 - a;
566    let m = month + 12 * a - 3;
567    let jdn = day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32_045;
568    let jd_whole = jdn as f64 - 0.5;
569    if scale == TimeScale::Utc && civil.second == 60 {
570        let remaining_s = 1.0 - civil.microsecond as f64 / 1_000_000.0;
571        return Ok(
572            JulianDateSplit::new(jd_whole + 1.0, -remaining_s / SECONDS_PER_DAY)
573                .expect("valid leap-second split Julian date"),
574        );
575    }
576
577    let day_seconds = civil.hour as f64 * 3600.0
578        + civil.minute as f64 * 60.0
579        + civil.second as f64
580        + civil.microsecond as f64 / 1_000_000.0;
581    Ok(
582        JulianDateSplit::new(jd_whole, day_seconds / SECONDS_PER_DAY)
583            .expect("valid split Julian date"),
584    )
585}
586
587fn civil_second_policy_for_time_scale(scale: TimeScale) -> validate::CivilSecondPolicy {
588    match scale {
589        TimeScale::Utc => validate::CivilSecondPolicy::UtcLike,
590        TimeScale::Tai
591        | TimeScale::Tt
592        | TimeScale::Tdb
593        | TimeScale::Gpst
594        | TimeScale::Gst
595        | TimeScale::Bdt => validate::CivilSecondPolicy::Continuous,
596    }
597}
598
599fn gps_seconds_from_civil(civil: validate::ValidCivilMicrosecond) -> Option<f64> {
600    if civil.year < 1 {
601        return None;
602    }
603
604    let days = days_since_gps_epoch(civil.year as i32, civil.month as u8, civil.day as u8);
605    let whole = days as f64 * SECONDS_PER_DAY
606        + (i64::from(civil.hour) * 3_600 + i64::from(civil.minute) * 60 + i64::from(civil.second))
607            as f64;
608    Some(whole + f64::from(civil.microsecond) / 1_000_000.0)
609}
610
611fn map_field_error(line_number: usize, error: FieldError, value: &str) -> RinexClockError {
612    RinexClockError::BadField {
613        line: line_number,
614        field: error.field(),
615        value: value.to_string(),
616    }
617}
618
619fn map_epoch_error(
620    line_number: usize,
621    error: FieldError,
622    epoch: ClockEpochFields<'_>,
623) -> RinexClockError {
624    match error {
625        FieldError::FloatParse { .. }
626        | FieldError::Missing { .. }
627        | FieldError::NonFinite { .. } => RinexClockError::BadField {
628            line: line_number,
629            field: "second",
630            value: epoch.second.to_string(),
631        },
632        _ => RinexClockError::BadField {
633            line: line_number,
634            field: "epoch",
635            value: format!(
636                "{} {} {} {} {} {}",
637                epoch.year,
638                epoch.month,
639                epoch.day,
640                epoch.hour,
641                epoch.minute,
642                normalized_second_text(epoch.second)
643            ),
644        },
645    }
646}
647
648fn normalized_second_text(second: &str) -> String {
649    validate::strict_f64(second, "second")
650        .map_or_else(|_| second.to_string(), |value| value.to_string())
651}
652
653fn build_series(
654    by_sat: BTreeMap<String, Vec<(ClockPoint, usize)>>,
655) -> BTreeMap<String, Vec<ClockPoint>> {
656    by_sat
657        .into_iter()
658        .map(|(sat, mut points)| {
659            points.sort_by(|(a, ai), (b, bi)| {
660                compare_instants(&a.epoch, &b.epoch).then_with(|| ai.cmp(bi))
661            });
662            (sat, dedup_by_time(points))
663        })
664        .collect()
665}
666
667fn dedup_by_time(points: Vec<(ClockPoint, usize)>) -> Vec<ClockPoint> {
668    let mut deduped = Vec::<ClockPoint>::new();
669    for (point, _) in points {
670        match deduped.last_mut() {
671            Some(prev) if prev.epoch == point.epoch => *prev = point,
672            _ => deduped.push(point),
673        }
674    }
675    deduped
676}
677
678fn interpolate(records: &[ClockPoint], epoch: Instant) -> Option<f64> {
679    let mut prev: Option<ClockPoint> = None;
680    for point in records {
681        match compare_instants_same_scale(&point.epoch, &epoch)? {
682            Ordering::Equal => return Some(point.bias_s),
683            Ordering::Greater => {
684                let p0 = prev?;
685                let p1 = *point;
686                let span_s = seconds_between(&p1.epoch, &p0.epoch)?;
687                if span_s <= 0.0 {
688                    return None;
689                }
690                let query_s = seconds_between(&epoch, &p0.epoch)?;
691                if query_s < 0.0 {
692                    return None;
693                }
694                return Some(p0.bias_s + (p1.bias_s - p0.bias_s) * query_s / span_s);
695            }
696            Ordering::Less => prev = Some(*point),
697        }
698    }
699    None
700}
701
702fn compare_instants(a: &Instant, b: &Instant) -> Ordering {
703    time_scale_rank(a.scale)
704        .cmp(&time_scale_rank(b.scale))
705        .then_with(|| match (a.julian_date(), b.julian_date()) {
706            (Some(a), Some(b)) => compare_julian_splits(a, b),
707            _ => Ordering::Equal,
708        })
709}
710
711fn compare_instants_same_scale(a: &Instant, b: &Instant) -> Option<Ordering> {
712    if a.scale != b.scale {
713        return None;
714    }
715    Some(compare_julian_splits(a.julian_date()?, b.julian_date()?))
716}
717
718fn compare_julian_splits(a: JulianDateSplit, b: JulianDateSplit) -> Ordering {
719    a.jd_whole
720        .partial_cmp(&b.jd_whole)
721        .unwrap_or(Ordering::Equal)
722        .then_with(|| {
723            a.fraction
724                .partial_cmp(&b.fraction)
725                .unwrap_or(Ordering::Equal)
726        })
727}
728
729fn seconds_between(later: &Instant, earlier: &Instant) -> Option<f64> {
730    if later.scale != earlier.scale {
731        return None;
732    }
733    let later = later.julian_date()?;
734    let earlier = earlier.julian_date()?;
735    let seconds = ((later.jd_whole - earlier.jd_whole) + (later.fraction - earlier.fraction))
736        * SECONDS_PER_DAY;
737    seconds.is_finite().then_some(seconds)
738}
739
740fn time_scale_rank(scale: TimeScale) -> u8 {
741    match scale {
742        TimeScale::Utc => 0,
743        TimeScale::Tai => 1,
744        TimeScale::Tt => 2,
745        TimeScale::Tdb => 3,
746        TimeScale::Gpst => 4,
747        TimeScale::Gst => 5,
748        TimeScale::Bdt => 6,
749    }
750}
751
752fn days_since_gps_epoch(year: i32, month: u8, day: u8) -> i64 {
753    days_before_date(year, month, day) - days_before_date(1980, 1, 6)
754}
755
756fn days_before_date(year: i32, month: u8, day: u8) -> i64 {
757    let mut days = days_before_year(year);
758    for m in 1..month {
759        days += i64::from(days_in_month(year, m));
760    }
761    days + i64::from(day - 1)
762}
763
764fn days_before_year(year: i32) -> i64 {
765    let y = i64::from(year - 1);
766    y * 365 + y / 4 - y / 100 + y / 400
767}
768
769fn days_in_month(year: i32, month: u8) -> u8 {
770    match month {
771        1 | 3 | 5 | 7 | 8 | 10 | 12 => 31,
772        4 | 6 | 9 | 11 => 30,
773        2 if is_leap_year(year) => 29,
774        2 => 28,
775        _ => 0,
776    }
777}
778
779fn is_leap_year(year: i32) -> bool {
780    (year % 4 == 0 && year % 100 != 0) || year % 400 == 0
781}
782
783#[cfg(test)]
784mod tests {
785    use super::*;
786
787    fn as_record(satellite: &str, bias: &str) -> String {
788        format!("AS {satellite} 2020 01 01 00 00 00.000000 1 {bias}")
789    }
790
791    #[test]
792    fn parse_rejects_non_finite_as_bias() {
793        let err = RinexClock::parse(&as_record("G01", "NaN")).unwrap_err();
794        assert_eq!(
795            err,
796            RinexClockError::BadField {
797                line: 1,
798                field: "bias",
799                value: "NaN".to_string(),
800            }
801        );
802    }
803
804    #[test]
805    fn parse_rejects_malformed_as_satellite_token() {
806        let err = RinexClock::parse(&as_record("X01", "1.0e-9")).unwrap_err();
807        assert_eq!(
808            err,
809            RinexClockError::BadField {
810                line: 1,
811                field: "satellite",
812                value: "X01".to_string(),
813            }
814        );
815    }
816
817    #[test]
818    fn explicit_utc_time_system_preserves_clock_epoch_scale() {
819        let text = " 3.00           C                                       RINEX VERSION / TYPE\n\
820                    UTC                                                     TIME SYSTEM ID\n\
821                                                                        END OF HEADER\n\
822                    AS G05  2017 01 01 00 00  0.000000  1   1.0e-04\n\
823                    AS G05  2017 01 01 00 00 30.000000  1   2.0e-04\n";
824        let clock = RinexClock::parse(text).expect("UTC RINEX clock");
825
826        assert_eq!(clock.time_scale, TimeScale::Utc);
827        assert_eq!(clock.series["G05"][0].epoch.scale, TimeScale::Utc);
828        let interpolated = clock
829            .clock_s(
830                "G05",
831                ClockEpoch {
832                    year: 2017,
833                    month: 1,
834                    day: 1,
835                    hour: 0,
836                    minute: 0,
837                    second: 15.0,
838                },
839            )
840            .expect("valid clock query")
841            .expect("UTC interpolated clock");
842        assert!((interpolated - 1.5e-4).abs() < 1.0e-18);
843
844        let gpst_query =
845            civil_to_clock_instant(TimeScale::Gpst, 2017, 1, 1, 0, 0, 15.0).expect("GPST instant");
846        assert_eq!(
847            clock
848                .clock_s_at_instant("G05", gpst_query)
849                .expect("valid clock query"),
850            None
851        );
852
853        let rows = clock.instant_series_rows();
854        assert_eq!(rows[0].1[0].0.scale, TimeScale::Utc);
855        let rebuilt = RinexClock::from_instant_series_rows(clock.time_scale, rows)
856            .expect("valid manual RINEX clock rows");
857        assert_eq!(rebuilt, clock);
858    }
859
860    #[test]
861    fn manual_series_rows_reject_non_finite_inputs() {
862        assert_eq!(
863            RinexClock::from_series_rows(vec![("G05".to_string(), vec![(f64::NAN, 1.0e-4)])])
864                .unwrap_err(),
865            RinexClockError::InvalidInput {
866                field: "gps_seconds",
867                reason: "must be finite",
868            }
869        );
870        assert_eq!(
871            RinexClock::from_series_rows(vec![(
872                "G05".to_string(),
873                vec![(1_463_904_000.0, f64::INFINITY)]
874            )])
875            .unwrap_err(),
876            RinexClockError::InvalidInput {
877                field: "bias_s",
878                reason: "must be finite",
879            }
880        );
881    }
882
883    #[test]
884    fn manual_series_rows_reject_unsorted_gps_seconds() {
885        assert_eq!(
886            RinexClock::from_series_rows(vec![(
887                "G05".to_string(),
888                vec![(1_463_904_030.0, 1.0e-4), (1_463_904_000.0, 2.0e-4)]
889            )])
890            .unwrap_err(),
891            RinexClockError::InvalidInput {
892                field: "gps_seconds",
893                reason: "must be strictly increasing",
894            }
895        );
896    }
897
898    #[test]
899    fn manual_instant_rows_reject_non_finite_inputs() {
900        let bad_epoch = Instant::from_julian_date(
901            TimeScale::Gpst,
902            JulianDateSplit {
903                jd_whole: f64::NAN,
904                fraction: 0.0,
905            },
906        );
907        assert_eq!(
908            RinexClock::from_instant_series_rows(
909                TimeScale::Gpst,
910                vec![("G05".to_string(), vec![(bad_epoch, 1.0e-4)])],
911            )
912            .unwrap_err(),
913            RinexClockError::InvalidInput {
914                field: "epoch",
915                reason: "must be finite",
916            }
917        );
918
919        let good_epoch =
920            civil_to_clock_instant(TimeScale::Gpst, 2026, 5, 13, 0, 0, 0.0).expect("GPST instant");
921        assert_eq!(
922            RinexClock::from_instant_series_rows(
923                TimeScale::Gpst,
924                vec![("G05".to_string(), vec![(good_epoch, f64::NAN)])],
925            )
926            .unwrap_err(),
927            RinexClockError::InvalidInput {
928                field: "bias_s",
929                reason: "must be finite",
930            }
931        );
932    }
933
934    #[test]
935    fn manual_instant_rows_reject_unsorted_epochs() {
936        let later =
937            civil_to_clock_instant(TimeScale::Gpst, 2026, 5, 13, 0, 0, 30.0).expect("later epoch");
938        let earlier =
939            civil_to_clock_instant(TimeScale::Gpst, 2026, 5, 13, 0, 0, 0.0).expect("earlier epoch");
940
941        assert_eq!(
942            RinexClock::from_instant_series_rows(
943                TimeScale::Gpst,
944                vec![("G05".to_string(), vec![(later, 1.0e-4), (earlier, 2.0e-4)])],
945            )
946            .unwrap_err(),
947            RinexClockError::InvalidInput {
948                field: "epoch",
949                reason: "must be strictly increasing",
950            }
951        );
952    }
953
954    #[test]
955    fn rinex_clock_queries_reject_non_finite_inputs() {
956        let clock = RinexClock::from_series_rows(vec![(
957            "G05".to_string(),
958            vec![(1_463_904_000.0, 1.0e-4)],
959        )])
960        .expect("valid manual RINEX clock rows");
961        let bad_epoch = Instant::from_julian_date(
962            TimeScale::Gpst,
963            JulianDateSplit {
964                jd_whole: f64::INFINITY,
965                fraction: 0.0,
966            },
967        );
968        assert_eq!(
969            clock.clock_s_at_instant("G05", bad_epoch).unwrap_err(),
970            RinexClockError::InvalidInput {
971                field: "epoch",
972                reason: "must be finite",
973            }
974        );
975        assert_eq!(
976            clock.clock_s_at_gps_seconds("G05", f64::NAN).unwrap_err(),
977            RinexClockError::InvalidInput {
978                field: "gps_seconds",
979                reason: "must be finite",
980            }
981        );
982        assert_eq!(
983            clock
984                .clock_s(
985                    "G05",
986                    ClockEpoch {
987                        year: 2026,
988                        month: 5,
989                        day: 13,
990                        hour: 0,
991                        minute: 0,
992                        second: f64::NAN,
993                    },
994                )
995                .unwrap_err(),
996            RinexClockError::InvalidInput {
997                field: "epoch",
998                reason: "invalid civil clock epoch",
999            }
1000        );
1001    }
1002
1003    #[test]
1004    fn interpolation_rejects_non_positive_bracket_span() {
1005        let day = 2_457_753.5;
1006        let p0 = Instant::from_julian_date(
1007            TimeScale::Utc,
1008            JulianDateSplit::new(day, 1.0).expect("valid split Julian date"),
1009        );
1010        let p1 = Instant::from_julian_date(
1011            TimeScale::Utc,
1012            JulianDateSplit::new(day + 1.0, 0.0).expect("valid split Julian date"),
1013        );
1014        let query = Instant::from_julian_date(
1015            TimeScale::Utc,
1016            JulianDateSplit::new(day + 1.0, 0.5 / SECONDS_PER_DAY)
1017                .expect("valid split Julian date"),
1018        );
1019        let records = [
1020            ClockPoint {
1021                epoch: p0,
1022                bias_s: 1.0e-4,
1023            },
1024            ClockPoint {
1025                epoch: p1,
1026                bias_s: 2.0e-4,
1027            },
1028        ];
1029
1030        assert_eq!(interpolate(&records, query), None);
1031    }
1032}
sidereon_core/rinex_clock.rs

sidereon_core/
rinex_clock.rs
