Merge pull request #458 from nyx-space/fix-gil-deadlock

Fix Python GIL deadlock in `many` functions

Author: ChristopherRabotin

Cargo.toml

@@ -3,7 +3,7 @@ resolver = "2"
 members = ["anise", "anise-cli", "anise-gui", "anise-py", "anise/fuzz"]
 
 [workspace.package]
-version = "0.6.1"
+version = "0.6.2"
 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."
@@ -46,7 +46,7 @@ numpy = "0.25"
 ndarray = ">= 0.15, < 0.17"
 rayon = "1.10.0"
 
-anise = { version = "0.6.1", path = "anise", default-features = false }
+anise = { version = "0.6.2", path = "anise", default-features = false }
 
 [profile.bench]
 debug = true

anise/src/almanac/python.rs

@@ -143,25 +143,28 @@ impl Almanac {
     ))]
     fn py_azimuth_elevation_range_sez_many(
         &self,
+        py: Python,
         rx_tx_states: Vec<(CartesianState, CartesianState)>,
         obstructing_body: Option<Frame>,
         ab_corr: Option<Aberration>,
     ) -> Vec<AzElRange> {
-        let mut rslt = rx_tx_states
-            .par_iter()
-            .filter_map(|(rx, tx)| {
-                self.azimuth_elevation_range_sez(*rx, *tx, obstructing_body, ab_corr)
-                    .map_or_else(
-                        |e| {
-                            println!("{e}");
-                            None
-                        },
-                        |aer| Some(aer),
-                    )
-            })
-            .collect::<Vec<AzElRange>>();
-        rslt.sort_by(|aer_a, aer_b| aer_a.epoch.cmp(&aer_b.epoch));
-        rslt
+        py.allow_threads(|| {
+            let mut rslt = rx_tx_states
+                .par_iter()
+                .filter_map(|(rx, tx)| {
+                    self.azimuth_elevation_range_sez(*rx, *tx, obstructing_body, ab_corr)
+                        .map_or_else(
+                            |e| {
+                                println!("{e}");
+                                None
+                            },
+                            |aer| Some(aer),
+                        )
+                })
+                .collect::<Vec<AzElRange>>();
+            rslt.sort_by(|aer_a, aer_b| aer_a.epoch.cmp(&aer_b.epoch));
+            rslt
+        })
     }
 
     /// Computes whether the line of sight between an observer and an observed Cartesian state is obstructed by the obstructing body.
@@ -283,25 +286,28 @@ impl Almanac {
     ))]
     fn py_solar_eclipsing_many(
         &self,
+        py: Python,
         eclipsing_frame: Frame,
         observers: Vec<Orbit>,
         ab_corr: Option<Aberration>,
     ) -> Vec<Occultation> {
-        let mut rslt = observers
-            .par_iter()
-            .filter_map(|observer| {
-                self.solar_eclipsing(eclipsing_frame, *observer, ab_corr)
-                    .map_or_else(
-                        |e| {
-                            println!("{e}");
-                            None
-                        },
-                        |aer| Some(aer),
-                    )
-            })
-            .collect::<Vec<Occultation>>();
-        rslt.sort_by(|aer_a, aer_b| aer_a.epoch.cmp(&aer_b.epoch));
-        rslt
+        py.allow_threads(|| {
+            let mut rslt = observers
+                .par_iter()
+                .filter_map(|observer| {
+                    self.solar_eclipsing(eclipsing_frame, *observer, ab_corr)
+                        .map_or_else(
+                            |e| {
+                                println!("{e}");
+                                None
+                            },
+                            |aer| Some(aer),
+                        )
+                })
+                .collect::<Vec<Occultation>>();
+            rslt.sort_by(|aer_a, aer_b| aer_a.epoch.cmp(&aer_b.epoch));
+            rslt
+        })
     }
 
     /// Computes the Beta angle (β) for a given orbital state, in degrees. A Beta angle of 0° indicates that the orbit plane is edge-on to the Sun, leading to maximum eclipse time. Conversely, a Beta angle of +90° or -90° means the orbit plane is face-on to the Sun, resulting in continuous sunlight exposure and no eclipses.
@@ -377,26 +383,29 @@ impl Almanac {
     ))]
     fn py_transform_many<'py>(
         &self,
+        py: Python,
         target_frame: Frame,
         observer_frame: Frame,
         time_series: TimeSeries,
         ab_corr: Option<Aberration>,
     ) -> Vec<CartesianState> {
-        let mut states = time_series
-            .par_bridge()
-            .filter_map(|epoch| {
-                self.transform(target_frame, observer_frame, epoch, ab_corr)
-                    .map_or_else(
-                        |e| {
-                            eprintln!("{e}");
-                            None
-                        },
-                        |state| Some(state),
-                    )
-            })
-            .collect::<Vec<CartesianState>>();
-        states.sort_by(|state_a, state_b| state_a.epoch.cmp(&state_b.epoch));
-        states
+        py.allow_threads(|| {
+            let mut states = time_series
+                .par_bridge()
+                .filter_map(|epoch| {
+                    self.transform(target_frame, observer_frame, epoch, ab_corr)
+                        .map_or_else(
+                            |e| {
+                                eprintln!("{e}");
+                                None
+                            },
+                            |state| Some(state),
+                        )
+                })
+                .collect::<Vec<CartesianState>>();
+            states.sort_by(|state_a, state_b| state_a.epoch.cmp(&state_b.epoch));
+            states
+        })
     }
 
     /// Returns the provided state as seen from the observer frame, given the aberration.
@@ -436,25 +445,28 @@ impl Almanac {
     ))]
     fn py_transform_many_to(
         &self,
+        py: Python,
         states: Vec<CartesianState>,
         observer_frame: Frame,
         ab_corr: Option<Aberration>,
     ) -> Vec<CartesianState> {
-        let mut rslt = states
-            .par_iter()
-            .filter_map(|state| {
-                self.transform_to(*state, observer_frame, ab_corr)
-                    .map_or_else(
-                        |e| {
-                            println!("{e}");
-                            None
-                        },
-                        |state| Some(state),
-                    )
-            })
-            .collect::<Vec<CartesianState>>();
-        rslt.sort_by(|state_a, state_b| state_a.epoch.cmp(&state_b.epoch));
-        rslt
+        py.allow_threads(|| {
+            let mut rslt = states
+                .par_iter()
+                .filter_map(|state| {
+                    self.transform_to(*state, observer_frame, ab_corr)
+                        .map_or_else(
+                            |e| {
+                                println!("{e}");
+                                None
+                            },
+                            |state| Some(state),
+                        )
+                })
+                .collect::<Vec<CartesianState>>();
+            rslt.sort_by(|state_a, state_b| state_a.epoch.cmp(&state_b.epoch));
+            rslt
+        })
     }
 
     /// Returns the Cartesian state of the object as seen from the provided observer frame (essentially `spkezr`).

anise/src/astro/orbit.rs

@@ -13,7 +13,7 @@ use super::PhysicsResult;
 
 use crate::{
     errors::{
-        HyperbolicTrueAnomalySnafu, InfiniteValueSnafu, ParabolicEccentricitySnafu,
+        HyperbolicTrueAnomalySnafu, InfiniteValueSnafu, MathError, ParabolicEccentricitySnafu,
         ParabolicSemiParamSnafu, PhysicsError, RadiusSnafu, VelocitySnafu,
     },
     math::{
@@ -25,7 +25,7 @@ use crate::{
     prelude::{uuid_from_epoch, Frame},
     NaifId,
 };
-use core::f64::consts::PI;
+use core::f64::consts::{PI, TAU};
 
 use core::fmt;
 use hifitime::{Duration, Epoch, TimeUnits, Unit};
@@ -37,7 +37,7 @@ use pyo3::prelude::*;
 #[cfg(feature = "python")]
 use pyo3::types::PyType;
 
-/// If an orbit has an eccentricity below the following value, it is considered circular (only affects warning messages)
+/// If an orbit has an eccentricity below the following value, it is considered circular.
 pub const ECC_EPSILON: f64 = 1e-11;
 
 /// A helper type alias, but no assumptions are made on the underlying validity of the frame.
@@ -787,8 +787,7 @@ impl Orbit {
     ///
     /// :rtype: Duration
     pub fn period(&self) -> PhysicsResult<Duration> {
-        Ok(2.0
-            * PI
+        Ok(TAU
             * (self.sma_km()?.powi(3) / self.frame.mu_km3_s2()?)
                 .sqrt()
                 .seconds())
@@ -904,7 +903,7 @@ impl Orbit {
                 Ok(0.0)
             }
         } else if self.evec()?[2] < 0.0 {
-            Ok((2.0 * PI - aop).to_degrees())
+            Ok((TAU - aop).to_degrees())
         } else {
             Ok(aop.to_degrees())
         }
@@ -965,7 +964,7 @@ impl Orbit {
                 Ok(0.0)
             }
         } else if n[1] < 0.0 {
-            Ok((2.0 * PI - raan).to_degrees())
+            Ok((TAU - raan).to_degrees())
         } else {
             Ok(raan.to_degrees())
         }
@@ -1039,7 +1038,7 @@ impl Orbit {
                 Ok(0.0)
             }
         } else if self.radius_km.dot(&self.velocity_km_s) < 0.0 {
-            Ok((2.0 * PI - ta).to_degrees())
+            Ok((TAU - ta).to_degrees())
         } else {
             Ok(ta.to_degrees())
         }
@@ -1301,6 +1300,30 @@ impl Orbit {
         Ok(-self.frame.mu_km3_s2()? / self.sma_km()?)
     }
 
+    /// Returns the Longitude of the Ascending Node (LTAN), or an error of equatorial orbits
+    ///
+    /// :rtype: float
+    pub fn ltan_deg(&self) -> PhysicsResult<f64> {
+        let n = Vector3::new(0.0, 0.0, 1.0).cross(&self.hvec()?);
+
+        if n.norm_squared() < f64::EPSILON {
+            return Err(PhysicsError::AppliedMath {
+                source: MathError::DomainError {
+                    value: n.norm_squared(),
+                    msg: "orbit is equatorial",
+                },
+            });
+        }
+
+        let mut ltan = n.y.atan2(n.x);
+
+        if ltan < 0.0 {
+            ltan += TAU;
+        }
+
+        Ok(ltan.to_degrees())
+    }
+
     /// Returns the radius of periapse in kilometers for the provided turn angle of this hyperbolic orbit.
     /// Returns an error if the orbit is not hyperbolic.
     ///

anise/tests/astro/orbit.rs

@@ -855,3 +855,21 @@ fn gh_regression_340(almanac: Almanac) {
         assert!(orbit.at_epoch(epoch).is_ok(), "error on {epoch}");
     }
 }
+
+#[rstest]
+fn misc_verif(almanac: Almanac) {
+    let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap();
+
+    // Blue Ghost landing day!
+    let epoch = Epoch::from_gregorian_utc_at_noon(2025, 3, 2);
+
+    // LTAN
+    let prograde = Orbit::new(0.0, 7000.0, 0.0, -1.0, 0.0, 7.5, epoch, eme2k);
+    f64_eq!(prograde.ltan_deg().unwrap(), 90.0, "prograde LTAN");
+
+    let retrograde = Orbit::new(0.0, -7000.0, 0.0, -1.0, 0.0, 7.5, epoch, eme2k);
+    f64_eq!(retrograde.ltan_deg().unwrap(), 270.0, "retrograde LTAN");
+
+    let equatorial = Orbit::keplerian(7000.0, 1e-4, 1e-12, 270.0, 45.0, 76.0, epoch, eme2k);
+    assert!(equatorial.ltan_deg().is_err())
+}