From d2f94978f98b51f8e92a0e8f596364f30bdb5b92 Mon Sep 17 00:00:00 2001
From: Christopher Rabotin <christopher.rabotin@gmail.com>
Date: Sat, 7 Sep 2024 01:36:29 -0600
Subject: [PATCH] Initial work upgrading to anise 0.4.3

Update the GEO test case because the eclipsing
is more thoroughly validated in ANISE
---
 examples/01_orbit_prop/main.rs             |   4 +-
 examples/02_jwst_covar_monte_carlo/main.rs |   4 +-
 examples/03_geo_analysis/drift.rs          |   2 +-
 examples/03_geo_analysis/raise.rs          |  10 +-
 examples/03_geo_analysis/stationkeeping.rs |   9 +-
 examples/04_lro_od/main.rs                 |   4 +-
 src/cosmic/eclipse.rs                      | 467 ++-------------------
 src/dynamics/guidance/ruggiero.rs          |  26 +-
 src/dynamics/solarpressure.rs              |   4 +-
 src/dynamics/sph_harmonics.rs              |   4 +-
 src/od/ground_station.rs                   | 101 ++---
 tests/cosmic/eclipse.rs                    |  33 +-
 tests/lib.rs                               |   2 +-
 tests/propagation/events.rs                |  19 +-
 14 files changed, 135 insertions(+), 554 deletions(-)

diff --git a/examples/01_orbit_prop/main.rs b/examples/01_orbit_prop/main.rs
index c2f46b0b..451f11c2 100644
--- a/examples/01_orbit_prop/main.rs
+++ b/examples/01_orbit_prop/main.rs
@@ -112,7 +112,7 @@ fn main() -> Result<(), Box<dyn Error>> {
         crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
     };
     // And let's download it if we don't have it yet.
-    jgm3_meta.process()?;
+    jgm3_meta.process(true)?;
 
     // Build the spherical harmonics.
     // The harmonics must be computed in the body fixed frame.
@@ -237,6 +237,8 @@ fn main() -> Result<(), Box<dyn Error>> {
                 let aer = almanac.azimuth_elevation_range_sez(
                     state.orbit,
                     boulder_station.to_orbit(this_epoch, &almanac)?,
+                    None,
+                    None,
                 )?;
                 azimuth_deg.push(aer.azimuth_deg);
                 elevation_deg.push(aer.elevation_deg);
diff --git a/examples/02_jwst_covar_monte_carlo/main.rs b/examples/02_jwst_covar_monte_carlo/main.rs
index e6b33a36..cf58f687 100644
--- a/examples/02_jwst_covar_monte_carlo/main.rs
+++ b/examples/02_jwst_covar_monte_carlo/main.rs
@@ -35,13 +35,13 @@ fn main() -> Result<(), Box<dyn Error>> {
         uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
         crc32: None,
     };
-    latest_jwst_ephem.process()?;
+    latest_jwst_ephem.process(true)?;
 
     // Load this ephem in the general Almanac we're using for this analysis.
     let almanac = Arc::new(
         MetaAlmanac::latest()
             .map_err(Box::new)?
-            .load_from_metafile(latest_jwst_ephem)?,
+            .load_from_metafile(latest_jwst_ephem, true)?,
     );
 
     // By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
diff --git a/examples/03_geo_analysis/drift.rs b/examples/03_geo_analysis/drift.rs
index ca590a62..82fa541f 100644
--- a/examples/03_geo_analysis/drift.rs
+++ b/examples/03_geo_analysis/drift.rs
@@ -82,7 +82,7 @@ fn main() -> Result<(), Box<dyn Error>> {
         crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
     };
     // And let's download it if we don't have it yet.
-    jgm3_meta.process()?;
+    jgm3_meta.process(true)?;
 
     // Build the spherical harmonics.
     // The harmonics must be computed in the body fixed frame.
diff --git a/examples/03_geo_analysis/raise.rs b/examples/03_geo_analysis/raise.rs
index 6d7750f7..649caf61 100644
--- a/examples/03_geo_analysis/raise.rs
+++ b/examples/03_geo_analysis/raise.rs
@@ -12,10 +12,7 @@ use anise::{
 };
 use hifitime::{Epoch, TimeUnits, Unit};
 use nyx::{
-    cosmic::{
-        eclipse::{EclipseLocator, EclipseState},
-        GuidanceMode, MetaAlmanac, Orbit, SrpConfig,
-    },
+    cosmic::{eclipse::EclipseLocator, GuidanceMode, MetaAlmanac, Orbit, SrpConfig},
     dynamics::{
         guidance::{GuidanceLaw, Ruggiero, Thruster},
         Harmonics, OrbitalDynamics, SolarPressure, SpacecraftDynamics,
@@ -75,8 +72,7 @@ fn main() -> Result<(), Box<dyn Error>> {
     ];
 
     // Ensure that we only thrust if we have more than 20% illumination.
-    let ruggiero_ctrl =
-        Ruggiero::from_max_eclipse(objectives, sc, EclipseState::Penumbra(0.2)).unwrap();
+    let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2).unwrap();
     println!("{ruggiero_ctrl}");
 
     // Define the high fidelity dynamics
@@ -94,7 +90,7 @@ fn main() -> Result<(), Box<dyn Error>> {
         crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
     };
     // And let's download it if we don't have it yet.
-    jgm3_meta.process()?;
+    jgm3_meta.process(true)?;
 
     // Build the spherical harmonics.
     // The harmonics must be computed in the body fixed frame.
diff --git a/examples/03_geo_analysis/stationkeeping.rs b/examples/03_geo_analysis/stationkeeping.rs
index ed8ecef1..afc6358c 100644
--- a/examples/03_geo_analysis/stationkeeping.rs
+++ b/examples/03_geo_analysis/stationkeeping.rs
@@ -12,10 +12,7 @@ use anise::{
 };
 use hifitime::{Epoch, TimeUnits, Unit};
 use nyx::{
-    cosmic::{
-        eclipse::{EclipseLocator, EclipseState},
-        GuidanceMode, MetaAlmanac, Orbit, SrpConfig,
-    },
+    cosmic::{eclipse::EclipseLocator, GuidanceMode, MetaAlmanac, Orbit, SrpConfig},
     dynamics::{
         guidance::{Ruggiero, Thruster},
         Harmonics, OrbitalDynamics, SolarPressure, SpacecraftDynamics,
@@ -63,7 +60,7 @@ fn main() -> Result<(), Box<dyn Error>> {
         Objective::within_tolerance(StateParameter::Inclination, 0.05, 1e-2),
     ];
 
-    let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, EclipseState::Penumbra(0.2))?;
+    let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2)?;
     println!("{ruggiero_ctrl}");
 
     let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
@@ -72,7 +69,7 @@ fn main() -> Result<(), Box<dyn Error>> {
         uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
         crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
     };
-    jgm3_meta.process()?;
+    jgm3_meta.process(true)?;
 
     let harmonics = Harmonics::from_stor(
         almanac.frame_from_uid(IAU_EARTH_FRAME)?,
diff --git a/examples/04_lro_od/main.rs b/examples/04_lro_od/main.rs
index ab926aa9..3a2a724f 100644
--- a/examples/04_lro_od/main.rs
+++ b/examples/04_lro_od/main.rs
@@ -50,7 +50,7 @@ fn main() -> Result<(), Box<dyn Error>> {
     // Load this ephem in the general Almanac we're using for this analysis.
     let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
         .map_err(Box::new)?
-        .process()
+        .process(true)
         .map_err(Box::new)?;
 
     let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
@@ -120,7 +120,7 @@ fn main() -> Result<(), Box<dyn Error>> {
         crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
     };
     // And let's download it if we don't have it yet.
-    jggrx_meta.process()?;
+    jggrx_meta.process(true)?;
 
     // Build the spherical harmonics.
     // The harmonics must be computed in the body fixed frame.
diff --git a/src/cosmic/eclipse.rs b/src/cosmic/eclipse.rs
index 1f5c1976..b4f1271e 100644
--- a/src/cosmic/eclipse.rs
+++ b/src/cosmic/eclipse.rs
@@ -17,107 +17,18 @@
 */
 
 use anise::almanac::Almanac;
+use anise::astro::Occultation;
 use anise::constants::frames::{EARTH_J2000, MOON_J2000, SUN_J2000};
-use anise::ephemerides::EphemerisPhysicsSnafu;
-use anise::errors::{AlmanacError, AlmanacResult, EphemerisSnafu};
+use anise::errors::AlmanacResult;
 use snafu::ResultExt;
 
 pub use super::{Frame, Orbit, Spacecraft};
 use crate::errors::{EventAlmanacSnafu, EventError};
 use crate::md::EventEvaluator;
 use crate::time::{Duration, Unit};
-use std::cmp::{Eq, Ord, Ordering, PartialOrd};
-use std::convert::Into;
 use std::fmt;
 use std::sync::Arc;
 
-/// Stores the eclipse state
-#[derive(Clone, Copy, Debug, PartialEq)]
-pub enum EclipseState {
-    Umbra,
-    /// The f64 is between ]0; 1[ and corresponds to the percentage of penumbra: the closer to 1, the more light is seen.
-    Penumbra(f64),
-    Visibilis,
-}
-
-#[allow(clippy::from_over_into)]
-impl Into<f64> for EclipseState {
-    fn into(self) -> f64 {
-        match self {
-            Self::Umbra => 0.0,
-            Self::Visibilis => 1.0,
-            Self::Penumbra(val) => val,
-        }
-    }
-}
-
-impl Eq for EclipseState {}
-
-impl Ord for EclipseState {
-    /// Orders eclipse states to the greatest eclipse.
-    ///
-    /// *Examples*
-    ///
-    /// ```
-    /// extern crate nyx_space as nyx;
-    /// use nyx::cosmic::eclipse::EclipseState;
-    /// assert!(EclipseState::Umbra == EclipseState::Umbra);
-    /// assert!(EclipseState::Visibilis == EclipseState::Visibilis);
-    /// assert!(EclipseState::Penumbra(0.5) == EclipseState::Penumbra(0.5));
-    /// assert!(EclipseState::Umbra > EclipseState::Penumbra(0.1));
-    /// assert!(EclipseState::Umbra > EclipseState::Penumbra(0.9));
-    /// assert!(EclipseState::Penumbra(0.1) < EclipseState::Umbra);
-    /// assert!(EclipseState::Penumbra(0.9) < EclipseState::Umbra);
-    /// assert!(EclipseState::Penumbra(0.9) > EclipseState::Visibilis);
-    /// assert!(EclipseState::Visibilis < EclipseState::Penumbra(0.9));
-    /// ```
-    fn cmp(&self, other: &Self) -> Ordering {
-        match *self {
-            EclipseState::Umbra => {
-                if *other == EclipseState::Umbra {
-                    Ordering::Equal
-                } else {
-                    Ordering::Greater
-                }
-            }
-            EclipseState::Visibilis => {
-                if *other == EclipseState::Visibilis {
-                    Ordering::Equal
-                } else {
-                    Ordering::Less
-                }
-            }
-            EclipseState::Penumbra(s) => match *other {
-                EclipseState::Penumbra(o) => {
-                    if s > o {
-                        Ordering::Greater
-                    } else {
-                        Ordering::Less
-                    }
-                }
-                EclipseState::Visibilis => Ordering::Greater,
-                EclipseState::Umbra => Ordering::Less,
-            },
-        }
-    }
-}
-
-impl PartialOrd for EclipseState {
-    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
-        Some(self.cmp(other))
-    }
-}
-
-impl fmt::Display for EclipseState {
-    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
-        match *self {
-            Self::Umbra => write!(f, "Umbra"),
-            Self::Visibilis => write!(f, "Visibilis"),
-            Self::Penumbra(v) => write!(f, "Penumbra {:.2}%", v * 100.0),
-        }
-    }
-}
-
 #[derive(Clone)]
 pub struct EclipseLocator {
     pub light_source: Frame,
@@ -153,11 +64,16 @@ impl EclipseLocator {
     }
 
     /// Compute the visibility/eclipse between an observer and an observed state
-    pub fn compute(&self, observer: Orbit, almanac: Arc<Almanac>) -> AlmanacResult<EclipseState> {
-        let mut state = EclipseState::Visibilis;
+    pub fn compute(&self, observer: Orbit, almanac: Arc<Almanac>) -> AlmanacResult<Occultation> {
+        let mut state = Occultation {
+            epoch: observer.epoch,
+            back_frame: SUN_J2000,
+            front_frame: observer.frame,
+            percentage: 0.0,
+        };
         for eclipsing_body in &self.shadow_bodies {
-            let this_state = eclipse_state(observer, self.light_source, *eclipsing_body, &almanac)?;
-            if this_state > state {
+            let this_state = almanac.solar_eclipsing(*eclipsing_body, observer, None)?;
+            if this_state.percentage > state.percentage {
                 state = this_state;
             }
         }
@@ -193,17 +109,22 @@ impl fmt::Display for UmbraEvent {
 }
 
 impl EventEvaluator<Spacecraft> for UmbraEvent {
-    // Evaluation of the event, returns 0.0 for umbra, 1.0 for visibility (no shadow) and some value in between for penumbra
+    // Evaluation of the event
     fn eval(&self, sc: &Spacecraft, almanac: Arc<Almanac>) -> Result<f64, EventError> {
-        match self
+        Ok(self
             .e_loc
             .compute(sc.orbit, almanac)
             .context(EventAlmanacSnafu)?
-        {
-            EclipseState::Umbra => Ok(0.0),
-            EclipseState::Visibilis => Ok(1.0),
-            EclipseState::Penumbra(val) => Ok(val),
-        }
+            .factor())
+        // match self
+        //     .e_loc
+        //     .compute(sc.orbit, almanac)
+        //     .context(EventAlmanacSnafu)?
+        // {
+        //     EclipseState::Umbra => Ok(0.0),
+        //     EclipseState::Visibilis => Ok(1.0),
+        //     EclipseState::Penumbra(val) => Ok(val),
+        // }
     }
 
     /// Stop searching when the time has converged to less than 0.1 seconds
@@ -237,15 +158,11 @@ impl fmt::Display for PenumbraEvent {
 
 impl EventEvaluator<Spacecraft> for PenumbraEvent {
     fn eval(&self, sc: &Spacecraft, almanac: Arc<Almanac>) -> Result<f64, EventError> {
-        match self
+        Ok(self
             .e_loc
             .compute(sc.orbit, almanac)
             .context(EventAlmanacSnafu)?
-        {
-            EclipseState::Umbra => Ok(0.0),
-            EclipseState::Visibilis => Ok(1.0),
-            EclipseState::Penumbra(val) => Ok(val),
-        }
+            .factor())
     }
 
     /// Stop searching when the time has converged to less than 0.1 seconds
@@ -266,337 +183,3 @@ impl EventEvaluator<Spacecraft> for PenumbraEvent {
         ))
     }
 }
-
-/// Computes the umbra/visibilis/penumbra state between between two states accounting for eclipsing of the providing geoid.
-pub fn eclipse_state(
-    observer: Orbit,
-    mut light_source: Frame,
-    mut eclipsing_body: Frame,
-    almanac: &Almanac,
-) -> AlmanacResult<EclipseState> {
-    if light_source.mean_equatorial_radius_km().is_err() {
-        light_source =
-            almanac
-                .frame_from_uid(light_source)
-                .map_err(|e| AlmanacError::GenericError {
-                    err: format!("{e} when fetching light source data ({light_source})"),
-                })?;
-    }
-
-    if eclipsing_body.mean_equatorial_radius_km().is_err() {
-        eclipsing_body =
-            almanac
-                .frame_from_uid(light_source)
-                .map_err(|e| AlmanacError::GenericError {
-                    err: format!("{e} when fetching eclipsing body data ({eclipsing_body})"),
-                })?;
-    }
-
-    let ls_mean_eq_radius_km = light_source
-        .mean_equatorial_radius_km()
-        .context(EphemerisPhysicsSnafu {
-            action: "fetching mean equatorial radius of light source",
-        })
-        .context(EphemerisSnafu {
-            action: "computing eclipse state",
-        })?;
-
-    // If the light source's radius is zero, just call the line of sight algorithm
-
-    if ls_mean_eq_radius_km < f64::EPSILON {
-        let observed = -almanac.transform_to(observer, light_source, None)?;
-        return line_of_sight(observer, observed, eclipsing_body, almanac);
-    }
-
-    // All of the computations happen with the observer as the center.
-    // `eb` stands for eclipsing body; `ls` stands for light source.
-    // Get the radius vector of the spacecraft to the eclipsing body
-
-    let r_eb = almanac
-        .transform_to(observer, eclipsing_body, None)?
-        .radius_km;
-
-    // Get the radius vector of the light source to the spacecraft
-    let r_ls = -almanac
-        .transform_to(observer, light_source, None)?
-        .radius_km;
-
-    // Compute the apparent radii of the light source and eclipsing body (preventing any NaN)
-    let r_ls_prime = if ls_mean_eq_radius_km >= r_ls.norm() {
-        ls_mean_eq_radius_km
-    } else {
-        (ls_mean_eq_radius_km / r_ls.norm()).asin()
-    };
-
-    let eb_mean_eq_radius_km = eclipsing_body
-        .mean_equatorial_radius_km()
-        .context(EphemerisPhysicsSnafu {
-            action: "fetching mean equatorial radius of eclipsing body",
-        })
-        .context(EphemerisSnafu {
-            action: "computing eclipse state",
-        })?;
-
-    let r_eb_prime = if eb_mean_eq_radius_km >= r_eb.norm() {
-        eb_mean_eq_radius_km
-    } else {
-        (eb_mean_eq_radius_km / r_eb.norm()).asin()
-    };
-
-    // Compute the apparent separation of both circles
-    let d_prime = (-(r_ls.dot(&r_eb)) / (r_eb.norm() * r_ls.norm())).acos();
-
-    if d_prime - r_ls_prime > r_eb_prime {
-        // If the closest point where the apparent radius of the light source _starts_ is further
-        // away than the furthest point where the eclipsing body's shadow can reach, then the light
-        // source is totally visible.
-        Ok(EclipseState::Visibilis)
-    } else if r_eb_prime > d_prime + r_ls_prime {
-        // The light source is fully hidden by the eclipsing body, hence we're in total eclipse.
-        Ok(EclipseState::Umbra)
-    } else if (r_ls_prime - r_eb_prime).abs() < d_prime && d_prime < r_ls_prime + r_eb_prime {
-        // If we have reached this point, we're in penumbra.
-        // Both circles, which represent the light source projected onto the plane and the eclipsing geoid,
-        // now overlap creating an asymmetrial lens.
-        // The following math comes from http://mathworld.wolfram.com/Circle-CircleIntersection.html
-        // and https://stackoverflow.com/questions/3349125/circle-circle-intersection-points .
-
-        // Compute the distances between the center of the eclipsing geoid and the line crossing the intersection
-        // points of both circles.
-        let d1 = (d_prime.powi(2) - r_ls_prime.powi(2) + r_eb_prime.powi(2)) / (2.0 * d_prime);
-        let d2 = (d_prime.powi(2) + r_ls_prime.powi(2) - r_eb_prime.powi(2)) / (2.0 * d_prime);
-
-        let shadow_area = circ_seg_area(r_eb_prime, d1) + circ_seg_area(r_ls_prime, d2);
-        if shadow_area.is_nan() {
-            error!(
-                "Shadow area is NaN! Please file a bug with initial states, eclipsing bodies, etc."
-            );
-            return Ok(EclipseState::Umbra);
-        }
-        // Compute the nominal area of the light source
-        let nominal_area = core::f64::consts::PI * r_ls_prime.powi(2);
-        // And return the percentage (between 0 and 1) of the eclipse.
-        Ok(EclipseState::Penumbra(1.0 - shadow_area / nominal_area))
-    } else {
-        // Annular eclipse.
-        // If r_eb_prime is very small, then the fraction is very small: however, we note a penumbra close to 1.0 as near full light source visibility, so let's subtract one from this.
-        Ok(EclipseState::Penumbra(
-            1.0 - r_eb_prime.powi(2) / r_ls_prime.powi(2),
-        ))
-    }
-}
-
-// Compute the area of the circular segment of radius r and chord length d
-fn circ_seg_area(r: f64, d: f64) -> f64 {
-    r.powi(2) * (d / r).acos() - d * (r.powi(2) - d.powi(2)).sqrt()
-}
-
-/// Computes the light of sight the provided time between two states accounting for eclipsing of the providing geoid.
-/// This works for visibility between spacecraft and a ground station. For eclipsing and penumbras, use `eclipse_state`.
-///
-/// Source: Algorithm 35 of Vallado, 4th edition, page 308.
-pub fn line_of_sight(
-    observer: Orbit,
-    observed: Orbit,
-    eclipsing_body: Frame,
-    almanac: &Almanac,
-) -> AlmanacResult<EclipseState> {
-    if observer == observed {
-        return Ok(EclipseState::Visibilis);
-    }
-
-    let eb_mean_eq_radius_km = eclipsing_body
-        .mean_equatorial_radius_km()
-        .context(EphemerisPhysicsSnafu {
-            action: "fetching mean equatorial radius of eclipsing body",
-        })
-        .context(EphemerisSnafu {
-            action: "computing eclipse state",
-        })?;
-
-    // Convert the states to the same frame as the eclipsing body (ensures we're in the same frame)
-    let r1 = almanac
-        .transform_to(observed, eclipsing_body, None)?
-        .radius_km;
-    let r2 = almanac
-        .transform_to(observer, eclipsing_body, None)?
-        .radius_km;
-
-    let r1sq = r1.dot(&r1);
-    let r2sq = r2.dot(&r2);
-    let r1dotr2 = r1.dot(&r2);
-
-    let tau = (r1sq - r1dotr2) / (r1sq + r2sq - 2.0 * r1dotr2);
-    if !(0.0..=1.0).contains(&tau)
-        || (1.0 - tau) * r1sq + r1dotr2 * tau > eb_mean_eq_radius_km.powi(2)
-    {
-        Ok(EclipseState::Visibilis)
-    } else {
-        Ok(EclipseState::Umbra)
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-    use hifitime::Epoch;
-    use rstest::*;
-
-    #[fixture]
-    pub fn almanac() -> Almanac {
-        use std::path::PathBuf;
-
-        let manifest_dir =
-            PathBuf::from(std::env::var("CARGO_MANIFEST_DIR").unwrap_or(".".to_string()));
-
-        Almanac::new(
-            &manifest_dir
-                .clone()
-                .join("data/de440s.bsp")
-                .to_string_lossy(),
-        )
-        .unwrap()
-        .load(
-            &manifest_dir
-                .clone()
-                .join("data/pck08.pca")
-                .to_string_lossy(),
-        )
-        .unwrap()
-        .load(
-            &manifest_dir
-                .join("data/earth_latest_high_prec.bpc")
-                .to_string_lossy(),
-        )
-        .unwrap()
-    }
-
-    #[rstest]
-    fn los_edge_case(almanac: Almanac) {
-        let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap();
-        let luna = almanac.frame_from_uid(MOON_J2000).unwrap();
-
-        let dt1 = Epoch::from_gregorian_tai_hms(2020, 1, 1, 6, 7, 40);
-        let dt2 = Epoch::from_gregorian_tai_hms(2020, 1, 1, 6, 7, 50);
-        let dt3 = Epoch::from_gregorian_tai_hms(2020, 1, 1, 6, 8, 0);
-
-        let xmtr1 = Orbit::new(
-            397_477.494_485,
-            -57_258.902_156,
-            -62_857.909_437,
-            0.230_482,
-            2.331_362,
-            0.615_501,
-            dt1,
-            eme2k,
-        );
-        let rcvr1 = Orbit::new(
-            338_335.467_589,
-            -55_439.526_977,
-            -13_327.354_273,
-            0.197_141,
-            0.944_261,
-            0.337_407,
-            dt1,
-            eme2k,
-        );
-        let xmtr2 = Orbit::new(
-            397_479.756_900,
-            -57_235.586_465,
-            -62_851.758_851,
-            0.222_000,
-            2.331_768,
-            0.614_614,
-            dt2,
-            eme2k,
-        );
-        let rcvr2 = Orbit::new(
-            338_337.438_860,
-            -55_430.084_340,
-            -13_323.980_229,
-            0.197_113,
-            0.944_266,
-            0.337_402,
-            dt2,
-            eme2k,
-        );
-        let xmtr3 = Orbit::new(
-            397_481.934_480,
-            -57_212.266_970,
-            -62_845.617_185,
-            0.213_516,
-            2.332_122,
-            0.613_717,
-            dt3,
-            eme2k,
-        );
-        let rcvr3 = Orbit::new(
-            338_339.409_858,
-            -55_420.641_651,
-            -13_320.606_228,
-            0.197_086,
-            0.944_272,
-            0.337_398,
-            dt3,
-            eme2k,
-        );
-
-        assert_eq!(
-            line_of_sight(xmtr1, rcvr1, luna, &almanac).unwrap(),
-            EclipseState::Umbra
-        );
-        assert_eq!(
-            line_of_sight(xmtr2, rcvr2, luna, &almanac).unwrap(),
-            EclipseState::Umbra
-        );
-        assert_eq!(
-            line_of_sight(xmtr3, rcvr3, luna, &almanac).unwrap(),
-            EclipseState::Umbra
-        );
-
-        // Test converse
-
-        assert_eq!(
-            line_of_sight(rcvr1, xmtr1, luna, &almanac).unwrap(),
-            EclipseState::Umbra
-        );
-        assert_eq!(
-            line_of_sight(rcvr2, xmtr2, luna, &almanac).unwrap(),
-            EclipseState::Umbra
-        );
-        assert_eq!(
-            line_of_sight(rcvr3, xmtr3, luna, &almanac).unwrap(),
-            EclipseState::Umbra
-        );
-    }
-
-    #[rstest]
-    fn los_earth_eclipse(almanac: Almanac) {
-        let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap();
-
-        let dt = Epoch::from_gregorian_tai_at_midnight(2020, 1, 1);
-
-        let sma = eme2k.mean_equatorial_radius_km().unwrap() + 300.0;
-
-        let sc1 = Orbit::keplerian(sma, 0.001, 0.1, 90.0, 75.0, 0.0, dt, eme2k);
-        let sc2 = Orbit::keplerian(sma + 1.0, 0.001, 0.1, 90.0, 75.0, 0.0, dt, eme2k);
-        let sc3 = Orbit::keplerian(sma, 0.001, 0.1, 90.0, 75.0, 180.0, dt, eme2k);
-
-        // Out of phase by pi.
-        assert_eq!(
-            line_of_sight(sc1, sc3, eme2k, &almanac).unwrap(),
-            EclipseState::Umbra
-        );
-
-        assert_eq!(
-            line_of_sight(sc2, sc1, eme2k, &almanac).unwrap(),
-            EclipseState::Visibilis
-        );
-
-        // Nearly identical orbits in the same phasing
-        assert_eq!(
-            line_of_sight(sc1, sc2, eme2k, &almanac).unwrap(),
-            EclipseState::Visibilis
-        );
-    }
-}
diff --git a/src/dynamics/guidance/ruggiero.rs b/src/dynamics/guidance/ruggiero.rs
index d953b0c6..c3eea4f7 100644
--- a/src/dynamics/guidance/ruggiero.rs
+++ b/src/dynamics/guidance/ruggiero.rs
@@ -23,7 +23,7 @@ use super::{
     unit_vector_from_plane_angles, GuidStateSnafu, GuidanceError, GuidanceLaw, GuidanceMode,
     GuidancePhysicsSnafu, NyxError, Orbit, Spacecraft, Vector3,
 };
-use crate::cosmic::eclipse::{EclipseLocator, EclipseState};
+use crate::cosmic::eclipse::EclipseLocator;
 pub use crate::md::objective::Objective;
 pub use crate::md::StateParameter;
 use crate::State;
@@ -32,14 +32,14 @@ use std::fmt;
 use std::sync::Arc;
 
 /// Ruggiero defines the closed loop guidance law from IEPC 2011-102
-#[derive(Copy, Clone, Default, Debug)]
+#[derive(Copy, Clone, Default)]
 pub struct Ruggiero {
     /// Stores the objectives
     pub objectives: [Option<Objective>; 5],
     /// Stores the minimum efficiency to correct a given orbital element, defaults to zero (i.e. always correct)
     pub ηthresholds: [f64; 5],
     /// If define, coast until vehicle is out of the provided eclipse state.
-    pub max_eclipse: Option<EclipseState>,
+    pub max_eclipse_prct: Option<f64>,
     init_state: Spacecraft,
 }
 
@@ -105,7 +105,7 @@ impl Ruggiero {
             objectives: objs,
             init_state: initial,
             ηthresholds: eff,
-            max_eclipse: None,
+            max_eclipse_prct: None,
         }))
     }
 
@@ -114,7 +114,7 @@ impl Ruggiero {
     pub fn from_max_eclipse(
         objectives: &[Objective],
         initial: Spacecraft,
-        max_eclipse: EclipseState,
+        max_eclipse: f64,
     ) -> Result<Arc<Self>, NyxError> {
         let mut objs: [Option<Objective>; 5] = [None, None, None, None, None];
         let eff: [f64; 5] = [0.0; 5];
@@ -151,13 +151,13 @@ impl Ruggiero {
             objectives: objs,
             init_state: initial,
             ηthresholds: eff,
-            max_eclipse: Some(max_eclipse),
+            max_eclipse_prct: Some(max_eclipse),
         }))
     }
 
     /// Sets the maximum eclipse during which we can thrust.
-    pub fn set_max_eclipse(&mut self, max_eclipse: EclipseState) {
-        self.max_eclipse = Some(max_eclipse);
+    pub fn set_max_eclipse(&mut self, max_eclipse: f64) {
+        self.max_eclipse_prct = Some(max_eclipse);
     }
 
     /// Returns the efficiency η ∈ [0; 1] of correcting a specific orbital element at the provided osculating orbit
@@ -270,8 +270,11 @@ impl fmt::Display for Ruggiero {
             .collect::<Vec<String>>();
         write!(
             f,
-            "Ruggiero Controller (max eclipse: {:?}): \n {}",
-            self.max_eclipse,
+            "Ruggiero Controller (max eclipse: {}): \n {}",
+            match self.max_eclipse_prct {
+                Some(eclp) => format!("{eclp}"),
+                None => "None".to_string(),
+            },
             obj_msg.join("\n")
         )
     }
@@ -422,11 +425,12 @@ impl GuidanceLaw for Ruggiero {
         if sc.mode() != GuidanceMode::Inhibit {
             if !self.achieved(sc).unwrap() {
                 // Check eclipse state if applicable.
-                if let Some(max_eclipse) = self.max_eclipse {
+                if let Some(max_eclipse) = self.max_eclipse_prct {
                     let locator = EclipseLocator::cislunar(almanac.clone());
                     if locator
                         .compute(sc.orbit, almanac)
                         .expect("cannot compute eclipse")
+                        .percentage
                         > max_eclipse
                     {
                         // Coast in eclipse
diff --git a/src/dynamics/solarpressure.rs b/src/dynamics/solarpressure.rs
index ca65761e..5e8a4ae6 100644
--- a/src/dynamics/solarpressure.rs
+++ b/src/dynamics/solarpressure.rs
@@ -135,7 +135,7 @@ impl ForceModel for SolarPressure {
             .context(DynamicsAlmanacSnafu {
                 action: "solar radiation pressure computation",
             })?
-            .into();
+            .factor();
 
         let r_sun_au = r_sun.norm() / AU;
         // in N/(m^2)
@@ -170,7 +170,7 @@ impl ForceModel for SolarPressure {
             .context(DynamicsAlmanacSnafu {
                 action: "solar radiation pressure computation",
             })?
-            .into();
+            .factor();
 
         let r_sun_au = norm(&r_sun_d) / AU;
         let inv_r_sun_au = OHyperdual::<f64, Const<9>>::from_real(1.0) / (r_sun_au);
diff --git a/src/dynamics/sph_harmonics.rs b/src/dynamics/sph_harmonics.rs
index 7d9c076d..df40c49c 100644
--- a/src/dynamics/sph_harmonics.rs
+++ b/src/dynamics/sph_harmonics.rs
@@ -250,7 +250,7 @@ impl AccelModel for Harmonics {
         // As discussed with Sai, if the Earth was spinning faster, would the acceleration due to the harmonics be any different?
         // No. Therefore, we do not need to account for the transport theorem here.
         let dcm = almanac
-            .rotate_from_to(self.compute_frame, osc.frame, osc.epoch)
+            .rotate(self.compute_frame, osc.frame, osc.epoch)
             .context(OrientationSnafu {
                 action: "transform state dcm",
             })
@@ -373,7 +373,7 @@ impl AccelModel for Harmonics {
         }
 
         let dcm = almanac
-            .rotate_from_to(self.compute_frame, osc.frame, osc.epoch)
+            .rotate(self.compute_frame, osc.frame, osc.epoch)
             .context(OrientationSnafu {
                 action: "transform state dcm",
             })
diff --git a/src/od/ground_station.rs b/src/od/ground_station.rs
index 1daeb146..3d8f753f 100644
--- a/src/od/ground_station.rs
+++ b/src/od/ground_station.rs
@@ -16,14 +16,13 @@
     along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
 
-use anise::astro::{AzElRange, PhysicsResult};
+use anise::astro::{Aberration, AzElRange, PhysicsResult};
 use anise::errors::AlmanacResult;
 use anise::prelude::{Almanac, Frame, Orbit};
 
 use super::msr::RangeDoppler;
 use super::noise::StochasticNoise;
-use super::{ODAlmanacSnafu, ODError, ODPlanetaryDataSnafu, ODTrajSnafu, TrackingDeviceSim};
-use crate::cosmic::eclipse::{line_of_sight, EclipseState};
+use super::{ODAlmanacSnafu, ODError, ODTrajSnafu, TrackingDeviceSim};
 use crate::errors::EventError;
 use crate::io::ConfigRepr;
 use crate::md::prelude::{Interpolatable, Traj};
@@ -159,8 +158,23 @@ impl GroundStation {
 
     /// Computes the azimuth and elevation of the provided object seen from this ground station, both in degrees.
     /// This is a shortcut to almanac.azimuth_elevation_range_sez.
-    pub fn azimuth_elevation_of(&self, rx: Orbit, almanac: &Almanac) -> AlmanacResult<AzElRange> {
-        almanac.azimuth_elevation_range_sez(rx, self.to_orbit(rx.epoch, almanac).unwrap())
+    pub fn azimuth_elevation_of(
+        &self,
+        rx: Orbit,
+        obstructing_body: Option<Frame>,
+        almanac: &Almanac,
+    ) -> AlmanacResult<AzElRange> {
+        let ab_corr = if self.light_time_correction {
+            Aberration::LT
+        } else {
+            Aberration::NONE
+        };
+        almanac.azimuth_elevation_range_sez(
+            rx,
+            self.to_orbit(rx.epoch, almanac).unwrap(),
+            obstructing_body,
+            ab_corr,
+        )
     }
 
     /// Return this ground station as an orbit in its current frame
@@ -234,16 +248,22 @@ impl TrackingDeviceSim<Spacecraft, RangeDoppler> for GroundStation {
                 let rx_0 = traj.at(epoch - integration_time).context(ODTrajSnafu)?;
                 let rx_1 = traj.at(epoch).context(ODTrajSnafu)?;
 
-                let aer_t0 =
-                    self.azimuth_elevation_of(rx_0.orbit, &almanac)
-                        .context(ODAlmanacSnafu {
-                            action: "computing AER",
-                        })?;
-                let aer_t1 =
-                    self.azimuth_elevation_of(rx_1.orbit, &almanac)
-                        .context(ODAlmanacSnafu {
-                            action: "computing AER",
-                        })?;
+                let obstructing_body = if !self.frame.ephem_origin_match(rx_0.frame()) {
+                    Some(rx_0.frame())
+                } else {
+                    None
+                };
+
+                let aer_t0 = self
+                    .azimuth_elevation_of(rx_0.orbit, obstructing_body, &almanac)
+                    .context(ODAlmanacSnafu {
+                        action: "computing AER",
+                    })?;
+                let aer_t1 = self
+                    .azimuth_elevation_of(rx_1.orbit, obstructing_body, &almanac)
+                    .context(ODAlmanacSnafu {
+                        action: "computing AER",
+                    })?;
 
                 if aer_t0.elevation_deg < self.elevation_mask_deg
                     || aer_t1.elevation_deg < self.elevation_mask_deg
@@ -255,29 +275,6 @@ impl TrackingDeviceSim<Spacecraft, RangeDoppler> for GroundStation {
                     return Ok(None);
                 }
 
-                // If the frame of the trajectory is different from that of the ground station, then check that there is no eclipse.
-                for rx in [rx_0, rx_1] {
-                    if !self.frame.ephem_origin_match(rx.frame()) {
-                        let observer = self.to_orbit(rx.orbit.epoch, &almanac).unwrap();
-                        if line_of_sight(
-                            observer,
-                            rx.orbit,
-                            almanac
-                                .frame_from_uid(rx.frame())
-                                .context(ODPlanetaryDataSnafu {
-                                    action: "computing line of sight",
-                                })?,
-                            &almanac,
-                        )
-                        .context(ODAlmanacSnafu {
-                            action: "computing line of sight",
-                        })? == EclipseState::Umbra
-                        {
-                            return Ok(None);
-                        }
-                    }
-                }
-
                 // Noises are computed at the midpoint of the integration time.
                 let (timestamp_noise_s, range_noise_km, doppler_noise_km_s) =
                     self.noises(epoch - integration_time * 0.5, rng)?;
@@ -308,32 +305,18 @@ impl TrackingDeviceSim<Spacecraft, RangeDoppler> for GroundStation {
         rng: Option<&mut Pcg64Mcg>,
         almanac: Arc<Almanac>,
     ) -> Result<Option<RangeDoppler>, ODError> {
+        let obstructing_body = if !self.frame.ephem_origin_match(rx.frame()) {
+            Some(rx.frame())
+        } else {
+            None
+        };
+
         let aer = self
-            .azimuth_elevation_of(rx.orbit, &almanac)
+            .azimuth_elevation_of(rx.orbit, obstructing_body, &almanac)
             .context(ODAlmanacSnafu {
                 action: "computing AER",
             })?;
 
-        if !self.frame.ephem_origin_match(rx.frame()) {
-            let observer = self.to_orbit(rx.orbit.epoch, &almanac).unwrap();
-            if line_of_sight(
-                observer,
-                rx.orbit,
-                almanac
-                    .frame_from_uid(rx.frame())
-                    .context(ODPlanetaryDataSnafu {
-                        action: "computing line of sight",
-                    })?,
-                &almanac,
-            )
-            .context(ODAlmanacSnafu {
-                action: "computing line of sight",
-            })? == EclipseState::Umbra
-            {
-                return Ok(None);
-            }
-        }
-
         if aer.elevation_deg >= self.elevation_mask_deg {
             // Only update the noises if the measurement is valid.
             let (timestamp_noise_s, range_noise_km, doppler_noise_km_s) =
diff --git a/tests/cosmic/eclipse.rs b/tests/cosmic/eclipse.rs
index fde889f4..e9cdb84c 100644
--- a/tests/cosmic/eclipse.rs
+++ b/tests/cosmic/eclipse.rs
@@ -1,8 +1,9 @@
 extern crate nyx_space as nyx;
 
+use anise::astro::Occultation;
 use anise::constants::celestial_objects::{JUPITER_BARYCENTER, SUN};
 use anise::constants::frames::SUN_J2000;
-use nyx::cosmic::eclipse::{EclipseLocator, EclipseState};
+use nyx::cosmic::eclipse::EclipseLocator;
 use nyx::cosmic::Orbit;
 use nyx::dynamics::orbital::OrbitalDynamics;
 use nyx::dynamics::SpacecraftDynamics;
@@ -52,14 +53,18 @@ fn leo_sun_earth_eclipses(almanac: Arc<Almanac>) {
     };
 
     // Receive the states on the main thread.
-    let mut prev_eclipse_state = EclipseState::Umbra;
+    let mut prev_eclipse_state: Option<Occultation> = None;
     let mut cnt_changes = 0;
     while let Ok(rx_state) = truth_rx.recv() {
         let new_eclipse_state = e_loc.compute(rx_state.orbit, almanac.clone()).unwrap();
-        if new_eclipse_state != prev_eclipse_state {
-            println!("{:.6} now in {:?}", rx_state.orbit.epoch, new_eclipse_state);
-            prev_eclipse_state = new_eclipse_state;
-            cnt_changes += 1;
+        if let Some(prev_state) = prev_eclipse_state {
+            if new_eclipse_state.percentage != prev_state.percentage {
+                println!("{:.6} now in {}", rx_state.orbit.epoch, new_eclipse_state);
+                prev_eclipse_state = Some(new_eclipse_state);
+                cnt_changes += 1;
+            }
+        } else {
+            prev_eclipse_state = Some(new_eclipse_state);
         }
     }
 
@@ -100,16 +105,20 @@ fn geo_sun_earth_eclipses(almanac: Arc<Almanac>) {
     };
 
     // Receive the states on the main thread.
-    let mut prev_eclipse_state = EclipseState::Umbra;
+    let mut prev_eclipse_state: Option<Occultation> = None;
     let mut cnt_changes = 0;
     while let Ok(rx_state) = truth_rx.recv() {
         let new_eclipse_state = e_loc.compute(rx_state.orbit, almanac.clone()).unwrap();
-        if new_eclipse_state != prev_eclipse_state {
-            println!("{:.6} now in {:?}", rx_state.orbit.epoch, new_eclipse_state);
-            prev_eclipse_state = new_eclipse_state;
-            cnt_changes += 1;
+        if let Some(prev_state) = prev_eclipse_state {
+            if new_eclipse_state.percentage != prev_state.percentage {
+                println!("{:.6} now in {}", rx_state.orbit.epoch, new_eclipse_state);
+                prev_eclipse_state = Some(new_eclipse_state);
+                cnt_changes += 1;
+            }
+        } else {
+            prev_eclipse_state = Some(new_eclipse_state);
         }
     }
 
-    assert_eq!(cnt_changes, 15, "wrong number of eclipse state changes");
+    assert_eq!(cnt_changes, 14, "wrong number of eclipse state changes");
 }
diff --git a/tests/lib.rs b/tests/lib.rs
index f41517b7..ceca0e65 100644
--- a/tests/lib.rs
+++ b/tests/lib.rs
@@ -56,7 +56,7 @@ pub fn test_almanac_gmat_arcd() -> Arc<Almanac> {
         uri: "http://public-data.nyxspace.com/anise/de438.bsp".to_string(),
         crc32: Some(1111895644),
     };
-    de438.process().unwrap();
+    de438.process(true).unwrap();
 
     let mut almanac = test_almanac();
     almanac = almanac.load(&de438.uri).unwrap();
diff --git a/tests/propagation/events.rs b/tests/propagation/events.rs
index 5507dab7..4de707b2 100644
--- a/tests/propagation/events.rs
+++ b/tests/propagation/events.rs
@@ -2,6 +2,7 @@ extern crate nyx_space as nyx;
 use std::{fmt::Write, sync::Arc};
 
 use anise::{
+    astro::Occultation,
     constants::frames::{EARTH_J2000, IAU_EARTH_FRAME, SUN_J2000},
     prelude::Almanac,
 };
@@ -15,7 +16,7 @@ fn almanac() -> Arc<Almanac> {
 
 #[rstest]
 fn event_tracker_true_anomaly(almanac: Arc<Almanac>) {
-    use nyx::cosmic::eclipse::{EclipseLocator, EclipseState};
+    use nyx::cosmic::eclipse::EclipseLocator;
     use nyx::md::prelude::*;
     use nyx::od::GroundStation;
 
@@ -64,7 +65,7 @@ fn event_tracker_true_anomaly(almanac: Arc<Almanac>) {
     };
 
     // Adding this print to confirm that the penumbra calculation continuously increases and then decreases.
-    let mut e_state = EclipseState::Umbra;
+    let mut e_state: Option<Occultation> = None;
     // Also see what is the max elevation of this spacecraft over the Grand Canyon
     let gc = GroundStation::from_point(
         "Grand Canyon".to_string(),
@@ -79,13 +80,19 @@ fn event_tracker_true_anomaly(almanac: Arc<Almanac>) {
     let mut max_dt = dt;
     for state in traj.every(10 * Unit::Second) {
         let new_e_state = e_loc.compute(state.orbit, almanac.clone()).unwrap();
-        if e_state != new_e_state {
-            println!("{:x}\t{}", state, new_e_state);
-            e_state = new_e_state;
+        if let Some(prev_state) = e_state {
+            if prev_state.percentage != new_e_state.percentage {
+                println!("{:x}\t{}", state, new_e_state);
+                e_state = Some(new_e_state);
+            }
+        } else {
+            e_state = Some(new_e_state);
         }
 
         // Compute the elevation
-        let aer = gc.azimuth_elevation_of(state.orbit, &almanac).unwrap();
+        let aer = gc
+            .azimuth_elevation_of(state.orbit, None, &almanac)
+            .unwrap();
         let elevation = aer.elevation_deg;
         if elevation > max_el {
             max_el = elevation;