feat[physics] :: improve numerical stability and refactor simulation loop

.github/workflows/pull_request_template.md

@@ -0,0 +1,41 @@
+## Summary
+<!-- Brief, high-level description of what this PR changes -->
+-
+
+<details>
+<summary>Why this change?</summary>
+
+<!-- Explain the motivation, bug, or problem this PR addresses.
+     Keep technical details here to avoid cluttering the summary. -->
+
+</details>
+
+---
+
+## What changed
+<!-- Optional: list concrete changes if the PR is non-trivial -->
+-
+
+---
+
+## Impact
+<!-- Optional: describe user-facing or system impact -->
+- [ ] Bug fix
+- [ ] Behavior change
+- [ ] Performance improvement
+- [ ] Refactor / cleanup
+- [ ] Documentation
+
+---
+
+## Testing
+<!-- Describe how this change was tested -->
+- [ ] Unit tests
+- [ ] Manual testing
+- [ ] Not tested (explain why)
+
+---
+
+## Notes for reviewers
+<!-- Optional: call out areas needing special attention -->
+-

.github/workflows/release.yml

@@ -1,3 +1,4 @@
+---
 name: Release
 
 on:

.github/workflows/rust.yml

@@ -125,7 +125,9 @@ jobs:
       - run: zip simulation-engine-${CARGO_RELEASE_TAG}-linux-x64.zip target/release/simulation_engine
       - run: |
           if ! gh release view $CARGO_RELEASE_TAG > /dev/null 2>&1; then
-            gh release create $CARGO_RELEASE_TAG simulation-engine-${CARGO_RELEASE_TAG}-linux-x64.zip --generate-notes --verify-tag
+            gh release create $CARGO_RELEASE_TAG \
+              simulation-engine-${CARGO_RELEASE_TAG}-linux-x64.zip \
+              --generate-notes --verify-tag
           fi
         env:
           GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}

.pre-commit-config.yaml

@@ -1,6 +1,6 @@
 repos:
 - repo: https://github.com/pre-commit/pre-commit-hooks
-  rev: v4.6.0
+  rev: v6.0.0
   hooks:
   - id: trailing-whitespace
   - id: end-of-file-fixer

engine/constraints/mod.rs

@@ -0,0 +1,2 @@
+/// Constraints module for physics engine constraints and limitations
+pub mod physics_constraints;

engine/constraints/physics_constraints.rs

@@ -0,0 +1,88 @@
+/// Physics constraints and boundary conditions
+/// Handles various types of constraints that can be applied to physics objects
+use nalgebra::Vector2;
+
+/// Represents different types of constraints that can be applied to physics objects
+#[derive(Debug, Clone)]
+pub enum Constraint {
+    /// Fixed position constraint - object cannot move from a specific point
+    FixedPosition(Vector2<f32>),
+    /// Distance constraint - maintains a fixed distance between two objects
+    Distance { target_distance: f32 },
+    /// Axis constraint - restricts movement to a specific axis
+    Axis { axis: Vector2<f32>, position: f32 },
+    /// Boundary constraint - keeps object within rectangular bounds
+    Boundary {
+        min: Vector2<f32>,
+        max: Vector2<f32>,
+    },
+}
+
+/// A constraint that can be applied to a physics object
+#[derive(Debug, Clone)]
+pub struct PhysicsConstraint {
+    pub constraint_type: Constraint,
+    pub stiffness: f32, // How strongly the constraint is enforced (0.0 to 1.0)
+    pub damping: f32,   // Damping factor for constraint forces
+}
+
+impl Default for PhysicsConstraint {
+    fn default() -> Self {
+        Self {
+            constraint_type: Constraint::Boundary {
+                min: Vector2::new(-100.0, -100.0),
+                max: Vector2::new(100.0, 100.0),
+            },
+            stiffness: 1.0,
+            damping: 0.1,
+        }
+    }
+}
+
+impl PhysicsConstraint {
+    /// Creates a new boundary constraint with the given bounds
+    pub fn boundary(min: Vector2<f32>, max: Vector2<f32>) -> Self {
+        Self {
+            constraint_type: Constraint::Boundary { min, max },
+            stiffness: 1.0,
+            damping: 0.1,
+        }
+    }
+
+    /// Applies the constraint to a position and velocity, returning the constrained values
+    pub fn apply(&self, position: &mut Vector2<f32>, velocity: &mut Vector2<f32>, delta_time: f32) {
+        match &self.constraint_type {
+            Constraint::Boundary { min, max } => {
+                self.apply_boundary_constraint(position, velocity, *min, *max, delta_time);
+            }
+            Constraint::FixedPosition(fixed_pos) => {
+                *position = *fixed_pos;
+                *velocity = Vector2::zeros();
+            }
+            _ => {
+                // Other constraint types not yet implemented
+            }
+        }
+    }
+
+    fn apply_boundary_constraint(
+        &self,
+        position: &mut Vector2<f32>,
+        velocity: &mut Vector2<f32>,
+        min: Vector2<f32>,
+        max: Vector2<f32>,
+        _delta_time: f32,
+    ) {
+        // Clamp position to bounds
+        position.x = position.x.max(min.x).min(max.x);
+        position.y = position.y.max(min.y).min(max.y);
+
+        // Apply damping when hitting boundaries
+        if position.x <= min.x || position.x >= max.x {
+            velocity.x *= 1.0 - self.damping;
+        }
+        if position.y <= min.y || position.y >= max.y {
+            velocity.y *= 1.0 - self.damping;
+        }
+    }
+}

engine/lib.rs

@@ -1,4 +1,6 @@
+pub mod constraints;
 pub mod managers;
+pub mod numerics;
 pub mod physics;
 
 #[cfg(feature = "python")]

engine/managers/simulation_manager.rs

@@ -26,7 +26,7 @@ impl Simulation {
             state: SimulationState::Stopped,
             time_step,
             duration,
-            physics_engine: PhysicsEngine::new(),
+            physics_engine: PhysicsEngine::new().expect("Failed to create physics engine"),
         }
     }
 
@@ -35,7 +35,9 @@ impl Simulation {
         println!("Simulation {} started", self.id);
         while self.state == SimulationState::Running {
             self.update();
-            self.physics_engine.simulate(self.time_step);
+            self.physics_engine
+                .simulate(self.time_step)
+                .expect("Simulation step failed");
             thread::sleep(Duration::from_millis((self.time_step * 1000.0) as u64));
         }
     }

engine/numerics/constants.rs

@@ -0,0 +1,21 @@
+//! Numerical constants for physics and general calculations
+//! Centralizes all magic numbers and thresholds used throughout the engine
+
+/// Minimum distance threshold to avoid division by very small values in collision calculations
+pub const MIN_DISTANCE_EPSILON: f32 = 1e-6;
+
+/// Maximum velocity magnitude to prevent numerical explosions
+pub const MAX_VELOCITY: f32 = 1000.0;
+
+/// Minimum mass threshold (objects below this are considered invalid)
+pub const MIN_MASS: f32 = 1e-3;
+
+/// Maximum mass threshold to prevent numerical issues
+pub const MAX_MASS: f32 = 1e6;
+
+/// Gravity constant (standard Earth gravity)
+pub const EARTH_GRAVITY: f32 = 9.81;
+
+/// Time step epsilon for floating point comparisons
+#[allow(dead_code)]
+pub const TIME_STEP_EPSILON: f32 = 1e-7;

engine/numerics/mod.rs

@@ -0,0 +1,7 @@
+/// Numerics module providing constants and utilities for numerical computations
+pub mod constants;
+pub mod utilities;
+
+// Re-export commonly used items for convenience
+pub use constants::*;
+pub use utilities::*;

engine/numerics/utilities.rs

@@ -0,0 +1,69 @@
+use crate::numerics::constants::*;
+/// Numerical utilities and helper functions
+/// Provides safe mathematical operations and validation functions
+use nalgebra::Vector2;
+
+/// Clamps a value between min and max
+#[inline]
+pub fn clamp<T: PartialOrd>(value: T, min: T, max: T) -> T {
+    if value < min {
+        min
+    } else if value > max {
+        max
+    } else {
+        value
+    }
+}
+
+/// Safely divides two floats, returning a default value if denominator is too small
+#[inline]
+pub fn safe_divide(numerator: f32, denominator: f32, default: f32) -> f32 {
+    if denominator.abs() < MIN_DISTANCE_EPSILON {
+        default
+    } else {
+        numerator / denominator
+    }
+}
+
+/// Clamps velocity components to prevent numerical instability
+#[inline]
+pub fn clamp_velocity(velocity: &mut Vector2<f32>) {
+    velocity.x = clamp(velocity.x, -MAX_VELOCITY, MAX_VELOCITY);
+    velocity.y = clamp(velocity.y, -MAX_VELOCITY, MAX_VELOCITY);
+}
+
+/// Validates that mass is within acceptable bounds
+#[inline]
+pub fn validate_mass(mass: f32) -> Result<f32, String> {
+    if mass < MIN_MASS {
+        Err(format!(
+            "Mass {} is below minimum threshold {}",
+            mass, MIN_MASS
+        ))
+    } else if mass > MAX_MASS {
+        Err(format!(
+            "Mass {} exceeds maximum threshold {}",
+            mass, MAX_MASS
+        ))
+    } else if !mass.is_finite() {
+        Err("Mass must be finite".to_string())
+    } else {
+        Ok(mass)
+    }
+}
+
+/// Validates that a vector contains only finite values
+#[inline]
+pub fn validate_vector(vec: &Vector2<f32>) -> Result<(), String> {
+    if !vec.x.is_finite() || !vec.y.is_finite() {
+        Err(format!("Vector contains non-finite values: {:?}", vec))
+    } else {
+        Ok(())
+    }
+}
+
+/// Checks if two floats are approximately equal within epsilon
+#[inline]
+pub fn approx_eq(a: f32, b: f32, epsilon: f32) -> bool {
+    (a - b).abs() < epsilon
+}

engine/physics/data_cleaner.rs

@@ -190,9 +190,13 @@ mod tests {
 
     #[test]
     fn test_cleanup_stats() {
-        let mut engine = PhysicsEngine::new();
-        engine.add_object(Vector2::new(0.0, 0.0), Vector2::new(1.0, 0.0), 1.0);
-        engine.add_object(Vector2::new(100.0, 100.0), Vector2::new(0.0, 0.0), 1.0); // Stationary, far away
+        let mut engine = PhysicsEngine::new().unwrap();
+        engine
+            .add_object(Vector2::new(0.0, 0.0), Vector2::new(1.0, 0.0), 1.0)
+            .unwrap();
+        engine
+            .add_object(Vector2::new(100.0, 100.0), Vector2::new(0.0, 0.0), 1.0)
+            .unwrap(); // Stationary, far away
 
         let cleaner = DataCleaner::new().with_cleanup_threshold(10);
         let stats = cleaner.get_stats(&engine);
@@ -206,9 +210,13 @@ mod tests {
 
     #[test]
     fn test_cleanup_inactive_objects() {
-        let mut engine = PhysicsEngine::new();
-        engine.add_object(Vector2::new(0.0, 0.0), Vector2::new(1.0, 0.0), 1.0); // Active
-        engine.add_object(Vector2::new(100.0, 100.0), Vector2::new(0.0, 0.0), 1.0); // Inactive, far away
+        let mut engine = PhysicsEngine::new().unwrap();
+        engine
+            .add_object(Vector2::new(0.0, 0.0), Vector2::new(1.0, 0.0), 1.0)
+            .unwrap(); // Active
+        engine
+            .add_object(Vector2::new(100.0, 100.0), Vector2::new(0.0, 0.0), 1.0)
+            .unwrap(); // Inactive, far away
 
         let cleaner = DataCleaner::new();
         let removed_count = cleaner.cleanup_inactive_objects(&mut engine, 0.0);
@@ -222,11 +230,13 @@ mod tests {
 
     #[test]
     fn test_needs_cleanup() {
-        let mut engine = PhysicsEngine::new();
+        let mut engine = PhysicsEngine::new().unwrap();
 
         // Add many objects to trigger cleanup threshold
         for i in 0..1500 {
-            engine.add_object(Vector2::new(i as f32, 0.0), Vector2::new(0.0, 0.0), 1.0);
+            engine
+                .add_object(Vector2::new(i as f32, 0.0), Vector2::new(0.0, 0.0), 1.0)
+                .unwrap();
         }
 
         let cleaner = DataCleaner::new().with_cleanup_threshold(1000);