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    /// Whether the parser encountered the mandatory terminal record.
339    had_eof: bool,
340    trailing_content_after_eof: bool,
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/// Incremental line-driven SP3 parser state machine.
591struct Parser {
592    version: Option<Sp3Version>,
593    data_type: Option<Sp3DataType>,
594    num_epochs: u64,
595    declared_start_j2000_s: Option<f64>,
596    coordinate_system: String,
597    orbit_type: String,
598    agency: String,
599    gnss_week: u32,
600    seconds_of_week: f64,
601    epoch_interval_s: f64,
602    mjd: u32,
603    mjd_fraction: f64,
604    time_system: Option<Sp3TimeSystem>,
605    /// `+`-line declared satellites, in file order.
606    sat_list: Vec<GnssSatelliteId>,
607    declared_satellite_count: Option<usize>,
608    declared_satellite_tokens: Vec<String>,
609    /// `++`-line per-satellite accuracy codes, in satellite-list order.
610    sat_accuracy_codes: Vec<u16>,
611    /// Number of real (non-padding) `+`-line satellite slots seen, including any
612    /// dropped because their token was unrepresentable. The `++` accuracy codes
613    /// are positionally aligned with these declaration slots, so this is the axis
614    /// the accuracy parser walks (not the filtered [`Self::sat_list`]).
615    declared_sat_slots: usize,
616    /// Declaration-slot indices (into the `declared_sat_slots` axis) whose token
617    /// was unrepresentable and dropped from [`Self::sat_list`]. Their `++`
618    /// accuracy columns must be skipped so the surviving satellites keep their own
619    /// codes. Empty for every well-formed file, making the accuracy parse a no-op
620    /// realignment in the common case.
621    dropped_sat_slots: Vec<usize>,
622    /// Cursor along the declaration-slot axis consumed by the `++` accuracy
623    /// parser across one or more `++` lines.
624    accuracy_slot_cursor: usize,
625    /// `%c` descriptor lines seen so far (the first carries the time system).
626    pc_count: u32,
627    satellite_header_lines: usize,
628    accuracy_header_lines: usize,
629    float_header_lines: usize,
630    integer_header_lines: usize,
631    header_comment_lines: usize,
632    /// Header line 1 parsed?
633    have_line1: bool,
634    /// Header line 2 parsed?
635    have_line2: bool,
636    /// Epoch currently being filled.
637    current_epoch: Option<Instant>,
638    epochs: Vec<Instant>,
639    epoch_j2000_s: Vec<f64>,
640    states: Vec<BTreeMap<GnssSatelliteId, Sp3State>>,
641    interp_raw: Vec<BTreeMap<GnssSatelliteId, RawNode>>,
642    epoch_position_tokens: Vec<Vec<String>>,
643    epoch_velocity_tokens: Vec<Vec<String>>,
644    epoch_state_record_sequence: Vec<Vec<(char, String)>>,
645    comments: Vec<String>,
646    diagnostics: Diagnostics,
647    done: bool,
648    had_eof: bool,
649    trailing_content_after_eof: bool,
650}
651
652impl Parser {
653    fn new() -> Self {
654        Self {
655            version: None,
656            data_type: None,
657            num_epochs: 0,
658            declared_start_j2000_s: None,
659            coordinate_system: String::new(),
660            orbit_type: String::new(),
661            agency: String::new(),
662            gnss_week: 0,
663            seconds_of_week: 0.0,
664            epoch_interval_s: 0.0,
665            mjd: 0,
666            mjd_fraction: 0.0,
667            time_system: None,
668            sat_list: Vec::new(),
669            declared_satellite_count: None,
670            declared_satellite_tokens: Vec::new(),
671            sat_accuracy_codes: Vec::new(),
672            declared_sat_slots: 0,
673            dropped_sat_slots: Vec::new(),
674            accuracy_slot_cursor: 0,
675            pc_count: 0,
676            satellite_header_lines: 0,
677            accuracy_header_lines: 0,
678            float_header_lines: 0,
679            integer_header_lines: 0,
680            header_comment_lines: 0,
681            have_line1: false,
682            have_line2: false,
683            current_epoch: None,
684            epochs: Vec::new(),
685            epoch_j2000_s: Vec::new(),
686            states: Vec::new(),
687            interp_raw: Vec::new(),
688            epoch_position_tokens: Vec::new(),
689            epoch_velocity_tokens: Vec::new(),
690            epoch_state_record_sequence: Vec::new(),
691            comments: Vec::new(),
692            diagnostics: Diagnostics::new(),
693            done: false,
694            had_eof: false,
695            trailing_content_after_eof: false,
696        }
697    }
698
699    fn feed(&mut self, raw: &str, line_number: usize) -> Result<()> {
700        if self.done {
701            if !raw.trim().is_empty() {
702                self.trailing_content_after_eof = true;
703            }
704            return Ok(());
705        }
706        // SP3 is fixed-column ASCII; trim only the trailing CR / newline noise,
707        // never leading spaces (columns are significant).
708        let line = raw.trim_end_matches(['\r', '\n']);
709
710        if line == "EOF" {
711            self.done = true;
712            self.had_eof = true;
713            return Ok(());
714        }
715        if line.starts_with("/*") {
716            if self.integer_header_lines >= 2 && self.epochs.is_empty() {
717                self.header_comment_lines += 1;
718            }
719            // Comment line; columns 4.. are the text.
720            // Blank records are mandatory structural padding, not semantic
721            // comments. Count them above but do not add empty strings to the
722            // public `comments` collection.
723            if line.len() > 3 {
724                let comment = line[3..].trim_end();
725                if !comment.is_empty() {
726                    self.comments.push(comment.to_string());
727                }
728            }
729            return Ok(());
730        }
731        // Header line 2 (`##`) must be tested before line 1 (`#`).
732        if line.starts_with("##") {
733            self.parse_line2(line)?;
734            return Ok(());
735        }
736        if line.starts_with('#') {
737            self.parse_line1(line)?;
738            return Ok(());
739        }
740        if line.starts_with('+') {
741            if line.starts_with("++") {
742                self.accuracy_header_lines += 1;
743            } else {
744                self.satellite_header_lines += 1;
745            }
746            self.parse_plus_line(line, line_number)?;
747            return Ok(());
748        }
749        if line.starts_with("%c") {
750            self.parse_pc_line(line)?;
751            return Ok(());
752        }
753        if line.starts_with("%f") {
754            self.float_header_lines += 1;
755            // Float accuracy descriptors are retained only as structural
756            // evidence for exact validation.
757            return Ok(());
758        }
759        if line.starts_with("%i") {
760            self.integer_header_lines += 1;
761            // Float/int accuracy descriptor lines - not needed for the typed
762            // state; skipped deterministically.
763            return Ok(());
764        }
765        if line.starts_with('*') {
766            self.parse_epoch_line(line)?;
767            return Ok(());
768        }
769        if line.starts_with('P') {
770            self.parse_position_line(line, line_number)?;
771            return Ok(());
772        }
773        if line.starts_with('V') {
774            self.parse_velocity_line(line, line_number)?;
775            return Ok(());
776        }
777        // Unknown / ignorable line (e.g. `%/`); skip without failing - SP3 has
778        // optional descriptor lines a parser must tolerate.
779        Ok(())
780    }
781
782    /// Header line 1: `#cP2020 ...` / `#dV...`.
783    fn parse_line1(&mut self, line: &str) -> Result<()> {
784        // Minimum well-formed line-1 length per the standard.
785        if line.len() < 55 {
786            return Err(Error::Parse(format!(
787                "SP3 header line 1 too short: {line:?}"
788            )));
789        }
790        let chars: Vec<char> = line.chars().collect();
791        let version = Sp3Version::from_char(chars[1])?;
792        self.version = Some(version);
793        self.data_type = Some(Sp3DataType::from_char(chars[2])?);
794        // SP3-a predates the %c time-system descriptor and is implicitly GPST.
795        // Set it here so a (correct) SP3-a file with no %c line still resolves,
796        // while SP3-b/c/d are left as None until a valid %c line proves the
797        // scale (a missing %c then becomes a hard error, not a GPST default).
798        if matches!(version, Sp3Version::A) {
799            self.time_system = Some(Sp3TimeSystem::Gps);
800        }
801
802        // Column layout per the SP3 standard, matching the (round-trip-tested)
803        // refs/sp3 line-1 reader: num_epochs 32..40, observables 40..45,
804        // coord_system 45..51, orbit_type 51..55, agency 55...
805        self.num_epochs = field(line, 32, 40)
806            .trim()
807            .parse::<u64>()
808            .map_err(|_| Error::Parse(format!("SP3 num_epochs unparsable in {line:?}")))?;
809        // Keep line-1 start metadata for exact-product validation. These fields
810        // were historically cosmetic to the base parser, so use a best-effort
811        // parse here: malformed values remain parse-compatible but are rejected
812        // by the exact validator as unavailable declared metadata.
813        self.declared_start_j2000_s = parse_declared_start_j2000_s(line);
814        self.coordinate_system = field(line, 45, 51).trim().to_string();
815        self.orbit_type = field(line, 51, 55).trim().to_string();
816        self.agency = field_from(line, 55).trim().to_string();
817        self.have_line1 = true;
818        Ok(())
819    }
820
821    /// Header line 2: `## 2276  21600.00000000   900.00000000 60176 0.25...`.
822    fn parse_line2(&mut self, line: &str) -> Result<()> {
823        self.gnss_week = field(line, 3, 7)
824            .trim()
825            .parse::<u32>()
826            .map_err(|_| Error::Parse(format!("SP3 GNSS week unparsable in {line:?}")))?;
827        self.seconds_of_week = field(line, 8, 23)
828            .trim()
829            .parse::<f64>()
830            .map_err(|_| Error::Parse(format!("SP3 seconds-of-week unparsable in {line:?}")))?;
831        self.epoch_interval_s = field(line, 24, 38)
832            .trim()
833            .parse::<f64>()
834            .map_err(|_| Error::Parse(format!("SP3 epoch interval unparsable in {line:?}")))?;
835        self.mjd = field(line, 39, 44)
836            .trim()
837            .parse::<u32>()
838            .map_err(|_| Error::Parse(format!("SP3 MJD unparsable in {line:?}")))?;
839        self.mjd_fraction = strict_f64(field_from(line, 45), "mjd_fraction")
840            .map_err(|error| map_field_error(error, line))?;
841        self.have_line2 = true;
842        Ok(())
843    }
844
845    /// `+` satellite-list line: `+   32   G01G02...` (3-char SV tokens from
846    /// column 9 in groups of 17). Continuation `+` lines append more tokens.
847    fn parse_plus_line(&mut self, line: &str, line_number: usize) -> Result<()> {
848        if line.starts_with("++") {
849            return self.parse_accuracy_line(line);
850        }
851        if self.satellite_header_lines == 1 {
852            self.declared_satellite_count = field(line, 3, 6).trim().parse::<usize>().ok();
853        }
854        // SV tokens start at column 9 (0-based), each 3 chars, up to 17 per line.
855        let mut col = 9;
856        while col + 3 <= line.len() {
857            let token = field(line, col, col + 3);
858            let trimmed = token.trim();
859            // Unused satellite slots are zero-filled, not a declaration. The
860            // SP3 zero-fill varies between producers (`  0`, ` 00`, `000`), so
861            // any all-zero (or blank) token is padding - never a satellite,
862            // whose token is a system letter + PRN (or, in SP3-a, a non-zero
863            // numeric PRN). Misreading ` 00` as an unrepresentable satellite
864            // inflates `skipped_records` and breaks the parse/write/parse
865            // round trip (the writer re-emits the canonical `  0`).
866            if trimmed.is_empty() || trimmed.bytes().all(|b| b == b'0') {
867                col += 3;
868                continue;
869            }
870            // This is a real declaration slot; the `++` accuracy codes are aligned
871            // to this axis, so track its index whether or not the token resolves.
872            let slot_index = self.declared_sat_slots;
873            self.declared_sat_slots += 1;
874            self.declared_satellite_tokens.push(trimmed.to_owned());
875            if let Some(id) = parse_sv_token(token, self.version) {
876                if !self.sat_list.contains(&id) {
877                    self.sat_list.push(id);
878                }
879            } else {
880                // A declared satellite whose token is not representable (e.g. an
881                // extended GLONASS slot R28 beyond the engine's PRN cap) is
882                // dropped from the satellite list, but counted rather than dropped
883                // silently - consistent with the position/velocity record paths
884                // (see `Sp3::skipped_records`). Record the slot so its accuracy
885                // column is skipped, keeping the surviving codes aligned.
886                self.push_unrepresentable_satellite_skip(line_number, token);
887                self.dropped_sat_slots.push(slot_index);
888            }
889            col += 3;
890        }
891        Ok(())
892    }
893
894    /// `++` per-satellite accuracy-code line: 3-char integer fields from column
895    /// 9, aligned with the `+` declaration slots.
896    ///
897    /// The columns track the `+` declaration order, so a column whose declaration
898    /// slot was dropped (an unrepresentable satellite) is read and discarded, not
899    /// pushed - otherwise the surviving satellites would inherit a neighbour's
900    /// accuracy code. With no dropped slots this is exactly the 1:1 push as before.
901    fn parse_accuracy_line(&mut self, line: &str) -> Result<()> {
902        let mut col = 9;
903        while col + 3 <= line.len() && self.accuracy_slot_cursor < self.declared_sat_slots {
904            let token = field(line, col, col + 3);
905            let trimmed = token.trim();
906            let code = if trimmed.is_empty() {
907                0
908            } else {
909                validate::strict_int::<u16>(trimmed, "satellite_accuracy_code")
910                    .map_err(|error| map_field_error(error, line))?
911            };
912            if !self.dropped_sat_slots.contains(&self.accuracy_slot_cursor) {
913                self.sat_accuracy_codes.push(code);
914            }
915            self.accuracy_slot_cursor += 1;
916            col += 3;
917        }
918        Ok(())
919    }
920
921    /// `%c` descriptor: the first one (chars 9-12) carries the time system.
922    fn parse_pc_line(&mut self, line: &str) -> Result<()> {
923        if self.pc_count == 0 {
924            // SP3-a is implicitly GPST regardless of descriptor content.
925            if matches!(self.version, Some(Sp3Version::A)) {
926                self.time_system = Some(Sp3TimeSystem::Gps);
927            } else if line.len() >= 12 {
928                let label = field(line, 9, 12);
929                let trimmed = label.trim();
930                // STRICT: a blank time-system field on the first %c is not GPST,
931                // it is malformed. Reject rather than silently defaulting so a
932                // precise pipeline never mis-attributes an epoch's scale.
933                if trimmed.is_empty() {
934                    return Err(Error::Parse(format!(
935                        "SP3 %c time system is blank in {line:?}"
936                    )));
937                }
938                self.time_system = Some(time_system_from_label(label)?);
939            } else {
940                // STRICT: a short %c line for SP3-b/c/d carries no time system
941                // we can trust. Reject rather than defaulting to GPST.
942                return Err(Error::Parse(format!(
943                    "SP3 %c descriptor too short to carry a time system: {line:?}"
944                )));
945            }
946        }
947        self.pc_count += 1;
948        Ok(())
949    }
950
951    /// Epoch line: `*  2020  6 24  0  0  0.00000000`.
952    fn parse_epoch_line(&mut self, line: &str) -> Result<()> {
953        // STRICT: by the time we reach data, the time system must be known -
954        // implicitly GPST for SP3-a (set at line 1), or from a valid first %c
955        // line for SP3-b/c/d. A missing/blank/short %c is an error, never GPST.
956        let time_system = self.time_system.ok_or_else(|| {
957            Error::Parse("SP3 epoch encountered with no time system (missing %c descriptor)".into())
958        })?;
959        let scale = time_system.time_scale();
960        // Fields after the leading `*  ` (3 chars), then space-delimited.
961        let body = &line[1..];
962        let mut it = body.split_whitespace();
963        let year: i64 = next_field(&mut it, "epoch year")?;
964        let month: i64 = next_field(&mut it, "epoch month")?;
965        let day: i64 = next_field(&mut it, "epoch day")?;
966        let hour: i64 = next_field(&mut it, "epoch hour")?;
967        let minute: i64 = next_field(&mut it, "epoch minute")?;
968        let seconds: f64 = next_field(&mut it, "epoch seconds")?;
969
970        let civil = validate::civil_datetime_with_second_policy(
971            year,
972            month,
973            day,
974            hour,
975            minute,
976            seconds,
977            time_system.civil_second_policy(),
978        )
979        .map_err(|error| map_field_error(error, line))?;
980        let split = civil_to_julian_split(civil)?;
981        let epoch_j2000_s = j2000_seconds(
982            civil.year as i32,
983            civil.month as i32,
984            civil.day as i32,
985            civil.hour as i32,
986            civil.minute as i32,
987            civil.second,
988        );
989        let epoch = Instant {
990            scale,
991            repr: InstantRepr::JulianDate(split),
992        };
993        self.epochs.push(epoch);
994        self.epoch_j2000_s.push(epoch_j2000_s);
995        self.states.push(BTreeMap::new());
996        self.interp_raw.push(BTreeMap::new());
997        self.epoch_position_tokens.push(Vec::new());
998        self.epoch_velocity_tokens.push(Vec::new());
999        self.epoch_state_record_sequence.push(Vec::new());
1000        self.current_epoch = Some(epoch);
1001        Ok(())
1002    }
1003
1004    /// Position+clock record: `PG01  x  y  z  clk  ...flags`.
1005    fn parse_position_line(&mut self, line: &str, line_number: usize) -> Result<()> {
1006        if self.current_epoch.is_none() {
1007            return Err(Error::Parse(
1008                "SP3 position record before any epoch line".into(),
1009            ));
1010        }
1011        if line.len() < 46 {
1012            return Err(Error::Parse(format!(
1013                "SP3 position record truncated before vector fields in {line:?}"
1014            )));
1015        }
1016        let token = field(line, 1, 4);
1017        self.epoch_position_tokens
1018            .last_mut()
1019            .expect("current epoch has a raw position-token list")
1020            .push(token.trim().to_owned());
1021        self.epoch_state_record_sequence
1022            .last_mut()
1023            .expect("current epoch has a raw state-record sequence")
1024            .push(('P', token.trim().to_owned()));
1025        let Some(sat) = parse_sv_token(token, self.version) else {
1026            // A token that does not parse to a representable `GnssSatelliteId`
1027            // (e.g. an extended GLONASS slot like R28 beyond the engine's PRN
1028            // cap) is an independent, unsupported record. One such record must
1029            // not reject the whole file - skip and count it, mirroring nav
1030            // `parse_glonass` and `parse_tle_file`.
1031            self.push_unrepresentable_satellite_skip(line_number, token);
1032            return Ok(());
1033        };
1034
1035        // The header `+` lines are the authoritative satellite declaration; a
1036        // position record for an undeclared satellite is malformed. Accepting it
1037        // would store a state the writer (which emits only declared satellites)
1038        // cannot reproduce, breaking parse/encode/parse round-tripping.
1039        if !self.sat_list.contains(&sat) {
1040            return Err(Error::Parse(format!(
1041                "SP3 position record for satellite {token:?} not in the header satellite list"
1042            )));
1043        }
1044
1045        let x_km = parse_coord(line, 4, 18)?;
1046        let y_km = parse_coord(line, 18, 32)?;
1047        let z_km = parse_coord(line, 32, 46)?;
1048
1049        // All-zero position is the missing-orbit sentinel: skip the record.
1050        if x_km == MISSING_POSITION_KM && y_km == MISSING_POSITION_KM && z_km == MISSING_POSITION_KM
1051        {
1052            return Ok(());
1053        }
1054
1055        let clock_us = parse_clock_us(line)?;
1056        let clock_s = clock_us.map(|us| us * US_TO_S);
1057
1058        let flags = parse_flags(line);
1059
1060        let position = ItrfPositionM::new(x_km * KM_TO_M, y_km * KM_TO_M, z_km * KM_TO_M)
1061            .map_err(|e| Error::Parse(format!("SP3 invalid position record: {e}")))?;
1062        let state = Sp3State {
1063            position,
1064            clock_s,
1065            velocity: None,
1066            clock_rate_s_s: None,
1067            flags,
1068        };
1069        let idx = self.states.len() - 1;
1070        self.states[idx].insert(sat, state);
1071        // Keep the native-unit node for the interpolation path (see RawNode):
1072        // the spline must fit the file's own km/us, not the km->m->km round trip.
1073        self.interp_raw[idx].insert(
1074            sat,
1075            RawNode {
1076                km: [x_km, y_km, z_km],
1077                clock_us,
1078                clock_event: flags.clock_event,
1079            },
1080        );
1081        Ok(())
1082    }
1083
1084    /// Velocity record: `VG01  vx  vy  vz  clkrate ...`. Augments the matching
1085    /// position record at the current epoch (must follow it).
1086    fn parse_velocity_line(&mut self, line: &str, line_number: usize) -> Result<()> {
1087        if self.current_epoch.is_none() {
1088            return Err(Error::Parse(
1089                "SP3 velocity record before any epoch line".into(),
1090            ));
1091        }
1092        if line.len() < 46 {
1093            return Err(Error::Parse(format!(
1094                "SP3 velocity record truncated before vector fields in {line:?}"
1095            )));
1096        }
1097        let token = field(line, 1, 4);
1098        self.epoch_velocity_tokens
1099            .last_mut()
1100            .expect("current epoch has a raw velocity-token list")
1101            .push(token.trim().to_owned());
1102        self.epoch_state_record_sequence
1103            .last_mut()
1104            .expect("current epoch has a raw state-record sequence")
1105            .push(('V', token.trim().to_owned()));
1106        let Some(sat) = parse_sv_token(token, self.version) else {
1107            // Unparsable / out-of-range satellite token: skip and count, same
1108            // as the position-record path above.
1109            self.push_unrepresentable_satellite_skip(line_number, token);
1110            return Ok(());
1111        };
1112
1113        // SP3 velocity is in dm/s; read each axis independently (the refs/sp3
1114        // crate has a bug here that reuses Y for X - we do not).
1115        let vx_dm_s = parse_coord(line, 4, 18)?;
1116        let vy_dm_s = parse_coord(line, 18, 32)?;
1117        let vz_dm_s = parse_coord(line, 32, 46)?;
1118
1119        let missing_velocity = vx_dm_s == MISSING_VELOCITY_DM_S
1120            && vy_dm_s == MISSING_VELOCITY_DM_S
1121            && vz_dm_s == MISSING_VELOCITY_DM_S;
1122        let velocity = ItrfVelocityMS::new(
1123            vx_dm_s * DM_S_TO_M_S,
1124            vy_dm_s * DM_S_TO_M_S,
1125            vz_dm_s * DM_S_TO_M_S,
1126        )
1127        .map_err(|e| Error::Parse(format!("SP3 invalid velocity record: {e}")))?;
1128
1129        // Clock-rate field shares the clock column; bad-clock sentinel applies.
1130        let clock_rate_s_s = parse_clock_us(line)?.map(|rate| rate * CLOCK_RATE_TO_S_PER_S);
1131
1132        let idx = self.states.len() - 1;
1133        match self.states[idx].get_mut(&sat) {
1134            Some(state) if !missing_velocity => {
1135                state.velocity = Some(velocity);
1136                state.clock_rate_s_s = clock_rate_s_s;
1137            }
1138            Some(_) => {}
1139            None => {
1140                // A V-record always follows its P-record for the same satellite
1141                // at the same epoch (SP3 format invariant). With no preceding
1142                // P-record this satellite has NO valid position at this epoch;
1143                // synthesizing one (e.g. the geocenter (0,0,0)) would fabricate
1144                // an orbit that the all-zero missing-orbit guard exists to
1145                // reject, and would leak through the public state()/states_at().
1146                // Treat it as malformed and skip - consistent with the parser's
1147                // tolerant skipping of other malformed records. No state is
1148                // inserted, so the satellite stays UnknownSatellite at this
1149                // epoch and no (0,0,0) position is ever exposed.
1150            }
1151        }
1152        Ok(())
1153    }
1154
1155    fn push_unrepresentable_satellite_skip(&mut self, line_number: usize, token: &str) {
1156        self.diagnostics.push_skip(Skip {
1157            at: RecordRef::at_line(line_number).with_satellite(token),
1158            reason: SkipReason::UnrepresentableSatellite,
1159        });
1160    }
1161
1162    fn finish(self) -> Result<Sp3> {
1163        if !self.have_line1 {
1164            return Err(Error::Parse("SP3 missing header line 1".into()));
1165        }
1166        if !self.have_line2 {
1167            return Err(Error::Parse("SP3 missing header line 2".into()));
1168        }
1169        let version = self
1170            .version
1171            .ok_or_else(|| Error::Parse("SP3 version not determined".into()))?;
1172        let data_type = self
1173            .data_type
1174            .ok_or_else(|| Error::Parse("SP3 data type not determined".into()))?;
1175        // STRICT: SP3-a is implicitly GPST (set at line 1); SP3-b/c/d must have
1176        // proved their scale from a valid first %c line. Never default here.
1177        let time_system = self.time_system.ok_or_else(|| {
1178            Error::Parse(
1179                "SP3 time system not determined (missing/short/blank %c descriptor)".into(),
1180            )
1181        })?;
1182        let time_scale = time_system.time_scale();
1183
1184        let mut satellite_accuracy_codes = self.sat_accuracy_codes;
1185        satellite_accuracy_codes.truncate(self.sat_list.len());
1186        satellite_accuracy_codes.resize(self.sat_list.len(), 0);
1187        let skipped_records = self.diagnostics.skips.len();
1188
1189        let header = Sp3Header {
1190            version,
1191            data_type,
1192            num_epochs: self.epochs.len() as u64,
1193            coordinate_system: self.coordinate_system,
1194            orbit_type: self.orbit_type,
1195            agency: self.agency,
1196            gnss_week: self.gnss_week,
1197            seconds_of_week: self.seconds_of_week,
1198            epoch_interval_s: self.epoch_interval_s,
1199            mjd: self.mjd,
1200            mjd_fraction: self.mjd_fraction,
1201            time_system,
1202            time_scale,
1203            satellites: self.sat_list,
1204            satellite_accuracy_codes,
1205        };
1206
1207        Ok(Sp3 {
1208            header,
1209            epochs: self.epochs,
1210            declared_num_epochs: self.num_epochs,
1211            declared_start_j2000_s: self.declared_start_j2000_s,
1212            had_eof: self.had_eof,
1213            trailing_content_after_eof: self.trailing_content_after_eof,
1214            satellite_header_lines: self.satellite_header_lines,
1215            accuracy_header_lines: self.accuracy_header_lines,
1216            time_system_header_lines: self.pc_count as usize,
1217            float_header_lines: self.float_header_lines,
1218            integer_header_lines: self.integer_header_lines,
1219            header_comment_lines: self.header_comment_lines,
1220            declared_satellite_count: self.declared_satellite_count,
1221            declared_satellite_tokens: self.declared_satellite_tokens,
1222            epoch_position_tokens: self.epoch_position_tokens,
1223            epoch_velocity_tokens: self.epoch_velocity_tokens,
1224            epoch_state_record_sequence: self.epoch_state_record_sequence,
1225            epoch_j2000_s: self.epoch_j2000_s,
1226            states: self.states,
1227            interp_raw: self.interp_raw,
1228            comments: self.comments,
1229            skipped_records,
1230        })
1231    }
1232}
1233
1234/// Best-effort parse of the civil start epoch carried on SP3 header line 1.
1235///
1236/// Exact validation treats `None` as an integrity failure. Keeping this helper
1237/// non-fallible preserves the long-standing permissive behavior of `Sp3::parse`
1238/// for callers that only consume epoch records.
1239fn parse_declared_start_j2000_s(line: &str) -> Option<f64> {
1240    let year = field(line, 3, 7).trim().parse::<i64>().ok()?;
1241    let month = field(line, 8, 10).trim().parse::<i64>().ok()?;
1242    let day = field(line, 11, 13).trim().parse::<i64>().ok()?;
1243    let hour = field(line, 14, 16).trim().parse::<i64>().ok()?;
1244    let minute = field(line, 17, 19).trim().parse::<i64>().ok()?;
1245    let second = field(line, 20, 31).trim().parse::<f64>().ok()?;
1246    let civil = validate::civil_datetime_with_second_policy(
1247        year,
1248        month,
1249        day,
1250        hour,
1251        minute,
1252        second,
1253        validate::CivilSecondPolicy::UtcLike,
1254    )
1255    .ok()?;
1256    Some(j2000_seconds(
1257        civil.year as i32,
1258        civil.month as i32,
1259        civil.day as i32,
1260        civil.hour as i32,
1261        civil.minute as i32,
1262        civil.second,
1263    ))
1264}
1265
1266/// Parse a fixed-column float coordinate, mapping failures to a parse error
1267/// that names the offending text.
1268fn parse_coord(line: &str, start: usize, end: usize) -> Result<f64> {
1269    let raw = field(line, start, end).trim();
1270    strict_f64(raw, "coordinate").map_err(|error| map_field_error(error, line))
1271}
1272
1273/// Parse the clock column (chars 46..60). Returns `None` for the bad-clock
1274/// sentinel `999999.999999` or an absent/blank field; `Some(us)` otherwise.
1275fn parse_clock_us(line: &str) -> Result<Option<f64>> {
1276    if line.len() <= 46 {
1277        return Ok(None);
1278    }
1279    let raw = field(line, 46, 60).trim();
1280    if raw.is_empty() {
1281        return Ok(None);
1282    }
1283    let value = strict_f64(raw, "clock").map_err(|error| map_field_error(error, line))?;
1284    // Sentinel: any value at or beyond the bad-clock magnitude is "no estimate".
1285    if value.abs() >= BAD_CLOCK_US {
1286        return Ok(None);
1287    }
1288    Ok(Some(value))
1289}
1290
1291fn map_field_error(error: validate::FieldError, line: &str) -> Error {
1292    Error::Parse(format!("SP3 {error} in {line:?}"))
1293}
1294
1295/// Parse the four status flags from their fixed columns (SP3-c/-d shared
1296/// layout): clock-event col 74 = `E`, clock-prediction col 75 = `P`,
1297/// maneuver col 78 = `M`, orbit-prediction col 79 = `P`.
1298fn parse_flags(line: &str) -> Sp3Flags {
1299    let at = |col: usize, want: char| -> bool { char_at(line, col) == Some(want) };
1300    Sp3Flags {
1301        clock_event: at(74, 'E'),
1302        clock_predicted: at(75, 'P'),
1303        maneuver: at(78, 'M'),
1304        orbit_predicted: at(79, 'P'),
1305    }
1306}
1307
1308/// Parse a 3-char SV token (e.g. `G01`, `C30`, or a bare `  1` in SP3-a) into a
1309/// [`GnssSatelliteId`]. Returns `None` on an unrecognized token.
1310fn parse_sv_token(token: &str, version: Option<Sp3Version>) -> Option<GnssSatelliteId> {
1311    let token = token.trim();
1312    if token.is_empty() {
1313        return None;
1314    }
1315    let first = token.chars().next()?;
1316    if first.is_ascii_digit() {
1317        // SP3-a GPS-only: bare numeric PRN, optionally space-padded.
1318        if matches!(version, Some(Sp3Version::A)) || version.is_none() {
1319            let prn = token.parse::<u8>().ok()?;
1320            if !is_valid_prn(GnssSystem::Gps, prn) {
1321                return None;
1322            }
1323            return GnssSatelliteId::new(GnssSystem::Gps, prn).ok();
1324        }
1325        return None;
1326    }
1327    token.parse::<GnssSatelliteId>().ok()
1328}
1329
1330/// Pull and parse the next whitespace-delimited field from an iterator.
1331fn next_field<T: std::str::FromStr>(
1332    it: &mut std::str::SplitWhitespace<'_>,
1333    what: &str,
1334) -> Result<T> {
1335    let tok = it
1336        .next()
1337        .ok_or_else(|| Error::Parse(format!("SP3 missing {what}")))?;
1338    tok.parse::<T>()
1339        .map_err(|_| Error::Parse(format!("SP3 {what} {tok:?} unparsable")))
1340}
1341
1342mod combine;
1343mod exact;
1344mod interp;
1345mod interpolant;
1346mod interpolant_store;
1347mod provenance;
1348mod samples;
1349mod verify;
1350mod write;
1351
1352pub use combine::{
1353    align_clock_reference, clock_reference_offset, merge, AgreementMetric, ClockReferenceOffset,
1354    EpochAgreement, MergeCombine, MergeFlag, MergeOptions, MergePrecedenceScope, MergeReport,
1355    OutlierRejectOptions, Sp3FrameLabelSet, Sp3FrameReconciliation, Sp3FrameReconciliationMethod,
1356    Sp3FrameReconciliationOptions,
1357};
1358pub use exact::{
1359    parse_exact_sp3, validate_exact_sp3, ExactSp3Coverage, ExactSp3Request, ExactSp3ValidationError,
1360};
1361pub use interpolant::{PreciseEphemerisInterpolant, PreciseInterpolantError};
1362pub use interpolant_store::{
1363    precise_interpolant_store_checksum64, MmapPreciseEphemerisInterpolant,
1364    PreciseInterpolantStoreError,
1365};
1366pub use provenance::{
1367    Sp3ArtifactIdentity, Sp3MergeInputIdentity, Sp3MergeInputIdentityError,
1368    SP3_MERGE_INPUT_ID_PREFIX, SP3_MERGE_INPUT_SCHEMA_VERSION,
1369};
1370pub use samples::{
1371    sp3_ecef_state_to_eci, PreciseEphemerisSample, PreciseEphemerisSamples,
1372    PreciseEphemerisStateSample, PreciseSamplesError,
1373};
1374pub use verify::{
1375    compare_position_series, InterpolationComparison, InterpolationDivergence, ReferenceState,
1376};
1377
1378#[cfg(all(test, sidereon_repo_tests))]
1379mod tests;