Merge pull request #326 from nyx-space/gh-314-line-of-sight

Fixing line of sight computation + SPICE validation of penumbra calculator

Author: ChristopherRabotin

.github/workflows/rust.yml

@@ -134,6 +134,9 @@ jobs:
           ./target/debug/anise-cli inspect data/gmat-hermite.bsp
           ./target/debug/anise-cli inspect data/de440.bsp
 
+      - name: Rust-SPICE occultation validation
+        run: cargo test validate_gh_283_multi_barycenter_and_los --release --workspace --exclude anise-gui --exclude anise-py -- --nocapture --include-ignored
+
       - name: Rust-SPICE JPL DE validation
         run: RUST_BACKTRACE=1 cargo test validate_jplde --features spkezr_validation --release --workspace --exclude anise-gui --exclude anise-py -- --nocapture --include-ignored --test-threads 1
 
@@ -213,6 +216,7 @@ jobs:
           cargo llvm-cov clean --workspace
           cargo llvm-cov test --no-report -- --test-threads=1
           cargo llvm-cov test --no-report --tests -- compile_fail
+          cargo llvm-cov test --no-report validate_gh_283_multi_barycenter_and_los -- --nocapture --ignored
           cargo llvm-cov test --no-report validate_iau_rotation_to_parent -- --nocapture --ignored
           cargo llvm-cov test --no-report validate_bpc_to_iau_rotations -- --nocapture --ignored
           cargo llvm-cov test --no-report validate_jplde_de440s_no_aberration --features spkezr_validation -- --nocapture --ignored

Cargo.toml

@@ -3,7 +3,7 @@ resolver = "2"
 members = ["anise", "anise-cli", "anise-gui", "anise-py"]
 
 [workspace.package]
-version = "0.4.3"
+version = "0.4.4"
 edition = "2021"
 authors = ["Christopher Rabotin <christopher.rabotin@gmail.com>"]
 description = "ANISE provides a toolkit and files for Attitude, Navigation, Instrument, Spacecraft, and Ephemeris data. It's a modern replacement of NAIF SPICE file."
@@ -51,7 +51,7 @@ serde = "1"
 serde_derive = "1"
 serde_dhall = "0.12"
 
-anise = { version = "0.4.3", path = "anise", default-features = false }
+anise = { version = "0.4.4", path = "anise", default-features = false }
 
 [profile.bench]
 debug = true

anise/src/almanac/aer.rs

@@ -10,6 +10,7 @@
 
 use crate::{
     astro::{Aberration, AzElRange},
+    constants::SPEED_OF_LIGHT_KM_S,
     ephemerides::{EphemerisError, EphemerisPhysicsSnafu},
     errors::{AlmanacError, EphemerisSnafu, PhysicsError},
     frames::Frame,
@@ -21,6 +22,7 @@ use crate::{
 use super::Almanac;
 use crate::errors::AlmanacResult;
 
+use hifitime::TimeUnits;
 use log::warn;
 
 use snafu::ResultExt;
@@ -33,6 +35,9 @@ impl Almanac {
     /// Computes the azimuth (in degrees), elevation (in degrees), and range (in kilometers) of the
     /// receiver state (`rx`) seen from the transmitter state (`tx`), once converted into the SEZ frame of the transmitter.
     ///
+    /// # Warning
+    /// The obstructing body _should_ be a tri-axial ellipsoid body, e.g. IAU_MOON_FRAME.
+    ///
     /// # Algorithm
     /// 1. If any obstructing_bodies are provided, ensure that none of these are obstructing the line of sight between the receiver and transmitter.
     /// 2. Compute the SEZ (South East Zenith) frame of the transmitter.
@@ -120,6 +125,7 @@ impl Almanac {
             range_km: rho_sez.norm(),
             range_rate_km_s,
             obstructed_by,
+            light_time: (rho_sez.norm() / SPEED_OF_LIGHT_KM_S).seconds(),
         })
     }
 }
@@ -239,7 +245,43 @@ mod ut_aer {
             ),
         ];
 
-        for (sno, state) in states.iter().enumerate() {
+        for (sno, state) in states.iter().copied().enumerate() {
+            // Rebuild the ground station at this new epoch
+            let madrid = Orbit::try_latlongalt(
+                latitude_deg,
+                longitude_deg,
+                height_km,
+                MEAN_EARTH_ANGULAR_VELOCITY_DEG_S,
+                state.epoch,
+                iau_earth,
+            )
+            .unwrap();
+
+            let aer = almanac
+                .azimuth_elevation_range_sez(state, madrid, None, None)
+                .unwrap();
+
+            if sno == 0 {
+                assert_eq!(
+                    format!("{aer}"),
+                    format!(
+                        "{}: az.: 133.599990 deg    el.: 7.237568 deg    range: 91457.271742 km    range-rate: -12.396849 km/s    obstruction: none",
+                        state.epoch
+                    )
+                );
+            }
+
+            let expect = gmat_ranges_km[sno];
+
+            // This only checks that our computation isn't total garbage.
+            assert!((aer.range_km - expect).abs() < 5.0);
+        }
+
+        // Ensure that if the state are in another frame, the results are identical.
+
+        let states = states.map(|state| almanac.transform_to(state, EARTH_ITRF93, None).unwrap());
+
+        for (sno, state) in states.iter().copied().enumerate() {
             // Rebuild the ground station at this new epoch
             let madrid = Orbit::try_latlongalt(
                 latitude_deg,
@@ -252,7 +294,7 @@ mod ut_aer {
             .unwrap();
 
             let aer = almanac
-                .azimuth_elevation_range_sez(*state, madrid, None, None)
+                .azimuth_elevation_range_sez(state, madrid, None, None)
                 .unwrap();
 
             if sno == 0 {

anise/src/almanac/transform.rs

@@ -12,6 +12,7 @@ use hifitime::{Epoch, Unit as TimeUnit};
 use snafu::ResultExt;
 
 use crate::{
+    constants::orientations::J2000,
     errors::{AlmanacResult, EphemerisSnafu, OrientationSnafu},
     math::{cartesian::CartesianState, units::LengthUnit, Vector3},
     orientations::OrientationPhysicsSnafu,
@@ -71,25 +72,29 @@ impl Almanac {
     #[allow(clippy::too_many_arguments)]
     pub fn transform_to(
         &self,
-        state: CartesianState,
+        mut state: CartesianState,
         observer_frame: Frame,
         ab_corr: Option<Aberration>,
     ) -> AlmanacResult<CartesianState> {
-        let state = self
+        // If the input and final rotations differ, rotate into J2000 first
+        state = if state.frame.orient_origin_match(observer_frame) {
+            state
+        } else {
+            self.rotate_to(state, state.frame.with_orient(J2000))
+                .context(OrientationSnafu {
+                    action: "transform state dcm",
+                })?
+        };
+
+        // Transform in the base frame (J2000) or the common frame
+        state = self
             .translate_to(state, observer_frame, ab_corr)
             .context(EphemerisSnafu {
                 action: "transform state",
             })?;
 
-        // Compute the frame rotation
-        let dcm = self
-            .rotate(state.frame, observer_frame, state.epoch)
-            .context(OrientationSnafu {
-                action: "transform state dcm",
-            })?;
-
-        (dcm * state)
-            .context(OrientationPhysicsSnafu {})
+        // Rotate into the observer frame
+        self.rotate_to(state, observer_frame)
             .context(OrientationSnafu {
                 action: "transform state",
             })

anise/src/astro/mod.rs

@@ -13,7 +13,7 @@ use std::fmt::Display;
 use crate::errors::PhysicsError;
 use crate::frames::Frame;
 
-use hifitime::Epoch;
+use hifitime::{Duration, Epoch};
 
 #[cfg(feature = "python")]
 use pyo3::exceptions::PyTypeError;
@@ -47,6 +47,7 @@ pub struct AzElRange {
     pub range_km: f64,
     pub range_rate_km_s: f64,
     pub obstructed_by: Option<Frame>,
+    pub light_time: Duration,
 }
 
 #[cfg_attr(feature = "python", pymethods)]
@@ -56,6 +57,11 @@ impl AzElRange {
         self.azimuth_deg.is_finite() && self.elevation_deg.is_finite() && self.range_km > 1e-6
     }
 
+    /// Returns whether there is an obstruction.
+    pub const fn is_obstructed(&self) -> bool {
+        self.obstructed_by.is_some()
+    }
+
     /// Initializes a new AzElRange instance
     #[cfg(feature = "python")]
     #[new]
@@ -67,13 +73,17 @@ impl AzElRange {
         range_rate_km_s: f64,
         obstructed_by: Option<Frame>,
     ) -> Self {
+        use crate::constants::SPEED_OF_LIGHT_KM_S;
+        use hifitime::TimeUnits;
+
         Self {
             epoch,
             azimuth_deg,
             elevation_deg,
             range_km,
             range_rate_km_s,
             obstructed_by,
+            light_time: (range_km / SPEED_OF_LIGHT_KM_S).seconds(),
         }
     }
 

anise/src/astro/occultation.rs

@@ -19,7 +19,7 @@ use pyo3::prelude::*;
 
 /// Stores the result of an occultation computation with the occulation percentage
 /// Refer to the [MathSpec](https://nyxspace.com/nyxspace/MathSpec/celestial/eclipse/) for modeling details.
-#[derive(Copy, Clone, PartialEq)]
+#[derive(Copy, Clone, Debug, PartialEq)]
 #[cfg_attr(feature = "python", pyclass)]
 #[cfg_attr(feature = "python", pyo3(module = "anise"))]
 #[cfg_attr(feature = "python", pyo3(get_all, set_all))]

anise/src/astro/orbit_geodetic.rs

@@ -227,11 +227,19 @@ impl CartesianState {
         Ok((lat_deg, long_deg, alt_km))
     }
 
-    /// Returns the geodetic longitude (λ) in degrees. Value is between 0 and 360 degrees.
+    /// Returns the geodetic longitude (λ) in degrees. Value is between -180 and 180 degrees.
     ///
     /// # Frame warning
     /// This state MUST be in the body fixed frame (e.g. ITRF93) prior to calling this function, or the computation is **invalid**.
     pub fn longitude_deg(&self) -> f64 {
+        between_pm_180(self.radius_km.y.atan2(self.radius_km.x).to_degrees())
+    }
+
+    /// Returns the geodetic longitude (λ) in degrees. Value is between 0 and 360 degrees.
+    ///
+    /// # Frame warning
+    /// This state MUST be in the body fixed frame (e.g. ITRF93) prior to calling this function, or the computation is **invalid**.
+    pub fn longitude_360_deg(&self) -> f64 {
         between_0_360(self.radius_km.y.atan2(self.radius_km.x).to_degrees())
     }
 

anise/src/frames/frame.rs

@@ -11,9 +11,11 @@
 use core::fmt;
 use core::fmt::Debug;
 use serde_derive::{Deserialize, Serialize};
-use serde_dhall::StaticType;
 use snafu::ResultExt;
 
+#[cfg(feature = "metaload")]
+use serde_dhall::StaticType;
+
 use crate::astro::PhysicsResult;
 use crate::constants::celestial_objects::{
     celestial_name_from_id, id_to_celestial_name, SOLAR_SYSTEM_BARYCENTER,
@@ -32,7 +34,8 @@ use pyo3::prelude::*;
 use pyo3::pyclass::CompareOp;
 
 /// A Frame uniquely defined by its ephemeris center and orientation. Refer to FrameDetail for frames combined with parameters.
-#[derive(Copy, Clone, Debug, Serialize, Deserialize, PartialEq, StaticType)]
+#[derive(Copy, Clone, Debug, Serialize, Deserialize, PartialEq)]
+#[cfg_attr(feature = "metaload", derive(StaticType))]
 #[cfg_attr(feature = "python", pyclass)]
 #[cfg_attr(feature = "python", pyo3(get_all, set_all))]
 #[cfg_attr(feature = "python", pyo3(module = "anise.astro"))]
@@ -317,6 +320,7 @@ mod frame_ut {
         assert_eq!(format!("{EME2000:e}"), "Earth");
     }
 
+    #[cfg(feature = "metaload")]
     #[test]
     fn dhall_serde() {
         let serialized = serde_dhall::serialize(&EME2000)

anise/src/structure/planetocentric/ellipsoid.rs

@@ -11,6 +11,8 @@
 use core::fmt;
 use der::{Decode, Encode, Reader, Writer};
 use serde_derive::{Deserialize, Serialize};
+
+#[cfg(feature = "metaload")]
 use serde_dhall::StaticType;
 
 #[cfg(feature = "python")]
@@ -30,7 +32,8 @@ use pyo3::pyclass::CompareOp;
 /// Example: Radii of the Earth.
 ///
 ///    BODY399_RADII     = ( 6378.1366   6378.1366   6356.7519 )
-#[derive(Copy, Clone, Debug, Serialize, Deserialize, PartialEq, StaticType)]
+#[derive(Copy, Clone, Debug, Serialize, Deserialize, PartialEq)]
+#[cfg_attr(feature = "metaload", derive(StaticType))]
 #[cfg_attr(feature = "python", pyclass)]
 #[cfg_attr(feature = "python", pyo3(get_all, set_all))]
 #[cfg_attr(feature = "python", pyo3(module = "anise.astro"))]

anise/tests/astro/orbit.rs

@@ -642,7 +642,8 @@ fn verif_geodetic_vallado(almanac: Almanac) {
     f64_eq!(r.radius_km.z, rk, "r_k");
     let r = Orbit::from_position(ri, rj, rk, epoch, eme2k);
     f64_eq!(r.latitude_deg().unwrap(), lat_val, "latitude (φ)");
-    f64_eq!(r.longitude_deg(), long, "longitude (λ)");
+    f64_eq!(r.longitude_deg(), long - 360.0, "longitude (λ)");
+    f64_eq!(r.longitude_360_deg(), long, "longitude (λ)");
     f64_eq!(r.height_km().unwrap(), height_val, "height");
 
     // Check reciprocity near poles