feat: implement constraint system and PPM rendering

2 files changed, 28 insertions, 8 deletions

diff --git a/src/solver/midpoint.rs b/src/solver/midpoint.rs
index c5a01c4..08a3e3c 100644
--- a/src/solver/midpoint.rs
+++ b/src/solver/midpoint.rs
@@ -6,7 +6,7 @@ impl Solver for ParticleSystem {
         let mut state = self.collect_phase_space();
 
         // Shift to the midpoint
-        self.scatter_phase_space(&&PhaseSpace {
+        self.scatter_phase_space(&PhaseSpace {
             0: state.0.clone() + self.compute_derivative().0 * dt / 2.0,
         });
 
diff --git a/src/solver/mod.rs b/src/solver/mod.rs
index 4a5fec5..4bb8747 100644
--- a/src/solver/mod.rs
+++ b/src/solver/mod.rs
@@ -7,10 +7,10 @@ mod midpoint;
 #[derive(Debug, Clone)]
 pub struct PhaseSpace(DVector<Scalar>);
 type ParticleView<'a> = Matrix<
-    f32,
+    Scalar,
     Const<{ PhaseSpace::PARTICLE_DIM }>,
     Const<1>,
-    ViewStorage<'a, f32, Const<{ PhaseSpace::PARTICLE_DIM }>, Const<1>, Const<1>, Dyn>,
+    ViewStorage<'a, Scalar, Const<{ PhaseSpace::PARTICLE_DIM }>, Const<1>, Const<1>, Dyn>,
 >;
 
 impl PhaseSpace {
@@ -63,7 +63,6 @@ impl ParticleSystem {
     fn scatter_phase_space(&mut self, phase_space: &PhaseSpace) {
         for (particle_index, particle) in &mut self.particles.iter_mut().enumerate() {
             let view = phase_space.particle_view(particle_index);
-
             for i in 0..N {
                 particle.position[i] = view[i];
                 particle.velocity[i] = view[i + N];
@@ -85,6 +84,7 @@ mod tests {
     fn test_collect_phase_space() {
         let system = ParticleSystem {
             particles: vec![Particle::new(Point::new(2.0, 3.0), 1.0)],
+            constraints: vec![],
             t: 0.0,
         };
         let phase_space = system.collect_phase_space();
@@ -105,6 +105,7 @@ mod tests {
                 Particle::new(Point::new(0.0, 0.0), 1.0),
                 Particle::new(Point::new(0.0, 0.0), 1.0),
             ],
+            constraints: vec![],
             t: 0.0,
         };
 
@@ -120,11 +121,16 @@ mod tests {
         );
     }
 
-    fn simulate_falling_ball(fall_time: Scalar, dt: Scalar) -> (Vector, Vector) {
+    fn simulate_falling_ball(
+        fall_time: Scalar,
+        dt: Scalar,
+        mass: Scalar,
+    ) -> (Vector, Vector, ParticleSystem) {
         let gravity = -9.8 * Vector::y();
 
         let mut system = ParticleSystem {
-            particles: vec![Particle::new(Point::origin(), 1.0)],
+            particles: vec![Particle::new(Point::origin(), mass)],
+            constraints: vec![],
             t: 0.0,
         };
 
@@ -133,7 +139,7 @@ mod tests {
         for _ in 0..iterations {
             for particle in &mut system.particles {
                 particle.reset_force();
-                particle.apply_force(gravity);
+                particle.apply_force(gravity * particle.mass);
             }
             system.step(dt);
         }
@@ -144,12 +150,13 @@ mod tests {
         (
             system.particles[0].position.coords - expected_position,
             system.particles[0].velocity - expected_velocity,
+            system,
         )
     }
 
     #[test]
     fn ball_should_fall() {
-        let (position_error, velocity_error) = simulate_falling_ball(10.0, 0.01);
+        let (position_error, velocity_error, _) = simulate_falling_ball(10.0, 0.01, 3.0);
         assert!(
             position_error.norm() < 0.01,
             "Position error is too high: {}",
@@ -161,4 +168,17 @@ mod tests {
             velocity_error,
         );
     }
+
+    #[test]
+    fn freefall_different_masses() {
+        let (_, _, system1) = simulate_falling_ball(10.0, 0.01, 2.0);
+        let (_, _, system2) = simulate_falling_ball(10.0, 0.01, 10.0);
+
+        let diff = system1.particles[0].position - system2.particles[0].position;
+        assert!(
+            diff.norm() < 0.01,
+            "Points with different masses should fall with the same speed, diff: {}",
+            diff
+        );
+    }
 }