Skip to main content

sidereon_core/sp3/
mod.rs

1//! SP3-c / SP3-d precise-ephemeris parser.
2//!
3//! Parses the IGS SP3 precise orbit/clock format, both **SP3-c** and **SP3-d**
4//! (Hilla 2016), into a typed [`Sp3`] product. The parser is multi-GNSS,
5//! handles position/clock records plus optional velocity records,
6//! missing-value sentinels, predicted / clock-event / maneuver flags, and a
7//! system-aware [`GnssSatelliteId`]; the product's time system is read from the
8//! header.
9//!
10//! # Build vs adopt
11//!
12//! The spec permits using the `sp3` crate (MPL-2.0) as a deterministic byte
13//! reader, OR hand-rolling the record parsing. **This module hand-rolls it**,
14//! deliberately:
15//!
16//! - The `refs/sp3` crate hard-depends on `hifitime` for its `Epoch`,
17//!   `TimeScale`, and `Duration`, and on `flate2`. `sidereon-core` models time
18//!   with the **core crate's own** [`Instant`] / [`TimeScale`]
19//!   family, which is hifitime-free; adopting the `sp3` crate would
20//!   invert that and pull a parallel time stack into the GNSS layer.
21//! - The `sp3` crate also carries its own `SV` / `Constellation` identifiers
22//!   that duplicate this crate's [`GnssSatelliteId`] / [`GnssSystem`].
23//! - The SP3 record grammar is small, fixed-column, and fully specified, so a
24//!   byte reader is low-risk. (Note: the `refs/sp3` velocity parser at
25//!   `parsing.rs:241-245` has an axis bug - it reuses the Y component for X;
26//!   this module reads each axis independently and is unit-tested for it.)
27//!
28//! Parsing only is adopted-grade work; it is **not** a contested float recipe.
29//! The interpolation that consumes this product is built
30//! separately to match the `scipy.interpolate` reference and is out of scope
31//! for this module.
32//!
33//! # Units
34//!
35//! SP3 stores positions in **kilometers** and clock offsets in **microseconds**
36//! (velocities in dm/s, clock-rate in 1e-4 us/s). This parser converts at parse
37//! time to the crate's internal SI base units - positions in **meters**
38//! (`km * 1000.0`), clocks in **seconds** (`us * 1e-6`), velocities in **m/s**
39//! (`(dm/s) * 1e-1`), clock-rate in **s/s** (`(1e-4 us/s) * 1e-10`). Each scale
40//! factor is applied as a single multiply so the operation order is fixed for
41//! the clock-unit-conversion golden test.
42//!
43//! # Frames
44//!
45//! Positions/velocities are returned as frame-tagged [`ItrfPositionM`] /
46//! [`ItrfVelocityMS`], never a bare `position_m`.
47
48use std::collections::BTreeMap;
49
50use crate::astro::time::civil::{j2000_seconds, split_julian_date};
51use crate::astro::time::model::{Instant, InstantRepr, JulianDateSplit, TimeScale};
52
53use crate::constants::{KM_TO_M, US_TO_S};
54use crate::format::columns::{
55    char_at, raw_field as field, raw_field_from as field_from, strict_f64,
56};
57use crate::format::{Diagnostics, RecordRef, Skip, SkipReason};
58use crate::frame::{ItrfPositionM, ItrfVelocityMS};
59use crate::id::{is_valid_prn, GnssSatelliteId, GnssSystem};
60use crate::validate;
61use crate::{Error, Result};
62
63/// SP3 missing/bad position component sentinel, in kilometers.
64///
65/// SP3 writes a satellite with no usable orbit as a position record of exactly
66/// `0.000000 0.000000 0.000000`. We treat an all-zero position as "missing"
67/// (matching the `refs/sp3` validity guard at `parsing.rs:186`): a satellite is
68/// never legitimately at the geocenter.
69const MISSING_POSITION_KM: f64 = 0.0;
70/// SP3 missing velocity component sentinel, in decimeters per second.
71///
72/// Velocity products still carry a `V` record for each `P` record. When no
73/// velocity estimate exists, the record uses the all-zero vector sentinel rather
74/// than being omitted; do not surface that as a fabricated stationary satellite.
75const MISSING_VELOCITY_DM_S: f64 = 0.0;
76
77/// SP3 bad-clock sentinel, in microseconds: `999999.999999`.
78///
79/// A clock value at or above this magnitude means "no clock estimate"; it is
80/// surfaced as `clock_s = None`, not converted.
81const BAD_CLOCK_US: f64 = 999_999.999_999;
82
83/// SP3 velocity records are in decimeters per second; dm/s -> m/s is `* 0.1`.
84const DM_S_TO_M_S: f64 = 1.0e-1;
85/// SP3 clock-rate is in 1e-4 microseconds/second; -> s/s is `* 1e-10`.
86const CLOCK_RATE_TO_S_PER_S: f64 = 1.0e-10;
87
88/// SP3 format version.
89#[derive(Debug, Clone, Copy, PartialEq, Eq)]
90pub enum Sp3Version {
91    /// SP3-a (legacy, GPS-only).
92    A,
93    /// SP3-b.
94    B,
95    /// SP3-c.
96    C,
97    /// SP3-d (multi-GNSS, Hilla 2016).
98    D,
99}
100
101impl Sp3Version {
102    fn from_char(c: char) -> Result<Self> {
103        match c {
104            'a' | 'A' => Ok(Sp3Version::A),
105            'b' | 'B' => Ok(Sp3Version::B),
106            'c' | 'C' => Ok(Sp3Version::C),
107            'd' | 'D' => Ok(Sp3Version::D),
108            other => Err(Error::Parse(format!("unknown SP3 version '{other}'"))),
109        }
110    }
111}
112
113/// What kind of records the file carries.
114#[derive(Debug, Clone, Copy, PartialEq, Eq)]
115pub enum Sp3DataType {
116    /// Position + clock records only (`#?P...`).
117    Position,
118    /// Position + velocity (+ clock + clock-rate) records (`#?V...`).
119    Velocity,
120}
121
122impl Sp3DataType {
123    fn from_char(c: char) -> Result<Self> {
124        match c {
125            'P' => Ok(Sp3DataType::Position),
126            'V' => Ok(Sp3DataType::Velocity),
127            other => Err(Error::Parse(format!("unknown SP3 data type '{other}'"))),
128        }
129    }
130}
131
132/// SP3 time-system labels from the `%c` descriptor.
133///
134/// The core [`TimeScale`] model does not distinguish every SP3 label as its own
135/// global scale. Keep the exact SP3 label here so products using GLONASS, QZSS,
136/// or IRNSS time are accepted and can be serialized without being silently
137/// relabeled.
138#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
139pub enum Sp3TimeSystem {
140    /// GPS time (`GPS`).
141    Gps,
142    /// GLONASS UTC time system (`GLO`).
143    Glonass,
144    /// Galileo system time (`GAL`).
145    Galileo,
146    /// International Atomic Time (`TAI`).
147    Tai,
148    /// Coordinated Universal Time (`UTC`).
149    Utc,
150    /// QZSS time (`QZS`).
151    Qzss,
152    /// BeiDou time (`BDT`).
153    Beidou,
154    /// IRNSS / NavIC time (`IRN`).
155    Irnss,
156}
157
158impl Sp3TimeSystem {
159    /// Canonical three-character SP3 label.
160    pub fn label(self) -> &'static str {
161        match self {
162            Sp3TimeSystem::Gps => "GPS",
163            Sp3TimeSystem::Glonass => "GLO",
164            Sp3TimeSystem::Galileo => "GAL",
165            Sp3TimeSystem::Tai => "TAI",
166            Sp3TimeSystem::Utc => "UTC",
167            Sp3TimeSystem::Qzss => "QZS",
168            Sp3TimeSystem::Beidou => "BDT",
169            Sp3TimeSystem::Irnss => "IRN",
170        }
171    }
172
173    /// Core time scale used to tag parsed [`Instant`] values.
174    ///
175    /// For labels the core model has exactly, this is the direct equivalent. For
176    /// SP3-only labels, the exact product label remains available through
177    /// [`Sp3Header::time_system`], and this value preserves the existing
178    /// interpolation-axis API until the global time model grows those scales.
179    pub fn time_scale(self) -> TimeScale {
180        match self {
181            Sp3TimeSystem::Gps | Sp3TimeSystem::Irnss => TimeScale::Gpst,
182            // QZSST is the exact core scale for the SP3 "QZS" label (nominally
183            // synchronous with GPST); IRNSS has no distinct core scale yet.
184            Sp3TimeSystem::Qzss => TimeScale::Qzsst,
185            Sp3TimeSystem::Glonass | Sp3TimeSystem::Utc => TimeScale::Utc,
186            Sp3TimeSystem::Galileo => TimeScale::Gst,
187            Sp3TimeSystem::Tai => TimeScale::Tai,
188            Sp3TimeSystem::Beidou => TimeScale::Bdt,
189        }
190    }
191
192    fn civil_second_policy(self) -> validate::CivilSecondPolicy {
193        match self {
194            Sp3TimeSystem::Glonass | Sp3TimeSystem::Utc => validate::CivilSecondPolicy::UtcLike,
195            Sp3TimeSystem::Gps
196            | Sp3TimeSystem::Galileo
197            | Sp3TimeSystem::Tai
198            | Sp3TimeSystem::Qzss
199            | Sp3TimeSystem::Beidou
200            | Sp3TimeSystem::Irnss => validate::CivilSecondPolicy::Continuous,
201        }
202    }
203}
204
205/// Per-record quality / status flags (SP3-c columns 75-80, SP3-d same layout).
206///
207/// All four flags are independent and any combination may appear (e.g. a
208/// predicted orbit during a maneuver). They are surfaced verbatim from the
209/// record and never alter the parsed numbers.
210#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
211pub struct Sp3Flags {
212    /// `E` in the clock-event column: a clock discontinuity occurred near this
213    /// epoch; clock interpolation across it is unsafe.
214    pub clock_event: bool,
215    /// `P` in the clock-prediction column: the clock is predicted, not fitted.
216    pub clock_predicted: bool,
217    /// `M` in the maneuver column: the satellite was being maneuvered; the
218    /// state is not suitable for precise navigation.
219    pub maneuver: bool,
220    /// `P` in the orbit-prediction column: the orbit is predicted, not fitted.
221    pub orbit_predicted: bool,
222}
223
224/// A single satellite state at one SP3 epoch.
225///
226/// This is the spec's `Sp3State { position: ItrfPositionM, clock_s, velocity?,
227/// clock_rate?, flags }`. The frame/units are encoded in the
228/// member types; missing optional values are `None` rather than sentinels.
229#[derive(Debug, Clone, Copy, PartialEq)]
230pub struct Sp3State {
231    /// Satellite position in the ITRF/IGS ECEF frame, meters.
232    pub position: ItrfPositionM,
233    /// Satellite clock offset in **seconds** (`None` if the bad-clock sentinel
234    /// `999999.999999` us was recorded).
235    pub clock_s: Option<f64>,
236    /// Satellite velocity in the ITRF/IGS ECEF frame, m/s (present only for
237    /// velocity products).
238    pub velocity: Option<ItrfVelocityMS>,
239    /// Satellite clock rate in **seconds per second** (present only for
240    /// velocity products that carry a clock-rate field).
241    pub clock_rate_s_s: Option<f64>,
242    /// Per-record status flags.
243    pub flags: Sp3Flags,
244}
245
246/// Prediction status aggregated over every satellite record at one SP3 epoch.
247#[derive(Debug, Clone, PartialEq)]
248pub struct Sp3EpochPrediction {
249    /// The parsed epoch.
250    pub epoch: Instant,
251    /// Satellites whose orbit record is marked predicted at this epoch.
252    pub orbit_predicted_satellites: Vec<GnssSatelliteId>,
253    /// Satellites whose clock record is marked predicted at this epoch.
254    pub clock_predicted_satellites: Vec<GnssSatelliteId>,
255}
256
257impl Sp3EpochPrediction {
258    /// True when no position or clock record at this epoch is marked predicted.
259    pub fn is_observed(&self) -> bool {
260        self.orbit_predicted_satellites.is_empty() && self.clock_predicted_satellites.is_empty()
261    }
262}
263
264/// Product-wide observed/predicted metadata derived from SP3 record flags.
265#[derive(Debug, Clone, PartialEq)]
266pub struct Sp3PredictionSummary {
267    /// Per-epoch prediction status in parsed epoch order.
268    pub epochs: Vec<Sp3EpochPrediction>,
269    /// Last epoch before the first epoch containing any predicted record. This
270    /// is the product's truthful contiguous observed-through boundary. It is
271    /// `None` when the first epoch is already predicted or the product is empty.
272    pub observed_through: Option<Instant>,
273}
274
275/// Parsed SP3 header.
276#[derive(Debug, Clone, PartialEq)]
277pub struct Sp3Header {
278    /// SP3 format version (`a`/`b`/`c`/`d`).
279    pub version: Sp3Version,
280    /// Whether the file carries velocity records.
281    pub data_type: Sp3DataType,
282    /// Number of parsed epochs in the canonical product.
283    pub num_epochs: u64,
284    /// Coordinate-system / IGS-realization label (e.g. `IGS14`, `ITRF2`).
285    pub coordinate_system: String,
286    /// Orbit-type label (e.g. `FIT`, `BHN`).
287    pub orbit_type: String,
288    /// Producing agency.
289    pub agency: String,
290    /// GNSS week number (in the file's time system).
291    pub gnss_week: u32,
292    /// Seconds of week of the first epoch.
293    pub seconds_of_week: f64,
294    /// Nominal epoch spacing in seconds.
295    pub epoch_interval_s: f64,
296    /// Modified Julian Day of the first epoch (integer part).
297    pub mjd: u32,
298    /// Fractional day of the first epoch.
299    pub mjd_fraction: f64,
300    /// Time system label the epochs are expressed in. For SP3-b/c/d this is read
301    /// strictly from the first `%c` descriptor (a missing/short/blank descriptor
302    /// is a parse error, never a silent GPST default); SP3-a is implicitly GPST.
303    pub time_system: Sp3TimeSystem,
304    /// Core [`TimeScale`] used to tag parsed [`Instant`] values. See
305    /// [`Sp3Header::time_system`] for the exact SP3 label when the product uses
306    /// a standard SP3 time system that is not modeled as a distinct core scale.
307    pub time_scale: TimeScale,
308    /// The satellite list declared in the `+` header lines.
309    pub satellites: Vec<GnssSatelliteId>,
310    /// Per-satellite accuracy exponent codes from the `++` header lines,
311    /// index-aligned with [`Sp3Header::satellites`].
312    pub satellite_accuracy_codes: Vec<u16>,
313}
314
315/// A parsed SP3 precise-ephemeris product.
316///
317/// Construct with [`Sp3::parse`]. Epochs are stored in file order; exact-product
318/// consumers can use [`validate_exact_sp3`] to require a strictly increasing,
319/// regular requested-cadence grid. Each epoch maps satellite -> [`Sp3State`].
320/// Per-satellite/per-epoch access is via [`Sp3::state`]; arbitrary-epoch
321/// interpolation is built separately to match the parity reference and is not
322/// part of this parser.
323#[derive(Debug, Clone)]
324pub struct Sp3 {
325    /// The parsed header.
326    pub header: Sp3Header,
327    /// Epochs in file order, tagged with the header time scale.
328    pub epochs: Vec<Instant>,
329    /// Epoch count declared on SP3 header line 1. Kept separately from
330    /// [`Sp3Header::num_epochs`], which intentionally remains the number of
331    /// epoch records actually parsed for backward compatibility.
332    declared_num_epochs: u64,
333    /// Start epoch declared on SP3 header line 1, expressed as seconds since
334    /// J2000 in the product time scale. `None` means the legacy permissive
335    /// parser could not interpret those otherwise-unused line-1 fields; exact
336    /// product validation rejects that condition.
337    declared_start_j2000_s: Option<f64>,
338    /// Exact-integrity facts about the logical terminal record. The general
339    /// parser remains permissive; exact validation interprets this state.
340    terminal_record: TerminalRecordState,
341    /// Raw mandatory header-record counts retained for exact validation.
342    satellite_header_lines: usize,
343    accuracy_header_lines: usize,
344    time_system_header_lines: usize,
345    float_header_lines: usize,
346    integer_header_lines: usize,
347    header_comment_lines: usize,
348    /// Raw line-3 satellite count and record sequences retained only for exact
349    /// acquisition validation. Tokens include declarations the typed parser
350    /// cannot represent, so exact validation can still prove count and order.
351    declared_satellite_count: Option<usize>,
352    declared_satellite_tokens: Vec<String>,
353    epoch_position_tokens: Vec<Vec<String>>,
354    epoch_velocity_tokens: Vec<Vec<String>>,
355    epoch_state_record_sequence: Vec<Vec<(char, String)>>,
356    /// Exact seconds since J2000 for each parsed epoch, in the product time
357    /// scale, formed from the epoch record's civil fields with integer
358    /// whole-second arithmetic.
359    epoch_j2000_s: Vec<f64>,
360    /// `epoch_index -> (satellite -> state)`. Parallel to [`Sp3::epochs`].
361    states: Vec<BTreeMap<GnssSatelliteId, Sp3State>>,
362    /// `epoch_index -> (satellite -> native-unit node)`. Parallel to
363    /// [`Sp3::epochs`]; populated **only** from genuine position records. The
364    /// interpolator fits its spline over these (km/us straight from the ASCII,
365    /// exactly as the `scipy`/`gnssanalysis` reference does); reconstructing km
366    /// from the public meters (`km->m->km`) drifts up to 1 ULP and breaks the
367    /// 0-ULP parity. See `sp3/interp.rs`.
368    interp_raw: Vec<BTreeMap<GnssSatelliteId, RawNode>>,
369    /// Free-form `/*` comment lines (notice retained for provenance).
370    pub comments: Vec<String>,
371    /// Count of entries skipped because their satellite token did not parse to a
372    /// representable [`GnssSatelliteId`] (e.g. an extended GLONASS slot like `R28`
373    /// beyond the engine's PRN cap): position/velocity records, plus `+`-header
374    /// satellite declarations. Lets callers tell a clean file
375    /// (`skipped_records == 0`) apart from one carrying unsupported satellites,
376    /// without aborting the whole parse on one such entry. Mirrors
377    /// [`crate::astro::sgp4::TleFile::skipped`].
378    pub skipped_records: usize,
379}
380
381// Preserve the parser's existing canonical-product equality contract. Raw
382// line-1 declarations are retained only as acquisition-integrity evidence and
383// are intentionally excluded: serialization canonicalizes them to the actual
384// body count and first epoch.
385impl PartialEq for Sp3 {
386    fn eq(&self, other: &Self) -> bool {
387        self.header == other.header
388            && self.epochs == other.epochs
389            && self.epoch_j2000_s == other.epoch_j2000_s
390            && self.states == other.states
391            && self.interp_raw == other.interp_raw
392            && self.comments == other.comments
393            && self.skipped_records == other.skipped_records
394    }
395}
396
397/// Native-unit interpolation node: the file's own km / microseconds, kept
398/// verbatim from the ASCII so the spline fit is bit-identical to the reference.
399/// Private - the public surface is meters/seconds via [`Sp3State`].
400#[derive(Debug, Clone, Copy, PartialEq)]
401struct RawNode {
402    /// ECEF position in native SP3 kilometers (X/Y/Z), exact ASCII->f64.
403    km: [f64; 3],
404    /// Clock offset in native SP3 microseconds (`None` for the bad-clock
405    /// sentinel), exact ASCII->f64.
406    clock_us: Option<f64>,
407    /// Whether this epoch carried the clock-event (`E`) flag (clock-arc split).
408    clock_event: bool,
409}
410
411impl Sp3 {
412    /// Parse an SP3-c or SP3-d byte buffer into a typed product.
413    ///
414    /// `bytes` is the full file content (already decompressed; this crate does
415    /// not do gzip - that is a caller-layer I/O concern). Returns
416    /// [`Error::Parse`] with a human-readable reason on malformed input.
417    pub fn parse(bytes: &[u8]) -> Result<Self> {
418        let text = std::str::from_utf8(bytes)
419            .map_err(|e| Error::Parse(format!("SP3 is not valid UTF-8: {e}")))?;
420        Self::parse_str(text)
421    }
422
423    /// Parse from a `&str` (the UTF-8 fast path used by [`Sp3::parse`]).
424    pub fn parse_str(text: &str) -> Result<Self> {
425        if !text.is_ascii() {
426            return Err(Error::Parse("SP3 product text must be ASCII".into()));
427        }
428        let mut parser = Parser::new();
429        for (index, raw) in text.lines().enumerate() {
430            parser.feed(raw, index + 1)?;
431        }
432        parser.finish()
433    }
434
435    /// The satellites present in this product (from the header satellite list).
436    pub fn satellites(&self) -> &[GnssSatelliteId] {
437        &self.header.satellites
438    }
439
440    /// Number of parsed epochs.
441    pub fn epoch_count(&self) -> usize {
442        self.epochs.len()
443    }
444
445    /// Epoch count written in SP3 header line 1.
446    ///
447    /// This can differ from [`Sp3::epoch_count`] when a truncated or otherwise
448    /// inconsistent file is parsed through the deliberately permissive base
449    /// parser. Use [`validate_exact_sp3`] when those fields must agree.
450    pub fn declared_epoch_count(&self) -> u64 {
451        self.declared_num_epochs
452    }
453
454    /// Start epoch written in SP3 header line 1, as seconds since J2000 in the
455    /// product time scale.
456    ///
457    /// The base parser historically ignored these civil fields, so malformed
458    /// values are represented as `None` rather than changing `Sp3::parse`
459    /// compatibility. Exact product validation requires `Some` and checks it
460    /// against both the request and the first parsed epoch.
461    pub fn declared_start_j2000_s(&self) -> Option<f64> {
462        self.declared_start_j2000_s
463    }
464
465    /// The state of `sat` at the parsed epoch with index `epoch_index`.
466    ///
467    /// Returns [`Error::EpochOutOfRange`] if the index is past the end, or
468    /// [`Error::UnknownSatellite`] if the satellite has no record at that epoch.
469    pub fn state(&self, sat: GnssSatelliteId, epoch_index: usize) -> Result<Sp3State> {
470        let per_epoch = self.states.get(epoch_index).ok_or(Error::EpochOutOfRange)?;
471        per_epoch
472            .get(&sat)
473            .copied()
474            .ok_or(Error::UnknownSatellite(sat))
475    }
476
477    /// All `(satellite, state)` pairs recorded at `epoch_index`, in ascending
478    /// satellite order.
479    pub fn states_at(&self, epoch_index: usize) -> Result<&BTreeMap<GnssSatelliteId, Sp3State>> {
480        self.states.get(epoch_index).ok_or(Error::EpochOutOfRange)
481    }
482
483    /// Aggregate the per-record SP3 orbit/clock prediction flags by epoch and
484    /// compute the contiguous observed-through boundary.
485    ///
486    /// This uses the actual `P` flags carried by position records; it never
487    /// assumes a fixed ultra-rapid observed duration. Individual cell flags
488    /// remain available through [`Sp3::state`] and [`Sp3::states_at`].
489    pub fn prediction_summary(&self) -> Sp3PredictionSummary {
490        let epochs: Vec<Sp3EpochPrediction> = self
491            .epochs
492            .iter()
493            .copied()
494            .zip(self.states.iter())
495            .map(|(epoch, states)| Sp3EpochPrediction {
496                epoch,
497                orbit_predicted_satellites: states
498                    .iter()
499                    .filter_map(|(satellite, state)| {
500                        state.flags.orbit_predicted.then_some(*satellite)
501                    })
502                    .collect(),
503                clock_predicted_satellites: states
504                    .iter()
505                    .filter_map(|(satellite, state)| {
506                        state.flags.clock_predicted.then_some(*satellite)
507                    })
508                    .collect(),
509            })
510            .collect();
511        let first_predicted = epochs.iter().position(|epoch| !epoch.is_observed());
512        let observed_through = match first_predicted {
513            Some(0) => None,
514            Some(index) => self.epochs.get(index - 1).copied(),
515            None => self.epochs.last().copied(),
516        };
517
518        Sp3PredictionSummary {
519            epochs,
520            observed_through,
521        }
522    }
523}
524
525impl core::str::FromStr for Sp3 {
526    type Err = Error;
527
528    fn from_str(s: &str) -> Result<Self> {
529        Self::parse_str(s)
530    }
531}
532
533#[cfg(test)]
534impl Sp3 {}
535
536/// Parse an SP3 time-system label.
537///
538/// SP3-c/-d encode the time system in the `%c` descriptor line (chars 9-12).
539/// SP3-a is implicitly GPST. Unknown labels error rather than silently
540/// defaulting, so a parity pipeline never mis-attributes an epoch's scale.
541fn time_system_from_label(label: &str) -> Result<Sp3TimeSystem> {
542    match label.trim() {
543        "GPS" => Ok(Sp3TimeSystem::Gps),
544        "GLO" => Ok(Sp3TimeSystem::Glonass),
545        "GAL" => Ok(Sp3TimeSystem::Galileo),
546        "TAI" => Ok(Sp3TimeSystem::Tai),
547        "UTC" => Ok(Sp3TimeSystem::Utc),
548        "QZS" => Ok(Sp3TimeSystem::Qzss),
549        "BDT" | "BDS" => Ok(Sp3TimeSystem::Beidou),
550        "IRN" => Ok(Sp3TimeSystem::Irnss),
551        trimmed => Err(Error::Parse(format!(
552            "unsupported SP3 time system '{trimmed}'"
553        ))),
554    }
555}
556
557/// Compute the integer-day / fraction split Julian date from a Gregorian UTC-ish
558/// civil epoch, with the day fraction carried separately (Skyfield split
559/// convention, matching [`JulianDateSplit`]).
560///
561/// SP3 epoch lines are civil dates in the file's *own* time system; we keep
562/// them in that scale (no leap-second shifting here - that is a conversion
563/// concern handled by the core `scales` machinery, not the parser). The
564/// algorithm is the standard Fliegel-Van Flandern Gregorian-to-JDN, then the
565/// time-of-day fraction. JDN is computed in integer arithmetic so the whole-day
566/// boundary is exact; only the sub-day fraction is floating point.
567fn civil_to_julian_split(civil: validate::ValidCivil) -> Result<JulianDateSplit> {
568    // Canonical civil-to-split conversion: the integer JDN places the `*.5`
569    // civil-midnight boundary and the within-day clock fields become the
570    // fraction. SP3 epochs are civil days in the file's own scale (no leap
571    // second). The carry below is retained for the rare epoch whose seconds
572    // overflow a day.
573    let (mut jd_whole, mut fraction) = split_julian_date(
574        civil.year as i32,
575        civil.month as i32,
576        civil.day as i32,
577        civil.hour as i32,
578        civil.minute as i32,
579        civil.second,
580    );
581    if fraction > 1.0 {
582        let carry = fraction.floor();
583        jd_whole += carry;
584        fraction -= carry;
585    }
586    JulianDateSplit::new(jd_whole, fraction)
587        .map_err(|error| Error::Parse(format!("invalid SP3 epoch Julian date: {error}")))
588}
589
590/// Conservative maximum for one logical SP3 record.
591///
592/// SP3-d defines the EOF field as `EOF` in columns 1-3 (`A3`) but does not
593/// explicitly prescribe padding after column 3. Official ESA products pad the
594/// record with ASCII spaces through column 80, official GFZ products have also
595/// been observed padding through column 40, and other official products stop
596/// at column 3. Sidereon accepts those interoperable forms while bounding the
597/// record at the conventional 80-column width.
598const SP3_RECORD_WIDTH: usize = 80;
599
600#[derive(Debug, Clone, Copy, PartialEq, Eq)]
601enum EofRecordKind {
602    Valid,
603    Malformed,
604    Other,
605}
606
607#[derive(Debug, Clone, Copy, PartialEq, Eq)]
608struct MalformedEofRecord {
609    line_number: usize,
610    record_length: usize,
611}
612
613#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
614struct TerminalRecordState {
615    had_valid_record: bool,
616    first_malformed_record: Option<MalformedEofRecord>,
617    had_trailing_content: bool,
618}
619
620impl TerminalRecordState {
621    const fn valid() -> Self {
622        Self {
623            had_valid_record: true,
624            first_malformed_record: None,
625            had_trailing_content: false,
626        }
627    }
628}
629
630/// Classify a complete logical record without performing a substring search.
631fn classify_eof_record(line: &str) -> EofRecordKind {
632    if let Some(padding) = line.strip_prefix("EOF") {
633        if line.len() <= SP3_RECORD_WIDTH && padding.bytes().all(|byte| byte == b' ') {
634            EofRecordKind::Valid
635        } else {
636            EofRecordKind::Malformed
637        }
638    } else if line
639        .trim_start_matches(|character: char| character.is_ascii_whitespace())
640        .starts_with("EOF")
641    {
642        EofRecordKind::Malformed
643    } else {
644        EofRecordKind::Other
645    }
646}
647
648/// Incremental line-driven SP3 parser state machine.
649struct Parser {
650    version: Option<Sp3Version>,
651    data_type: Option<Sp3DataType>,
652    num_epochs: u64,
653    declared_start_j2000_s: Option<f64>,
654    coordinate_system: String,
655    orbit_type: String,
656    agency: String,
657    gnss_week: u32,
658    seconds_of_week: f64,
659    epoch_interval_s: f64,
660    mjd: u32,
661    mjd_fraction: f64,
662    time_system: Option<Sp3TimeSystem>,
663    /// `+`-line declared satellites, in file order.
664    sat_list: Vec<GnssSatelliteId>,
665    declared_satellite_count: Option<usize>,
666    declared_satellite_tokens: Vec<String>,
667    /// `++`-line per-satellite accuracy codes, in satellite-list order.
668    sat_accuracy_codes: Vec<u16>,
669    /// Number of real (non-padding) `+`-line satellite slots seen, including any
670    /// dropped because their token was unrepresentable. The `++` accuracy codes
671    /// are positionally aligned with these declaration slots, so this is the axis
672    /// the accuracy parser walks (not the filtered [`Self::sat_list`]).
673    declared_sat_slots: usize,
674    /// Declaration-slot indices (into the `declared_sat_slots` axis) whose token
675    /// was unrepresentable and dropped from [`Self::sat_list`]. Their `++`
676    /// accuracy columns must be skipped so the surviving satellites keep their own
677    /// codes. Empty for every well-formed file, making the accuracy parse a no-op
678    /// realignment in the common case.
679    dropped_sat_slots: Vec<usize>,
680    /// Cursor along the declaration-slot axis consumed by the `++` accuracy
681    /// parser across one or more `++` lines.
682    accuracy_slot_cursor: usize,
683    /// `%c` descriptor lines seen so far (the first carries the time system).
684    pc_count: u32,
685    satellite_header_lines: usize,
686    accuracy_header_lines: usize,
687    float_header_lines: usize,
688    integer_header_lines: usize,
689    header_comment_lines: usize,
690    /// Header line 1 parsed?
691    have_line1: bool,
692    /// Header line 2 parsed?
693    have_line2: bool,
694    /// Epoch currently being filled.
695    current_epoch: Option<Instant>,
696    epochs: Vec<Instant>,
697    epoch_j2000_s: Vec<f64>,
698    states: Vec<BTreeMap<GnssSatelliteId, Sp3State>>,
699    interp_raw: Vec<BTreeMap<GnssSatelliteId, RawNode>>,
700    epoch_position_tokens: Vec<Vec<String>>,
701    epoch_velocity_tokens: Vec<Vec<String>>,
702    epoch_state_record_sequence: Vec<Vec<(char, String)>>,
703    comments: Vec<String>,
704    diagnostics: Diagnostics,
705    terminal_record: TerminalRecordState,
706}
707
708impl Parser {
709    fn new() -> Self {
710        Self {
711            version: None,
712            data_type: None,
713            num_epochs: 0,
714            declared_start_j2000_s: None,
715            coordinate_system: String::new(),
716            orbit_type: String::new(),
717            agency: String::new(),
718            gnss_week: 0,
719            seconds_of_week: 0.0,
720            epoch_interval_s: 0.0,
721            mjd: 0,
722            mjd_fraction: 0.0,
723            time_system: None,
724            sat_list: Vec::new(),
725            declared_satellite_count: None,
726            declared_satellite_tokens: Vec::new(),
727            sat_accuracy_codes: Vec::new(),
728            declared_sat_slots: 0,
729            dropped_sat_slots: Vec::new(),
730            accuracy_slot_cursor: 0,
731            pc_count: 0,
732            satellite_header_lines: 0,
733            accuracy_header_lines: 0,
734            float_header_lines: 0,
735            integer_header_lines: 0,
736            header_comment_lines: 0,
737            have_line1: false,
738            have_line2: false,
739            current_epoch: None,
740            epochs: Vec::new(),
741            epoch_j2000_s: Vec::new(),
742            states: Vec::new(),
743            interp_raw: Vec::new(),
744            epoch_position_tokens: Vec::new(),
745            epoch_velocity_tokens: Vec::new(),
746            epoch_state_record_sequence: Vec::new(),
747            comments: Vec::new(),
748            diagnostics: Diagnostics::new(),
749            terminal_record: TerminalRecordState::default(),
750        }
751    }
752
753    fn feed(&mut self, raw: &str, line_number: usize) -> Result<()> {
754        if self.terminal_record.had_valid_record {
755            if !raw.bytes().all(|byte| byte == b' ') {
756                self.terminal_record.had_trailing_content = true;
757            }
758            return Ok(());
759        }
760        // `str::lines` removes one LF or one CRLF separator. Any remaining CR,
761        // LF, leading whitespace, or other byte is record content; SP3 columns
762        // are significant, so it must not be trimmed here.
763        let line = raw;
764
765        match classify_eof_record(line) {
766            EofRecordKind::Valid => {
767                self.terminal_record.had_valid_record = true;
768                return Ok(());
769            }
770            EofRecordKind::Malformed => {
771                self.terminal_record
772                    .first_malformed_record
773                    .get_or_insert(MalformedEofRecord {
774                        line_number,
775                        record_length: line.len(),
776                    });
777                return Ok(());
778            }
779            EofRecordKind::Other => {}
780        }
781        if line.starts_with("/*") {
782            if self.integer_header_lines >= 2 && self.epochs.is_empty() {
783                self.header_comment_lines += 1;
784            }
785            // Comment line; columns 4.. are the text.
786            // Blank records are mandatory structural padding, not semantic
787            // comments. Count them above but do not add empty strings to the
788            // public `comments` collection.
789            if line.len() > 3 {
790                let comment = line[3..].trim_end();
791                if !comment.is_empty() {
792                    self.comments.push(comment.to_string());
793                }
794            }
795            return Ok(());
796        }
797        // Header line 2 (`##`) must be tested before line 1 (`#`).
798        if line.starts_with("##") {
799            self.parse_line2(line)?;
800            return Ok(());
801        }
802        if line.starts_with('#') {
803            self.parse_line1(line)?;
804            return Ok(());
805        }
806        if line.starts_with('+') {
807            if line.starts_with("++") {
808                self.accuracy_header_lines += 1;
809            } else {
810                self.satellite_header_lines += 1;
811            }
812            self.parse_plus_line(line, line_number)?;
813            return Ok(());
814        }
815        if line.starts_with("%c") {
816            self.parse_pc_line(line)?;
817            return Ok(());
818        }
819        if line.starts_with("%f") {
820            self.float_header_lines += 1;
821            // Float accuracy descriptors are retained only as structural
822            // evidence for exact validation.
823            return Ok(());
824        }
825        if line.starts_with("%i") {
826            self.integer_header_lines += 1;
827            // Float/int accuracy descriptor lines - not needed for the typed
828            // state; skipped deterministically.
829            return Ok(());
830        }
831        if line.starts_with('*') {
832            self.parse_epoch_line(line)?;
833            return Ok(());
834        }
835        if line.starts_with('P') {
836            self.parse_position_line(line, line_number)?;
837            return Ok(());
838        }
839        if line.starts_with('V') {
840            self.parse_velocity_line(line, line_number)?;
841            return Ok(());
842        }
843        // Unknown / ignorable line (e.g. `%/`); skip without failing - SP3 has
844        // optional descriptor lines a parser must tolerate.
845        Ok(())
846    }
847
848    /// Header line 1: `#cP2020 ...` / `#dV...`.
849    fn parse_line1(&mut self, line: &str) -> Result<()> {
850        // Minimum well-formed line-1 length per the standard.
851        if line.len() < 55 {
852            return Err(Error::Parse(format!(
853                "SP3 header line 1 too short: {line:?}"
854            )));
855        }
856        let chars: Vec<char> = line.chars().collect();
857        let version = Sp3Version::from_char(chars[1])?;
858        self.version = Some(version);
859        self.data_type = Some(Sp3DataType::from_char(chars[2])?);
860        // SP3-a predates the %c time-system descriptor and is implicitly GPST.
861        // Set it here so a (correct) SP3-a file with no %c line still resolves,
862        // while SP3-b/c/d are left as None until a valid %c line proves the
863        // scale (a missing %c then becomes a hard error, not a GPST default).
864        if matches!(version, Sp3Version::A) {
865            self.time_system = Some(Sp3TimeSystem::Gps);
866        }
867
868        // Column layout per the SP3 standard, matching the (round-trip-tested)
869        // refs/sp3 line-1 reader: num_epochs 32..40, observables 40..45,
870        // coord_system 45..51, orbit_type 51..55, agency 55...
871        self.num_epochs = field(line, 32, 40)
872            .trim()
873            .parse::<u64>()
874            .map_err(|_| Error::Parse(format!("SP3 num_epochs unparsable in {line:?}")))?;
875        // Keep line-1 start metadata for exact-product validation. These fields
876        // were historically cosmetic to the base parser, so use a best-effort
877        // parse here: malformed values remain parse-compatible but are rejected
878        // by the exact validator as unavailable declared metadata.
879        self.declared_start_j2000_s = parse_declared_start_j2000_s(line);
880        self.coordinate_system = field(line, 45, 51).trim().to_string();
881        self.orbit_type = field(line, 51, 55).trim().to_string();
882        self.agency = field_from(line, 55).trim().to_string();
883        self.have_line1 = true;
884        Ok(())
885    }
886
887    /// Header line 2: `## 2276  21600.00000000   900.00000000 60176 0.25...`.
888    fn parse_line2(&mut self, line: &str) -> Result<()> {
889        self.gnss_week = field(line, 3, 7)
890            .trim()
891            .parse::<u32>()
892            .map_err(|_| Error::Parse(format!("SP3 GNSS week unparsable in {line:?}")))?;
893        self.seconds_of_week = field(line, 8, 23)
894            .trim()
895            .parse::<f64>()
896            .map_err(|_| Error::Parse(format!("SP3 seconds-of-week unparsable in {line:?}")))?;
897        self.epoch_interval_s = field(line, 24, 38)
898            .trim()
899            .parse::<f64>()
900            .map_err(|_| Error::Parse(format!("SP3 epoch interval unparsable in {line:?}")))?;
901        self.mjd = field(line, 39, 44)
902            .trim()
903            .parse::<u32>()
904            .map_err(|_| Error::Parse(format!("SP3 MJD unparsable in {line:?}")))?;
905        self.mjd_fraction = strict_f64(field_from(line, 45), "mjd_fraction")
906            .map_err(|error| map_field_error(error, line))?;
907        self.have_line2 = true;
908        Ok(())
909    }
910
911    /// `+` satellite-list line: `+   32   G01G02...` (3-char SV tokens from
912    /// column 9 in groups of 17). Continuation `+` lines append more tokens.
913    fn parse_plus_line(&mut self, line: &str, line_number: usize) -> Result<()> {
914        if line.starts_with("++") {
915            return self.parse_accuracy_line(line);
916        }
917        if self.satellite_header_lines == 1 {
918            self.declared_satellite_count = field(line, 3, 6).trim().parse::<usize>().ok();
919        }
920        // SV tokens start at column 9 (0-based), each 3 chars, up to 17 per line.
921        let mut col = 9;
922        while col + 3 <= line.len() {
923            let token = field(line, col, col + 3);
924            let trimmed = token.trim();
925            // Unused satellite slots are zero-filled, not a declaration. The
926            // SP3 zero-fill varies between producers (`  0`, ` 00`, `000`), so
927            // any all-zero (or blank) token is padding - never a satellite,
928            // whose token is a system letter + PRN (or, in SP3-a, a non-zero
929            // numeric PRN). Misreading ` 00` as an unrepresentable satellite
930            // inflates `skipped_records` and breaks the parse/write/parse
931            // round trip (the writer re-emits the canonical `  0`).
932            if trimmed.is_empty() || trimmed.bytes().all(|b| b == b'0') {
933                col += 3;
934                continue;
935            }
936            // This is a real declaration slot; the `++` accuracy codes are aligned
937            // to this axis, so track its index whether or not the token resolves.
938            let slot_index = self.declared_sat_slots;
939            self.declared_sat_slots += 1;
940            self.declared_satellite_tokens.push(trimmed.to_owned());
941            if let Some(id) = parse_sv_token(token, self.version) {
942                if !self.sat_list.contains(&id) {
943                    self.sat_list.push(id);
944                }
945            } else {
946                // A declared satellite whose token is not representable (e.g. an
947                // extended GLONASS slot R28 beyond the engine's PRN cap) is
948                // dropped from the satellite list, but counted rather than dropped
949                // silently - consistent with the position/velocity record paths
950                // (see `Sp3::skipped_records`). Record the slot so its accuracy
951                // column is skipped, keeping the surviving codes aligned.
952                self.push_unrepresentable_satellite_skip(line_number, token);
953                self.dropped_sat_slots.push(slot_index);
954            }
955            col += 3;
956        }
957        Ok(())
958    }
959
960    /// `++` per-satellite accuracy-code line: 3-char integer fields from column
961    /// 9, aligned with the `+` declaration slots.
962    ///
963    /// The columns track the `+` declaration order, so a column whose declaration
964    /// slot was dropped (an unrepresentable satellite) is read and discarded, not
965    /// pushed - otherwise the surviving satellites would inherit a neighbour's
966    /// accuracy code. With no dropped slots this is exactly the 1:1 push as before.
967    fn parse_accuracy_line(&mut self, line: &str) -> Result<()> {
968        let mut col = 9;
969        while col + 3 <= line.len() && self.accuracy_slot_cursor < self.declared_sat_slots {
970            let token = field(line, col, col + 3);
971            let trimmed = token.trim();
972            let code = if trimmed.is_empty() {
973                0
974            } else {
975                validate::strict_int::<u16>(trimmed, "satellite_accuracy_code")
976                    .map_err(|error| map_field_error(error, line))?
977            };
978            if !self.dropped_sat_slots.contains(&self.accuracy_slot_cursor) {
979                self.sat_accuracy_codes.push(code);
980            }
981            self.accuracy_slot_cursor += 1;
982            col += 3;
983        }
984        Ok(())
985    }
986
987    /// `%c` descriptor: the first one (chars 9-12) carries the time system.
988    fn parse_pc_line(&mut self, line: &str) -> Result<()> {
989        if self.pc_count == 0 {
990            // SP3-a is implicitly GPST regardless of descriptor content.
991            if matches!(self.version, Some(Sp3Version::A)) {
992                self.time_system = Some(Sp3TimeSystem::Gps);
993            } else if line.len() >= 12 {
994                let label = field(line, 9, 12);
995                let trimmed = label.trim();
996                // STRICT: a blank time-system field on the first %c is not GPST,
997                // it is malformed. Reject rather than silently defaulting so a
998                // precise pipeline never mis-attributes an epoch's scale.
999                if trimmed.is_empty() {
1000                    return Err(Error::Parse(format!(
1001                        "SP3 %c time system is blank in {line:?}"
1002                    )));
1003                }
1004                self.time_system = Some(time_system_from_label(label)?);
1005            } else {
1006                // STRICT: a short %c line for SP3-b/c/d carries no time system
1007                // we can trust. Reject rather than defaulting to GPST.
1008                return Err(Error::Parse(format!(
1009                    "SP3 %c descriptor too short to carry a time system: {line:?}"
1010                )));
1011            }
1012        }
1013        self.pc_count += 1;
1014        Ok(())
1015    }
1016
1017    /// Epoch line: `*  2020  6 24  0  0  0.00000000`.
1018    fn parse_epoch_line(&mut self, line: &str) -> Result<()> {
1019        // STRICT: by the time we reach data, the time system must be known -
1020        // implicitly GPST for SP3-a (set at line 1), or from a valid first %c
1021        // line for SP3-b/c/d. A missing/blank/short %c is an error, never GPST.
1022        let time_system = self.time_system.ok_or_else(|| {
1023            Error::Parse("SP3 epoch encountered with no time system (missing %c descriptor)".into())
1024        })?;
1025        let scale = time_system.time_scale();
1026        // Fields after the leading `*  ` (3 chars), then space-delimited.
1027        let body = &line[1..];
1028        let mut it = body.split_whitespace();
1029        let year: i64 = next_field(&mut it, "epoch year")?;
1030        let month: i64 = next_field(&mut it, "epoch month")?;
1031        let day: i64 = next_field(&mut it, "epoch day")?;
1032        let hour: i64 = next_field(&mut it, "epoch hour")?;
1033        let minute: i64 = next_field(&mut it, "epoch minute")?;
1034        let seconds: f64 = next_field(&mut it, "epoch seconds")?;
1035
1036        let civil = validate::civil_datetime_with_second_policy(
1037            year,
1038            month,
1039            day,
1040            hour,
1041            minute,
1042            seconds,
1043            time_system.civil_second_policy(),
1044        )
1045        .map_err(|error| map_field_error(error, line))?;
1046        let split = civil_to_julian_split(civil)?;
1047        let epoch_j2000_s = j2000_seconds(
1048            civil.year as i32,
1049            civil.month as i32,
1050            civil.day as i32,
1051            civil.hour as i32,
1052            civil.minute as i32,
1053            civil.second,
1054        );
1055        let epoch = Instant {
1056            scale,
1057            repr: InstantRepr::JulianDate(split),
1058        };
1059        self.epochs.push(epoch);
1060        self.epoch_j2000_s.push(epoch_j2000_s);
1061        self.states.push(BTreeMap::new());
1062        self.interp_raw.push(BTreeMap::new());
1063        self.epoch_position_tokens.push(Vec::new());
1064        self.epoch_velocity_tokens.push(Vec::new());
1065        self.epoch_state_record_sequence.push(Vec::new());
1066        self.current_epoch = Some(epoch);
1067        Ok(())
1068    }
1069
1070    /// Position+clock record: `PG01  x  y  z  clk  ...flags`.
1071    fn parse_position_line(&mut self, line: &str, line_number: usize) -> Result<()> {
1072        if self.current_epoch.is_none() {
1073            return Err(Error::Parse(
1074                "SP3 position record before any epoch line".into(),
1075            ));
1076        }
1077        if line.len() < 46 {
1078            return Err(Error::Parse(format!(
1079                "SP3 position record truncated before vector fields in {line:?}"
1080            )));
1081        }
1082        let token = field(line, 1, 4);
1083        self.epoch_position_tokens
1084            .last_mut()
1085            .expect("current epoch has a raw position-token list")
1086            .push(token.trim().to_owned());
1087        self.epoch_state_record_sequence
1088            .last_mut()
1089            .expect("current epoch has a raw state-record sequence")
1090            .push(('P', token.trim().to_owned()));
1091        let Some(sat) = parse_sv_token(token, self.version) else {
1092            // A token that does not parse to a representable `GnssSatelliteId`
1093            // (e.g. an extended GLONASS slot like R28 beyond the engine's PRN
1094            // cap) is an independent, unsupported record. One such record must
1095            // not reject the whole file - skip and count it, mirroring nav
1096            // `parse_glonass` and `parse_tle_file`.
1097            self.push_unrepresentable_satellite_skip(line_number, token);
1098            return Ok(());
1099        };
1100
1101        // The header `+` lines are the authoritative satellite declaration; a
1102        // position record for an undeclared satellite is malformed. Accepting it
1103        // would store a state the writer (which emits only declared satellites)
1104        // cannot reproduce, breaking parse/encode/parse round-tripping.
1105        if !self.sat_list.contains(&sat) {
1106            return Err(Error::Parse(format!(
1107                "SP3 position record for satellite {token:?} not in the header satellite list"
1108            )));
1109        }
1110
1111        let x_km = parse_coord(line, 4, 18)?;
1112        let y_km = parse_coord(line, 18, 32)?;
1113        let z_km = parse_coord(line, 32, 46)?;
1114
1115        // All-zero position is the missing-orbit sentinel: skip the record.
1116        if x_km == MISSING_POSITION_KM && y_km == MISSING_POSITION_KM && z_km == MISSING_POSITION_KM
1117        {
1118            return Ok(());
1119        }
1120
1121        let clock_us = parse_clock_us(line)?;
1122        let clock_s = clock_us.map(|us| us * US_TO_S);
1123
1124        let flags = parse_flags(line);
1125
1126        let position = ItrfPositionM::new(x_km * KM_TO_M, y_km * KM_TO_M, z_km * KM_TO_M)
1127            .map_err(|e| Error::Parse(format!("SP3 invalid position record: {e}")))?;
1128        let state = Sp3State {
1129            position,
1130            clock_s,
1131            velocity: None,
1132            clock_rate_s_s: None,
1133            flags,
1134        };
1135        let idx = self.states.len() - 1;
1136        self.states[idx].insert(sat, state);
1137        // Keep the native-unit node for the interpolation path (see RawNode):
1138        // the spline must fit the file's own km/us, not the km->m->km round trip.
1139        self.interp_raw[idx].insert(
1140            sat,
1141            RawNode {
1142                km: [x_km, y_km, z_km],
1143                clock_us,
1144                clock_event: flags.clock_event,
1145            },
1146        );
1147        Ok(())
1148    }
1149
1150    /// Velocity record: `VG01  vx  vy  vz  clkrate ...`. Augments the matching
1151    /// position record at the current epoch (must follow it).
1152    fn parse_velocity_line(&mut self, line: &str, line_number: usize) -> Result<()> {
1153        if self.current_epoch.is_none() {
1154            return Err(Error::Parse(
1155                "SP3 velocity record before any epoch line".into(),
1156            ));
1157        }
1158        if line.len() < 46 {
1159            return Err(Error::Parse(format!(
1160                "SP3 velocity record truncated before vector fields in {line:?}"
1161            )));
1162        }
1163        let token = field(line, 1, 4);
1164        self.epoch_velocity_tokens
1165            .last_mut()
1166            .expect("current epoch has a raw velocity-token list")
1167            .push(token.trim().to_owned());
1168        self.epoch_state_record_sequence
1169            .last_mut()
1170            .expect("current epoch has a raw state-record sequence")
1171            .push(('V', token.trim().to_owned()));
1172        let Some(sat) = parse_sv_token(token, self.version) else {
1173            // Unparsable / out-of-range satellite token: skip and count, same
1174            // as the position-record path above.
1175            self.push_unrepresentable_satellite_skip(line_number, token);
1176            return Ok(());
1177        };
1178
1179        // SP3 velocity is in dm/s; read each axis independently (the refs/sp3
1180        // crate has a bug here that reuses Y for X - we do not).
1181        let vx_dm_s = parse_coord(line, 4, 18)?;
1182        let vy_dm_s = parse_coord(line, 18, 32)?;
1183        let vz_dm_s = parse_coord(line, 32, 46)?;
1184
1185        let missing_velocity = vx_dm_s == MISSING_VELOCITY_DM_S
1186            && vy_dm_s == MISSING_VELOCITY_DM_S
1187            && vz_dm_s == MISSING_VELOCITY_DM_S;
1188        let velocity = ItrfVelocityMS::new(
1189            vx_dm_s * DM_S_TO_M_S,
1190            vy_dm_s * DM_S_TO_M_S,
1191            vz_dm_s * DM_S_TO_M_S,
1192        )
1193        .map_err(|e| Error::Parse(format!("SP3 invalid velocity record: {e}")))?;
1194
1195        // Clock-rate field shares the clock column; bad-clock sentinel applies.
1196        let clock_rate_s_s = parse_clock_us(line)?.map(|rate| rate * CLOCK_RATE_TO_S_PER_S);
1197
1198        let idx = self.states.len() - 1;
1199        match self.states[idx].get_mut(&sat) {
1200            Some(state) if !missing_velocity => {
1201                state.velocity = Some(velocity);
1202                state.clock_rate_s_s = clock_rate_s_s;
1203            }
1204            Some(_) => {}
1205            None => {
1206                // A V-record always follows its P-record for the same satellite
1207                // at the same epoch (SP3 format invariant). With no preceding
1208                // P-record this satellite has NO valid position at this epoch;
1209                // synthesizing one (e.g. the geocenter (0,0,0)) would fabricate
1210                // an orbit that the all-zero missing-orbit guard exists to
1211                // reject, and would leak through the public state()/states_at().
1212                // Treat it as malformed and skip - consistent with the parser's
1213                // tolerant skipping of other malformed records. No state is
1214                // inserted, so the satellite stays UnknownSatellite at this
1215                // epoch and no (0,0,0) position is ever exposed.
1216            }
1217        }
1218        Ok(())
1219    }
1220
1221    fn push_unrepresentable_satellite_skip(&mut self, line_number: usize, token: &str) {
1222        self.diagnostics.push_skip(Skip {
1223            at: RecordRef::at_line(line_number).with_satellite(token),
1224            reason: SkipReason::UnrepresentableSatellite,
1225        });
1226    }
1227
1228    fn finish(self) -> Result<Sp3> {
1229        if !self.have_line1 {
1230            return Err(Error::Parse("SP3 missing header line 1".into()));
1231        }
1232        if !self.have_line2 {
1233            return Err(Error::Parse("SP3 missing header line 2".into()));
1234        }
1235        let version = self
1236            .version
1237            .ok_or_else(|| Error::Parse("SP3 version not determined".into()))?;
1238        let data_type = self
1239            .data_type
1240            .ok_or_else(|| Error::Parse("SP3 data type not determined".into()))?;
1241        // STRICT: SP3-a is implicitly GPST (set at line 1); SP3-b/c/d must have
1242        // proved their scale from a valid first %c line. Never default here.
1243        let time_system = self.time_system.ok_or_else(|| {
1244            Error::Parse(
1245                "SP3 time system not determined (missing/short/blank %c descriptor)".into(),
1246            )
1247        })?;
1248        let time_scale = time_system.time_scale();
1249
1250        let mut satellite_accuracy_codes = self.sat_accuracy_codes;
1251        satellite_accuracy_codes.truncate(self.sat_list.len());
1252        satellite_accuracy_codes.resize(self.sat_list.len(), 0);
1253        let skipped_records = self.diagnostics.skips.len();
1254
1255        let header = Sp3Header {
1256            version,
1257            data_type,
1258            num_epochs: self.epochs.len() as u64,
1259            coordinate_system: self.coordinate_system,
1260            orbit_type: self.orbit_type,
1261            agency: self.agency,
1262            gnss_week: self.gnss_week,
1263            seconds_of_week: self.seconds_of_week,
1264            epoch_interval_s: self.epoch_interval_s,
1265            mjd: self.mjd,
1266            mjd_fraction: self.mjd_fraction,
1267            time_system,
1268            time_scale,
1269            satellites: self.sat_list,
1270            satellite_accuracy_codes,
1271        };
1272
1273        Ok(Sp3 {
1274            header,
1275            epochs: self.epochs,
1276            declared_num_epochs: self.num_epochs,
1277            declared_start_j2000_s: self.declared_start_j2000_s,
1278            terminal_record: self.terminal_record,
1279            satellite_header_lines: self.satellite_header_lines,
1280            accuracy_header_lines: self.accuracy_header_lines,
1281            time_system_header_lines: self.pc_count as usize,
1282            float_header_lines: self.float_header_lines,
1283            integer_header_lines: self.integer_header_lines,
1284            header_comment_lines: self.header_comment_lines,
1285            declared_satellite_count: self.declared_satellite_count,
1286            declared_satellite_tokens: self.declared_satellite_tokens,
1287            epoch_position_tokens: self.epoch_position_tokens,
1288            epoch_velocity_tokens: self.epoch_velocity_tokens,
1289            epoch_state_record_sequence: self.epoch_state_record_sequence,
1290            epoch_j2000_s: self.epoch_j2000_s,
1291            states: self.states,
1292            interp_raw: self.interp_raw,
1293            comments: self.comments,
1294            skipped_records,
1295        })
1296    }
1297}
1298
1299/// Best-effort parse of the civil start epoch carried on SP3 header line 1.
1300///
1301/// Exact validation treats `None` as an integrity failure. Keeping this helper
1302/// non-fallible preserves the long-standing permissive behavior of `Sp3::parse`
1303/// for callers that only consume epoch records.
1304fn parse_declared_start_j2000_s(line: &str) -> Option<f64> {
1305    let year = field(line, 3, 7).trim().parse::<i64>().ok()?;
1306    let month = field(line, 8, 10).trim().parse::<i64>().ok()?;
1307    let day = field(line, 11, 13).trim().parse::<i64>().ok()?;
1308    let hour = field(line, 14, 16).trim().parse::<i64>().ok()?;
1309    let minute = field(line, 17, 19).trim().parse::<i64>().ok()?;
1310    let second = field(line, 20, 31).trim().parse::<f64>().ok()?;
1311    let civil = validate::civil_datetime_with_second_policy(
1312        year,
1313        month,
1314        day,
1315        hour,
1316        minute,
1317        second,
1318        validate::CivilSecondPolicy::UtcLike,
1319    )
1320    .ok()?;
1321    Some(j2000_seconds(
1322        civil.year as i32,
1323        civil.month as i32,
1324        civil.day as i32,
1325        civil.hour as i32,
1326        civil.minute as i32,
1327        civil.second,
1328    ))
1329}
1330
1331/// Parse a fixed-column float coordinate, mapping failures to a parse error
1332/// that names the offending text.
1333fn parse_coord(line: &str, start: usize, end: usize) -> Result<f64> {
1334    let raw = field(line, start, end).trim();
1335    strict_f64(raw, "coordinate").map_err(|error| map_field_error(error, line))
1336}
1337
1338/// Parse the clock column (chars 46..60). Returns `None` for the bad-clock
1339/// sentinel `999999.999999` or an absent/blank field; `Some(us)` otherwise.
1340fn parse_clock_us(line: &str) -> Result<Option<f64>> {
1341    if line.len() <= 46 {
1342        return Ok(None);
1343    }
1344    let raw = field(line, 46, 60).trim();
1345    if raw.is_empty() {
1346        return Ok(None);
1347    }
1348    let value = strict_f64(raw, "clock").map_err(|error| map_field_error(error, line))?;
1349    // Sentinel: any value at or beyond the bad-clock magnitude is "no estimate".
1350    if value.abs() >= BAD_CLOCK_US {
1351        return Ok(None);
1352    }
1353    Ok(Some(value))
1354}
1355
1356fn map_field_error(error: validate::FieldError, line: &str) -> Error {
1357    Error::Parse(format!("SP3 {error} in {line:?}"))
1358}
1359
1360/// Parse the four status flags from their fixed columns (SP3-c/-d shared
1361/// layout): clock-event col 74 = `E`, clock-prediction col 75 = `P`,
1362/// maneuver col 78 = `M`, orbit-prediction col 79 = `P`.
1363fn parse_flags(line: &str) -> Sp3Flags {
1364    let at = |col: usize, want: char| -> bool { char_at(line, col) == Some(want) };
1365    Sp3Flags {
1366        clock_event: at(74, 'E'),
1367        clock_predicted: at(75, 'P'),
1368        maneuver: at(78, 'M'),
1369        orbit_predicted: at(79, 'P'),
1370    }
1371}
1372
1373/// Parse a 3-char SV token (e.g. `G01`, `C30`, or a bare `  1` in SP3-a) into a
1374/// [`GnssSatelliteId`]. Returns `None` on an unrecognized token.
1375fn parse_sv_token(token: &str, version: Option<Sp3Version>) -> Option<GnssSatelliteId> {
1376    let token = token.trim();
1377    if token.is_empty() {
1378        return None;
1379    }
1380    let first = token.chars().next()?;
1381    if first.is_ascii_digit() {
1382        // SP3-a GPS-only: bare numeric PRN, optionally space-padded.
1383        if matches!(version, Some(Sp3Version::A)) || version.is_none() {
1384            let prn = token.parse::<u8>().ok()?;
1385            if !is_valid_prn(GnssSystem::Gps, prn) {
1386                return None;
1387            }
1388            return GnssSatelliteId::new(GnssSystem::Gps, prn).ok();
1389        }
1390        return None;
1391    }
1392    token.parse::<GnssSatelliteId>().ok()
1393}
1394
1395/// Pull and parse the next whitespace-delimited field from an iterator.
1396fn next_field<T: std::str::FromStr>(
1397    it: &mut std::str::SplitWhitespace<'_>,
1398    what: &str,
1399) -> Result<T> {
1400    let tok = it
1401        .next()
1402        .ok_or_else(|| Error::Parse(format!("SP3 missing {what}")))?;
1403    tok.parse::<T>()
1404        .map_err(|_| Error::Parse(format!("SP3 {what} {tok:?} unparsable")))
1405}
1406
1407mod combine;
1408mod exact;
1409mod interp;
1410mod interpolant;
1411mod interpolant_store;
1412mod provenance;
1413mod samples;
1414mod verify;
1415mod write;
1416
1417pub use combine::{
1418    align_clock_reference, clock_reference_offset, merge, AgreementMetric, ClockReferenceOffset,
1419    EpochAgreement, MergeCombine, MergeFlag, MergeOptions, MergePrecedenceScope, MergeReport,
1420    OutlierRejectOptions, Sp3FrameLabelSet, Sp3FrameReconciliation, Sp3FrameReconciliationMethod,
1421    Sp3FrameReconciliationOptions,
1422};
1423pub use exact::{
1424    parse_exact_sp3, validate_exact_sp3, ExactSp3Coverage, ExactSp3Request, ExactSp3ValidationError,
1425};
1426pub use interpolant::{PreciseEphemerisInterpolant, PreciseInterpolantError};
1427pub use interpolant_store::{
1428    precise_interpolant_store_checksum64, MmapPreciseEphemerisInterpolant,
1429    PreciseInterpolantStoreError,
1430};
1431pub use provenance::{
1432    Sp3ArtifactIdentity, Sp3MergeInputIdentity, Sp3MergeInputIdentityError,
1433    SP3_MERGE_INPUT_ID_PREFIX, SP3_MERGE_INPUT_SCHEMA_VERSION,
1434};
1435pub use samples::{
1436    sp3_ecef_state_to_eci, PreciseEphemerisSample, PreciseEphemerisSamples,
1437    PreciseEphemerisStateSample, PreciseSamplesError,
1438};
1439pub use verify::{
1440    compare_position_series, InterpolationComparison, InterpolationDivergence, ReferenceState,
1441};
1442
1443#[cfg(all(test, sidereon_repo_tests))]
1444mod tests;