siderust 0.9.0

// SPDX-License-Identifier: AGPL-3.0-or-later
// Copyright (C) 2026 Vallés Puig, Ramon

//! Frame-rotation providers for **inertial and quasi-inertial** frames.
//!
//! This module covers rotations between ICRS (the hub), the J2000 ecliptic
//! and equatorial frames, the mean/true equators of date, and the
//! identity-equivalent frames ICRF and GCRS.
//!
//! Key rotation chains:
//!
//! - **Frame bias**: ICRS ↔ EquatorialMeanJ2000 (constant ~80 µas matrix).
//! - **Obliquity**: EquatorialMeanJ2000 ↔ EclipticMeanJ2000 (Rx by −ε₀).
//! - **Precession**: EquatorialMeanJ2000 → EquatorialMeanOfDate (IAU 2006).
//! - **Nutation**: EquatorialMeanOfDate → EquatorialTrueOfDate (IAU 2000B).
//! - **EME2000 ≡ EquatorialMeanJ2000**: identity alias.
//! - **ICRF ≡ ICRS**: identity alias.
//! - **GCRS ≈ ICRS**: treated as identity (no aberration modelling).

use super::*;
use crate::coordinates::transform::frames::bias;

macro_rules! impl_alias_frame_rotations {
    ($alias:ty, $base:ty; $($target:ty),+ $(,)?) => {
        $(
            impl FrameRotationProvider<$alias, $target> for () {
                #[inline]
                fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
                    jd: JulianDate,
                    ctx: &AstroContext<Eph, Eop>,
                ) -> Rotation3 {
                    <() as FrameRotationProvider<$base, $target>>::rotation::<Eph, Eop, Nut>(jd, ctx)
                }
            }

            impl FrameRotationProvider<$target, $alias> for () {
                #[inline]
                fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
                    jd: JulianDate,
                    ctx: &AstroContext<Eph, Eop>,
                ) -> Rotation3 {
                    inverse_rotation::<$target, $alias, Eph, Eop, Nut>(jd, ctx)
                }
            }
        )+
    };
}

/// ICRS → EclipticMeanJ2000 rotation (J2000 mean ecliptic).
impl FrameRotationProvider<ICRS, EclipticMeanJ2000> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        _jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        bias::icrs_to_ecliptic_j2000()
    }
}

impl FrameRotationProvider<EclipticMeanJ2000, ICRS> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        inverse_rotation::<EclipticMeanJ2000, ICRS, Eph, Eop, Nut>(jd, ctx)
    }
}

impl FrameRotationProvider<ICRS, EquatorialMeanJ2000> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        _jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        bias::frame_bias_icrs_to_j2000()
    }
}

impl FrameRotationProvider<EquatorialMeanJ2000, ICRS> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        _jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        bias::frame_bias_j2000_to_icrs()
    }
}

impl FrameRotationProvider<EME2000, EquatorialMeanJ2000> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        _jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        Rotation3::IDENTITY
    }
}

impl FrameRotationProvider<EquatorialMeanJ2000, EME2000> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        inverse_rotation::<EquatorialMeanJ2000, EME2000, Eph, Eop, Nut>(jd, ctx)
    }
}

impl_alias_frame_rotations!(
    EME2000,
    EquatorialMeanJ2000;
    ICRS,
    EclipticMeanJ2000,
    EquatorialMeanOfDate,
    EquatorialTrueOfDate,
    ICRF,
    GCRSFrame,
);

impl FrameRotationProvider<EquatorialMeanJ2000, EclipticMeanJ2000> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        _jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        bias::obliquity_eq_to_ecl()
    }
}

impl FrameRotationProvider<EclipticMeanJ2000, EquatorialMeanJ2000> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        _jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        bias::obliquity_ecl_to_eq()
    }
}

impl FrameRotationProvider<EquatorialMeanJ2000, EquatorialMeanOfDate> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        // precession_matrix_iau2006 returns the full Fukushima-Williams matrix
        // (ICRS → MeanOfDate, including frame bias).  Strip the bias to get
        // pure precession: P = FW(date) · rb⁻¹.
        precession::precession_matrix_iau2006(jd) * bias::frame_bias_j2000_to_icrs()
    }
}

impl FrameRotationProvider<EquatorialMeanOfDate, EquatorialMeanJ2000> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        (precession::precession_matrix_iau2006(jd) * bias::frame_bias_j2000_to_icrs()).inverse()
    }
}

impl FrameRotationProvider<EquatorialMeanOfDate, EquatorialTrueOfDate> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        let nut = Nut::nutation(jd);
        Rotation3::fused_rx_rz_rx(nut.mean_obliquity + nut.deps, nut.dpsi, -nut.mean_obliquity)
    }
}

impl FrameRotationProvider<EquatorialTrueOfDate, EquatorialMeanOfDate> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        inverse_rotation::<EquatorialTrueOfDate, EquatorialMeanOfDate, Eph, Eop, Nut>(jd, ctx)
    }
}

impl FrameRotationProvider<EquatorialMeanJ2000, EquatorialTrueOfDate> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        let nut = Nut::nutation(jd);
        precession::precession_nutation_matrix(jd, nut.dpsi, nut.deps)
            * bias::frame_bias_j2000_to_icrs()
    }
}

impl FrameRotationProvider<EquatorialTrueOfDate, EquatorialMeanJ2000> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        let nut = Nut::nutation(jd);
        (precession::precession_nutation_matrix(jd, nut.dpsi, nut.deps)
            * bias::frame_bias_j2000_to_icrs())
        .inverse()
    }
}

impl FrameRotationProvider<ICRS, EquatorialMeanOfDate> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        precession::precession_matrix_iau2006(jd)
    }
}

impl FrameRotationProvider<EquatorialMeanOfDate, ICRS> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        precession::precession_matrix_iau2006(jd).inverse()
    }
}

/// ICRS → EquatorialTrueOfDate rotation.
///
/// Combines nutation (dispatched via `Nut`) with IERS celestial-pole
/// corrections (dX, dY) from the EOP provider into a full
/// precession–nutation matrix.
///
/// The first-order correction is:
///   dψ_eop = dX / sin(εA),  dε_eop = dY
impl FrameRotationProvider<ICRS, EquatorialTrueOfDate> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        let nut = Nut::nutation(jd);
        let mut dpsi = nut.dpsi;
        let mut deps = nut.deps;
        let eop = ctx.eop_at_tt(jd);
        let dx_rad = crate::qtty::Radians::from(eop.dx);
        let dy_rad = crate::qtty::Radians::from(eop.dy);

        let zero = crate::qtty::Radians::new(0.0);
        if dx_rad != zero || dy_rad != zero {
            let sin_eps = nut.mean_obliquity.sin();
            if sin_eps.abs() > 1e-15 {
                dpsi += dx_rad / sin_eps;
            }
            deps += dy_rad;
        }

        precession::precession_nutation_matrix(jd, dpsi, deps)
    }
}

impl FrameRotationProvider<EquatorialTrueOfDate, ICRS> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        inverse_rotation::<EquatorialTrueOfDate, ICRS, Eph, Eop, Nut>(jd, ctx)
    }
}

impl FrameRotationProvider<ICRF, ICRS> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        _jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        Rotation3::IDENTITY
    }
}

impl FrameRotationProvider<ICRS, ICRF> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        jd: JulianDate,
        ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        inverse_rotation::<ICRS, ICRF, Eph, Eop, Nut>(jd, ctx)
    }
}

impl_alias_frame_rotations!(
    ICRF,
    ICRS;
    EquatorialMeanJ2000,
    EclipticMeanJ2000,
    EquatorialMeanOfDate,
    EquatorialTrueOfDate,
);

impl FrameRotationProvider<GCRSFrame, ICRS> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        _jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        Rotation3::IDENTITY
    }
}

impl FrameRotationProvider<ICRS, GCRSFrame> for () {
    #[inline]
    fn rotation<Eph, Eop: EopProvider, Nut: NutationModel>(
        _jd: JulianDate,
        _ctx: &AstroContext<Eph, Eop>,
    ) -> Rotation3 {
        Rotation3::IDENTITY
    }
}

impl_alias_frame_rotations!(
    GCRSFrame,
    ICRS;
    EquatorialMeanJ2000,
    EquatorialTrueOfDate,
    EclipticMeanJ2000,
);