astrodynamics 0.13.0

//! Conical Earth-shadow (eclipse) geometry.
//!
//! Determines whether a satellite is in full sunlight, penumbra, or umbra
//! using a conical shadow model: the penumbra and umbra cones cast by Earth
//! given the Sun's position, and the satellite's perpendicular distance from
//! the Earth-Sun shadow axis. This is the authoritative implementation; the
//! Elixir binding is a thin marshaling layer over it.

use crate::constants::{astro::SOLAR_RADIUS_KM, earth::MEAN_EARTH_RADIUS_KM};
use crate::math::vec3;

/// Illumination state of a satellite relative to Earth's conical shadow.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum EclipseStatus {
    /// Full sunlight.
    Sunlit,
    /// Partial shadow between the umbra and penumbra cones.
    Penumbra,
    /// Full shadow inside the umbra cone.
    Umbra,
}

/// Shadow fraction in `[0.0, 1.0]`: `0.0` is full sunlight, `1.0` is full
/// umbra, intermediate values are penumbra.
///
/// `sat_pos` is the satellite GCRS position (km); `sun_pos` is the vector from
/// Earth center to the Sun (km).
pub fn shadow_fraction(sat_pos: [f64; 3], sun_pos: [f64; 3]) -> f64 {
    let d_sun = vec3::norm3(sun_pos);

    // Unit vector from Earth to Sun. The division form (not reciprocal-multiply
    // via `vec3::unit3`) is required to stay bit-exact with the reference.
    let sun_u = [sun_pos[0] / d_sun, sun_pos[1] / d_sun, sun_pos[2] / d_sun];

    // Project the satellite onto the Earth-Sun line. Positive points toward the
    // Sun; a satellite on the sunlit side cannot be in shadow.
    let proj = vec3::dot3(sat_pos, sun_u);
    if proj >= 0.0 {
        return 0.0;
    }

    // Perpendicular distance from the satellite to the shadow axis.
    let perp = vec3::sub3(sat_pos, vec3::scale3(sun_u, proj));
    let rho = vec3::norm3(perp);

    // Distance behind Earth along the shadow axis (positive).
    let dist_behind = -proj;

    // Umbra cone (converging): sin(alpha) = (R_sun - R_earth) / d_sun.
    let alpha_umbra = ((SOLAR_RADIUS_KM - MEAN_EARTH_RADIUS_KM) / d_sun).asin();
    // Penumbra cone (diverging): sin(alpha) = (R_sun + R_earth) / d_sun.
    let alpha_penumbra = ((SOLAR_RADIUS_KM + MEAN_EARTH_RADIUS_KM) / d_sun).asin();

    // Cone radii at the satellite's distance behind Earth.
    let r_umbra = MEAN_EARTH_RADIUS_KM - dist_behind * alpha_umbra.tan();
    let r_penumbra = MEAN_EARTH_RADIUS_KM + dist_behind * alpha_penumbra.tan();

    if rho >= r_penumbra {
        // Beyond the penumbra cone: full sunlight.
        0.0
    } else if r_umbra > 0.0 && rho <= r_umbra {
        // Inside the umbra cone: full shadow.
        1.0
    } else if r_umbra > 0.0 {
        // Penumbra region: linear interpolation between the cone radii.
        (r_penumbra - rho) / (r_penumbra - r_umbra)
    } else {
        // Past the umbra tip (antumbra): penumbra still applies, fraction
        // decreasing from the axis.
        0.0_f64.max((r_penumbra - rho) / r_penumbra)
    }
}

/// Classify the satellite's eclipse status from its [`shadow_fraction`].
pub fn status(sat_pos: [f64; 3], sun_pos: [f64; 3]) -> EclipseStatus {
    let fraction = shadow_fraction(sat_pos, sun_pos);
    if fraction >= 1.0 {
        EclipseStatus::Umbra
    } else if fraction > 0.0 {
        EclipseStatus::Penumbra
    } else {
        EclipseStatus::Sunlit
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    // Approximate Sun-Earth distance, ~1 AU in km, matching the reference
    // oracle inputs.
    const AU_KM: f64 = 149_597_870.7;

    #[test]
    fn shadow_fraction_matches_reference_bits() {
        // Frozen bits captured from the reference (Elixir) conical-shadow
        // implementation. Cross-language 0-ULP equality.
        let sun = [AU_KM, 0.0, 0.0];

        assert_eq!(
            shadow_fraction([7000.0, 0.0, 0.0], sun).to_bits(),
            0.0_f64.to_bits()
        );
        assert_eq!(
            shadow_fraction([0.0, 7000.0, 0.0], sun).to_bits(),
            0.0_f64.to_bits()
        );
        assert_eq!(
            shadow_fraction([-7000.0, 0.0, 0.0], sun).to_bits(),
            1.0_f64.to_bits()
        );
        assert_eq!(
            shadow_fraction([-7000.0, 6370.0, 0.0], sun).to_bits(),
            0x3fe0_a32f_08e7_fb1f
        );
        assert_eq!(
            shadow_fraction([-7000.0, 6410.0, 0.0], sun).to_bits(),
            0.0_f64.to_bits()
        );
        assert_eq!(
            shadow_fraction([-7000.0, 6330.0, 0.0], sun).to_bits(),
            1.0_f64.to_bits()
        );
    }

    #[test]
    fn status_classifies_cones() {
        let sun = [AU_KM, 0.0, 0.0];

        assert_eq!(status([7000.0, 0.0, 0.0], sun), EclipseStatus::Sunlit);
        assert_eq!(status([0.0, 7000.0, 0.0], sun), EclipseStatus::Sunlit);
        assert_eq!(status([-7000.0, 0.0, 0.0], sun), EclipseStatus::Umbra);
        assert_eq!(status([-7000.0, 6370.0, 0.0], sun), EclipseStatus::Penumbra);
    }

    #[test]
    fn shadow_fraction_increases_toward_axis() {
        let sun = [AU_KM, 0.0, 0.0];
        let outside = shadow_fraction([-7000.0, 6410.0, 0.0], sun);
        let penumbra = shadow_fraction([-7000.0, 6370.0, 0.0], sun);
        let umbra = shadow_fraction([-7000.0, 6330.0, 0.0], sun);
        let center = shadow_fraction([-7000.0, 0.0, 0.0], sun);

        assert_eq!(outside, 0.0);
        assert!(penumbra > 0.0 && penumbra < 1.0);
        assert_eq!(umbra, 1.0);
        assert_eq!(center, 1.0);
    }
}