space_dust/
tle.rs

1//! Two-Line Element Set (TLE) parsing and SGP4 propagation.
2//!
3//! This module provides:
4//! - TLE parsing from standard two-line format
5//! - SGP4 satellite propagation using the `sgp4` crate
6//! - Conversion to state vectors in TEME frame
7//!
8//! ## Example
9//!
10//! ```rust,no_run
11//! use space_dust::tle::Tle;
12//! use chrono::Utc;
13//!
14//! let line1 = "1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9002";
15//! let line2 = "2 25544  51.6400 208.1200 0001234  85.0000 275.0000 15.48919100123456";
16//!
17//! let tle = Tle::parse(line1, line2).unwrap();
18//! println!("Satellite: {}", tle.catalog_number());
19//! println!("Epoch: {:?}", tle.epoch());
20//!
21//! // Propagate to current time
22//! let state = tle.propagate(&Utc::now()).unwrap();
23//! println!("Position: {:?}", state.position);
24//! ```
25
26use crate::math::Vector3;
27use crate::state::TEMEState;
28use chrono::{DateTime, TimeZone, Utc};
29use sgp4::{Constants as Sgp4Constants, Elements, Prediction};
30use std::error::Error;
31use std::fmt;
32
33// ============================================================================
34// Error Types
35// ============================================================================
36
37/// Errors that can occur during TLE parsing or propagation.
38#[derive(Debug)]
39pub enum TleError {
40    /// TLE line has incorrect length
41    InvalidLineLength {
42        line: u8,
43        expected: usize,
44        actual: usize,
45    },
46    /// Failed to parse a field in the TLE
47    ParseError { field: String, message: String },
48    /// SGP4 propagation error
49    PropagationError { message: String },
50    /// Elements are invalid for SGP4
51    InvalidElements { message: String },
52}
53
54impl fmt::Display for TleError {
55    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
56        match self {
57            TleError::InvalidLineLength {
58                line,
59                expected,
60                actual,
61            } => {
62                write!(
63                    f,
64                    "TLE line {} has invalid length: expected {}, got {}",
65                    line, expected, actual
66                )
67            }
68            TleError::ParseError { field, message } => {
69                write!(f, "Failed to parse TLE field '{}': {}", field, message)
70            }
71            TleError::PropagationError { message } => {
72                write!(f, "SGP4 propagation error: {}", message)
73            }
74            TleError::InvalidElements { message } => {
75                write!(f, "Invalid TLE elements: {}", message)
76            }
77        }
78    }
79}
80
81impl Error for TleError {}
82
83/// Result type for TLE operations.
84pub type TleResult<T> = Result<T, TleError>;
85
86// ============================================================================
87// TLE Structure
88// ============================================================================
89
90/// Two-Line Element Set (TLE) data structure.
91///
92/// A TLE is a standardized format for describing the orbit of an Earth-orbiting
93/// object. TLEs are used with the SGP4/SDP4 propagators to predict satellite
94/// positions at future times.
95#[derive(Debug, Clone)]
96pub struct Tle {
97    /// Raw TLE line 1
98    line1: String,
99    /// Raw TLE line 2
100    line2: String,
101    /// NORAD catalog number
102    catalog_number: String,
103    /// Security classification (U = Unclassified)
104    classification: char,
105    /// International designator
106    international_designator: String,
107    /// TLE epoch as `DateTime<Utc>`
108    epoch: DateTime<Utc>,
109    /// First derivative of mean motion (rev/day²)
110    mean_motion_dot: f64,
111    /// Second derivative of mean motion (rev/day³)
112    mean_motion_double_dot: f64,
113    /// BSTAR drag coefficient (1/Earth radii)
114    bstar: f64,
115    /// Ephemeris type (usually 0)
116    ephemeris_type: u8,
117    /// Element set number
118    element_set_number: u32,
119    /// Orbital inclination (degrees)
120    inclination_deg: f64,
121    /// Right Ascension of Ascending Node (degrees)
122    raan_deg: f64,
123    /// Orbital eccentricity
124    eccentricity: f64,
125    /// Argument of perigee (degrees)
126    arg_perigee_deg: f64,
127    /// Mean anomaly at epoch (degrees)
128    mean_anomaly_deg: f64,
129    /// Mean motion (revolutions per day)
130    mean_motion: f64,
131    /// Revolution number at epoch
132    rev_number: u32,
133    /// Parsed SGP4 elements (kept for potential future use/debugging)
134    #[allow(dead_code)]
135    sgp4_elements: Elements,
136    /// SGP4 propagator constants
137    sgp4_constants: Sgp4Constants,
138}
139
140impl Tle {
141    /// Parse a TLE from two lines.
142    ///
143    /// # Arguments
144    /// * `line1` - First line of the TLE (69 characters)
145    /// * `line2` - Second line of the TLE (69 characters)
146    ///
147    /// # Returns
148    /// * `Ok(Tle)` - Successfully parsed TLE
149    /// * `Err(TleError)` - Parsing failed
150    pub fn parse(line1: &str, line2: &str) -> TleResult<Self> {
151        // Validate line lengths
152        if line1.len() != 69 {
153            return Err(TleError::InvalidLineLength {
154                line: 1,
155                expected: 69,
156                actual: line1.len(),
157            });
158        }
159        if line2.len() != 69 {
160            return Err(TleError::InvalidLineLength {
161                line: 2,
162                expected: 69,
163                actual: line2.len(),
164            });
165        }
166
167        // Parse epoch
168        let epoch_year = parse_field::<i32>(&line1[18..20], "epoch_year")?;
169        let epoch_day = parse_field::<f64>(&line1[20..32], "epoch_day")?;
170
171        let full_year = if epoch_year > 56 {
172            1900 + epoch_year
173        } else {
174            2000 + epoch_year
175        };
176
177        let epoch = epoch_from_year_and_day(full_year, epoch_day)?;
178
179        // Parse mean motion derivatives
180        let mean_motion_dot = parse_mean_motion_dot(&line1[33..43])?;
181        let mean_motion_double_dot = parse_exponential(&line1[44..52])?;
182        let bstar = parse_exponential(&line1[53..61])?;
183
184        // Parse other line 1 fields
185        let catalog_number = line1[2..7].trim().to_string();
186        let classification = line1.chars().nth(7).unwrap_or('U');
187        let international_designator = line1[9..17].trim().to_string();
188        let ephemeris_type = parse_field::<u8>(&line1[62..63], "ephemeris_type").unwrap_or(0);
189        let element_set_number = parse_field::<u32>(line1[64..68].trim(), "element_set_number")?;
190
191        // Parse line 2 fields
192        let inclination_deg = parse_field::<f64>(line2[8..16].trim(), "inclination")?;
193        let raan_deg = parse_field::<f64>(line2[17..25].trim(), "raan")?;
194        let eccentricity =
195            parse_field::<f64>(&format!("0.{}", &line2[26..33].trim()), "eccentricity")?;
196        let arg_perigee_deg = parse_field::<f64>(line2[34..42].trim(), "arg_perigee")?;
197        let mean_anomaly_deg = parse_field::<f64>(line2[43..51].trim(), "mean_anomaly")?;
198        let mean_motion = parse_field::<f64>(line2[52..63].trim(), "mean_motion")?;
199        let rev_number = parse_field::<u32>(line2[63..68].trim(), "rev_number")?;
200
201        // Parse using sgp4 crate for propagation
202        let sgp4_elements = sgp4::Elements::from_tle(None, line1.as_bytes(), line2.as_bytes())
203            .map_err(|e| TleError::ParseError {
204                field: "sgp4_elements".to_string(),
205                message: format!("{:?}", e),
206            })?;
207
208        // Create SGP4 constants
209        let sgp4_constants = Sgp4Constants::from_elements(&sgp4_elements).map_err(|e| {
210            TleError::InvalidElements {
211                message: format!("{:?}", e),
212            }
213        })?;
214
215        Ok(Self {
216            line1: line1.to_string(),
217            line2: line2.to_string(),
218            catalog_number,
219            classification,
220            international_designator,
221            epoch,
222            mean_motion_dot,
223            mean_motion_double_dot,
224            bstar,
225            ephemeris_type,
226            element_set_number,
227            inclination_deg,
228            raan_deg,
229            eccentricity,
230            arg_perigee_deg,
231            mean_anomaly_deg,
232            mean_motion,
233            rev_number,
234            sgp4_elements,
235            sgp4_constants,
236        })
237    }
238
239    /// Parse a TLE from three lines (name + two TLE lines).
240    ///
241    /// # Arguments
242    /// * `name` - Satellite name (ignored)
243    /// * `line1` - First line of the TLE
244    /// * `line2` - Second line of the TLE
245    pub fn parse_3le(_name: &str, line1: &str, line2: &str) -> TleResult<Self> {
246        Self::parse(line1, line2)
247    }
248
249    // ========================================================================
250    // Accessors
251    // ========================================================================
252
253    /// Get the raw TLE line 1.
254    pub fn line1(&self) -> &str {
255        &self.line1
256    }
257
258    /// Get the raw TLE line 2.
259    pub fn line2(&self) -> &str {
260        &self.line2
261    }
262
263    /// Get the NORAD catalog number.
264    pub fn catalog_number(&self) -> &str {
265        &self.catalog_number
266    }
267
268    /// Get the security classification.
269    pub fn classification(&self) -> char {
270        self.classification
271    }
272
273    /// Get the international designator.
274    pub fn international_designator(&self) -> &str {
275        &self.international_designator
276    }
277
278    /// Get the TLE epoch.
279    pub fn epoch(&self) -> DateTime<Utc> {
280        self.epoch
281    }
282
283    /// Get the first derivative of mean motion (rev/day²).
284    pub fn mean_motion_dot(&self) -> f64 {
285        self.mean_motion_dot
286    }
287
288    /// Get the second derivative of mean motion (rev/day³).
289    pub fn mean_motion_double_dot(&self) -> f64 {
290        self.mean_motion_double_dot
291    }
292
293    /// Get the BSTAR drag coefficient.
294    pub fn bstar(&self) -> f64 {
295        self.bstar
296    }
297
298    /// Get the ephemeris type.
299    pub fn ephemeris_type(&self) -> u8 {
300        self.ephemeris_type
301    }
302
303    /// Get the element set number.
304    pub fn element_set_number(&self) -> u32 {
305        self.element_set_number
306    }
307
308    /// Get the orbital inclination in degrees.
309    pub fn inclination_deg(&self) -> f64 {
310        self.inclination_deg
311    }
312
313    /// Get the Right Ascension of Ascending Node in degrees.
314    pub fn raan_deg(&self) -> f64 {
315        self.raan_deg
316    }
317
318    /// Get the orbital eccentricity.
319    pub fn eccentricity(&self) -> f64 {
320        self.eccentricity
321    }
322
323    /// Get the argument of perigee in degrees.
324    pub fn arg_perigee_deg(&self) -> f64 {
325        self.arg_perigee_deg
326    }
327
328    /// Get the mean anomaly at epoch in degrees.
329    pub fn mean_anomaly_deg(&self) -> f64 {
330        self.mean_anomaly_deg
331    }
332
333    /// Get the mean motion in revolutions per day.
334    pub fn mean_motion(&self) -> f64 {
335        self.mean_motion
336    }
337
338    /// Get the revolution number at epoch.
339    pub fn rev_number(&self) -> u32 {
340        self.rev_number
341    }
342
343    /// Get the orbital period in minutes.
344    pub fn period_minutes(&self) -> f64 {
345        1440.0 / self.mean_motion
346    }
347
348    /// Get the orbital period in seconds.
349    pub fn period_seconds(&self) -> f64 {
350        86400.0 / self.mean_motion
351    }
352
353    // ========================================================================
354    // Propagation
355    // ========================================================================
356
357    /// Propagate the TLE to a specific epoch using SGP4.
358    ///
359    /// Returns the satellite state in the TEME (True Equator Mean Equinox) frame.
360    ///
361    /// # Arguments
362    /// * `epoch` - Target epoch for propagation
363    ///
364    /// # Returns
365    /// * `Ok(TEMEState)` - Position and velocity in TEME frame
366    /// * `Err(TleError)` - Propagation failed
367    pub fn propagate(&self, epoch: &DateTime<Utc>) -> TleResult<TEMEState> {
368        // Calculate minutes since TLE epoch
369        let minutes_since_epoch = self.minutes_since_epoch(epoch);
370
371        // Propagate using SGP4
372        let prediction = self
373            .sgp4_constants
374            .propagate(sgp4::MinutesSinceEpoch(minutes_since_epoch))
375            .map_err(|e| TleError::PropagationError {
376                message: format!("{:?}", e),
377            })?;
378
379        // Convert to TEMEState
380        Ok(prediction_to_teme_state(&prediction, *epoch))
381    }
382
383    /// Propagate the TLE by a number of minutes from the TLE epoch.
384    ///
385    /// # Arguments
386    /// * `minutes` - Minutes since TLE epoch (can be negative)
387    pub fn propagate_minutes(&self, minutes: f64) -> TleResult<TEMEState> {
388        let prediction = self
389            .sgp4_constants
390            .propagate(sgp4::MinutesSinceEpoch(minutes))
391            .map_err(|e| TleError::PropagationError {
392                message: format!("{:?}", e),
393            })?;
394
395        // Calculate the target epoch
396        let target_epoch =
397            self.epoch + chrono::Duration::milliseconds((minutes * 60.0 * 1000.0) as i64);
398
399        Ok(prediction_to_teme_state(&prediction, target_epoch))
400    }
401
402    /// Get position and velocity at a specific epoch.
403    ///
404    /// Returns ((x, y, z), (vx, vy, vz)) in km and km/s.
405    pub fn get_rv_at_time(
406        &self,
407        epoch: &DateTime<Utc>,
408    ) -> TleResult<((f64, f64, f64), (f64, f64, f64))> {
409        let state = self.propagate(epoch)?;
410        Ok((state.position.to_tuple(), state.velocity.to_tuple()))
411    }
412
413    /// Calculate minutes since TLE epoch.
414    pub fn minutes_since_epoch(&self, target: &DateTime<Utc>) -> f64 {
415        let duration = target.signed_duration_since(self.epoch);
416        duration.num_milliseconds() as f64 / 60_000.0
417    }
418
419    /// Check if propagation to a given epoch is likely to be accurate.
420    ///
421    /// TLEs are typically valid for a few days before/after their epoch.
422    /// This function returns a warning level based on the time difference.
423    ///
424    /// # Returns
425    /// * `0` - Within ±1 day (good accuracy)
426    /// * `1` - Within ±3 days (acceptable)
427    /// * `2` - Within ±7 days (degraded)
428    /// * `3` - Beyond ±7 days (poor accuracy)
429    pub fn accuracy_warning(&self, epoch: &DateTime<Utc>) -> u8 {
430        let days = self.minutes_since_epoch(epoch).abs() / 1440.0;
431        if days <= 1.0 {
432            0
433        } else if days <= 3.0 {
434            1
435        } else if days <= 7.0 {
436            2
437        } else {
438            3
439        }
440    }
441}
442
443// ============================================================================
444// Helper Functions
445// ============================================================================
446
447/// Parse a field from a TLE string.
448fn parse_field<T: std::str::FromStr>(s: &str, field_name: &str) -> TleResult<T> {
449    s.trim().parse().map_err(|_| TleError::ParseError {
450        field: field_name.to_string(),
451        message: format!("Could not parse '{}'", s.trim()),
452    })
453}
454
455/// Parse the mean motion derivative field (special format).
456fn parse_mean_motion_dot(s: &str) -> TleResult<f64> {
457    let trimmed = s.trim();
458    if trimmed.is_empty() || trimmed == "." {
459        return Ok(0.0);
460    }
461
462    // Handle sign
463    let (sign, value_str) = if trimmed.starts_with('-') {
464        (-1.0, &trimmed[1..])
465    } else if trimmed.starts_with('+') {
466        (1.0, &trimmed[1..])
467    } else {
468        (1.0, trimmed)
469    };
470
471    // Parse as decimal
472    let value: f64 = if value_str.starts_with('.') {
473        format!("0{}", value_str)
474    } else {
475        value_str.to_string()
476    }
477    .trim()
478    .parse()
479    .map_err(|_| TleError::ParseError {
480        field: "mean_motion_dot".to_string(),
481        message: format!("Could not parse '{}'", s),
482    })?;
483
484    Ok(sign * value)
485}
486
487/// Parse an exponential notation field (like BSTAR).
488fn parse_exponential(s: &str) -> TleResult<f64> {
489    let trimmed = s.trim();
490    if trimmed.is_empty() {
491        return Ok(0.0);
492    }
493
494    // Handle sign
495    let (sign, rest) = if trimmed.starts_with('-') {
496        (-1.0, &trimmed[1..])
497    } else if trimmed.starts_with('+') {
498        (1.0, &trimmed[1..])
499    } else if trimmed.starts_with(' ') {
500        (1.0, &trimmed[1..])
501    } else {
502        (1.0, trimmed)
503    };
504
505    // TLE format: mantissa followed by exponent (e.g., "12345-4" means 0.12345e-4)
506    if rest.len() < 2 {
507        return Ok(0.0);
508    }
509
510    // Find the exponent (last 2 characters including sign)
511    let exp_start = rest.len() - 2;
512    let mantissa_str = &rest[..exp_start];
513    let exp_str = &rest[exp_start..];
514
515    // Parse mantissa as 0.xxxxx
516    let mantissa: f64 = format!("0.{}", mantissa_str.trim()).parse().unwrap_or(0.0);
517
518    // Parse exponent
519    let exponent: i32 = exp_str.trim().parse().unwrap_or(0);
520
521    Ok(sign * mantissa * 10.0_f64.powi(exponent))
522}
523
524/// Convert year and day of year to DateTime<Utc>.
525fn epoch_from_year_and_day(year: i32, day_of_year: f64) -> TleResult<DateTime<Utc>> {
526    // Day 1 is January 1st, so we subtract 1 to get days to add
527    let days_to_add = (day_of_year - 1.0).floor() as i64;
528    let fraction = day_of_year - day_of_year.floor();
529    let microseconds = (fraction * 86_400_000_000.0) as i64;
530
531    let start_of_year = Utc
532        .with_ymd_and_hms(year, 1, 1, 0, 0, 0)
533        .single()
534        .ok_or_else(|| TleError::ParseError {
535            field: "epoch".to_string(),
536            message: format!("Invalid year: {}", year),
537        })?;
538
539    let epoch = start_of_year
540        + chrono::Duration::days(days_to_add)
541        + chrono::Duration::microseconds(microseconds);
542
543    Ok(epoch)
544}
545
546/// Convert SGP4 prediction to TEMEState.
547fn prediction_to_teme_state(prediction: &Prediction, epoch: DateTime<Utc>) -> TEMEState {
548    let position = Vector3::new(
549        prediction.position[0],
550        prediction.position[1],
551        prediction.position[2],
552    );
553    let velocity = Vector3::new(
554        prediction.velocity[0],
555        prediction.velocity[1],
556        prediction.velocity[2],
557    );
558    TEMEState::new(epoch, position, velocity)
559}
560
561#[cfg(test)]
562mod tests {
563    use super::*;
564    use chrono::Datelike;
565
566    // Real ISS TLE from Celestrak with correct checksums
567    const ISS_LINE1: &str = "1 25544U 98067A   26014.62805721  .00006818  00000+0  13044-3 0  9991";
568    const ISS_LINE2: &str = "2 25544  51.6333 339.6562 0007763  17.9854 342.1408 15.49289811547943";
569
570    #[test]
571    fn test_tle_parse() {
572        let tle = Tle::parse(ISS_LINE1, ISS_LINE2);
573        assert!(tle.is_ok());
574
575        let tle = tle.unwrap();
576        assert_eq!(tle.catalog_number(), "25544");
577        assert_eq!(tle.classification(), 'U');
578        assert!(tle.inclination_deg() > 51.0 && tle.inclination_deg() < 52.0);
579        assert!(tle.eccentricity() < 0.01); // Nearly circular
580    }
581
582    #[test]
583    fn test_tle_invalid_length() {
584        let short_line = "1 25544U";
585        let result = Tle::parse(short_line, ISS_LINE2);
586        assert!(result.is_err());
587
588        if let Err(TleError::InvalidLineLength { line, .. }) = result {
589            assert_eq!(line, 1);
590        } else {
591            panic!("Expected InvalidLineLength error");
592        }
593    }
594
595    #[test]
596    fn test_tle_epoch() {
597        let tle = Tle::parse(ISS_LINE1, ISS_LINE2).unwrap();
598        let epoch = tle.epoch();
599
600        // Should be in 2026
601        assert_eq!(epoch.year(), 2026);
602        assert_eq!(epoch.month(), 1);
603    }
604
605    #[test]
606    fn test_tle_period() {
607        let tle = Tle::parse(ISS_LINE1, ISS_LINE2).unwrap();
608
609        // ISS period is about 92 minutes
610        let period = tle.period_minutes();
611        assert!(period > 90.0 && period < 95.0);
612    }
613
614    #[test]
615    fn test_tle_propagate() {
616        let tle = Tle::parse(ISS_LINE1, ISS_LINE2).unwrap();
617
618        // Propagate to TLE epoch (0 minutes)
619        let state = tle.propagate_minutes(0.0);
620        assert!(state.is_ok());
621
622        let state = state.unwrap();
623
624        // Position should be roughly in LEO (300-450 km altitude, so ~6700-6800 km from center)
625        let radius = state.position.magnitude();
626        assert!(radius > 6500.0 && radius < 7000.0);
627
628        // Velocity should be roughly 7.5 km/s for LEO
629        let speed = state.velocity.magnitude();
630        assert!(speed > 7.0 && speed < 8.0);
631    }
632
633    #[test]
634    fn test_tle_propagate_forward() {
635        let tle = Tle::parse(ISS_LINE1, ISS_LINE2).unwrap();
636
637        // Propagate forward 1 orbit (about 92 minutes)
638        let state0 = tle.propagate_minutes(0.0).unwrap();
639        let state1 = tle.propagate_minutes(92.0).unwrap();
640
641        // After one orbit, should be approximately back to same position
642        // (not exact due to J2 perturbations and precession)
643        let diff = state0.position - state1.position;
644        let distance = diff.magnitude();
645
646        // Should be within ~200 km after one orbit
647        assert!(distance < 500.0);
648    }
649
650    #[test]
651    fn test_tle_accuracy_warning() {
652        let tle = Tle::parse(ISS_LINE1, ISS_LINE2).unwrap();
653
654        // At TLE epoch
655        assert_eq!(tle.accuracy_warning(&tle.epoch()), 0);
656
657        // 2 days later
658        let later = tle.epoch() + chrono::Duration::days(2);
659        assert_eq!(tle.accuracy_warning(&later), 1);
660
661        // 5 days later
662        let much_later = tle.epoch() + chrono::Duration::days(5);
663        assert_eq!(tle.accuracy_warning(&much_later), 2);
664
665        // 10 days later
666        let way_later = tle.epoch() + chrono::Duration::days(10);
667        assert_eq!(tle.accuracy_warning(&way_later), 3);
668    }
669
670    #[test]
671    fn test_parse_exponential() {
672        // Test BSTAR-style exponential parsing
673        assert!((parse_exponential(" 10270-3").unwrap() - 0.0001027).abs() < 1e-10);
674        assert!((parse_exponential("-10270-3").unwrap() - (-0.0001027)).abs() < 1e-10);
675        assert!((parse_exponential(" 00000+0").unwrap()).abs() < 1e-15);
676    }
677
678    #[test]
679    fn test_get_rv_at_time() {
680        let tle = Tle::parse(ISS_LINE1, ISS_LINE2).unwrap();
681        let result = tle.get_rv_at_time(&tle.epoch());
682
683        assert!(result.is_ok());
684        let ((x, y, z), (vx, vy, vz)) = result.unwrap();
685
686        // Position magnitude check
687        let r = (x * x + y * y + z * z).sqrt();
688        assert!(r > 6500.0 && r < 7000.0);
689
690        // Velocity magnitude check
691        let v = (vx * vx + vy * vy + vz * vz).sqrt();
692        assert!(v > 7.0 && v < 8.0);
693    }
694}
space_dust/tle.rs

space_dust/
tle.rs
