//! SGP4 near-earth orbit propagator.
//!
//! Implements the SGP4 algorithm from Vallado, Crawford, Hujsak, Kelso (2006),
//! "Revisiting Spacetrack Report Number 3."
//!
//! All computation is pure and immutable.  The Kepler equation is solved via
//! bounded recursion rather than loops.

use crate::coordinate::vec3::Vec3;
use crate::error::Error;
use crate::propagation::constants::{
    A30, EARTH_RADIUS_KM, J2, J3, K2, K4, KE, MINUTES_PER_DAY, QOMS2T, SO, TWO_PI, X2O3,
};
use crate::tle::Tle;
use crate::tle::newtypes::MinutesSinceEpoch;

// ---------------------------------------------------------------------------
// ECI output
// ---------------------------------------------------------------------------

/// Position and velocity in Earth-Centered Inertial (ECI) coordinates.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct EciOutput {
    /// Position vector [km].
    position_km: Vec3,
    /// Velocity vector [km/s].
    velocity_km_s: Vec3,
}

impl EciOutput {
    #[must_use]
    pub fn position_km(&self) -> Vec3 {
        self.position_km
    }

    #[must_use]
    pub fn velocity_km_s(&self) -> Vec3 {
        self.velocity_km_s
    }
}

// ---------------------------------------------------------------------------
// Helper structs for initialization sub-computations
// ---------------------------------------------------------------------------

/// Trigonometric quantities derived from the inclination.
struct InclinationTrig {
    cosio: f64,
    sinio: f64,
    x1mth2: f64,
    x3thm1: f64,
    x7thm1: f64,
}

/// Atmospheric density parameters adjusted for perigee height.
struct AtmosphericParams {
    ss_param: f64,
    qzms24: f64,
}

/// Secular rate coefficients.
struct SecularRates {
    mdot: f64,
    argpdot: f64,
    nodedot: f64,
    xhdot1: f64,
}

/// Higher-order drag terms (zero when `isimp` is true).
struct HigherOrderTerms {
    d2: f64,
    d3: f64,
    d4: f64,
    t3cof: f64,
    t4cof: f64,
    t5cof: f64,
}

/// C-coefficients and related drag parameters.
struct DragCoefficients {
    cc1: f64,
    cc4: f64,
    cc5: f64,
    omgcof: f64,
    sinmao: f64,
    mean_mot_drag_coef: f64,
    long_per_l_coef: f64,
    aycof: f64,
    delmo: f64,
}

// ---------------------------------------------------------------------------
// Helper functions for initialization
// ---------------------------------------------------------------------------

/// Compute trigonometric quantities from the inclination angle.
fn compute_inclination_trig(inclo: f64) -> InclinationTrig {
    let cosio = f64::cos(inclo);
    let sinio = f64::sin(inclo);
    let cosio_sq = cosio * cosio;
    InclinationTrig {
        cosio,
        sinio,
        x1mth2: 1.0 - cosio_sq,
        x3thm1: 3.0 * cosio_sq - 1.0,
        x7thm1: 7.0 * cosio_sq - 1.0,
    }
}

/// Adjust atmospheric density reference for low-perigee orbits.
///
/// Returns the (ss, qzms24) pair for the C-coefficient computation.
fn compute_atmospheric_params(perigee_km: f64) -> AtmosphericParams {
    if perigee_km < 156.0 {
        let ss_val = if perigee_km < 98.0 {
            20.0 / EARTH_RADIUS_KM + 1.0
        } else {
            perigee_km / EARTH_RADIUS_KM + 1.0
        };
        let qzms24_base = (120.0 - (ss_val - 1.0) * EARTH_RADIUS_KM) / EARTH_RADIUS_KM;
        AtmosphericParams {
            ss_param: ss_val,
            qzms24: qzms24_base.powi(4),
        }
    } else {
        AtmosphericParams {
            ss_param: SO,
            qzms24: QOMS2T,
        }
    }
}

/// Compute the secular rate coefficients for mean anomaly, argument of
/// perigee, and right ascension of the ascending node.
fn compute_secular_rates(
    no_unkozai: f64,
    po_sq: f64,
    etasq: f64,
    cosio: f64,
    x3thm1: f64,
    con42: f64,
) -> SecularRates {
    let cosio_sq = cosio * cosio;
    let rate1 = 1.5 * K2 * no_unkozai / po_sq;
    let rate2 = 0.5 * rate1 * K2 / f64::sqrt(po_sq);
    let rate3 = -0.46875 * K4 * no_unkozai / (po_sq * po_sq);

    let mdot = no_unkozai
        + 0.5 * rate1 * (1.0 - etasq) * x3thm1
        + 0.0625 * rate2 * (1.0 - etasq) * (13.0 - 78.0 * cosio_sq + 137.0 * cosio_sq * cosio_sq);

    let argpdot = -0.5 * rate1 * con42
        + 0.0625 * rate2 * (7.0 - 114.0 * cosio_sq + 395.0 * cosio_sq * cosio_sq)
        + rate3 * (3.0 - 36.0 * cosio_sq + 49.0 * cosio_sq * cosio_sq);

    let xhdot1 = -rate1 * cosio;
    let nodedot = xhdot1
        + (0.5 * rate2 * (4.0 - 19.0 * cosio_sq) + 2.0 * rate3 * (3.0 - 7.0 * cosio_sq)) * cosio;

    SecularRates {
        mdot,
        argpdot,
        nodedot,
        xhdot1,
    }
}

/// Compute higher-order drag terms for non-simplified orbits.
///
/// When `isimp` is true, all terms are zero.
fn compute_higher_order_terms(
    isimp: bool,
    ao: f64,
    tsi: f64,
    cc1: f64,
    ss_param: f64,
) -> HigherOrderTerms {
    if isimp {
        HigherOrderTerms {
            d2: 0.0,
            d3: 0.0,
            d4: 0.0,
            t3cof: 0.0,
            t4cof: 0.0,
            t5cof: 0.0,
        }
    } else {
        let cc1sq = cc1 * cc1;
        let d2 = 4.0 * ao * tsi * cc1sq;
        let scratch = d2 * tsi * cc1 / 3.0;
        let d3 = (17.0 * ao + ss_param) * scratch;
        let d4 = 0.5 * scratch * ao * tsi * (221.0 * ao + 31.0 * ss_param) * cc1;
        let t3cof = d2 + 2.0 * cc1sq;
        let t4cof = 0.25 * (3.0 * d3 + cc1 * (12.0 * d2 + 10.0 * cc1sq));
        let t5cof = 0.2
            * (3.0 * d4
                + 12.0 * cc1 * d3
                + 6.0 * d2 * d2
                + 15.0 * cc1sq * (2.0 * d2 + cc1sq));
        HigherOrderTerms {
            d2,
            d3,
            d4,
            t3cof,
            t4cof,
            t5cof,
        }
    }
}

/// Compute the C-coefficients and related drag parameters from the
/// intermediate orbital quantities.
#[allow(clippy::too_many_arguments)]
fn compute_drag_coefficients(
    trig: &InclinationTrig,
    ecco: f64,
    argpo: f64,
    mo: f64,
    bstar: f64,
    no_unkozai: f64,
    ao: f64,
    tsi: f64,
    eta: f64,
    con41: f64,
    coef: f64,
    coef1: f64,
    etasq: f64,
    ecc_eta: f64,
    psisq: f64,
) -> DragCoefficients {
    let cc2 = coef1
        * no_unkozai
        * (ao * (1.0 + 1.5 * etasq + ecc_eta * (4.0 + etasq))
            + 0.375 * K2 * tsi / psisq * trig.x3thm1
                * (8.0 + 3.0 * etasq * (8.0 + etasq)));
    let cc1 = bstar * cc2;
    let cc3 = if ecco > 1.0e-4 {
        -2.0 * coef * tsi * A30 / KE * no_unkozai * trig.sinio / ecco
    } else {
        0.0
    };
    let cc4 = 2.0
        * no_unkozai
        * coef1
        * ao
        * (1.0 - etasq)
        * (eta * (2.0 + 0.5 * etasq)
            + ecco * (0.5 + 2.0 * etasq)
            - K2 * tsi / (ao * psisq)
                * (-3.0 * con41 * (1.0 - 2.0 * ecc_eta + etasq * (1.5 - 0.5 * ecc_eta))
                    + 0.75
                        * trig.x1mth2
                        * (2.0 * etasq - ecc_eta * (1.0 + etasq))
                        * f64::cos(2.0 * argpo)));
    let cc5 = 2.0
        * coef1
        * ao
        * (1.0 - etasq)
        * (1.0 + 2.75 * (etasq + ecc_eta) + ecc_eta * etasq);

    let sinmao = f64::sin(mo);
    let omgcof = bstar * cc3 * f64::cos(argpo);
    let mean_mot_drag_coef = if ecco > 1.0e-4 {
        -X2O3 * coef * bstar / ecc_eta
    } else {
        0.0
    };

    let long_per_l_coef = if f64::abs(trig.cosio + 1.0) > 1.5e-12 {
        -0.25 * J3 / J2 * trig.sinio * (3.0 + 5.0 * trig.cosio) / (1.0 + trig.cosio)
    } else {
        -0.25 * J3 / J2 * trig.sinio * (3.0 + 5.0 * trig.cosio) / 1.5e-12
    };
    let aycof = -0.5 * J3 / J2 * trig.sinio;
    let delmo = (1.0 + eta * f64::cos(mo)).powi(3);

    DragCoefficients {
        cc1,
        cc4,
        cc5,
        omgcof,
        sinmao,
        mean_mot_drag_coef,
        long_per_l_coef,
        aycof,
        delmo,
    }
}

// ---------------------------------------------------------------------------
// Initialization record (computed once per TLE)
// ---------------------------------------------------------------------------

/// Pre-computed constants from a TLE, used for repeated propagation at
/// different time offsets.
#[derive(Debug, Clone)]
pub struct SgpInit {
    // Keplerian elements in radians / internal units
    epoch_jd: f64,
    bstar: f64,
    inclo: f64,
    nodeo: f64,
    ecco: f64,
    argpo: f64,
    mo: f64,
    no_unkozai: f64, // mean motion in rad/min (un-Kozai'd / Brouwer)

    // Initialization constants
    cosio: f64,
    sinio: f64,
    x1mth2: f64,
    x3thm1: f64,
    x7thm1: f64,
    long_per_l_coef: f64,
    aycof: f64,
    cc1: f64,
    cc4: f64,
    cc5: f64,
    d2: f64,
    d3: f64,
    d4: f64,
    eta: f64,
    argpdot: f64,
    omgcof: f64,
    sinmao: f64,
    t2cof: f64,
    t3cof: f64,
    t4cof: f64,
    t5cof: f64,
    mean_mot_drag_coef: f64,
    mdot: f64,
    nodedot: f64,
    nodecf: f64,
    delmo: f64,
    aodp: f64,
    isimp: bool, // true = simplified (low perigee, skip higher-order terms)
}

impl SgpInit {
    /// Initialize SGP4 constants from a TLE.
    ///
    /// # Errors
    ///
    /// Returns an error if the TLE's eccentricity or other elements
    /// produce singular conditions.
    pub fn from_tle(tle: &Tle) -> Result<Self, Error> {
        let ecco = tle.eccentricity().value();
        let inclo = tle.inclination().radians();
        let nodeo = tle.raan().radians();
        let argpo = tle.arg_perigee().radians();
        let mo = tle.mean_anomaly().radians();
        let bstar = tle.bstar().value();
        let no_kozai = tle.mean_motion().revs_per_day() * TWO_PI / MINUTES_PER_DAY;
        let epoch_jd = crate::time::epoch::JulianDate::from_epoch(
            tle.epoch_year(),
            tle.epoch_day(),
        )
        .value();

        let trig = compute_inclination_trig(inclo);
        let con42 = 1.0 - 5.0 * trig.cosio * trig.cosio;

        // Un-Kozai the mean motion
        let a1_val = f64::powf(KE / no_kozai, X2O3);
        let d1 = 1.5 * K2 * trig.x3thm1 / (a1_val * a1_val);
        let adel = a1_val
            * (1.0 - d1 * d1 - d1 * (1.0 / 3.0 + 134.0 * d1 * d1 / 81.0));
        let del_ = 1.5 * K2 * trig.x3thm1 / (adel * adel);
        let no_unkozai = no_kozai / (1.0 + del_);

        let ao = f64::powf(KE / no_unkozai, X2O3);
        let po = ao * (1.0 - ecco * ecco);
        let con41 = -con42 - 2.0 * trig.cosio * trig.cosio;
        let perigee_km = (ao * (1.0 - ecco) - 1.0) * EARTH_RADIUS_KM;
        let atm = compute_atmospheric_params(perigee_km);

        let tsi = 1.0 / (ao - atm.ss_param);
        let eta = ao * ecco * tsi;
        let etasq = eta * eta;
        let ecc_eta = ecco * eta;
        let psisq = f64::abs(1.0 - etasq);
        let coef = atm.qzms24 * tsi.powi(4);
        let coef1 = coef / f64::powf(psisq, 3.5);

        let dc = compute_drag_coefficients(
            &trig, ecco, argpo, mo, bstar, no_unkozai, ao, tsi,
            eta, con41, coef, coef1, etasq, ecc_eta, psisq,
        );
        let rates = compute_secular_rates(no_unkozai, po * po, etasq, trig.cosio, trig.x3thm1, con42);
        let nodecf = 3.5 * (1.0 - etasq) * rates.xhdot1 * dc.cc1;
        let isimp = perigee_km < 220.0;
        let ho = compute_higher_order_terms(isimp, ao, tsi, dc.cc1, atm.ss_param);

        Ok(Self {
            epoch_jd, bstar, inclo, nodeo, ecco, argpo, mo, no_unkozai,
            cosio: trig.cosio, sinio: trig.sinio, x1mth2: trig.x1mth2,
            x3thm1: trig.x3thm1, x7thm1: trig.x7thm1,
            long_per_l_coef: dc.long_per_l_coef, aycof: dc.aycof,
            cc1: dc.cc1, cc4: dc.cc4, cc5: dc.cc5,
            d2: ho.d2, d3: ho.d3, d4: ho.d4,
            eta, argpdot: rates.argpdot, omgcof: dc.omgcof,
            sinmao: dc.sinmao, t2cof: 1.5 * dc.cc1,
            t3cof: ho.t3cof, t4cof: ho.t4cof, t5cof: ho.t5cof,
            mean_mot_drag_coef: dc.mean_mot_drag_coef,
            mdot: rates.mdot, nodedot: rates.nodedot, nodecf,
            delmo: dc.delmo, aodp: ao, isimp,
        })
    }

    /// Epoch as Julian Date.
    #[must_use]
    pub fn epoch_jd(&self) -> f64 {
        self.epoch_jd
    }
}

// ---------------------------------------------------------------------------
// Kepler equation solver (bounded recursion, no mutation)
// ---------------------------------------------------------------------------

/// Solve Kepler's equation `E - e*sin(E) = M` via Newton-Raphson.
///
/// Returns the eccentric anomaly `E` in radians.
///
/// # Errors
///
/// Returns an error if convergence is not reached within `depth` iterations.
fn solve_kepler(mean_anom: f64, ecc: f64, depth: u8) -> Result<f64, Error> {
    solve_kepler_inner(mean_anom, ecc, mean_anom, depth)
}

fn solve_kepler_inner(mean_anom: f64, ecc: f64, current: f64, depth: u8) -> Result<f64, Error> {
    if depth == 0 {
        Err(Error::PropagationDiverged {
            context: "Kepler equation did not converge",
        })
    } else {
        let sin_e = f64::sin(current);
        let cos_e = f64::cos(current);
        let delta = (mean_anom - current + ecc * sin_e) / (1.0 - ecc * cos_e);
        let next = current + delta;
        if f64::abs(delta) < 1.0e-12 {
            Ok(next)
        } else {
            solve_kepler_inner(mean_anom, ecc, next, depth - 1)
        }
    }
}

// ---------------------------------------------------------------------------
// Propagation helper: secular drag and higher-order terms
// ---------------------------------------------------------------------------

/// Result of applying secular drag corrections.
struct SecularCorrection {
    xmp: f64,
    argpm: f64,
    a_factor: f64,
    e_factor: f64,
    l_factor: f64,
}

/// Apply secular effects of atmospheric drag including higher-order
/// terms when applicable.
fn apply_secular_drag(init: &SgpInit, t: f64, tsq: f64) -> SecularCorrection {
    let xmdf = init.mo + init.mdot * t;
    let argpdf = init.argpo + init.argpdot * t;
    let a_base = 1.0 - init.cc1 * t;
    let e_base = init.bstar * init.cc4 * t;
    let l_base = init.t2cof * tsq;

    if init.isimp {
        SecularCorrection {
            xmp: xmdf,
            argpm: argpdf,
            a_factor: a_base,
            e_factor: e_base,
            l_factor: l_base,
        }
    } else {
        let delomg = init.omgcof * t;
        let delmf = (1.0 + init.eta * f64::cos(xmdf)).powi(3) - init.delmo;
        let correction = delomg + init.mean_mot_drag_coef * delmf;
        let xmp_val = xmdf + correction;
        let tcube = tsq * t;
        let tfour = t * tcube;
        SecularCorrection {
            xmp: xmp_val,
            argpm: argpdf - correction,
            a_factor: a_base - init.d2 * tsq - init.d3 * tcube - init.d4 * tfour,
            e_factor: e_base + init.bstar * init.cc5 * (f64::sin(xmp_val) - init.sinmao),
            l_factor: l_base + init.t3cof * tcube + tfour * (init.t4cof + t * init.t5cof),
        }
    }
}

// ---------------------------------------------------------------------------
// Propagation helper: short-period periodics and ECI output
// ---------------------------------------------------------------------------

/// Compute ECI position and velocity from orbital elements after
/// short-period periodic corrections.
///
/// # Errors
///
/// Returns an error if the Kepler equation fails to converge or the
/// orbit produces a negative semi-latus rectum.
#[allow(clippy::similar_names)]
fn compute_eci_output(
    init: &SgpInit,
    a: f64,
    ecc: f64,
    xmp: f64,
    argpm: f64,
    xnode: f64,
    n: f64,
) -> Result<EciOutput, Error> {
    // Long-period periodics
    let axn = ecc * f64::cos(argpm);
    let xl = xmp + argpm + xnode + init.no_unkozai * n;
    let ayn = ecc * f64::sin(argpm) + init.aycof * (1.0 + init.long_per_l_coef * (xl - xnode));
    let xlt = xl + init.long_per_l_coef * axn;

    // Solve Kepler's equation
    let capu = (xlt - xnode) % TWO_PI;
    let eo1 = solve_kepler(capu, (axn * axn + ayn * ayn).sqrt(), 30)?;

    let sin_eo1 = f64::sin(eo1);
    let cos_eo1 = f64::cos(eo1);

    let ecose = axn * cos_eo1 + ayn * sin_eo1;
    let esine = axn * sin_eo1 - ayn * cos_eo1;

    let el2 = axn * axn + ayn * ayn;
    let pl = a * (1.0 - el2);
    if pl < 0.0 {
        Err(Error::PropagationDiverged {
            context: "negative semi-latus rectum",
        })
    } else {
        let freq = KE / f64::powf(a, 1.5);
        let rl = a * (1.0 - ecose);
        let rdot_raw = (a * esine) / rl * (KE / a.sqrt());
        let rvdot_raw = pl.sqrt() * KE / rl;

        let betal = (1.0 - el2).sqrt();
        let sin_u = a / rl * (sin_eo1 - ayn - axn * esine / (1.0 + betal));
        let cos_u = a / rl * (cos_eo1 - axn + ayn * esine / (1.0 + betal));
        let su = f64::atan2(sin_u, cos_u);

        let sin_2u = 2.0 * cos_u * sin_u;
        let cos_2u = 2.0 * cos_u * cos_u - 1.0;

        let inv_pl = 1.0 / pl;
        let inv_pl_sq = inv_pl * inv_pl;

        // Short-period periodics
        let rk = rl * (1.0 - 1.5 * K2 * inv_pl_sq * init.x3thm1 * (0.5 * cos_2u))
            + 0.5 * K2 * inv_pl * init.x1mth2 * cos_2u;
        let uk = su - 0.25 * K2 * inv_pl_sq * init.x7thm1 * sin_2u;
        let xnodek = xnode + 1.5 * K2 * inv_pl_sq * init.cosio * sin_2u;
        let xinck = init.inclo + 1.5 * K2 * inv_pl_sq * init.sinio * init.cosio * cos_2u;
        let rdot_corr = rdot_raw - freq * K2 * inv_pl * init.x1mth2 * sin_2u;
        let rfdot_corr = rvdot_raw + freq * K2 * inv_pl * (init.x1mth2 * cos_2u + 1.5 * init.x3thm1);

        // Orientation vectors
        let sin_uk = f64::sin(uk);
        let cos_uk = f64::cos(uk);
        let sin_ik = f64::sin(xinck);
        let cos_ik = f64::cos(xinck);
        let sin_nk = f64::sin(xnodek);
        let cos_nk = f64::cos(xnodek);

        let xmx = -sin_nk * cos_ik;
        let xmy = cos_nk * cos_ik;

        let ux = xmx * sin_uk + cos_nk * cos_uk;
        let uy = xmy * sin_uk + sin_nk * cos_uk;
        let uz = sin_ik * sin_uk;

        let vx = xmx * cos_uk - cos_nk * sin_uk;
        let vy = xmy * cos_uk - sin_nk * sin_uk;
        let vz = sin_ik * cos_uk;

        // Position (earth radii) and velocity (earth radii / min)
        let r_er = Vec3::new(rk * ux, rk * uy, rk * uz);
        let v_er_min = Vec3::new(
            rdot_corr * ux + rfdot_corr * vx,
            rdot_corr * uy + rfdot_corr * vy,
            rdot_corr * uz + rfdot_corr * vz,
        );

        Ok(EciOutput {
            position_km: r_er.scale(EARTH_RADIUS_KM),
            velocity_km_s: v_er_min.scale(EARTH_RADIUS_KM / 60.0),
        })
    }
}

// ---------------------------------------------------------------------------
// Propagation
// ---------------------------------------------------------------------------

/// Propagate a near-earth satellite to `tsince` minutes after epoch.
///
/// # Errors
///
/// Returns an error if the Kepler equation fails to converge or
/// the orbit decays (eccentricity diverges).
pub fn propagate(init: &SgpInit, tsince: MinutesSinceEpoch) -> Result<EciOutput, Error> {
    let t = tsince.value();
    let tsq = t * t;

    let nodedf = init.nodeo + init.nodedot * t;
    let xnode = nodedf + init.nodecf * tsq;

    let sec = apply_secular_drag(init, t, tsq);

    let a = init.aodp * sec.a_factor * sec.a_factor;
    let ecc_raw = init.ecco - sec.e_factor;

    // Clamp eccentricity
    let ecc = if ecc_raw < 1.0e-6 { 1.0e-6 } else { ecc_raw };
    if ecc >= 1.0 {
        Err(Error::PropagationDiverged {
            context: "eccentricity >= 1.0 during propagation",
        })
    } else {
        compute_eci_output(init, a, ecc, sec.xmp, sec.argpm, xnode, sec.l_factor)
    }
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;
    use crate::tle::parse::parse_tle;

    // ISS TLE
    const ISS_L0: &str = "ISS (ZARYA)";
    const ISS_L1: &str =
        "1 25544U 98067A   24045.54339622  .00016717  00000-0  10270-3 0  9118";
    const ISS_L2: &str =
        "2 25544  51.6400 208.2546 0006703  44.2624  67.4637 15.50056823 23455";

    #[test]
    fn sgp4_init_succeeds_for_iss() -> Result<(), Error> {
        let tle = parse_tle(ISS_L0, ISS_L1, ISS_L2).run()?;
        SgpInit::from_tle(&tle).map(|_| ())
    }

    #[test]
    fn sgp4_propagate_at_epoch_returns_reasonable_position() -> Result<(), Error> {
        let tle = parse_tle(ISS_L0, ISS_L1, ISS_L2).run()?;
        let init = SgpInit::from_tle(&tle)?;
        let result = propagate(&init, MinutesSinceEpoch::new(0.0))?;

        // ISS orbits at ~408 km altitude => radius ~6786 km
        let r = result.position_km().magnitude();
        assert!(
            r > 6400.0 && r < 7200.0,
            "ISS radius {r} km out of expected range"
        );
        Ok(())
    }

    #[test]
    fn sgp4_propagate_at_60_min_differs_from_epoch() -> Result<(), Error> {
        let tle = parse_tle(ISS_L0, ISS_L1, ISS_L2).run()?;
        let init = SgpInit::from_tle(&tle)?;
        let r0 = propagate(&init, MinutesSinceEpoch::new(0.0))?;
        let r60 = propagate(&init, MinutesSinceEpoch::new(60.0))?;

        // Positions should differ (ISS completes ~2/3 of an orbit in 60 min)
        let dist = (r0.position_km() - r60.position_km()).magnitude();
        assert!(
            dist > 1000.0,
            "expected significant position change after 60 min, got {dist} km"
        );
        Ok(())
    }

    #[test]
    fn sgp4_velocity_is_reasonable_for_iss() -> Result<(), Error> {
        let tle = parse_tle(ISS_L0, ISS_L1, ISS_L2).run()?;
        let init = SgpInit::from_tle(&tle)?;
        let result = propagate(&init, MinutesSinceEpoch::new(0.0))?;

        // ISS velocity is ~7.66 km/s
        let v = result.velocity_km_s().magnitude();
        assert!(
            v > 7.0 && v < 8.5,
            "ISS velocity {v} km/s out of expected range"
        );
        Ok(())
    }
}
