astrodynamics 0.6.0

//! SGP4 satellite propagation with 0 ULP parity to the Vallado C++ reference
//! implementation (v2020-07-13).
//!
//! This module is a faithful pure-Rust port of Vallado's SGP4 verified bit-for-bit
//! against the canonical 33-satellite / 198-propagation-point Vallado verification
//! suite. TLEs are curve-fit parameters generated by Vallado's SGP4; using a
//! different propagator introduces errors. This module preserves the exact
//! floating-point computation order of the C++ reference so output matches
//! Python's `sgp4` C extension (which compiles the same source) bit-for-bit.
//!
//! ## Quick start
//!
//! ```
//! use astrodynamics::sgp4::{Satellite, MinutesSinceEpoch};
//!
//! let line1 = "1 25544U 98067A   18184.80969102  .00001614  00000-0  31745-4 0  9993";
//! let line2 = "2 25544  51.6414 295.8524 0003435 262.6267 204.2868 15.54005638121106";
//!
//! let sat = Satellite::from_tle(line1, line2).unwrap();
//! let pred = sat.propagate(MinutesSinceEpoch(0.0)).unwrap();
//!
//! // position in km, velocity in km/s, TEME frame
//! let _ = pred.position;
//! let _ = pred.velocity;
//! ```

#[allow(
    dead_code,
    unused_variables,
    unused_assignments,
    unused_mut,
    non_snake_case,
    non_camel_case_types,
    clippy::approx_constant,
    clippy::excessive_precision,
    clippy::too_many_arguments,
    clippy::needless_return,
    clippy::assign_op_pattern,
    clippy::manual_range_contains,
    clippy::collapsible_if,
    clippy::collapsible_else_if,
    clippy::float_cmp,
    clippy::needless_late_init,
    clippy::field_reassign_with_default
)]
mod vallado;

use thiserror::Error;

// ── Error ────────────────────────────────────────────────────────────

/// Error from SGP4 initialization or propagation.
#[derive(Error, Debug, Clone, PartialEq)]
pub enum Error {
    /// TLE line has invalid format.
    #[error("invalid TLE: {0}")]
    InvalidTle(String),
    /// SGP4 returned a non-zero error code.
    ///
    /// Codes: 1 = mean elements, 2 = mean motion, 3 = perturbed elements,
    /// 4 = semi-latus rectum, 5 = epoch elements sub-orbital,
    /// 6 = satellite decayed.
    #[error("SGP4 error code {code}")]
    Sgp4 { code: i32 },
}

// ── Types ────────────────────────────────────────────────────────────

/// Minutes since the TLE epoch. Newtype to prevent mixing with raw `f64`.
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct MinutesSinceEpoch(pub f64);

/// Propagation result in the TEME (True Equator, Mean Equinox) frame.
#[derive(Debug, Clone, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Prediction {
    /// Position in km, TEME frame.
    pub position: [f64; 3],
    /// Velocity in km/s, TEME frame.
    pub velocity: [f64; 3],
}

/// Julian date split as `(whole, fraction)` for high-precision time input.
///
/// Skyfield convention: `whole = floor(JD)`, `fraction = remainder`.
/// For example, 2018-07-04 00:00:00 UTC = `JulianDate(2458303.0, 0.5)`.
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct JulianDate(pub f64, pub f64);

// ── Satellite ────────────────────────────────────────────────────────

/// A parsed TLE ready for propagation.
///
/// Holds the raw TLE lines plus the initialized SGP4 satellite record. The
/// TLE is parsed and `sgp4init` is run exactly once during `from_tle`, so
/// subsequent `propagate` calls just invoke the propagation kernel directly
/// — fast, and crucially, with no precision loss from JD round-tripping.
#[derive(Clone)]
pub struct Satellite {
    line1: String,
    line2: String,
    /// Pre-initialized satellite record. Boxed to keep `Satellite` small on
    /// the stack — `ElsetRec` has ~150 fields.
    satrec: Box<vallado::ElsetRec>,
}

impl std::fmt::Debug for Satellite {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("Satellite")
            .field("line1", &self.line1)
            .field("line2", &self.line2)
            .finish_non_exhaustive()
    }
}

impl Satellite {
    /// Parse a two-line element set.
    ///
    /// Lines should be the standard 69-character TLE format. Leading and
    /// trailing whitespace is trimmed before length validation. Runs the full
    /// SGP4 initialization (`sgp4init`) once and caches the resulting state
    /// so propagation calls are pure step kernels.
    pub fn from_tle(line1: &str, line2: &str) -> Result<Self, Error> {
        let l1 = line1.trim();
        let l2 = line2.trim();

        if l1.len() < 69 || !l1.starts_with('1') {
            return Err(Error::InvalidTle(
                "line 1 too short or doesn't start with '1'".into(),
            ));
        }
        if l2.len() < 69 || !l2.starts_with('2') {
            return Err(Error::InvalidTle(
                "line 2 too short or doesn't start with '2'".into(),
            ));
        }

        let satrec = init_satrec_from_tle(l1, l2)?;

        Ok(Satellite {
            line1: l1.to_string(),
            line2: l2.to_string(),
            satrec: Box::new(satrec),
        })
    }

    /// Propagate to a time given as minutes since the TLE epoch.
    ///
    /// Calls the SGP4 step kernel directly with the supplied tsince — no JD
    /// round-trip, no precision loss.
    pub fn propagate(&self, t: MinutesSinceEpoch) -> Result<Prediction, Error> {
        // Clone the satrec so propagation doesn't mutate the cached state
        // (sgp4 writes back into the satrec — error code, atime, etc.).
        let mut working = (*self.satrec).clone();
        let mut r = [0.0_f64; 3];
        let mut v = [0.0_f64; 3];
        let ok = vallado::sgp4(&mut working, t.0, &mut r, &mut v);
        if !ok || working.error != 0 {
            return Err(Error::Sgp4 {
                code: working.error,
            });
        }
        Ok(Prediction {
            position: r,
            velocity: v,
        })
    }

    /// Propagate to a Julian date, split as `(whole, fraction)`.
    ///
    /// Computes the tsince from the cached epoch via the same subtraction
    /// the C++ wrapper uses:
    ///
    /// ```text
    /// tsince = (jd - jdsatepoch) * 1440 + (fr - jdsatepochF) * 1440
    /// ```
    pub fn propagate_jd(&self, jd: JulianDate) -> Result<Prediction, Error> {
        let tsince = (jd.0 - self.satrec.jdsatepoch) * 1440.0
            + (jd.1 - self.satrec.jdsatepochF) * 1440.0;
        self.propagate(MinutesSinceEpoch(tsince))
    }

    /// Raw TLE line 1.
    pub fn line1(&self) -> &str {
        &self.line1
    }

    /// Raw TLE line 2.
    pub fn line2(&self) -> &str {
        &self.line2
    }
}

#[cfg(feature = "serde")]
impl serde::Serialize for Satellite {
    fn serialize<S: serde::Serializer>(&self, s: S) -> Result<S::Ok, S::Error> {
        use serde::ser::SerializeStruct;
        let mut st = s.serialize_struct("Satellite", 2)?;
        st.serialize_field("line1", &self.line1)?;
        st.serialize_field("line2", &self.line2)?;
        st.end()
    }
}

#[cfg(feature = "serde")]
impl<'de> serde::Deserialize<'de> for Satellite {
    fn deserialize<D: serde::Deserializer<'de>>(d: D) -> Result<Self, D::Error> {
        #[derive(serde::Deserialize)]
        struct Wire {
            line1: String,
            line2: String,
        }
        let w = Wire::deserialize(d)?;
        Satellite::from_tle(&w.line1, &w.line2).map_err(serde::de::Error::custom)
    }
}

// ── Internal helpers ─────────────────────────────────────────────────

/// Parse a TLE and run `sgp4init`, returning the initialized satellite
/// record. Mirrors the parse logic in `vallado::twoline2rv_propagate` exactly.
fn init_satrec_from_tle(line1: &str, line2: &str) -> Result<vallado::ElsetRec, Error> {
    let l1 = line1.trim_end();
    let l2 = line2.trim_end();
    if l1.len() < 64 || l2.len() < 68 {
        return Err(Error::InvalidTle("TLE lines too short".into()));
    }

    let deg2rad = std::f64::consts::PI / 180.0;
    let xpdotp = 1440.0 / (2.0 * std::f64::consts::PI);

    // Line 1 fields
    let two_digit_year: i32 = l1[18..20]
        .trim()
        .parse()
        .map_err(|_| Error::InvalidTle("bad epochyr".into()))?;
    let epochdays: f64 = l1[20..32]
        .trim()
        .parse()
        .map_err(|_| Error::InvalidTle("bad epochdays".into()))?;
    let ndot_raw: f64 = l1[33..43]
        .trim()
        .parse()
        .map_err(|_| Error::InvalidTle("bad ndot".into()))?;
    let nddot_str = format!("{}.{}", &l1[44..45], &l1[45..50]);
    let nddot_mantissa: f64 = nddot_str.trim().parse().unwrap_or(0.0);
    let nexp: i32 = l1[50..52].trim().parse().unwrap_or(0);
    let bstar_str = format!("{}.{}", &l1[53..54], &l1[54..59]);
    let bstar_mantissa: f64 = bstar_str.trim().parse().unwrap_or(0.0);
    let ibexp: i32 = l1[59..61].trim().parse().unwrap_or(0);

    // Line 2 fields
    let inclo_deg: f64 = l2[8..16]
        .trim()
        .parse()
        .map_err(|_| Error::InvalidTle("bad inclo".into()))?;
    let nodeo_deg: f64 = l2[17..25]
        .trim()
        .parse()
        .map_err(|_| Error::InvalidTle("bad nodeo".into()))?;
    let ecco_str = format!("0.{}", l2[26..33].replace(' ', "0"));
    let ecco: f64 = ecco_str
        .parse()
        .map_err(|_| Error::InvalidTle("bad ecco".into()))?;
    let argpo_deg: f64 = l2[34..42]
        .trim()
        .parse()
        .map_err(|_| Error::InvalidTle("bad argpo".into()))?;
    let mo_deg: f64 = l2[43..51]
        .trim()
        .parse()
        .map_err(|_| Error::InvalidTle("bad mo".into()))?;
    let no_kozai_revday: f64 = l2[52..63]
        .trim()
        .parse()
        .map_err(|_| Error::InvalidTle("bad no_kozai".into()))?;

    // Convert to SGP4 internal units (matches Python sgp4's twoline2rv).
    let no_kozai = no_kozai_revday / xpdotp;
    let nddot = nddot_mantissa * 10.0_f64.powi(nexp) / (xpdotp * 1440.0 * 1440.0);
    let bstar = bstar_mantissa * 10.0_f64.powi(ibexp);
    let ndot = ndot_raw / (xpdotp * 1440.0);
    let inclo = inclo_deg * deg2rad;
    let nodeo = nodeo_deg * deg2rad;
    let argpo = argpo_deg * deg2rad;
    let mo = mo_deg * deg2rad;

    // Epoch JD with Python sgp4 v2.20+ rounding.
    let year_full = if two_digit_year < 57 {
        two_digit_year + 2000
    } else {
        two_digit_year + 1900
    };
    let (mon, day, hr, minute, sec) = vallado::days2mdhms_SGP4(year_full, epochdays);
    let (jd, jdfrac_raw) = vallado::jday_SGP4(year_full, mon, day, hr, minute, sec);
    let jdfrac = (jdfrac_raw * 100_000_000.0).round() / 100_000_000.0;
    let epoch_sgp4 = jd + jdfrac - 2433281.5;

    let satnum = l1[2..7].trim();

    let mut satrec = vallado::ElsetRec::default();
    satrec.epochyr = two_digit_year;
    satrec.epochdays = epochdays;
    satrec.jdsatepoch = jd;
    satrec.jdsatepochF = jdfrac;

    vallado::sgp4init(
        vallado::GravConstType::Wgs72,
        'i',
        satnum,
        epoch_sgp4,
        bstar,
        ndot,
        nddot,
        ecco,
        argpo,
        inclo,
        mo,
        no_kozai,
        nodeo,
        &mut satrec,
    );

    // Note: a non-zero `satrec.error` is NOT a hard failure. Vallado's
    // sgp4init populates the satrec even on errors (1..=6) and the original
    // C++ / Python API surfaces these via per-propagate calls. We mirror that
    // here so callers can construct a Satellite from any well-formatted TLE
    // and discover propagation errors when they actually try to propagate.

    // Restore the split epoch — sgp4init/initl may have rewritten it.
    satrec.jdsatepoch = jd;
    satrec.jdsatepochF = jdfrac;

    Ok(satrec)
}