Skip to main content

sidereon_core/astro/sgp4/
mod.rs

1//! SGP4 satellite propagation with 0 ULP parity to the Vallado C++ reference
2//! implementation (v2020-07-13).
3//!
4//! This module is a faithful pure-Rust port of Vallado's SGP4 verified bit-for-bit
5//! against the canonical 33-satellite / 198-propagation-point Vallado verification
6//! suite. TLEs are curve-fit parameters generated by Vallado's SGP4; using a
7//! different propagator introduces errors. This module preserves the exact
8//! floating-point computation order of the C++ reference so output matches
9//! Python's `sgp4` C extension (which compiles the same source) bit-for-bit.
10//!
11//! ## Quick start
12//!
13//! ```
14//! use sidereon_core::astro::sgp4::{Satellite, MinutesSinceEpoch};
15//!
16//! let line1 = "1 25544U 98067A   18184.80969102  .00001614  00000-0  31745-4 0  9993";
17//! let line2 = "2 25544  51.6414 295.8524 0003435 262.6267 204.2868 15.54005638121106";
18//!
19//! let sat = Satellite::from_tle(line1, line2).unwrap();
20//! let pred = sat.propagate(MinutesSinceEpoch(0.0)).unwrap();
21//!
22//! // position in km, velocity in km/s, TEME frame
23//! let _ = pred.position;
24//! let _ = pred.velocity;
25//! ```
26
27#[allow(
28    dead_code,
29    unused_variables,
30    unused_assignments,
31    unused_mut,
32    non_snake_case,
33    non_camel_case_types,
34    clippy::approx_constant,
35    clippy::excessive_precision,
36    clippy::too_many_arguments,
37    clippy::needless_return,
38    clippy::assign_op_pattern,
39    clippy::manual_range_contains,
40    clippy::collapsible_if,
41    clippy::collapsible_else_if,
42    clippy::float_cmp,
43    clippy::needless_late_init,
44    clippy::field_reassign_with_default
45)]
46mod vallado;
47
48use crate::astro::tle;
49use crate::validate::{self, FieldError};
50use thiserror::Error;
51
52const MAX_VALLADO_SATNUM: u32 = 99_999;
53
54// ── Error ────────────────────────────────────────────────────────────
55
56/// Validation failure category for public SGP4 inputs.
57#[derive(Debug, Clone, Copy, PartialEq, Eq)]
58pub enum Sgp4InputErrorKind {
59    /// A numeric input was NaN or infinite.
60    NonFinite,
61    /// A positive physical input was zero or negative.
62    NotPositive,
63    /// A non-negative physical input was negative.
64    Negative,
65    /// A numeric input was finite but outside the SGP4 domain.
66    OutOfRange,
67    /// A required input field was absent.
68    Missing,
69    /// A text field could not be parsed as a float.
70    FloatParse,
71    /// A text field could not be parsed as an integer.
72    IntParse,
73    /// A civil date field was out of range.
74    InvalidCivilDate,
75    /// A civil time field was out of range.
76    InvalidCivilTime,
77}
78
79impl core::fmt::Display for Sgp4InputErrorKind {
80    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
81        let label = match self {
82            Self::NonFinite => "not finite",
83            Self::NotPositive => "not positive",
84            Self::Negative => "negative",
85            Self::OutOfRange => "out of range",
86            Self::Missing => "missing",
87            Self::FloatParse => "invalid float",
88            Self::IntParse => "invalid integer",
89            Self::InvalidCivilDate => "invalid civil date",
90            Self::InvalidCivilTime => "invalid civil time",
91        };
92        f.write_str(label)
93    }
94}
95
96impl From<&FieldError> for Sgp4InputErrorKind {
97    fn from(error: &FieldError) -> Self {
98        match error {
99            FieldError::Missing { .. } => Self::Missing,
100            FieldError::NonFinite { .. } => Self::NonFinite,
101            FieldError::NotPositive { .. } => Self::NotPositive,
102            FieldError::Negative { .. } => Self::Negative,
103            FieldError::OutOfRange { .. } => Self::OutOfRange,
104            FieldError::FloatParse { .. } => Self::FloatParse,
105            FieldError::IntParse { .. } => Self::IntParse,
106            FieldError::InvalidCivilDate { .. } => Self::InvalidCivilDate,
107            FieldError::InvalidCivilTime { .. } => Self::InvalidCivilTime,
108        }
109    }
110}
111
112/// Error from SGP4 initialization or propagation.
113#[derive(Error, Debug, Clone, PartialEq)]
114pub enum Error {
115    /// A public SGP4 input was malformed, non-finite, or outside the model
116    /// domain. Boundary validation rejects this before the Vallado kernel runs.
117    #[error("invalid SGP4 input {field}: {kind}")]
118    InvalidInput {
119        /// The invalid input field.
120        field: &'static str,
121        /// The validation failure category.
122        kind: Sgp4InputErrorKind,
123    },
124    /// The Vallado step kernel returned a non-finite state vector.
125    #[error("SGP4 returned non-finite {field}")]
126    NonFiniteOutput {
127        /// The output vector that was non-finite.
128        field: &'static str,
129    },
130    /// TLE line has invalid format.
131    #[error("invalid TLE: {0}")]
132    InvalidTle(String),
133    /// SGP4 returned a non-zero error code.
134    ///
135    /// Codes: 1 = mean elements, 2 = mean motion, 3 = perturbed elements,
136    /// 4 = semi-latus rectum, 5 = epoch elements sub-orbital,
137    /// 6 = satellite decayed.
138    #[error("SGP4 error code {code}")]
139    Sgp4 { code: i32 },
140}
141
142const MAX_MINUTES_SINCE_EPOCH: f64 = 10_000_000.0;
143
144// ── Types ────────────────────────────────────────────────────────────
145
146/// Minutes since the TLE epoch. Newtype to prevent mixing with raw `f64`.
147#[derive(Debug, Clone, Copy, PartialEq)]
148#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
149pub struct MinutesSinceEpoch(pub f64);
150
151/// Propagation result in the TEME (True Equator, Mean Equinox) frame.
152#[derive(Debug, Clone, PartialEq)]
153#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
154pub struct Prediction {
155    /// Position in km, TEME frame.
156    pub position: [f64; 3],
157    /// Velocity in km/s, TEME frame.
158    pub velocity: [f64; 3],
159}
160
161/// Julian date split as `(whole, fraction)` for high-precision time input.
162///
163/// Skyfield convention: `whole = floor(JD)`, `fraction = remainder`.
164/// For example, 2018-07-04 00:00:00 UTC = `JulianDate(2458303.0, 0.5)`.
165#[derive(Debug, Clone, Copy, PartialEq)]
166#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
167pub struct JulianDate(pub f64, pub f64);
168
169pub(crate) fn sgp4_julian_date_from_calendar(
170    year: i32,
171    mon: i32,
172    day: i32,
173    hr: i32,
174    minute: i32,
175    sec: f64,
176) -> JulianDate {
177    let (jd, jdfrac) = vallado::jday_SGP4(year, mon, day, hr, minute, sec);
178    JulianDate(jd, jdfrac)
179}
180
181pub(crate) fn sgp4_julian_date_from_day_of_year(year: i32, days: f64) -> JulianDate {
182    let (mon, day, hr, minute, sec) = vallado::days2mdhms_SGP4(year, days);
183    let JulianDate(jd, jdfrac_raw) =
184        sgp4_julian_date_from_calendar(year, mon, day, hr, minute, sec);
185    let jdfrac = (jdfrac_raw * 100_000_000.0).round() / 100_000_000.0;
186    JulianDate(jd, jdfrac)
187}
188
189/// Vallado SGP4 operation mode. Controls which initialization branch
190/// `sgp4init` follows. The two modes produce subtly different results;
191/// the divergence grows with propagation time and can reach hundreds of
192/// millimeters over a few orbits at LEO.
193#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
194#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
195pub enum OpsMode {
196    /// Improved formulation (Vallado `'i'`). Default for new work and
197    /// matches the default of Python's `sgp4` package. Recommended unless
198    /// you specifically need AFSPC operational parity.
199    #[default]
200    Improved,
201    /// AFSPC operational compatibility mode (Vallado `'a'`). Use this
202    /// when reproducing outputs from operational AFSPC systems or matching
203    /// reference values from older crates / catalogs that ran in AFSPC mode.
204    Afspc,
205}
206
207impl OpsMode {
208    fn as_char(self) -> char {
209        match self {
210            OpsMode::Improved => 'i',
211            OpsMode::Afspc => 'a',
212        }
213    }
214}
215
216/// Pre-parsed Vallado SGP4 element set.
217///
218/// Use this when the TLE has already been parsed externally (e.g. from an
219/// OMM message, JSON catalog, or another system) and you want to feed the
220/// element values directly into the SGP4 initializer instead of going through
221/// the TLE string parser.
222///
223/// Field units match the Vallado SGP4 reference inputs:
224/// angles in **degrees**, mean motion in **revolutions per day**, drag term
225/// in the dimensionless TLE convention.
226#[derive(Debug, Clone, PartialEq)]
227#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
228pub struct ElementSet {
229    /// Epoch as a split Julian date `(whole, fraction)`.
230    ///
231    /// This is format-agnostic and loss-free for the SGP4 initializer: TLE
232    /// conversion stores the exact split JD produced by Vallado's
233    /// `days2mdhms`/`jday` path with the legacy 8-decimal fraction rounding,
234    /// while OMM conversion stores the split JD from its full calendar
235    /// timestamp directly.
236    pub epoch: JulianDate,
237    /// SGP4 drag term (Vallado B\*). Dimensionless TLE convention.
238    pub bstar: f64,
239    /// First derivative of mean motion in rev/day². TLE "ndot".
240    pub mean_motion_dot: f64,
241    /// Second derivative of mean motion in rev/day³. TLE "nddot".
242    pub mean_motion_double_dot: f64,
243    /// Eccentricity, dimensionless, in [0, 1).
244    pub eccentricity: f64,
245    /// Argument of perigee, degrees.
246    pub argument_of_perigee_deg: f64,
247    /// Inclination, degrees.
248    pub inclination_deg: f64,
249    /// Mean anomaly, degrees.
250    pub mean_anomaly_deg: f64,
251    /// Mean motion, revolutions per day.
252    pub mean_motion_rev_per_day: f64,
253    /// Right ascension of ascending node (RAAN), degrees.
254    pub right_ascension_deg: f64,
255    /// Catalog (NORAD) number for this object. Used only for diagnostic
256    /// reporting inside SGP4 - propagation results do not depend on it.
257    /// Pass `0` if unknown.
258    pub catalog_number: u32,
259}
260
261// ── Satellite ────────────────────────────────────────────────────────
262
263/// A parsed TLE ready for propagation.
264///
265/// Holds the raw TLE lines plus the initialized SGP4 satellite record. The
266/// TLE is parsed and `sgp4init` is run exactly once during `from_tle`, so
267/// subsequent `propagate` calls just invoke the propagation kernel directly:
268/// fast, and crucially, with no precision loss from JD round-tripping.
269#[derive(Clone)]
270pub struct Satellite {
271    line1: String,
272    line2: String,
273    /// Source elements used to initialize the cached SGP4 record. Kept so
274    /// element-built satellites can serialize without inventing TLE text.
275    elements: ElementSet,
276    opsmode: OpsMode,
277    /// Pre-initialized satellite record. Boxed to keep `Satellite` small on
278    /// the stack - `ElsetRec` has ~150 fields.
279    satrec: Box<vallado::ElsetRec>,
280}
281
282impl std::fmt::Debug for Satellite {
283    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
284        f.debug_struct("Satellite")
285            .field("line1", &self.line1)
286            .field("line2", &self.line2)
287            .field("elements", &self.elements)
288            .field("opsmode", &self.opsmode)
289            .finish_non_exhaustive()
290    }
291}
292
293impl Satellite {
294    /// Parse a two-line element set using the default `Improved` opsmode.
295    ///
296    /// Lines should be the standard 69-character TLE format. Leading and
297    /// trailing whitespace is trimmed before length validation. Runs the full
298    /// SGP4 initialization (`sgp4init`) once and caches the resulting state
299    /// so propagation calls are pure step kernels.
300    pub fn from_tle(line1: &str, line2: &str) -> Result<Self, Error> {
301        Self::from_tle_with_opsmode(line1, line2, OpsMode::Improved)
302    }
303
304    /// Parse a two-line element set with an explicit `OpsMode`.
305    ///
306    /// Use `OpsMode::Afspc` to reproduce results from operational AFSPC
307    /// systems or older crates that ran in AFSPC compatibility mode.
308    ///
309    /// The TLE flows through the canonical element-set IR: it is parsed by the
310    /// forgiving [`crate::astro::tle`] grammar into [`crate::astro::tle::TleElements`],
311    /// converted to an [`ElementSet`], and initialized through the single
312    /// [`Satellite::from_elements`] path - the same path OMM and any other input
313    /// format use. There is no separate TLE-direct initialization. The parse is
314    /// lenient on cosmetics (trailing content past column 69, leading-dot and
315    /// assumed-decimal field forms, an advisory checksum) but still rejects
316    /// genuinely corrupt input (non-ASCII, wrong line structure, mismatched
317    /// satellite numbers).
318    pub fn from_tle_with_opsmode(
319        line1: &str,
320        line2: &str,
321        opsmode: OpsMode,
322    ) -> Result<Self, Error> {
323        let l1 = line1.trim();
324        let l2 = line2.trim();
325
326        let parsed = tle::parse(l1, l2).map_err(|e| Error::InvalidTle(e.to_string()))?;
327        let elements = parsed
328            .elements
329            .to_element_set()
330            .map_err(map_tle_bridge_error)?;
331        let satrec = init_satrec_from_elements(&elements, opsmode)?;
332
333        Ok(Satellite {
334            line1: l1.to_string(),
335            line2: l2.to_string(),
336            elements,
337            opsmode,
338            satrec: Box::new(satrec),
339        })
340    }
341
342    /// Parse a satellite from a TLE block that may carry a leading name line.
343    ///
344    /// CelesTrak and Space-Track serve TLEs in the common "3-line" (TLE/3LE)
345    /// form: an object-name line followed by the two element lines. This accepts
346    /// that block, strips an optional leading name line, and parses the first
347    /// `1 `/`2 ` element-line pair. A plain 2-line block (no name line) also
348    /// parses. Uses the default `Improved` opsmode.
349    pub fn from_3line(block: &str) -> Result<Self, Error> {
350        Self::from_3line_with_opsmode(block, OpsMode::Improved)
351    }
352
353    /// Parse a name-line-prefixed TLE block with an explicit `OpsMode`. See
354    /// [`Satellite::from_3line`].
355    pub fn from_3line_with_opsmode(block: &str, opsmode: OpsMode) -> Result<Self, Error> {
356        let mut l1 = None;
357        let mut l2 = None;
358        for line in block.lines() {
359            let line = line.trim();
360            if l1.is_none() && line.starts_with("1 ") {
361                l1 = Some(line.to_string());
362            } else if l2.is_none() && line.starts_with("2 ") {
363                l2 = Some(line.to_string());
364            }
365        }
366        let l1 = l1.ok_or_else(|| Error::InvalidTle("no line 1 in TLE block".into()))?;
367        let l2 = l2.ok_or_else(|| Error::InvalidTle("no line 2 in TLE block".into()))?;
368        Self::from_tle_with_opsmode(&l1, &l2, opsmode)
369    }
370
371    /// Construct a `Satellite` from pre-parsed Vallado SGP4 elements using
372    /// the default `Improved` opsmode.
373    ///
374    /// Useful when TLE data has already been parsed externally (OMM, JSON
375    /// catalog, another system) and you only need the element values to flow
376    /// into SGP4 initialization. Equivalent to `from_tle` for propagation
377    /// behavior, but bypasses the TLE string parser.
378    ///
379    /// Note: a `Satellite` constructed this way has empty `line1()` and
380    /// `line2()` accessors since there is no source TLE to return.
381    pub fn from_elements(elements: &ElementSet) -> Result<Self, Error> {
382        Self::from_elements_with_opsmode(elements, OpsMode::Improved)
383    }
384
385    /// Construct a `Satellite` from pre-parsed elements with an explicit
386    /// `OpsMode`. See `from_tle_with_opsmode` for the rationale.
387    pub fn from_elements_with_opsmode(
388        elements: &ElementSet,
389        opsmode: OpsMode,
390    ) -> Result<Self, Error> {
391        let satrec = init_satrec_from_elements(elements, opsmode)?;
392        Ok(Satellite {
393            line1: String::new(),
394            line2: String::new(),
395            elements: elements.clone(),
396            opsmode,
397            satrec: Box::new(satrec),
398        })
399    }
400
401    /// Propagate to a time given as minutes since the TLE epoch.
402    ///
403    /// Calls the SGP4 step kernel directly with the supplied tsince - no JD
404    /// round-trip, no precision loss.
405    pub fn propagate(&self, t: MinutesSinceEpoch) -> Result<Prediction, Error> {
406        // Clone the satrec so propagation doesn't mutate the cached state
407        // (sgp4 writes back into the satrec - error code, atime, etc.).
408        propagate_satrec((*self.satrec).clone(), t)
409    }
410
411    /// Propagate to a Julian date, split as `(whole, fraction)`.
412    ///
413    /// Computes the tsince from the cached epoch via the same subtraction
414    /// the C++ wrapper uses:
415    ///
416    /// ```text
417    /// tsince = (jd - jdsatepoch) * 1440 + (fr - jdsatepochF) * 1440
418    /// ```
419    pub fn propagate_jd(&self, jd: JulianDate) -> Result<Prediction, Error> {
420        validate::finite(jd.0, "julian_date.whole").map_err(map_input_error)?;
421        validate::finite_in_range_exclusive_upper(jd.1, 0.0, 1.0, "julian_date.fraction")
422            .map_err(map_input_error)?;
423        let tsince =
424            (jd.0 - self.satrec.jdsatepoch) * 1440.0 + (jd.1 - self.satrec.jdsatepochF) * 1440.0;
425        validate::finite(tsince, "minutes_since_epoch").map_err(map_input_error)?;
426        self.propagate(MinutesSinceEpoch(tsince))
427    }
428
429    pub(crate) fn mean_motion_rad_per_min(&self) -> f64 {
430        self.satrec.no_kozai
431    }
432
433    pub(crate) fn eccentricity(&self) -> f64 {
434        self.satrec.ecco
435    }
436
437    /// Raw TLE line 1. Returns an empty string when this `Satellite` was
438    /// constructed via `from_elements` (no source TLE).
439    pub fn line1(&self) -> &str {
440        &self.line1
441    }
442
443    /// Raw TLE line 2. Returns an empty string when this `Satellite` was
444    /// constructed via `from_elements` (no source TLE).
445    pub fn line2(&self) -> &str {
446        &self.line2
447    }
448
449    /// Cached TLE epoch as a split Julian date `(jdsatepoch, jdsatepochF)`.
450    ///
451    /// Useful for computing time offsets between two TLEs without losing
452    /// precision through floating-point round-trips.
453    pub fn epoch_jd(&self) -> JulianDate {
454        JulianDate(self.satrec.jdsatepoch, self.satrec.jdsatepochF)
455    }
456
457    #[cfg(feature = "serde")]
458    fn has_source_tle(&self) -> bool {
459        !self.line1.is_empty() && !self.line2.is_empty()
460    }
461}
462
463/// One-shot SGP4 propagation from pre-parsed elements using the default
464/// `Improved` opsmode.
465///
466/// Equivalent to `Satellite::from_elements(&e)?.propagate(t)` but without
467/// allocating a cached `Satellite`. Suitable for one-call use cases where
468/// the satellite record is not reused (e.g. NIF entry points that get
469/// elements + a single time per call).
470pub fn propagate_elements(
471    elements: &ElementSet,
472    t: MinutesSinceEpoch,
473) -> Result<Prediction, Error> {
474    propagate_elements_with_opsmode(elements, t, OpsMode::Improved)
475}
476
477/// One-shot SGP4 propagation with an explicit `OpsMode`.
478pub fn propagate_elements_with_opsmode(
479    elements: &ElementSet,
480    t: MinutesSinceEpoch,
481    opsmode: OpsMode,
482) -> Result<Prediction, Error> {
483    let satrec = init_satrec_from_elements(elements, opsmode)?;
484    propagate_satrec(satrec, t)
485}
486
487#[cfg(feature = "serde")]
488impl serde::Serialize for Satellite {
489    fn serialize<S: serde::Serializer>(&self, s: S) -> Result<S::Ok, S::Error> {
490        use serde::ser::SerializeStruct;
491        let mut st = s.serialize_struct("Satellite", 2)?;
492        if self.has_source_tle() {
493            st.serialize_field("line1", &self.line1)?;
494            st.serialize_field("line2", &self.line2)?;
495        } else {
496            st.serialize_field("elements", &self.elements)?;
497            st.serialize_field("opsmode", &self.opsmode)?;
498        }
499        st.end()
500    }
501}
502
503#[cfg(feature = "serde")]
504impl<'de> serde::Deserialize<'de> for Satellite {
505    fn deserialize<D: serde::Deserializer<'de>>(d: D) -> Result<Self, D::Error> {
506        #[derive(serde::Deserialize)]
507        struct Wire {
508            line1: Option<String>,
509            line2: Option<String>,
510            elements: Option<ElementSet>,
511            opsmode: Option<OpsMode>,
512        }
513        let w = Wire::deserialize(d)?;
514        let opsmode = w.opsmode.unwrap_or_default();
515        let has_tle_line = w
516            .line1
517            .as_deref()
518            .is_some_and(|line| !line.trim().is_empty())
519            || w.line2
520                .as_deref()
521                .is_some_and(|line| !line.trim().is_empty());
522        if let Some(elements) = w.elements {
523            if has_tle_line {
524                Err(serde::de::Error::custom(
525                    "ambiguous Satellite wire format: use either TLE lines or elements",
526                ))
527            } else {
528                Satellite::from_elements_with_opsmode(&elements, opsmode)
529                    .map_err(serde::de::Error::custom)
530            }
531        } else if let (Some(line1), Some(line2)) = (w.line1, w.line2) {
532            if line1.trim().is_empty() || line2.trim().is_empty() {
533                Err(serde::de::Error::custom(
534                    "Satellite wire format requires non-empty line1/line2 or elements",
535                ))
536            } else {
537                Satellite::from_tle_with_opsmode(&line1, &line2, opsmode)
538                    .map_err(serde::de::Error::custom)
539            }
540        } else {
541            Err(serde::de::Error::custom(
542                "Satellite wire format requires non-empty line1/line2 or elements",
543            ))
544        }
545    }
546}
547
548// ── Internal helpers ─────────────────────────────────────────────────
549
550fn propagate_satrec(
551    mut satrec: vallado::ElsetRec,
552    t: MinutesSinceEpoch,
553) -> Result<Prediction, Error> {
554    validate::finite(t.0, "minutes_since_epoch").map_err(map_input_error)?;
555    if t.0.abs() > MAX_MINUTES_SINCE_EPOCH {
556        return Err(invalid_domain("minutes_since_epoch"));
557    }
558
559    let mut r = [0.0_f64; 3];
560    let mut v = [0.0_f64; 3];
561    let ok = vallado::sgp4(&mut satrec, t.0, &mut r, &mut v);
562    if !ok || satrec.error != 0 {
563        return Err(Error::Sgp4 { code: satrec.error });
564    }
565    validate_prediction(r, v)?;
566    Ok(Prediction {
567        position: r,
568        velocity: v,
569    })
570}
571
572fn validate_prediction(position: [f64; 3], velocity: [f64; 3]) -> Result<(), Error> {
573    validate::finite_vec3(position, "position_km").map_err(map_output_error)?;
574    validate::finite_vec3(velocity, "velocity_km_s").map_err(map_output_error)?;
575    Ok(())
576}
577
578/// Run `sgp4init` from a pre-parsed element set, returning the initialized
579/// satellite record. Performs the same angle/units conversion as
580/// `vallado::twoline2rv_propagate` so a `Satellite` constructed from
581/// elements is equivalent to one constructed from the matching TLE.
582fn init_satrec_from_elements(
583    elements: &ElementSet,
584    opsmode: OpsMode,
585) -> Result<vallado::ElsetRec, Error> {
586    validate_elements(elements)?;
587
588    let deg2rad = std::f64::consts::PI / 180.0;
589    let xpdotp = 1440.0 / (2.0 * std::f64::consts::PI);
590
591    let inclo = elements.inclination_deg * deg2rad;
592    let nodeo = elements.right_ascension_deg * deg2rad;
593    let argpo = elements.argument_of_perigee_deg * deg2rad;
594    let mo = elements.mean_anomaly_deg * deg2rad;
595    let no_kozai = elements.mean_motion_rev_per_day / xpdotp;
596    // ndot rev/day² → rad/min², nddot rev/day³ → rad/min³.
597    // Matches the conversion in `vallado::twoline2rv_propagate`.
598    let ndot = elements.mean_motion_dot / (xpdotp * 1440.0);
599    let nddot = elements.mean_motion_double_dot / (xpdotp * 1440.0 * 1440.0);
600
601    let JulianDate(jd, jdfrac) = elements.epoch;
602    let epoch_sgp4 = jd + jdfrac - 2433281.5;
603
604    let satnum_str = format!("{:>5}", elements.catalog_number);
605
606    let mut satrec = vallado::ElsetRec {
607        jdsatepoch: jd,
608        jdsatepochF: jdfrac,
609        ..vallado::ElsetRec::default()
610    };
611
612    vallado::sgp4init(
613        vallado::GravConstType::Wgs72,
614        opsmode.as_char(),
615        &satnum_str,
616        epoch_sgp4,
617        elements.bstar,
618        ndot,
619        nddot,
620        elements.eccentricity,
621        argpo,
622        inclo,
623        mo,
624        no_kozai,
625        nodeo,
626        &mut satrec,
627    );
628
629    // sgp4init may have rewritten jdsatepoch via initl - restore the split.
630    satrec.jdsatepoch = jd;
631    satrec.jdsatepochF = jdfrac;
632
633    Ok(satrec)
634}
635
636fn validate_elements(elements: &ElementSet) -> Result<(), Error> {
637    if elements.catalog_number > MAX_VALLADO_SATNUM {
638        return Err(invalid_domain("element.catalog_number"));
639    }
640    validate_epoch(elements.epoch)?;
641    validate::finite(elements.bstar, "element.bstar").map_err(map_input_error)?;
642    validate::finite(elements.mean_motion_dot, "element.mean_motion_dot")
643        .map_err(map_input_error)?;
644    validate::finite(
645        elements.mean_motion_double_dot,
646        "element.mean_motion_double_dot",
647    )
648    .map_err(map_input_error)?;
649    validate::finite_in_range_exclusive_upper(
650        elements.eccentricity,
651        0.0,
652        1.0,
653        "element.eccentricity",
654    )
655    .map_err(map_input_error)?;
656    validate::finite(
657        elements.argument_of_perigee_deg,
658        "element.argument_of_perigee_deg",
659    )
660    .map_err(map_input_error)?;
661    validate::finite(elements.inclination_deg, "element.inclination_deg")
662        .map_err(map_input_error)?;
663    validate::finite(elements.mean_anomaly_deg, "element.mean_anomaly_deg")
664        .map_err(map_input_error)?;
665    validate::finite_positive(
666        elements.mean_motion_rev_per_day,
667        "element.mean_motion_rev_per_day",
668    )
669    .map_err(map_input_error)?;
670    validate::finite(elements.right_ascension_deg, "element.right_ascension_deg")
671        .map_err(map_input_error)?;
672
673    Ok(())
674}
675
676fn validate_epoch(epoch: JulianDate) -> Result<(), Error> {
677    validate::finite(epoch.0, "element.epoch.whole").map_err(map_input_error)?;
678    validate::finite(epoch.1, "element.epoch.fraction").map_err(map_input_error)?;
679
680    let total = epoch.0 + epoch.1;
681    validate::finite(total, "element.epoch").map_err(map_input_error)?;
682    if !(0.0..=5_000_000.0).contains(&total) {
683        return Err(invalid_domain("element.epoch"));
684    }
685    Ok(())
686}
687
688fn map_input_error(error: FieldError) -> Error {
689    Error::InvalidInput {
690        field: error.field(),
691        kind: Sgp4InputErrorKind::from(&error),
692    }
693}
694
695fn invalid_domain(field: &'static str) -> Error {
696    Error::InvalidInput {
697        field,
698        kind: Sgp4InputErrorKind::OutOfRange,
699    }
700}
701
702fn map_tle_bridge_error(error: tle::TleError) -> Error {
703    match error {
704        tle::TleError::InvalidField { field, reason } => Error::InvalidInput {
705            field,
706            kind: match reason {
707                "not finite" => Sgp4InputErrorKind::NonFinite,
708                "not positive" => Sgp4InputErrorKind::NotPositive,
709                "negative" => Sgp4InputErrorKind::Negative,
710                "out of range" => Sgp4InputErrorKind::OutOfRange,
711                _ => Sgp4InputErrorKind::OutOfRange,
712            },
713        },
714        other => Error::InvalidTle(other.to_string()),
715    }
716}
717
718fn map_output_error(error: FieldError) -> Error {
719    Error::NonFiniteOutput {
720        field: error.field(),
721    }
722}
723
724#[cfg(test)]
725mod tests {
726    use super::{
727        propagate_elements, ElementSet, Error, JulianDate, MinutesSinceEpoch, Satellite,
728        Sgp4InputErrorKind, MAX_MINUTES_SINCE_EPOCH,
729    };
730
731    /// A TLE carrying a multibyte character inside a fixed-width column must
732    /// return a typed [`Error::InvalidTle`] rather than panicking. The field
733    /// extractor slices by byte column (`l1[18..20]`, ...); a non-ASCII byte
734    /// inside such a window used to panic on a non-char-boundary slice. The
735    /// ASCII guard now rejects the line with a typed error; valid input parses.
736    #[test]
737    fn non_ascii_tle_returns_invalid_tle_not_panic() {
738        let line1 = "1 25544U 98067A   18184.80969102  .00001614  00000-0  31745-4 0  9993";
739        let line2 = "2 25544  51.6414 295.8524 0003435 262.6267 204.2868 15.54005638121106";
740        assert!(
741            Satellite::from_tle(line1, line2).is_ok(),
742            "clean ASCII TLE must still parse"
743        );
744
745        // Drop a 3-byte character into the epoch-year column (bytes 18..20),
746        // straddling byte 20 so the pre-guard `l1[18..20]` slice would panic.
747        let mut bad1 = String::from(&line1[..18]);
748        bad1.push('\u{20ac}');
749        bad1.push_str(&line1[19..]);
750        assert!(
751            !bad1.is_char_boundary(20),
752            "corruption must straddle byte 20"
753        );
754
755        let err = Satellite::from_tle(&bad1, line2).expect_err("non-ASCII TLE must not parse");
756        assert!(
757            matches!(err, Error::InvalidTle(_)),
758            "expected a typed InvalidTle error, got: {err:?}"
759        );
760    }
761
762    // ── Forgiving-inbound leniency (now that `from_tle` flows through the
763    // canonical IR via the lenient `tle` parser) ───────────────────────────
764
765    const ISS_L1: &str = "1 25544U 98067A   18184.80969102  .00001614  00000-0  31745-4 0  9993";
766    const ISS_L2: &str = "2 25544  51.6414 295.8524 0003435 262.6267 204.2868 15.54005638121106";
767
768    fn iss_elements() -> ElementSet {
769        crate::astro::tle::parse(ISS_L1, ISS_L2)
770            .unwrap()
771            .elements
772            .to_element_set()
773            .expect("valid TLE bridge")
774    }
775
776    fn assert_invalid_input<T>(
777        result: Result<T, Error>,
778        field: &'static str,
779        kind: Sgp4InputErrorKind,
780    ) {
781        match result {
782            Err(Error::InvalidInput {
783                field: actual_field,
784                kind: actual_kind,
785            }) => {
786                assert_eq!(actual_field, field);
787                assert_eq!(actual_kind, kind);
788            }
789            Err(err) => panic!("expected InvalidInput({field}, {kind}), got {err:?}"),
790            Ok(_) => panic!("expected InvalidInput({field}, {kind}), got Ok"),
791        }
792    }
793
794    /// Assert two satellites are bit-identical (same cached epoch and same
795    /// propagated state at several offsets).
796    fn assert_same(a: &Satellite, b: &Satellite) {
797        let (ea, eb) = (a.epoch_jd(), b.epoch_jd());
798        assert_eq!(
799            (ea.0.to_bits(), ea.1.to_bits()),
800            (eb.0.to_bits(), eb.1.to_bits()),
801            "epoch JD differs"
802        );
803        for &t in &[0.0, 100.0, 1440.0] {
804            let pa = a.propagate(MinutesSinceEpoch(t)).unwrap();
805            let pb = b.propagate(MinutesSinceEpoch(t)).unwrap();
806            for axis in 0..3 {
807                assert_eq!(
808                    pa.position[axis].to_bits(),
809                    pb.position[axis].to_bits(),
810                    "position[{axis}] differs at t={t}"
811                );
812                assert_eq!(
813                    pa.velocity[axis].to_bits(),
814                    pb.velocity[axis].to_bits(),
815                    "velocity[{axis}] differs at t={t}"
816                );
817            }
818        }
819    }
820
821    #[cfg(feature = "serde")]
822    #[test]
823    fn serde_round_trips_tle_satellites() {
824        let sat = Satellite::from_tle(ISS_L1, ISS_L2).unwrap();
825        let encoded = serde_json::to_string(&sat).unwrap();
826        assert!(encoded.contains("\"line1\""));
827        assert!(encoded.contains("\"line2\""));
828        assert!(!encoded.contains("\"elements\""));
829
830        let decoded: Satellite = serde_json::from_str(&encoded).unwrap();
831        assert_eq!(decoded.line1(), ISS_L1);
832        assert_eq!(decoded.line2(), ISS_L2);
833        assert_same(&sat, &decoded);
834    }
835
836    #[cfg(feature = "serde")]
837    #[test]
838    fn serde_round_trips_element_built_satellites() {
839        let elements = iss_elements();
840        let sat = Satellite::from_elements(&elements).unwrap();
841        let encoded = serde_json::to_string(&sat).unwrap();
842        assert!(encoded.contains("\"elements\""));
843        assert!(encoded.contains("\"opsmode\""));
844        assert!(!encoded.contains("\"line1\""));
845        assert!(!encoded.contains("\"line2\""));
846
847        let decoded: Satellite = serde_json::from_str(&encoded).unwrap();
848        assert!(decoded.line1().is_empty());
849        assert!(decoded.line2().is_empty());
850        assert_same(&sat, &decoded);
851    }
852
853    #[test]
854    fn from_elements_rejects_non_finite_fields_before_sgp4init() {
855        let mut elements = iss_elements();
856        elements.bstar = f64::NAN;
857
858        assert_invalid_input(
859            Satellite::from_elements(&elements),
860            "element.bstar",
861            Sgp4InputErrorKind::NonFinite,
862        );
863    }
864
865    #[test]
866    fn from_elements_rejects_sgp4_domain_before_sgp4init() {
867        let mut elements = iss_elements();
868        elements.mean_motion_rev_per_day = 0.0;
869        assert_invalid_input(
870            Satellite::from_elements(&elements),
871            "element.mean_motion_rev_per_day",
872            Sgp4InputErrorKind::NotPositive,
873        );
874
875        let mut elements = iss_elements();
876        elements.eccentricity = -0.1;
877        assert_invalid_input(
878            Satellite::from_elements(&elements),
879            "element.eccentricity",
880            Sgp4InputErrorKind::OutOfRange,
881        );
882
883        let mut elements = iss_elements();
884        elements.eccentricity = 1.0;
885        assert_invalid_input(
886            Satellite::from_elements(&elements),
887            "element.eccentricity",
888            Sgp4InputErrorKind::OutOfRange,
889        );
890
891        let mut elements = iss_elements();
892        elements.catalog_number = 100_000;
893        assert_invalid_input(
894            Satellite::from_elements(&elements),
895            "element.catalog_number",
896            Sgp4InputErrorKind::OutOfRange,
897        );
898    }
899
900    #[test]
901    fn from_elements_rejects_invalid_epoch() {
902        let mut elements = iss_elements();
903        elements.epoch = JulianDate(f64::NAN, 0.0);
904        assert_invalid_input(
905            Satellite::from_elements(&elements),
906            "element.epoch.whole",
907            Sgp4InputErrorKind::NonFinite,
908        );
909
910        let mut elements = iss_elements();
911        elements.epoch = JulianDate(9_000_000.0, 0.0);
912        assert_invalid_input(
913            Satellite::from_elements(&elements),
914            "element.epoch",
915            Sgp4InputErrorKind::OutOfRange,
916        );
917    }
918
919    #[test]
920    fn from_elements_accepts_full_julian_epoch() {
921        let mut elements = iss_elements();
922        elements.epoch = super::sgp4_julian_date_from_calendar(2057, 1, 1, 0, 0, 0.0);
923        Satellite::from_elements(&elements).expect("full 2057 epoch is valid");
924    }
925
926    #[test]
927    fn from_tle_accepts_epoch_after_parser_conversion_to_full_jd() {
928        let mut line1 = ISS_L1.to_string();
929        line1.replace_range(18..32, "19366.00000000");
930
931        Satellite::from_tle(&line1, ISS_L2).expect("TLE epoch is converted to full JD");
932    }
933
934    #[test]
935    fn propagation_rejects_non_finite_time_inputs() {
936        let sat = Satellite::from_tle(ISS_L1, ISS_L2).unwrap();
937        assert_invalid_input(
938            sat.propagate(MinutesSinceEpoch(f64::NAN)),
939            "minutes_since_epoch",
940            Sgp4InputErrorKind::NonFinite,
941        );
942        assert_invalid_input(
943            sat.propagate_jd(JulianDate(f64::INFINITY, 0.0)),
944            "julian_date.whole",
945            Sgp4InputErrorKind::NonFinite,
946        );
947
948        let elements = iss_elements();
949        assert_invalid_input(
950            propagate_elements(&elements, MinutesSinceEpoch(f64::INFINITY)),
951            "minutes_since_epoch",
952            Sgp4InputErrorKind::NonFinite,
953        );
954    }
955
956    #[test]
957    fn propagation_rejects_out_of_domain_time_inputs() {
958        let sat = Satellite::from_tle(ISS_L1, ISS_L2).unwrap();
959        assert_invalid_input(
960            sat.propagate(MinutesSinceEpoch(MAX_MINUTES_SINCE_EPOCH.next_up())),
961            "minutes_since_epoch",
962            Sgp4InputErrorKind::OutOfRange,
963        );
964        assert_invalid_input(
965            sat.propagate_jd(JulianDate(2_458_304.0, 1.0)),
966            "julian_date.fraction",
967            Sgp4InputErrorKind::OutOfRange,
968        );
969    }
970
971    #[test]
972    fn lenient_trailing_whitespace_and_content_past_col_69() {
973        let clean = Satellite::from_tle(ISS_L1, ISS_L2).unwrap();
974
975        // Trailing whitespace on both lines.
976        let pad = Satellite::from_tle(&format!("{ISS_L1}   "), &format!("{ISS_L2}\t ")).unwrap();
977        assert_same(&clean, &pad);
978
979        // Extra content past the 69-column record (CelesTrak/Space-Track blobs
980        // sometimes carry it); it is trimmed before parsing.
981        let extra =
982            Satellite::from_tle(&format!("{ISS_L1} EXTRA-JUNK"), &format!("{ISS_L2} 999999"))
983                .unwrap();
984        assert_same(&clean, &extra);
985    }
986
987    #[test]
988    fn lenient_leading_dot_and_assumed_decimal_fields() {
989        // The ISS TLE already carries a leading-dot first derivative
990        // (` .00001614`), an assumed-decimal B\* (`31745-4`), and the implicit
991        // `0.` eccentricity (`0003435`). A successful parse + finite LEO state
992        // proves those normalizations run through the public entry point.
993        let sat = Satellite::from_tle(ISS_L1, ISS_L2).unwrap();
994        let p = sat.propagate(MinutesSinceEpoch(0.0)).unwrap();
995        let r = (p.position[0].powi(2) + p.position[1].powi(2) + p.position[2].powi(2)).sqrt();
996        assert!(
997            (6500.0..=7200.0).contains(&r),
998            "ISS radius {r} km outside LEO"
999        );
1000    }
1001
1002    #[test]
1003    fn lenient_missing_optional_bookkeeping_fields() {
1004        // Blank element-set number, ephemeris type, and revolution number are
1005        // cosmetic for propagation and must default rather than reject. Build
1006        // such a TLE by blanking those columns (cols 63, 65-68 on line 1 and
1007        // 64-68 on line 2) while leaving the orbital fields intact.
1008        let l1: String = ISS_L1
1009            .char_indices()
1010            .map(|(i, c)| {
1011                if i == 62 || (64..=67).contains(&i) {
1012                    ' '
1013                } else {
1014                    c
1015                }
1016            })
1017            .collect();
1018        let l2: String = ISS_L2
1019            .char_indices()
1020            .map(|(i, c)| if (63..=67).contains(&i) { ' ' } else { c })
1021            .collect();
1022        // Same orbital elements as the clean TLE → bit-identical propagation
1023        // (the blanked fields do not feed SGP4).
1024        let clean = Satellite::from_tle(ISS_L1, ISS_L2).unwrap();
1025        let blanked = Satellite::from_tle(&l1, &l2).unwrap();
1026        assert_same(&clean, &blanked);
1027    }
1028
1029    #[test]
1030    fn three_line_form_strips_name_line() {
1031        let clean = Satellite::from_tle(ISS_L1, ISS_L2).unwrap();
1032
1033        let block = format!("ISS (ZARYA)\n{ISS_L1}\n{ISS_L2}\n");
1034        let three = Satellite::from_3line(&block).unwrap();
1035        assert_same(&clean, &three);
1036
1037        // A plain two-line block (no name line) also parses.
1038        let two = Satellite::from_3line(&format!("{ISS_L1}\n{ISS_L2}")).unwrap();
1039        assert_same(&clean, &two);
1040    }
1041
1042    #[test]
1043    fn three_line_form_rejects_block_without_element_lines() {
1044        assert!(Satellite::from_3line("just a name\nand some text").is_err());
1045        assert!(Satellite::from_3line("").is_err());
1046    }
1047
1048    #[test]
1049    fn rejects_genuine_corruption() {
1050        // Empty.
1051        assert!(Satellite::from_tle("", "").is_err());
1052        // Non-TLE text.
1053        assert!(Satellite::from_tle("hello world", "goodbye world").is_err());
1054        // Swapped lines (line 2 first).
1055        assert!(Satellite::from_tle(ISS_L2, ISS_L1).is_err());
1056        // Mismatched satellite numbers between line 1 and line 2.
1057        let l2_wrong = "2 25545  51.6414 295.8524 0003435 262.6267 204.2868 15.54005638121106";
1058        assert!(matches!(
1059            Satellite::from_tle(ISS_L1, l2_wrong),
1060            Err(Error::InvalidTle(_))
1061        ));
1062    }
1063}