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