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> { 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> { 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> { 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> { 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> { ]; // 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> { 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> { 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> { 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> { // 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> { 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..05686a98 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 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 { - 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) -> AlmanacResult { - let mut state = EclipseState::Visibilis; + pub fn compute(&self, observer: Orbit, almanac: Arc) -> AlmanacResult { + 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,24 @@ impl fmt::Display for UmbraEvent { } impl EventEvaluator 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) -> Result { - match self + let occult = 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(); + + Ok((occult - 1.0).abs()) + // 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 +160,13 @@ impl fmt::Display for PenumbraEvent { impl EventEvaluator for PenumbraEvent { fn eval(&self, sc: &Spacecraft, almanac: Arc) -> Result { - match self + let occult = 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(); + + Ok((occult - 1.0).abs()) } /// Stop searching when the time has converged to less than 0.1 seconds @@ -266,337 +187,3 @@ impl EventEvaluator 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 { - 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 { - 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; 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, + pub max_eclipse_prct: Option, 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, NyxError> { let mut objs: [Option; 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::>(); 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..05139654 100644 --- a/src/dynamics/solarpressure.rs +++ b/src/dynamics/solarpressure.rs @@ -128,14 +128,17 @@ impl ForceModel for SolarPressure { let r_sun_unit = r_sun / r_sun.norm(); - // Compute the shaddowing factor. - let k: f64 = self + // ANISE returns the occultation percentage (or factor), which is the opposite as the illumination factor. + let occult = self .e_loc .compute(osc, almanac) .context(DynamicsAlmanacSnafu { action: "solar radiation pressure computation", })? - .into(); + .factor(); + + // Compute the illumination factor. + let k: f64 = (occult - 1.0).abs(); let r_sun_au = r_sun.norm() / AU; // in N/(m^2) @@ -163,14 +166,17 @@ impl ForceModel for SolarPressure { let r_sun_d: Vector3>> = hyperspace_from_vector(&r_sun); let r_sun_unit = r_sun_d / norm(&r_sun_d); - // Compute the shadowing factor. - let k: f64 = self + // ANISE returns the occultation percentage (or factor), which is the opposite as the illumination factor. + let occult = self .e_loc .compute(osc, almanac.clone()) .context(DynamicsAlmanacSnafu { action: "solar radiation pressure computation", })? - .into(); + .factor(); + + // Compute the illumination factor. + let k: f64 = (occult - 1.0).abs(); let r_sun_au = norm(&r_sun_d) / AU; let inv_r_sun_au = OHyperdual::>::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 . */ -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 { - almanac.azimuth_elevation_range_sez(rx, self.to_orbit(rx.epoch, almanac).unwrap()) + pub fn azimuth_elevation_of( + &self, + rx: Orbit, + obstructing_body: Option, + almanac: &Almanac, + ) -> AlmanacResult { + 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 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 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 for GroundStation { rng: Option<&mut Pcg64Mcg>, almanac: Arc, ) -> Result, 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) { }; // Receive the states on the main thread. - let mut prev_eclipse_state = EclipseState::Umbra; + let mut prev_eclipse_state: Option = 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) { }; // Receive the states on the main thread. - let mut prev_eclipse_state = EclipseState::Umbra; + let mut prev_eclipse_state: Option = 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 { 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/mission_design/force_models.rs b/tests/mission_design/force_models.rs index ab3ddde6..41fcf65c 100644 --- a/tests/mission_design/force_models.rs +++ b/tests/mission_design/force_models.rs @@ -67,8 +67,8 @@ fn srp_earth_full_vis(almanac: Arc) { err_r * 1e3, err_v * 1e3 ); - assert!(err_r < 5e-3, "position error too large for SRP"); - assert!(err_v < 1e-7, "velocity error too large for SRP"); + assert!(err_r < 5e-4, "position error too large for SRP"); + assert!(err_v < 9e-8, "velocity error too large for SRP"); } #[rstest] @@ -119,8 +119,8 @@ fn srp_earth_leo(almanac: Arc) { err_r * 1e3, err_v * 1e3 ); - assert!(err_r < 3e-1, "position error too large for SRP"); - assert!(err_v < 5e-4, "velocity error too large for SRP"); + assert!(err_r < 7e-3, "position error too large for SRP"); + assert!(err_v < 8e-6, "velocity error too large for SRP"); } #[rstest] @@ -173,8 +173,8 @@ fn srp_earth_meo_ecc_inc(almanac: Arc) { err_r * 1e3, err_v * 1e3 ); - assert!(err_r < 5e-2, "position error too large for SRP"); - assert!(err_v < 2e-5, "velocity error too large for SRP"); + assert!(err_r < 2e-3, "position error too large for SRP"); + assert!(err_v < 1e-6, "velocity error too large for SRP"); // Compare the case with the hyperdual EOMs (computation uses another part of the code) let mut prop = setup.with(sc.with_stm(), almanac); diff --git a/tests/propagation/events.rs b/tests/propagation/events.rs index 5507dab7..b48f0bad 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 { #[rstest] fn event_tracker_true_anomaly(almanac: Arc) { - use nyx::cosmic::eclipse::{EclipseLocator, EclipseState}; + use nyx::cosmic::eclipse::EclipseLocator; use nyx::md::prelude::*; use nyx::od::GroundStation; @@ -41,6 +42,10 @@ fn event_tracker_true_anomaly(almanac: Arc) { let mut prop = setup.with(state.into(), almanac.clone()); let (_, traj) = prop.for_duration_with_traj(prop_time).unwrap(); + println!("Initial state {state:x}"); + + println!("{traj}"); + // Find all of the events for e in &events { let found_events = traj.find(e, almanac.clone()).unwrap(); @@ -64,7 +69,7 @@ fn event_tracker_true_anomaly(almanac: Arc) { }; // 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 = 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 +84,19 @@ fn event_tracker_true_anomaly(almanac: Arc) { 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; @@ -106,6 +117,7 @@ fn event_tracker_true_anomaly(almanac: Arc) { let pretty = umbra_events .iter() + .skip(1) .fold(String::new(), |mut output, orbit_event| { let orbit = orbit_event.state.orbit; let _ = writeln!(