astrodynamics 0.12.1

//! Analytic low-precision Sun and Moon positions in ECI and ECEF.
//!
//! The Earth-centered-inertial series below are the standard low-precision
//! analytic solar/lunar formulae (Montenbruck & Gill, "Satellite Orbits",
//! sec. 3.3.2 / 3.3.3) as carried by common open GNSS code. They take the time
//! argument `t` in Julian centuries of Terrestrial Time since the J2000.0 epoch
//! (2000-01-01T12:00:00 TT). The five fundamental (Delaunay) arguments
//! `{l, l', F, D, Omega}` are evaluated by [`ast_args`] from the IAU 1980
//! nutation series coefficients.
//!
//! The series are referred to the mean equator and equinox of date (their mean
//! longitude and obliquity are of-date), so the resulting vectors are rotated
//! into the Earth-fixed (ITRS) frame with
//! [`crate::frames::transforms::mean_of_date_to_itrs_matrix`] (IAU 2000A
//! nutation + GAST). Precession is already implicit in the of-date series and is
//! deliberately NOT re-applied; the full GCRS->ITRS transform would double-count
//! it (about 0.37 deg of Sun-direction error at epoch 2026). This mirrors the
//! GMST/GAST + nutation rotation the analytic references consume the series with.

use crate::constants::astro::MONTENBRUCK_AU_M;
use crate::constants::earth::WGS84_A_M;
use crate::constants::time::{DAYS_PER_JULIAN_CENTURY, J2000_JD};
use crate::constants::units::{ARCSEC_TO_RAD, DEG_TO_RAD};
use crate::frames::transforms::{mat3_vec3_mul, mean_of_date_to_itrs_matrix};
use crate::time::scales::TimeScales;

/// Sun and Moon positions (metres) in a single frame.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct SunMoon {
    /// Sun position vector (m).
    pub sun: [f64; 3],
    /// Moon position vector (m).
    pub moon: [f64; 3],
}

/// Fundamental (Delaunay) arguments `f = {l, l', F, D, Omega}` in radians for
/// the given time `t` (Julian centuries TT since J2000). IAU 1980 series.
fn ast_args(t: f64) -> [f64; 5] {
    // Coefficients for IAU 1980 nutation: l, l', F, D, OMG.
    const FC: [[f64; 5]; 5] = [
        [
            134.96340251,
            1717915923.2178,
            31.8792,
            0.051635,
            -0.00024470,
        ],
        [357.52910918, 129596581.0481, -0.5532, 0.000136, -0.00001149],
        [
            93.27209062,
            1739527262.8478,
            -12.7512,
            -0.001037,
            0.00000417,
        ],
        [
            297.85019547,
            1602961601.2090,
            -6.3706,
            0.006593,
            -0.00003169,
        ],
        [125.04455501, -6962890.2665, 7.4722, 0.007702, -0.00005939],
    ];

    let mut tt = [0.0_f64; 4];
    tt[0] = t;
    for i in 1..4 {
        tt[i] = tt[i - 1] * t;
    }

    let mut f = [0.0_f64; 5];
    for i in 0..5 {
        let mut v = FC[i][0] * 3600.0;
        for j in 0..4 {
            v += FC[i][j + 1] * tt[j];
        }
        f[i] = (v * ARCSEC_TO_RAD).rem_euclid(2.0 * std::f64::consts::PI);
    }
    f
}

/// Analytic Sun and Moon position in ECI (metres) for time `t` in Julian
/// centuries of TT since J2000.0.
pub fn sun_moon_eci(t: f64) -> SunMoon {
    let f = ast_args(t);

    // Obliquity of the ecliptic (deg -> rad).
    let eps = 23.439291 - 0.0130042 * t;
    let sine = (eps * DEG_TO_RAD).sin();
    let cose = (eps * DEG_TO_RAD).cos();

    // Sun position in ECI.
    let ms = 357.5277233 + 35999.05034 * t;
    let ls = 280.460
        + 36000.770 * t
        + 1.914666471 * (ms * DEG_TO_RAD).sin()
        + 0.019994643 * (2.0 * ms * DEG_TO_RAD).sin();
    let rs = MONTENBRUCK_AU_M
        * (1.000140612
            - 0.016708617 * (ms * DEG_TO_RAD).cos()
            - 0.000139589 * (2.0 * ms * DEG_TO_RAD).cos());
    let sinl = (ls * DEG_TO_RAD).sin();
    let cosl = (ls * DEG_TO_RAD).cos();
    let sun = [rs * cosl, rs * cose * sinl, rs * sine * sinl];

    // Moon position in ECI.
    let lm = 218.32 + 481267.883 * t + 6.29 * f[0].sin() - 1.27 * (f[0] - 2.0 * f[3]).sin()
        + 0.66 * (2.0 * f[3]).sin()
        + 0.21 * (2.0 * f[0]).sin()
        - 0.19 * f[1].sin()
        - 0.11 * (2.0 * f[2]).sin();
    let pm = 5.13 * f[2].sin() + 0.28 * (f[0] + f[2]).sin()
        - 0.28 * (f[2] - f[0]).sin()
        - 0.17 * (f[2] - 2.0 * f[3]).sin();
    let rm = WGS84_A_M
        / ((0.9508
            + 0.0518 * f[0].cos()
            + 0.0095 * (f[0] - 2.0 * f[3]).cos()
            + 0.0078 * (2.0 * f[3]).cos()
            + 0.0028 * (2.0 * f[0]).cos())
            * DEG_TO_RAD)
            .sin();
    let sinlm = (lm * DEG_TO_RAD).sin();
    let coslm = (lm * DEG_TO_RAD).cos();
    let sinp = (pm * DEG_TO_RAD).sin();
    let cosp = (pm * DEG_TO_RAD).cos();
    let moon = [
        rm * cosp * coslm,
        rm * (cose * cosp * sinlm - sine * sinp),
        rm * (sine * cosp * sinlm + cose * sinp),
    ];

    SunMoon { sun, moon }
}

/// Analytic Sun and Moon geocentric positions in the Earth-fixed (ITRS) frame
/// (metres) for the given UTC instant.
///
/// The analytic series ([`sun_moon_eci`]) are referred to the **mean equator and
/// equinox of date** (of-date mean longitude and obliquity), so they are rotated
/// to ITRS with [`mean_of_date_to_itrs_matrix`] (nutation + GAST). Precession is
/// already implicit in the of-date series and must NOT be applied a second time;
/// using the full GCRS->ITRS transform here would double-count precession (about
/// 0.37 deg of Sun-direction error at epoch 2026).
pub fn sun_moon_ecef(ts: &TimeScales) -> SunMoon {
    let jd_tt = ts.jd_tt;
    let t = (jd_tt - J2000_JD) / DAYS_PER_JULIAN_CENTURY;

    let eci = sun_moon_eci(t);
    let r = mean_of_date_to_itrs_matrix(ts);
    let sun = mat3_vec3_mul(&r, &eci.sun);
    let moon = mat3_vec3_mul(&r, &eci.moon);
    SunMoon { sun, moon }
}