sidereon_core/astro/events/
eclipse.rs

1//! Conical Earth-shadow (eclipse) geometry.
2//!
3//! Determines whether a satellite is in full sunlight, penumbra, or umbra
4//! using a conical shadow model: the penumbra and umbra cones cast by Earth
5//! given the Sun's position, and the satellite's perpendicular distance from
6//! the Earth-Sun shadow axis. This is the authoritative implementation; the
7//! Elixir binding is a thin marshaling layer over it.
8
9use crate::astro::constants::{astro::SOLAR_RADIUS_KM, earth::MEAN_EARTH_RADIUS_KM};
10use crate::astro::math::vec3;
11
12/// Illumination state of a satellite relative to Earth's conical shadow.
13#[derive(Debug, Clone, Copy, PartialEq, Eq)]
14pub enum EclipseStatus {
15    /// Full sunlight.
16    Sunlit,
17    /// Partial shadow between the umbra and penumbra cones.
18    Penumbra,
19    /// Full shadow inside the umbra cone.
20    Umbra,
21}
22
23/// Error while computing conical eclipse geometry.
24#[derive(Debug, Clone, Copy, PartialEq, Eq, thiserror::Error)]
25pub enum EclipseError {
26    #[error("invalid eclipse input {field}: {reason}")]
27    InvalidInput {
28        field: &'static str,
29        reason: &'static str,
30    },
31}
32
33/// Shadow fraction in `[0.0, 1.0]`: `0.0` is full sunlight, `1.0` is full
34/// umbra, intermediate values are penumbra.
35///
36/// `sat_pos` is the satellite GCRS position (km); `sun_pos` is the vector from
37/// Earth center to the Sun (km).
38pub fn shadow_fraction(sat_pos: [f64; 3], sun_pos: [f64; 3]) -> Result<f64, EclipseError> {
39    validate_nonzero_vec3(sat_pos, "sat_pos")?;
40    validate_nonzero_vec3(sun_pos, "sun_pos")?;
41
42    let d_sun = vec3::norm3(sun_pos);
43
44    // Unit vector from Earth to Sun. The division form (not reciprocal-multiply
45    // via `vec3::unit3`) is required to stay bit-exact with the reference.
46    let sun_u = [sun_pos[0] / d_sun, sun_pos[1] / d_sun, sun_pos[2] / d_sun];
47
48    // Project the satellite onto the Earth-Sun line. Positive points toward the
49    // Sun; a satellite on the sunlit side cannot be in shadow.
50    let proj = vec3::dot3(sat_pos, sun_u);
51    if proj >= 0.0 {
52        return Ok(0.0);
53    }
54
55    // Perpendicular distance from the satellite to the shadow axis.
56    let perp = vec3::sub3(sat_pos, vec3::scale3(sun_u, proj));
57    let rho = vec3::norm3(perp);
58
59    // Distance behind Earth along the shadow axis (positive).
60    let dist_behind = -proj;
61
62    // Umbra cone (converging): sin(alpha) = (R_sun - R_earth) / d_sun.
63    let alpha_umbra = ((SOLAR_RADIUS_KM - MEAN_EARTH_RADIUS_KM) / d_sun).asin();
64    // Penumbra cone (diverging): sin(alpha) = (R_sun + R_earth) / d_sun.
65    let alpha_penumbra = ((SOLAR_RADIUS_KM + MEAN_EARTH_RADIUS_KM) / d_sun).asin();
66
67    // Cone radii at the satellite's distance behind Earth.
68    let r_umbra = MEAN_EARTH_RADIUS_KM - dist_behind * alpha_umbra.tan();
69    let r_penumbra = MEAN_EARTH_RADIUS_KM + dist_behind * alpha_penumbra.tan();
70
71    if rho >= r_penumbra {
72        // Beyond the penumbra cone: full sunlight.
73        Ok(0.0)
74    } else if r_umbra > 0.0 && rho <= r_umbra {
75        // Inside the umbra cone: full shadow.
76        Ok(1.0)
77    } else if r_umbra > 0.0 {
78        // Penumbra region: linear interpolation between the cone radii.
79        Ok((r_penumbra - rho) / (r_penumbra - r_umbra))
80    } else {
81        // Past the umbra tip (antumbra): penumbra still applies, fraction
82        // decreasing from the axis.
83        Ok(0.0_f64.max((r_penumbra - rho) / r_penumbra))
84    }
85}
86
87/// Classify the satellite's eclipse status from its [`shadow_fraction`].
88pub fn status(sat_pos: [f64; 3], sun_pos: [f64; 3]) -> Result<EclipseStatus, EclipseError> {
89    let fraction = shadow_fraction(sat_pos, sun_pos)?;
90    if fraction >= 1.0 {
91        Ok(EclipseStatus::Umbra)
92    } else if fraction > 0.0 {
93        Ok(EclipseStatus::Penumbra)
94    } else {
95        Ok(EclipseStatus::Sunlit)
96    }
97}
98
99fn validate_nonzero_vec3(v: [f64; 3], field: &'static str) -> Result<(), EclipseError> {
100    if !v.iter().all(|value| value.is_finite()) {
101        return Err(invalid_eclipse_input(field, "not finite"));
102    }
103    let norm = vec3::norm3(v);
104    if norm == 0.0 {
105        return Err(invalid_eclipse_input(field, "zero vector"));
106    }
107    if !norm.is_finite() {
108        return Err(invalid_eclipse_input(field, "out of range"));
109    }
110    Ok(())
111}
112
113fn invalid_eclipse_input(field: &'static str, reason: &'static str) -> EclipseError {
114    EclipseError::InvalidInput { field, reason }
115}
116
117#[cfg(test)]
118mod tests {
119    use crate::astro::events::{CrossingDirection, EventFinder};
120
121    use super::*;
122
123    // Approximate Sun-Earth distance, ~1 AU in km, matching the reference
124    // oracle inputs.
125    const AU_KM: f64 = 149_597_870.7;
126
127    #[test]
128    fn shadow_fraction_matches_reference_bits() {
129        // Frozen bits captured from the reference (Elixir) conical-shadow
130        // implementation. Cross-language 0-ULP equality.
131        let sun = [AU_KM, 0.0, 0.0];
132
133        assert_eq!(
134            shadow_fraction([7000.0, 0.0, 0.0], sun)
135                .expect("valid eclipse geometry")
136                .to_bits(),
137            0.0_f64.to_bits()
138        );
139        assert_eq!(
140            shadow_fraction([0.0, 7000.0, 0.0], sun)
141                .expect("valid eclipse geometry")
142                .to_bits(),
143            0.0_f64.to_bits()
144        );
145        assert_eq!(
146            shadow_fraction([-7000.0, 0.0, 0.0], sun)
147                .expect("valid eclipse geometry")
148                .to_bits(),
149            1.0_f64.to_bits()
150        );
151        assert_eq!(
152            shadow_fraction([-7000.0, 6370.0, 0.0], sun)
153                .expect("valid eclipse geometry")
154                .to_bits(),
155            0x3fe0_a32f_08e7_fb1f
156        );
157        assert_eq!(
158            shadow_fraction([-7000.0, 6410.0, 0.0], sun)
159                .expect("valid eclipse geometry")
160                .to_bits(),
161            0.0_f64.to_bits()
162        );
163        assert_eq!(
164            shadow_fraction([-7000.0, 6330.0, 0.0], sun)
165                .expect("valid eclipse geometry")
166                .to_bits(),
167            1.0_f64.to_bits()
168        );
169    }
170
171    #[test]
172    fn status_classifies_cones() {
173        let sun = [AU_KM, 0.0, 0.0];
174
175        assert_eq!(
176            status([7000.0, 0.0, 0.0], sun).expect("valid eclipse geometry"),
177            EclipseStatus::Sunlit
178        );
179        assert_eq!(
180            status([0.0, 7000.0, 0.0], sun).expect("valid eclipse geometry"),
181            EclipseStatus::Sunlit
182        );
183        assert_eq!(
184            status([-7000.0, 0.0, 0.0], sun).expect("valid eclipse geometry"),
185            EclipseStatus::Umbra
186        );
187        assert_eq!(
188            status([-7000.0, 6370.0, 0.0], sun).expect("valid eclipse geometry"),
189            EclipseStatus::Penumbra
190        );
191    }
192
193    #[test]
194    fn shadow_fraction_increases_toward_axis() {
195        let sun = [AU_KM, 0.0, 0.0];
196        let outside = shadow_fraction([-7000.0, 6410.0, 0.0], sun).expect("valid eclipse geometry");
197        let penumbra =
198            shadow_fraction([-7000.0, 6370.0, 0.0], sun).expect("valid eclipse geometry");
199        let umbra = shadow_fraction([-7000.0, 6330.0, 0.0], sun).expect("valid eclipse geometry");
200        let center = shadow_fraction([-7000.0, 0.0, 0.0], sun).expect("valid eclipse geometry");
201
202        assert_eq!(outside, 0.0);
203        assert!(penumbra > 0.0 && penumbra < 1.0);
204        assert_eq!(umbra, 1.0);
205        assert_eq!(center, 1.0);
206    }
207
208    #[test]
209    fn event_finder_refines_eclipse_fraction_crossings() {
210        let sun = [AU_KM, 0.0, 0.0];
211        let x_km = -7000.0;
212        let start_y_km = -7000.0;
213        let end_y_km = 7000.0;
214        let duration_seconds = 1200.0;
215        let threshold = 0.5;
216        let y_at = |time_seconds: f64| {
217            start_y_km + (end_y_km - start_y_km) * time_seconds / duration_seconds
218        };
219        let fraction_at = |time_seconds: f64| {
220            shadow_fraction([x_km, y_at(time_seconds), 0.0], sun)
221                .expect("valid synthetic eclipse geometry")
222        };
223
224        let events = EventFinder::new(0.0, duration_seconds, 60.0, 1.0e-7)
225            .expect("valid eclipse finder")
226            .find_crossings(fraction_at, threshold)
227            .expect("finite eclipse predicate");
228
229        assert_eq!(events.len(), 2);
230        assert_eq!(events[0].direction, CrossingDirection::Rising);
231        assert_eq!(events[1].direction, CrossingDirection::Falling);
232        assert_eq!(events[0].threshold, threshold);
233        assert_eq!(events[1].threshold, threshold);
234
235        let half_shadow_radius = shadow_radius_for_fraction(threshold, -x_km);
236        let expected_enter =
237            time_for_y(-half_shadow_radius, start_y_km, end_y_km, duration_seconds);
238        let expected_exit = time_for_y(half_shadow_radius, start_y_km, end_y_km, duration_seconds);
239
240        assert_close(events[0].time_seconds, expected_enter, 1.0e-6);
241        assert_close(events[1].time_seconds, expected_exit, 1.0e-6);
242        assert_close(events[0].value, threshold, 1.0e-7);
243        assert_close(events[1].value, threshold, 1.0e-7);
244        assert!(fraction_at(0.0) < threshold);
245        assert!(fraction_at(duration_seconds * 0.5) > threshold);
246        assert!(fraction_at(duration_seconds) < threshold);
247    }
248
249    fn shadow_radius_for_fraction(fraction: f64, dist_behind_km: f64) -> f64 {
250        let alpha_umbra = ((SOLAR_RADIUS_KM - MEAN_EARTH_RADIUS_KM) / AU_KM).asin();
251        let alpha_penumbra = ((SOLAR_RADIUS_KM + MEAN_EARTH_RADIUS_KM) / AU_KM).asin();
252        let r_umbra = MEAN_EARTH_RADIUS_KM - dist_behind_km * alpha_umbra.tan();
253        let r_penumbra = MEAN_EARTH_RADIUS_KM + dist_behind_km * alpha_penumbra.tan();
254        r_penumbra - fraction * (r_penumbra - r_umbra)
255    }
256
257    fn time_for_y(y_km: f64, start_y_km: f64, end_y_km: f64, duration_seconds: f64) -> f64 {
258        (y_km - start_y_km) * duration_seconds / (end_y_km - start_y_km)
259    }
260
261    fn assert_close(actual: f64, expected: f64, tolerance: f64) {
262        assert!(
263            (actual - expected).abs() <= tolerance,
264            "{actual} differs from {expected} by more than {tolerance}"
265        );
266    }
267
268    #[test]
269    fn eclipse_helpers_reject_invalid_vectors() {
270        assert_invalid_eclipse_field(
271            shadow_fraction([0.0, 0.0, 0.0], [AU_KM, 0.0, 0.0]).unwrap_err(),
272            "sat_pos",
273            "zero vector",
274        );
275        assert_invalid_eclipse_field(
276            status([7000.0, 0.0, 0.0], [f64::NAN, 0.0, 0.0]).unwrap_err(),
277            "sun_pos",
278            "not finite",
279        );
280    }
281
282    fn assert_invalid_eclipse_field(
283        error: EclipseError,
284        expected: &'static str,
285        expected_reason: &'static str,
286    ) {
287        let EclipseError::InvalidInput { field, reason } = error;
288        assert_eq!(field, expected);
289        assert_eq!(reason, expected_reason);
290    }
291}
sidereon_core/astro/events/eclipse.rs

sidereon_core/astro/events/
eclipse.rs
