fix(ode): validate ExponentialMap::integrate_attitude's timestep and …

CHANGELOG.md

@@ -7,6 +7,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Fixed
+
+- **`ExponentialMap::integrate_attitude` input validation.** A non-finite timestep, a
+non-positive timestep, or a non-finite rate from the caller's `angular_rate_at` callback
+used to silently produce a NaN orientation or integrate backwards without comment.
+`integrate_attitude` now returns the new `IntegrateError::NonPositiveTimestep` or the
+existing `IntegrateError::NonFinite` instead. `attitude_step` and
+`attitude_step_with_angular_acceleration` stay infallible, as they sit on
+`RigidBody::stepped`'s panic-free per-tick path; their behavior with non-finite or
+negative input is now documented instead of silent. (#302)
+
 ## [0.10.0] - 2026-08-09
 
 A feature release adding signal processing, polynomials and minimum-snap trajectories, LQR and

crates/multicalc/src/error.rs

@@ -65,6 +65,8 @@ pub enum IntegrateError {
     NonFinite,
     /// A variable index was `>=` the number of variables in the point.
     IndexOutOfRange,
+    /// A fixed integration timestep was zero or negative.
+    NonPositiveTimestep,
 }
 
 /// Errors from the solver modules (root finding, Gauss-Newton, Levenberg-Marquardt).
@@ -532,6 +534,9 @@ impl core::fmt::Display for IntegrateError {
                 f.write_str("integrand or state contained a non-finite value")
             }
             IntegrateError::IndexOutOfRange => f.write_str("variable index out of range"),
+            IntegrateError::NonPositiveTimestep => {
+            f.write_str("integration timestep must be strictly positive")
+            }
         }
     }
 }

crates/multicalc/src/ode/exponential_map.rs

@@ -1,5 +1,6 @@
 //! Turning an orientation forward in time by the turn it makes over the step.
 
+use crate::error::IntegrateError;
 use crate::linear_algebra::Vector3D;
 use crate::scalar::Numeric;
 use crate::spatial::SO3;
@@ -22,6 +23,11 @@ impl ExponentialMap {
     /// the fourth power. Use [`ExponentialMap::attitude_step_with_angular_acceleration`] when the
     /// rate is changing and the extra accuracy is wanted for the same one exponential.
     ///
+    /// Behavior: this is infallible and does not validate its input. A non-finite `angular_rate`
+    /// or a non-finite/negative `dt` produces a non-finite result rather than an error; callers
+    /// that need validated input can use [`ExponentialMap::integrate_attitude`] instead, or
+    /// validate upstream - this stays a raw, panic-free primitive for hot per-tick call sites.
+    ///
     /// ```
     /// use multicalc::ode::ExponentialMap;
     /// use multicalc::spatial::SO3;
@@ -59,6 +65,9 @@ impl ExponentialMap {
     /// the error with the square of the step size instead of in proportion to it. The result is
     /// still a true rotation to within rounding.
     ///
+    /// Behavior: this is infallible and does not validate its input; see the "Behavior" note on
+    /// [`ExponentialMap::attitude_step`] for how non-finite or negative input is handled.
+    ///
     /// ```
     /// use multicalc::ode::ExponentialMap;
     /// use multicalc::spatial::SO3;
@@ -95,6 +104,12 @@ impl ExponentialMap {
     /// of the step size, without the caller having to work out how fast the turn rate is changing.
     /// The orientation stays a true rotation to within rounding the whole way through.
     ///
+    /// # Errors
+    ///
+    /// [`IntegrateError::NonPositiveTimestep`] if `dt` is not strictly positive, or
+    /// [`IntegrateError::NonFinite`] if `dt` or a rate returned by `angular_rate_at` is not
+    /// finite.
+    ///
     /// ```
     /// use multicalc::ode::ExponentialMap;
     /// use multicalc::spatial::SO3;
@@ -107,7 +122,7 @@ impl ExponentialMap {
     /// let mut nodes = 0;
     /// let facing = ExponentialMap::integrate_attitude(
     ///     &steady, 0.0, SO3::identity(), timestep, steps, |_time, _orientation| nodes += 1,
-    /// );
+    /// ).unwrap();
     /// assert_eq!(nodes, steps + 1);
     ///
     /// let swung: Vector3D<f64> = facing.act(Vector::new([1.0, 0.0, 0.0]));
@@ -121,24 +136,36 @@ impl ExponentialMap {
         dt: T,
         steps: usize,
         mut observer: O,
-    ) -> SO3<T>
+    ) -> Result<SO3<T>, IntegrateError>
     where
         T: Numeric,
         F: Fn(T, SO3<T>) -> Vector3D<T>,
         O: FnMut(T, SO3<T>),
     {
+        if !dt.is_finite() {
+            return Err(IntegrateError::NonFinite);
+        }
+        if dt <= T::ZERO {
+            return Err(IntegrateError::NonPositiveTimestep);
+        }
         let half = dt * T::HALF;
         let mut time = t0;
         let mut orientation = start_orientation;
         observer(time, orientation);
         for _ in 0..steps {
             let rate_at_start = angular_rate_at(time, orientation);
+            if !rate_at_start.is_finite() {
+                return Err(IntegrateError::NonFinite);
+            }
             let half_way = Self::attitude_step(orientation, rate_at_start, half);
             let half_way_rate = angular_rate_at(time + half, half_way);
+            if !half_way_rate.is_finite() {
+                return Err(IntegrateError::NonFinite);
+            }
             orientation = Self::attitude_step(orientation, half_way_rate, dt);
             time += dt;
             observer(time, orientation);
         }
-        orientation
+        Ok(orientation)
     }
 }

crates/multicalc/tests/suite/ode/exponential_map.rs

@@ -1,5 +1,6 @@
 use multicalc::Dual;
 use multicalc::linear_algebra::{Vector, Vector3D};
+use multicalc::error::IntegrateError;
 use multicalc::ode::ExponentialMap;
 use multicalc::spatial::SO3;
 
@@ -23,7 +24,7 @@ fn reference_orientation(steps: usize, final_time: f64) -> SO3<f64> {
         final_time / steps as f64,
         steps,
         |_, _| {},
-    )
+    ).unwrap()
 }
 
 // How far apart two orientations are, in radians.
@@ -63,7 +64,7 @@ fn unit_length_holds_over_a_long_run() {
     // Twenty seconds of tumbling at a tenth of a millisecond a step: the length has to stay put.
     let rate = |time: f64, _orientation: SO3<f64>| prescribed_rate(time);
     let facing =
-        ExponentialMap::integrate_attitude(&rate, 0.0, SO3::identity(), 1e-4, 200_000, |_, _| {});
+        ExponentialMap::integrate_attitude(&rate, 0.0, SO3::identity(), 1e-4, 200_000, |_, _| {}).unwrap();
     assert!((facing.quaternion().norm() - 1.0).abs() < 1e-12);
 }
 
@@ -120,7 +121,7 @@ fn integrate_attitude_converges_second_order() {
             1.0 / steps as f64,
             steps,
             |_, _| {},
-        );
+        ).unwrap();
         angle_between(facing, reference)
     };
     let ratio = endpoint_error(200) / endpoint_error(400);
@@ -141,7 +142,7 @@ fn observer_sees_every_node_starting_with_the_first() {
         timestep,
         steps,
         |time, orientation| nodes.push((time, orientation.log()[2])),
-    );
+    ).unwrap();
 
     assert_eq!(nodes.len(), steps + 1);
     assert_eq!(nodes[0].0, 0.0);
@@ -195,10 +196,32 @@ fn f32_holds_unit_length_and_round_trips() {
         1e-4_f32,
         100_000,
         |_, _| {},
-    );
+    ).unwrap();
     assert!((facing.quaternion().norm() - 1.0).abs() < 1e-3);
 
     let turn = Vector::new([0.2_f32, -0.1, 0.4]);
     let round_tripped = SO3::<f32>::exp(turn).log();
     assert!((round_tripped - turn).norm() < 1e-5);
 }
+
+#[test]
+fn integrate_attitude_rejects_non_positive_timestep() {
+    let rate = |_time: f64, _orientation: SO3<f64>| Vector::new([0.0, 0.0, 1.0]);
+    assert_eq!(
+        ExponentialMap::integrate_attitude(&rate, 0.0, SO3::identity(), 0.0, 10, |_, _| {}),
+        Err(IntegrateError::NonPositiveTimestep)
+    );
+    assert_eq!(
+        ExponentialMap::integrate_attitude(&rate, 0.0, SO3::identity(), -0.1, 10, |_, _| {}),
+        Err(IntegrateError::NonPositiveTimestep)
+    );
+}
+
+#[test]
+fn integrate_attitude_rejects_non_finite_rate() {
+    let rate = |_time: f64, _orientation: SO3<f64>| Vector::new([f64::NAN, 0.0, 0.0]);
+    assert_eq!(
+        ExponentialMap::integrate_attitude(&rate, 0.0, SO3::identity(), 0.01, 10, |_, _| {}),
+        Err(IntegrateError::NonFinite)
+    );
+}