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;