Merge pull request #382 from nyx-space/fix/gh-381

Fix range-rate calculation in azimuth_elevation_range_sez

Author: ChristopherRabotin

.github/workflows/benchmarks.yml

@@ -49,7 +49,7 @@ jobs:
         run: cargo bench --bench "crit_planetary_data" --workspace --exclude anise-py
 
       - name: Save benchmark artifacts
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: jpl-development-ephemerides-benchmark
-          path: target/criterion/**/report/*
+          path: target/criterion/**/report/

.github/workflows/gui.yaml

@@ -26,7 +26,7 @@ jobs:
 
             - name: Install packages (Linux)
               if: runner.os == 'Linux'
-              uses: rerun-io/cache-apt-pkgs-action@59534850182063abf1b2c11bb3686722a12a8397
+              uses: awalsh128/cache-apt-pkgs-action@latest
               with:
                 packages: libxcb-render0-dev libxcb-shape0-dev libxcb-xfixes0-dev libxkbcommon-dev libssl-dev libgtk-3-dev # libgtk-3-dev is used by rfd
                 version: 1.0
@@ -36,7 +36,7 @@ jobs:
               run: cargo build --release --bin anise-gui --workspace --exclude anise-py
 
             - name: Save executable
-              uses: actions/upload-artifact@v3
+              uses: actions/upload-artifact@v4
               with:
                 name: anise-gui-linux
                 path: target/release/anise-gui
@@ -50,13 +50,13 @@ jobs:
             - uses: dtolnay/rust-toolchain@stable
 
             - name: Set up cargo cache
-              uses: Swatinem/rust-cache@v2
+              uses: Swatinem/rust-cache@v2.7.7
 
             - name: Build Windows executable
               run: cargo build --release --bin anise-gui --workspace --exclude anise-py
 
             - name: Save executable
-              uses: actions/upload-artifact@v3
+              uses: actions/upload-artifact@v4
               with:
                 name: anise-gui-windows
                 path: target\release\anise-gui.exe

.github/workflows/python.yml

@@ -59,9 +59,9 @@ jobs:
             fi
 
       - name: Upload wheels
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
-          name: wheels
+          name: wheels-linux-${{ matrix.target }}
           path: anise-py/dist
 
       - name: pytest
@@ -126,9 +126,9 @@ jobs:
           working-directory: anise-py
 
       - name: Upload wheels
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
-          name: wheels
+          name: wheels-windows-${{ matrix.target }}
           path: anise-py/dist
 
       - name: pytest
@@ -161,9 +161,9 @@ jobs:
           working-directory: anise-py
 
       - name: Upload wheels
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
-          name: wheels
+          name: wheels-macos-13
           path: anise-py/dist
 
       - name: pytest
@@ -196,9 +196,9 @@ jobs:
           working-directory: anise-py
 
       - name: Upload wheels
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
-          name: wheels
+          name: wheels-macos-14
           path: anise-py/dist
 
       - name: pytest
@@ -233,9 +233,9 @@ jobs:
           working-directory: anise-py
 
       - name: Upload sdist
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
-          name: wheels
+          name: wheels-sdist
           path: anise-py/dist
 
   release:
@@ -244,9 +244,7 @@ jobs:
     if: github.ref_type == 'tag'
     needs: [linux, windows, macos-13, macos-14, sdist]
     steps:
-      - uses: actions/download-artifact@v3
-        with:
-          name: wheels
+      - uses: actions/download-artifact@v4 # No `name` to download all artifacts.
 
       - name: Publish to PyPI
         uses: PyO3/maturin-action@v1

.github/workflows/rust.yml

@@ -166,7 +166,7 @@ jobs:
           python spk_validation_plots.py
 
       - name: Save validation artifacts
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: validation-artifacts
           path: target/*.html

Cargo.toml

@@ -3,7 +3,7 @@ resolver = "2"
 members = ["anise", "anise-cli", "anise-gui", "anise-py"]
 
 [workspace.package]
-version = "0.5.2"
+version = "0.5.3"
 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."
@@ -45,7 +45,7 @@ pyo3-log = "0.12"
 numpy = "0.23"
 ndarray = ">= 0.15, < 0.17"
 
-anise = { version = "0.5.2", path = "anise", default-features = false }
+anise = { version = "0.5.3", path = "anise", default-features = false }
 
 [profile.bench]
 debug = true

anise-py/tests/test_almanac.py

@@ -79,7 +79,7 @@ def test_state_transformation():
     topo_dcm = orig_state.dcm_from_topocentric_to_body_fixed(123)
     assert topo_dcm.get_state_dcm().shape == (6, 6)
     assert topo_dcm.rot_mat.shape == (3, 3)
-    assert topo_dcm.rot_mat_dt is None
+    assert (topo_dcm.rot_mat_dt is not None and topo_dcm.rot_mat_dt.shape == (3, 3)) or topo_dcm.rot_mat_dt is None
 
     # In Python, we can set the aberration to None
     aberration = None

anise/src/almanac/aer.rs

@@ -16,7 +16,6 @@ use crate::{
     frames::Frame,
     math::angles::{between_0_360, between_pm_180},
     prelude::Orbit,
-    time::uuid_from_epoch,
 };
 
 use super::Almanac;
@@ -81,10 +80,9 @@ impl Almanac {
         }
 
         // Compute the SEZ DCM
-        let from = uuid_from_epoch(tx.frame.orientation_id, rx.epoch);
         // SEZ DCM is topo to fixed
         let sez_dcm = tx
-            .dcm_from_topocentric_to_body_fixed(from)
+            .dcm_from_topocentric_to_body_fixed(-1)
             .context(EphemerisPhysicsSnafu { action: "" })
             .context(EphemerisSnafu {
                 action: "computing SEZ DCM for AER",
@@ -96,7 +94,7 @@ impl Almanac {
                 action: "transforming transmitter to SEZ",
             })?;
 
-        // Convert the receiver into the transmitter frame.
+        // Convert the receiver into the body fixed transmitter frame.
         let rx_in_tx_frame = self.transform_to(rx, tx.frame, ab_corr)?;
         // Convert into SEZ frame
         let rx_sez = (sez_dcm.transpose() * rx_in_tx_frame)
@@ -105,12 +103,17 @@ impl Almanac {
                 action: "transforming received to SEZ",
             })?;
 
-        // Compute the range ρ.
+        // Convert receiver into the transmitter frame
+        let rx_in_tx_frame = self.transform_to(rx, tx.frame, ab_corr)?;
+
+        // Compute the range ρ in the SEZ frame for az/el
         let rho_sez = rx_sez.radius_km - tx_sez.radius_km;
+        // And in the body-fixed transmitter frame for range and range-rate.
+        // While the norms of these vectors are identical, we need the exact vectors themselves for the range rate calculation.
+        let rho_tx_frame = rx_in_tx_frame.radius_km - tx.radius_km;
 
-        // Compute the range-rate \dot ρ
-        let range_rate_km_s =
-            rho_sez.dot(&(rx_sez.velocity_km_s - tx_sez.velocity_km_s)) / rho_sez.norm();
+        // Compute the range-rate \dot ρ. Note that rx_in_tx_frame is already the relative velocity of rx wrt tx!
+        let range_rate_km_s = rho_tx_frame.dot(&rx_in_tx_frame.velocity_km_s) / rho_tx_frame.norm();
 
         // Finally, compute the elevation (math is the same as declination)
         // Source: Vallado, section 4.4.3
@@ -270,7 +273,7 @@ mod ut_aer {
                 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",
+                        "{}: az.: 133.599990 deg    el.: 7.237568 deg    range: 91457.271742 km    range-rate: 2.198786 km/s    obstruction: none",
                         state.epoch
                     )
                 );
@@ -306,7 +309,7 @@ mod ut_aer {
                 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",
+                        "{}: az.: 133.599990 deg    el.: 7.237568 deg    range: 91457.271742 km    range-rate: 2.198786 km/s    obstruction: none",
                         state.epoch
                     )
                 );

anise/src/astro/orbit.rs

@@ -320,20 +320,52 @@ impl Orbit {
 #[allow(clippy::too_many_arguments)]
 #[cfg_attr(feature = "python", pymethods)]
 impl Orbit {
+    /// Builds the rotation matrix that rotates from the topocentric frame (SEZ) into the body fixed frame of this state.
+    ///
+    /// # Frame warnings
+    /// + If the state is NOT in a body fixed frame (i.e. ITRF93), then this computation is INVALID.
+    /// + (Usually) no time derivative can be computed: the orbit is expected to be a body fixed frame where the `at_epoch` function will fail. Exceptions for Moon body fixed frames.
+    ///
+    /// # UNUSED Arguments
+    /// + `from`: ID of this new frame. Only used to set the "from" frame of the DCM. -- No longer used since 0.5.3
+    ///
+    /// # Source
+    /// From the GMAT MathSpec, page 30 section 2.6.9 and from `Calculate_RFT` in `TopocentricAxes.cpp`, this returns the
+    /// rotation matrix from the topocentric frame (SEZ) to body fixed frame.
+    /// In the GMAT MathSpec notation, R_{IF} is the DCM from body fixed to inertial. Similarly, R{FT} is from topocentric
+    /// to body fixed.
+    pub fn dcm_from_topocentric_to_body_fixed(&self, _from: NaifId) -> PhysicsResult<DCM> {
+        let rot_mat_dt = if let Ok(pre) = self.at_epoch(self.epoch - Unit::Second * 1) {
+            if let Ok(post) = self.at_epoch(self.epoch + Unit::Second * 1) {
+                let dcm_pre = pre.dcm3x3_from_topocentric_to_body_fixed()?;
+                let dcm_post = post.dcm3x3_from_topocentric_to_body_fixed()?;
+                Some(0.5 * dcm_post.rot_mat - 0.5 * dcm_pre.rot_mat)
+            } else {
+                None
+            }
+        } else {
+            None
+        };
+
+        Ok(DCM {
+            rot_mat: self.dcm3x3_from_topocentric_to_body_fixed()?.rot_mat,
+            rot_mat_dt,
+            from: uuid_from_epoch(self.frame.orientation_id, self.epoch),
+            to: self.frame.orientation_id,
+        })
+    }
+
     /// Builds the rotation matrix that rotates from the topocentric frame (SEZ) into the body fixed frame of this state.
     ///
     /// # Frame warning
     /// If the state is NOT in a body fixed frame (i.e. ITRF93), then this computation is INVALID.
     ///
-    /// # Arguments
-    /// + `from`: ID of this new frame, must be unique if it'll be added to the Almanac. Only used to set the "from" frame of the DCM.
-    ///
     /// # Source
     /// From the GMAT MathSpec, page 30 section 2.6.9 and from `Calculate_RFT` in `TopocentricAxes.cpp`, this returns the
     /// rotation matrix from the topocentric frame (SEZ) to body fixed frame.
     /// In the GMAT MathSpec notation, R_{IF} is the DCM from body fixed to inertial. Similarly, R{FT} is from topocentric
     /// to body fixed.
-    pub fn dcm_from_topocentric_to_body_fixed(&self, from: NaifId) -> PhysicsResult<DCM> {
+    pub fn dcm3x3_from_topocentric_to_body_fixed(&self) -> PhysicsResult<DCM> {
         if (self.radius_km.x.powi(2) + self.radius_km.y.powi(2)).sqrt() < 1e-3 {
             warn!("SEZ frame ill-defined when close to the poles");
         }
@@ -357,7 +389,7 @@ impl Orbit {
         Ok(DCM {
             rot_mat,
             rot_mat_dt: None,
-            from,
+            from: uuid_from_epoch(self.frame.orientation_id, self.epoch),
             to: self.frame.orientation_id,
         })
     }

anise/src/frames/frame.rs

@@ -91,6 +91,15 @@ impl Frame {
         self.shape = Some(shape);
         self
     }
+
+    /// Returns a copy of this frame with the graviational parameter and the shape information from this frame.
+    /// Use this to prevent astrodynamical computations.
+    ///
+    /// :rtype: None
+    pub fn stripped(mut self) -> Self {
+        self.strip();
+        self
+    }
 }
 
 #[cfg(feature = "python")]