//! SDP4 deep-space orbital propagator extension.
//!
//! Extends SGP4 for satellites with orbital periods >= 225 minutes
//! (MEO, GEO, HEO). Handles third-body perturbations from the Sun
//! and Moon, and resonance effects for 12-hour and 24-hour orbits.
//!
//! Based on Vallado et al., "Revisiting Spacetrack Report #3" (2006).

use serde::{Deserialize, Serialize};

use crate::constants::*;

/// Deep-space resonance type.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
pub(crate) enum ResonanceType {
    /// No resonance.
    None,
    /// Synchronous (24-hour) resonance.
    Synchronous,
    /// Half-day (12-hour) resonance.
    HalfDay,
}

/// Deep-space state for the SDP4 propagator extension.
///
/// Stores lunar-solar perturbation terms, resonance coefficients,
/// and integration state for deep-space secular effects.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub(crate) struct DeepSpaceState {
    // ── Lunar-solar perturbation constants ────────────────────
    pub(crate) se2: f64,
    pub(crate) se3: f64,
    pub(crate) si2: f64,
    pub(crate) si3: f64,
    pub(crate) sl2: f64,
    pub(crate) sl3: f64,
    pub(crate) sl4: f64,
    pub(crate) sgh2: f64,
    pub(crate) sgh3: f64,
    pub(crate) sgh4: f64,
    pub(crate) sh2: f64,
    pub(crate) sh3: f64,

    pub(crate) ee2: f64,
    pub(crate) e3: f64,
    pub(crate) xi2: f64,
    pub(crate) xi3: f64,
    pub(crate) xl2: f64,
    pub(crate) xl3: f64,
    pub(crate) xl4: f64,
    pub(crate) xgh2: f64,
    pub(crate) xgh3: f64,
    pub(crate) xgh4: f64,
    pub(crate) xh2: f64,
    pub(crate) xh3: f64,

    // ── Solar and lunar mean elements ────────────────────────
    pub(crate) zmol: f64,
    pub(crate) zmos: f64,

    // ── Resonance state ──────────────────────────────────────
    pub(crate) resonance: ResonanceType,

    // Synchronous resonance terms
    pub(crate) del1: f64,
    pub(crate) del2: f64,
    pub(crate) del3: f64,
    pub(crate) xfact: f64,
    pub(crate) xlamo: f64,

    // Half-day resonance terms
    pub(crate) d2201: f64,
    pub(crate) d2211: f64,
    pub(crate) d3210: f64,
    pub(crate) d3222: f64,
    pub(crate) d4410: f64,
    pub(crate) d4422: f64,
    pub(crate) d5220: f64,
    pub(crate) d5232: f64,
    pub(crate) d5421: f64,
    pub(crate) d5433: f64,

    // Common integration terms
    pub(crate) dedt: f64,
    pub(crate) didt: f64,
    pub(crate) dmdt: f64,
    pub(crate) dnodt: f64,
    pub(crate) domdt: f64,

    // Saved epoch elements for integration
    pub(crate) atime: f64,
    pub(crate) xli: f64,
    pub(crate) xni: f64,

    // Secular rates from deep space
    pub(crate) gsto: f64,
}

impl DeepSpaceState {
    /// Initialize deep-space terms from the orbital elements.
    ///
    /// This computes lunar-solar perturbation coefficients and
    /// identifies any orbital resonance conditions.
    #[allow(clippy::too_many_arguments)] // SGP4/SDP4 deep-space state init
    pub(crate) fn initialize(
        epoch_jd: f64,
        ecco: f64,
        inclo: f64,
        nodeo: f64,
        argpo: f64,
        mo: f64,
        no_kozai: f64,
        _bstar: f64,
    ) -> Self {
        let cosio = inclo.cos();
        let sinio = inclo.sin();

        // Julian centuries from J2000
        let tc = (epoch_jd - 2451545.0) / 36525.0;

        // ── Solar terms ──────────────────────────────────────
        let _zns = 1.19459e-5;
        let zes = 0.01675;
        let _znl = 1.5835218e-4;
        let zel = 0.05490;
        let c1ss = 2.9864797e-6;
        let _c1l = 4.7968065e-7;

        // Mean anomaly of the sun
        let zmos = (6.2565837 + 0.017201977 * (epoch_jd - 2451545.0))
            .rem_euclid(TWOPI);
        // Mean anomaly of the moon
        let zmol = (2.3555557 + 0.228737 * (epoch_jd - 2451545.0) - 0.0003029 * tc * tc)
            .rem_euclid(TWOPI);

        let theta2 = cosio * cosio;

        // Eccentricity functions
        let eosq = ecco * ecco;

        // ── Compute solar perturbation coefficients ──────────
        let (se2, se3, si2, si3, sl2, sl3, sl4, sgh2, sgh3, sgh4, sh2, sh3) =
            Self::solar_lunar_coefficients(
                inclo, nodeo, argpo, ecco, true, zmos, zes,
            );

        // ── Compute lunar perturbation coefficients ──────────
        let (ee2, e3, xi2, xi3, xl2, xl3, xl4, xgh2, xgh3, xgh4, xh2, xh3) =
            Self::solar_lunar_coefficients(
                inclo, nodeo, argpo, ecco, false, zmol, zel,
            );

        // ── Resonance detection ──────────────────────────────
        let period = TWOPI / no_kozai;
        let resonance;
        let mut del1 = 0.0;
        let mut del2 = 0.0;
        let mut del3 = 0.0;
        let mut xfact = 0.0;
        let mut xlamo = 0.0;
        let mut d2201 = 0.0;
        let mut d2211 = 0.0;
        let mut d3210 = 0.0;
        let mut d3222 = 0.0;
        let mut d4410 = 0.0;
        let mut d4422 = 0.0;
        let mut d5220 = 0.0;
        let mut d5232 = 0.0;
        let mut d5421 = 0.0;
        let mut d5433 = 0.0;

        // Compute GMST at epoch
        let gsto = gmst_from_jd(epoch_jd);

        // Check for synchronous resonance (period ~ 1436 min)
        if period > 1200.0 && period < 1800.0 {
            resonance = ResonanceType::Synchronous;

            // Synchronous resonance coefficients
            let g200 = 1.0 + eosq * (-2.5 + 0.8125 * eosq);
            let g310 = 1.0 + 2.0 * eosq;
            let g300 = 1.0 + eosq * (-6.0 + 6.60937 * eosq);

            let f220 = 0.75 * (1.0 + cosio) * (1.0 + cosio);
            let f311 = 0.9375 * sinio * sinio * (1.0 + 3.0 * cosio) - 0.75 * (1.0 + cosio);
            let f330 = 1.0 + cosio;
            let f330 = 1.875 * f330 * f330 * f330;

            del1 = 3.0 * no_kozai * no_kozai * (1.0 / (ao_from_no(no_kozai) * ao_from_no(no_kozai)));
            del1 = del1 * f220 * g200 * CK2;
            del2 = 2.0 * del1 * f311 * g310 * CK2 / (ao_from_no(no_kozai));
            del3 = 3.0 * del1 * f330 * g300 * CK2 / (ao_from_no(no_kozai) * ao_from_no(no_kozai));

            xlamo = (mo + nodeo + argpo - gsto).rem_euclid(TWOPI);
            xfact = argpo + mo + nodeo + no_kozai - gsto - TWOPI / MINUTES_PER_DAY;
        } else if period > 680.0 && period < 760.0 {
            // Half-day resonance (period ~ 720 min)
            resonance = ResonanceType::HalfDay;

            // Inclination functions for half-day resonance
            let g201 = -0.306 - (ecco - 0.64) * 0.440;

            let (g211, g310, g322, g410, g422, g520, g521, g532, g533) = if ecco <= 0.65 {
                let g211 = 3.616 - 13.247 * ecco + 16.290 * eosq;
                let g310 = -19.302 + 117.390 * ecco - 228.419 * eosq + 156.591 * eosq * ecco;
                let g322 = -18.9068 + 109.7927 * ecco - 214.6334 * eosq + 146.5816 * eosq * ecco;
                let g410 = -41.122 + 242.694 * ecco - 471.094 * eosq + 313.953 * eosq * ecco;
                let g422 = -146.407 + 841.880 * ecco - 1629.014 * eosq + 1083.435 * eosq * ecco;
                let g520 = -532.114 + 3017.977 * ecco - 5740.032 * eosq + 3708.276 * eosq * ecco;
                let g521 = -822.71072 + 4568.6173 * ecco - 8491.4146 * eosq + 5337.524 * eosq * ecco;
                let g532 = -853.666 + 4690.25 * ecco - 8624.77 * eosq + 5341.4 * eosq * ecco;
                let g533 = -919.2277 + 4988.61 * ecco - 9064.77 * eosq + 5542.21 * eosq * ecco;
                (g211, g310, g322, g410, g422, g520, g521, g532, g533)
            } else {
                let g211 = -72.099 + 331.819 * ecco - 508.738 * eosq + 266.724 * eosq * ecco;
                let g310 = -346.844 + 1582.851 * ecco - 2415.925 * eosq + 1246.113 * eosq * ecco;
                let g322 = -342.585 + 1554.908 * ecco - 2366.899 * eosq + 1215.972 * eosq * ecco;
                let g410 = -1052.797 + 4758.686 * ecco - 7193.992 * eosq + 3651.957 * eosq * ecco;
                let g422 = -3581.69 + 16178.11 * ecco - 24462.77 * eosq + 12422.52 * eosq * ecco;
                let g520 = if ecco > 0.715 {
                    -5149.66 + 29936.92 * ecco - 54087.36 * eosq + 31324.56 * eosq * ecco
                } else {
                    1464.74 - 4664.75 * ecco + 3763.64 * eosq
                };
                let g521 = -3.908 + 155.237 * ecco - 498.867 * eosq;
                let g532 = -37.995 + 379.826 * ecco - 1241.41 * eosq + 1263.28 * eosq * ecco;
                let g533 = -113.479 + 691.178 * ecco - 1769.05 * eosq + 1571.84 * eosq * ecco;
                (g211, g310, g322, g410, g422, g520, g521, g532, g533)
            };

            let sini2o = sinio * sinio;
            let f220 = 0.75 * (1.0 + 2.0 * cosio + theta2);
            let f221 = 1.5 * sini2o;
            let f321 = 1.875 * sinio * (1.0 - 2.0 * cosio - 3.0 * theta2);
            let f322 = -1.875 * sinio * (1.0 + 2.0 * cosio - 3.0 * theta2);
            let f441 = 35.0 * sini2o * f220;
            let f442 = 39.375 * sini2o * sini2o;
            let f522 = 9.84375
                * sinio
                * (sini2o * (1.0 - 2.0 * cosio - 5.0 * theta2)
                    + 0.33333333 * (-2.0 + 4.0 * cosio + 6.0 * theta2));
            let f523 = sinio
                * (4.92187512 * sini2o * (-2.0 - 4.0 * cosio + 10.0 * theta2)
                    + 6.56250012 * (1.0 + 2.0 * cosio - 3.0 * theta2));
            let f542 = 29.53125
                * sinio
                * (2.0 - 8.0 * cosio + theta2 * (-12.0 + 8.0 * cosio + 10.0 * theta2));
            let f543 = 29.53125
                * sinio
                * (-2.0 - 8.0 * cosio + theta2 * (12.0 + 8.0 * cosio - 10.0 * theta2));

            let _ao = ao_from_no(no_kozai);

            d2201 = c1ss * f220 * g201;
            d2211 = c1ss * f221 * g211;
            d3210 = c1ss * f321 * g310;
            d3222 = c1ss * f322 * g322;
            d4410 = c1ss * f441 * g410;
            d4422 = c1ss * f442 * g422;
            d5220 = c1ss * f522 * g520;
            d5232 = c1ss * f523 * g532;
            d5421 = c1ss * f542 * g521;
            d5433 = c1ss * f543 * g533;

            xlamo = (mo + nodeo + nodeo - gsto - gsto).rem_euclid(TWOPI);
            xfact = mo + argpo + nodeo + nodeo - TWOPI / MINUTES_PER_DAY * 2.0 - gsto - gsto;
        } else {
            resonance = ResonanceType::None;
        }

        // Compute secular rates
        let dedt = se2 * zmos.sin() + sl2 * (2.0 * zmos).sin();
        let didt = si2 * zmos.sin();
        let dmdt = sl2 * zmos.cos();
        let dnodt = sh2 * zmos.sin();
        let domdt = sgh2 * zmos.sin();

        DeepSpaceState {
            se2, se3, si2, si3, sl2, sl3, sl4,
            sgh2, sgh3, sgh4, sh2, sh3,
            ee2, e3, xi2, xi3, xl2, xl3, xl4,
            xgh2, xgh3, xgh4, xh2, xh3,
            zmol, zmos,
            resonance, del1, del2, del3,
            xfact, xlamo,
            d2201, d2211, d3210, d3222,
            d4410, d4422, d5220, d5232, d5421, d5433,
            dedt, didt, dmdt, dnodt, domdt,
            atime: 0.0,
            xli: if resonance != ResonanceType::None { xlamo } else { 0.0 },
            xni: no_kozai,
            gsto,
        }
    }

    /// Compute solar or lunar perturbation coefficients.
    ///
    /// `is_solar`: true for solar, false for lunar.
    fn solar_lunar_coefficients(
        inclo: f64,
        nodeo: f64,
        argpo: f64,
        ecco: f64,
        is_solar: bool,
        _zm: f64,
        _ze: f64,
    ) -> (f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64) {
        let cosio = inclo.cos();
        let sinio = inclo.sin();

        // Third-body inclination and node
        let (sinx, _cosx) = if is_solar {
            (0.39785416_f64.sin(), 0.39785416_f64.cos())
        } else {
            (0.08979719_f64.sin(), 0.08979719_f64.cos())
        };

        let cc = if is_solar { 2.9864797e-6 } else { 4.7968065e-7 };

        let a9 = if is_solar { 1.19459e-5 } else { 1.5835218e-4 };
        let a10 = a9;

        // Simplified perturbation coefficients
        let sinzf = argpo.sin();
        let coszf = argpo.cos();
        let sinzy = nodeo.sin();
        let coszy = nodeo.cos();

        let ses = cc * sinx * (a9 + a10 * (sinio * sinzy * ecco));
        let sis = cc * (a9 * cosio * sinzy);
        let sls = cc * (a9 * (sinio * coszf - cosio * sinzf * coszy)
            + a10 * sinio * ecco * coszy);

        let sghs = cc * cosio * (a9 * sinzf + a10 * ecco * coszf);
        let shs = cc * sinio * (a9 * coszy - a10 * ecco * sinzy);

        let se2 = 2.0 * ses;
        let se3 = 0.0;
        let si2 = 2.0 * sis;
        let si3 = 0.0;
        let sl2 = 2.0 * sls;
        let sl3 = 0.0;
        let sl4 = 0.0;
        let sgh2 = 2.0 * sghs;
        let sgh3 = 0.0;
        let sgh4 = 0.0;
        let sh2 = 2.0 * shs;
        let sh3 = 0.0;

        (se2, se3, si2, si3, sl2, sl3, sl4, sgh2, sgh3, sgh4, sh2, sh3)
    }

    /// Apply deep-space secular perturbations.
    ///
    /// Returns updated (mp, node, argp, eccentricity, inclination, mean_motion).
    #[allow(clippy::too_many_arguments)] // SGP4/SDP4 deep-space secular update
    pub(crate) fn secular(
        &self,
        tsince: f64,
        mp: f64,
        node: f64,
        argp: f64,
        ecco: f64,
        inclo: f64,
        no_kozai: f64,
    ) -> (f64, f64, f64, f64, f64, f64) {
        // Solar and lunar mean anomalies at current time
        let zms = (self.zmos + 1.19459e-5 * tsince).rem_euclid(TWOPI);
        let zml = (self.zmol + 1.5835218e-4 * tsince).rem_euclid(TWOPI);

        // Solar perturbations
        let mut em = ecco
            + self.se2 * zms.sin()
            + self.se3 * (2.0 * zms).sin();
        let mut xinc = inclo
            + self.si2 * zms.sin()
            + self.si3 * (2.0 * zms).sin();
        let mut xm = mp
            + self.sl2 * zms.sin()
            + self.sl3 * (2.0 * zms).sin()
            + self.sl4 * (3.0 * zms).sin();
        let mut xn = node
            + self.sh2 * zms.sin()
            + self.sh3 * (2.0 * zms).sin();
        let mut xg = argp
            + self.sgh2 * zms.sin()
            + self.sgh3 * (2.0 * zms).sin()
            + self.sgh4 * (3.0 * zms).sin();

        // Lunar perturbations
        em += self.ee2 * zml.sin() + self.e3 * (2.0 * zml).sin();
        xinc += self.xi2 * zml.sin() + self.xi3 * (2.0 * zml).sin();
        xm += self.xl2 * zml.sin()
            + self.xl3 * (2.0 * zml).sin()
            + self.xl4 * (3.0 * zml).sin();
        xn += self.xh2 * zml.sin() + self.xh3 * (2.0 * zml).sin();
        xg += self.xgh2 * zml.sin()
            + self.xgh3 * (2.0 * zml).sin()
            + self.xgh4 * (3.0 * zml).sin();

        // Resonance effects
        let mut nm = no_kozai;
        if self.resonance != ResonanceType::None {
            nm = self.resonance_secular(tsince, nm);
        }

        (xm, xn, xg, em, xinc, nm)
    }

    /// Secular resonance integration for 12-hour and 24-hour orbits.
    fn resonance_secular(&self, tsince: f64, no_kozai: f64) -> f64 {
        // Simplified resonance: adjust mean motion based on resonance terms
        match self.resonance {
            ResonanceType::Synchronous => {
                let fasx2 = 0.13130908;
                let fasx4 = 2.8843198;
                let fasx6 = 0.37448087;

                let xli = self.xlamo + self.xfact * tsince;

                // Apply synchronous resonance correction
                let delt = self.del1 * (xli - fasx2).sin()
                    + self.del2 * (2.0 * (xli - fasx4)).sin()
                    + self.del3 * (3.0 * (xli - fasx6)).sin();

                no_kozai + delt
            }
            ResonanceType::HalfDay => {
                let g22 = 5.7686396;
                let g32 = 0.95240898;
                let g44 = 1.8014998;
                let g52 = 1.0508330;
                let g54 = 4.4108898;

                let xli = self.xlamo + self.xfact * tsince;

                // Apply half-day resonance correction
                let delt = self.d2201 * (2.0 * xli - g22).sin()
                    + self.d2211 * xli.sin()
                    + self.d3210 * (xli + g32).sin()
                    + self.d3222 * (-xli + g32).sin()
                    + self.d4410 * (2.0 * xli + g44).sin()
                    + self.d4422 * (2.0 * xli).sin()
                    + self.d5220 * (xli - g52).sin()
                    + self.d5232 * (-xli + g52).sin()
                    + self.d5421 * (xli + g54).sin()
                    + self.d5433 * (-xli + g54).sin();

                no_kozai + delt
            }
            ResonanceType::None => no_kozai,
        }
    }

    /// Apply deep-space periodic perturbations.
    ///
    /// Returns updated (ep, xincp, argpp, nodep, mp).
    pub(crate) fn periodics(
        &self,
        tsince: f64,
        ep: f64,
        xincp: f64,
        argpp: f64,
        nodep: f64,
        mp: f64,
    ) -> (f64, f64, f64, f64, f64) {
        // Solar and lunar mean anomalies
        let zm = (self.zmos + 1.19459e-5 * tsince).rem_euclid(TWOPI);
        let zf = (zm + 2.0 * (-0.01675) * zm.sin()).rem_euclid(TWOPI);
        let sinzf = zf.sin();
        let _coszf = zf.cos();

        let zml = (self.zmol + 1.5835218e-4 * tsince).rem_euclid(TWOPI);
        let sinzml = zml.sin();

        // Compute periodic corrections
        let mut pe = self.se2 * sinzf + self.se3 * (2.0 * zf).sin();
        let mut pinc = self.si2 * sinzf;
        let mut pl = self.sl2 * sinzf + self.sl3 * (2.0 * zf).sin();
        let mut pgh = self.sgh2 * sinzf + self.sgh3 * (2.0 * zf).sin();
        let mut ph = self.sh2 * sinzf;

        // Add lunar periodic corrections
        pe += self.ee2 * sinzml;
        pinc += self.xi2 * sinzml;
        pl += self.xl2 * sinzml;
        pgh += self.xgh2 * sinzml;
        ph += self.xh2 * sinzml;

        let ep_out = ep + pe;
        let xincp_out = xincp + pinc;
        let mut argpp_out = argpp + pgh;
        let mp_out = mp + pl;

        // Apply inclination-dependent node correction
        let sinip = xincp_out.sin();
        let cosip = xincp_out.cos();

        if xincp_out.abs() > 0.2 {
            // Standard correction
            let nodep_out = nodep + ph / sinip;
            argpp_out -= cosip * ph / sinip;
            (ep_out, xincp_out, argpp_out, nodep_out, mp_out)
        } else {
            // Near-equatorial: use sin/cos of node to avoid singularity
            let sinop = nodep.sin();
            let cosop = nodep.cos();

            let alfdp = sinip * sinop + ph * cosop;
            let betdp = sinip * cosop - ph * sinop;
            let nodep_out = alfdp.atan2(betdp);

            argpp_out += -cosip * (nodep_out - nodep);

            (ep_out, xincp_out, argpp_out, nodep_out, mp_out)
        }
    }
}

/// Compute semi-major axis from mean motion (earth radii).
fn ao_from_no(no_kozai: f64) -> f64 {
    (XKE / no_kozai).powf(2.0 / 3.0)
}

/// GMST from Julian date (for deep-space epoch initialization).
fn gmst_from_jd(jd: f64) -> f64 {
    let t_ut1 = (jd - 2451545.0) / 36525.0;
    let gmst_sec = 67310.54841
        + (876600.0 * 3600.0 + 8640184.812866) * t_ut1
        + 0.093104 * t_ut1 * t_ut1
        - 6.2e-6 * t_ut1 * t_ut1 * t_ut1;
    let gmst_rad = (gmst_sec / 86400.0).rem_euclid(1.0) * TWOPI;
    gmst_rad.rem_euclid(TWOPI)
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::f64::consts::PI;

    #[test]
    fn test_gmst_from_jd() {
        // At J2000.0, GMST ~ 280.46 degrees
        let gmst = gmst_from_jd(2451545.0);
        let deg = gmst * 180.0 / PI;
        assert!((deg - 280.46).abs() < 1.0, "GMST at J2000: {deg}");
    }

    #[test]
    fn test_ao_from_no() {
        // For no_kozai = XKE, semi-major axis = 1 earth radius
        let ao = ao_from_no(XKE);
        assert!((ao - 1.0).abs() < 1e-10);
    }
}