astrodynamics 0.6.3

//! 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);

/// Vallado SGP4 operation mode. Controls which initialization branch
/// `sgp4init` follows. The two modes produce subtly different results;
/// the divergence grows with propagation time and can reach hundreds of
/// millimeters over a few orbits at LEO.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub enum OpsMode {
    /// Improved formulation (Vallado `'i'`). Default for new work and
    /// matches the default of Python's `sgp4` package. Recommended unless
    /// you specifically need AFSPC operational parity.
    #[default]
    Improved,
    /// AFSPC operational compatibility mode (Vallado `'a'`). Use this
    /// when reproducing outputs from operational AFSPC systems or matching
    /// reference values from older crates / catalogs that ran in AFSPC mode.
    Afspc,
}

impl OpsMode {
    fn as_char(self) -> char {
        match self {
            OpsMode::Improved => 'i',
            OpsMode::Afspc => 'a',
        }
    }
}

/// Pre-parsed Vallado SGP4 element set.
///
/// Use this when the TLE has already been parsed externally (e.g. from an
/// OMM message, JSON catalog, or another system) and you want to feed the
/// element values directly into the SGP4 initializer instead of going through
/// the TLE string parser.
///
/// Field units match the Vallado SGP4 reference inputs:
/// angles in **degrees**, mean motion in **revolutions per day**, drag term
/// in the dimensionless TLE convention.
#[derive(Debug, Clone, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct ElementSet {
    /// Two-digit epoch year as it appears in the TLE. Years 0..=56 are
    /// interpreted as 2000..=2056; 57..=99 as 1957..=1999.
    pub epoch_year_two_digit: i32,
    /// Day of year, fractional. e.g. `184.80969102` for 2018-07-04 19:25 UTC.
    pub epoch_days: f64,
    /// SGP4 drag term (Vallado B\*). Dimensionless TLE convention.
    pub bstar: f64,
    /// First derivative of mean motion in rev/day². TLE "ndot".
    pub mean_motion_dot: f64,
    /// Second derivative of mean motion in rev/day³. TLE "nddot".
    pub mean_motion_double_dot: f64,
    /// Eccentricity, dimensionless, in [0, 1).
    pub eccentricity: f64,
    /// Argument of perigee, degrees.
    pub argument_of_perigee_deg: f64,
    /// Inclination, degrees.
    pub inclination_deg: f64,
    /// Mean anomaly, degrees.
    pub mean_anomaly_deg: f64,
    /// Mean motion, revolutions per day.
    pub mean_motion_rev_per_day: f64,
    /// Right ascension of ascending node (RAAN), degrees.
    pub right_ascension_deg: f64,
    /// Catalog (NORAD) number for this object. Used only for diagnostic
    /// reporting inside SGP4 — propagation results do not depend on it.
    /// Pass `0` if unknown.
    pub catalog_number: u32,
}

// ── 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 using the default `Improved` opsmode.
    ///
    /// 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> {
        Self::from_tle_with_opsmode(line1, line2, OpsMode::Improved)
    }

    /// Parse a two-line element set with an explicit `OpsMode`.
    ///
    /// Use `OpsMode::Afspc` to reproduce results from operational AFSPC
    /// systems or older crates that ran in AFSPC compatibility mode.
    pub fn from_tle_with_opsmode(
        line1: &str,
        line2: &str,
        opsmode: OpsMode,
    ) -> 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, opsmode)?;

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

    /// Construct a `Satellite` from pre-parsed Vallado SGP4 elements using
    /// the default `Improved` opsmode.
    ///
    /// Useful when TLE data has already been parsed externally (OMM, JSON
    /// catalog, another system) and you only need the element values to flow
    /// into SGP4 initialization. Equivalent to `from_tle` for propagation
    /// behavior, but bypasses the TLE string parser.
    ///
    /// Note: a `Satellite` constructed this way has empty `line1()` and
    /// `line2()` accessors since there is no source TLE to return.
    pub fn from_elements(elements: &ElementSet) -> Result<Self, Error> {
        Self::from_elements_with_opsmode(elements, OpsMode::Improved)
    }

    /// Construct a `Satellite` from pre-parsed elements with an explicit
    /// `OpsMode`. See `from_tle_with_opsmode` for the rationale.
    pub fn from_elements_with_opsmode(
        elements: &ElementSet,
        opsmode: OpsMode,
    ) -> Result<Self, Error> {
        let satrec = init_satrec_from_elements(elements, opsmode)?;
        Ok(Satellite {
            line1: String::new(),
            line2: String::new(),
            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. Returns an empty string when this `Satellite` was
    /// constructed via `from_elements` (no source TLE).
    pub fn line1(&self) -> &str {
        &self.line1
    }

    /// Raw TLE line 2. Returns an empty string when this `Satellite` was
    /// constructed via `from_elements` (no source TLE).
    pub fn line2(&self) -> &str {
        &self.line2
    }

    /// Cached TLE epoch as a split Julian date `(jdsatepoch, jdsatepochF)`.
    ///
    /// Useful for computing time offsets between two TLEs without losing
    /// precision through floating-point round-trips.
    pub fn epoch_jd(&self) -> JulianDate {
        JulianDate(self.satrec.jdsatepoch, self.satrec.jdsatepochF)
    }
}

/// One-shot SGP4 propagation from pre-parsed elements using the default
/// `Improved` opsmode.
///
/// Equivalent to `Satellite::from_elements(&e)?.propagate(t)` but without
/// allocating a cached `Satellite`. Suitable for one-call use cases where
/// the satellite record is not reused (e.g. NIF entry points that get
/// elements + a single time per call).
pub fn propagate_elements(
    elements: &ElementSet,
    t: MinutesSinceEpoch,
) -> Result<Prediction, Error> {
    propagate_elements_with_opsmode(elements, t, OpsMode::Improved)
}

/// One-shot SGP4 propagation with an explicit `OpsMode`.
pub fn propagate_elements_with_opsmode(
    elements: &ElementSet,
    t: MinutesSinceEpoch,
    opsmode: OpsMode,
) -> Result<Prediction, Error> {
    let mut satrec = init_satrec_from_elements(elements, opsmode)?;
    let mut r = [0.0_f64; 3];
    let mut v = [0.0_f64; 3];
    let ok = vallado::sgp4(&mut satrec, t.0, &mut r, &mut v);
    if !ok || satrec.error != 0 {
        return Err(Error::Sgp4 {
            code: satrec.error,
        });
    }
    Ok(Prediction {
        position: r,
        velocity: v,
    })
}

#[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 ─────────────────────────────────────────────────

/// Run `sgp4init` from a pre-parsed element set, returning the initialized
/// satellite record. Performs the same angle/units conversion as
/// `vallado::twoline2rv_propagate` so a `Satellite` constructed from
/// elements is equivalent to one constructed from the matching TLE.
fn init_satrec_from_elements(
    elements: &ElementSet,
    opsmode: OpsMode,
) -> Result<vallado::ElsetRec, Error> {
    let deg2rad = std::f64::consts::PI / 180.0;
    let xpdotp = 1440.0 / (2.0 * std::f64::consts::PI);

    let inclo = elements.inclination_deg * deg2rad;
    let nodeo = elements.right_ascension_deg * deg2rad;
    let argpo = elements.argument_of_perigee_deg * deg2rad;
    let mo = elements.mean_anomaly_deg * deg2rad;
    let no_kozai = elements.mean_motion_rev_per_day / xpdotp;
    // ndot rev/day² → rad/min², nddot rev/day³ → rad/min³.
    // Matches the conversion in `vallado::twoline2rv_propagate`.
    let ndot = elements.mean_motion_dot / (xpdotp * 1440.0);
    let nddot = elements.mean_motion_double_dot / (xpdotp * 1440.0 * 1440.0);

    // Compute jdsatepoch / jdsatepochF the same way twoline2rv_propagate does
    // (so a Satellite created via from_elements has the same epoch precision
    // as one created via from_tle for the equivalent TLE).
    let year_full = if elements.epoch_year_two_digit < 57 {
        elements.epoch_year_two_digit + 2000
    } else {
        elements.epoch_year_two_digit + 1900
    };
    let (mon, day, hr, minute, sec) =
        vallado::days2mdhms_SGP4(year_full, elements.epoch_days);
    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_str = format!("{:>5}", elements.catalog_number);

    let mut satrec = vallado::ElsetRec {
        epochyr: elements.epoch_year_two_digit,
        epochdays: elements.epoch_days,
        jdsatepoch: jd,
        jdsatepochF: jdfrac,
        ..vallado::ElsetRec::default()
    };

    vallado::sgp4init(
        vallado::GravConstType::Wgs72,
        opsmode.as_char(),
        &satnum_str,
        epoch_sgp4,
        elements.bstar,
        ndot,
        nddot,
        elements.eccentricity,
        argpo,
        inclo,
        mo,
        no_kozai,
        nodeo,
        &mut satrec,
    );

    // sgp4init may have rewritten jdsatepoch via initl — restore the split.
    satrec.jdsatepoch = jd;
    satrec.jdsatepochF = jdfrac;

    Ok(satrec)
}

/// 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,
    opsmode: OpsMode,
) -> 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 {
        epochyr: two_digit_year,
        epochdays,
        jdsatepoch: jd,
        jdsatepochF: jdfrac,
        ..vallado::ElsetRec::default()
    };

    vallado::sgp4init(
        vallado::GravConstType::Wgs72,
        opsmode.as_char(),
        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)
}