feat: create basic particle-spring system

1 files changed, 167 insertions, 0 deletions

diff --git a/src/main.rs b/src/main.rs
new file mode 100644
index 0000000..ab2e15f
--- /dev/null
+++ b/src/main.rs
@@ -0,0 +1,167 @@
+use std::io::Write;
+use std::{fs::File, path::PathBuf};
+use cgmath::{Angle, ElementWise, InnerSpace, Matrix3, MetricSpace, Rad, Vector3, Zero};
+
+type Vector = Vector3<f32>;
+
+#[derive(Debug, Clone, Copy)]
+struct Particle {
+    position: Vector,
+    position_old: Vector,
+    // velocity: Vector,
+    acceleration: Vector,
+    mass: f32,
+    fixed: bool,
+}
+
+impl Particle {
+    pub fn new(position: Vector, mass: f32, fixed: bool) -> Self {
+        Self {
+            position,
+            position_old: position,
+            acceleration: Vector::zero(),
+            mass,
+            fixed,
+        }
+
+    }
+
+    fn apply_force(&mut self, force: Vector) {
+        self.acceleration += force / self.mass;
+    }
+
+    fn tick(&mut self, dt: f32) {
+        if !self.fixed {
+            // Verlet method
+            let velocity = self.position - self.position_old;
+            self.position_old = self.position;
+            self.position += velocity + self.acceleration * dt * dt;
+        }
+
+        self.acceleration = Vector::zero();
+    }
+}
+
+#[derive(Debug, Clone, Copy)]
+struct Spring {
+    rest_length: f32,
+    stiffness: f32,
+    particle_ids: [usize; 2],
+}
+
+#[derive(Debug, Clone)]
+struct Simulation {
+    particles: Vec<Particle>,
+    springs: Vec<Spring>,
+}
+
+impl Simulation {
+    fn tick(&mut self, dt: f32) {
+        for spring in self.springs.iter() {
+            let [a, b] = spring.particle_ids;
+            let dist_vec = self.particles[b].position - self.particles[a].position;
+            let force_magnitude = spring.stiffness * (dist_vec.magnitude() - spring.rest_length);
+            let normalized = dist_vec.normalize();
+
+            self.particles[a].apply_force(normalized.mul_element_wise(force_magnitude));
+            self.particles[b].apply_force(normalized.mul_element_wise(force_magnitude * -1.0));
+        }
+        for p in &mut self.particles {
+            p.tick(dt);
+        }
+    }
+}
+
+
+struct PPM {
+    prefix: PathBuf,
+    width: usize,
+    height: usize,
+    // pixels_per_unit: usize,
+}
+
+impl PPM {
+    fn render_particles(&self, particles: &Vec<Particle>) -> String {
+        let mut s = format!("P3\n{} {}\n255\n", self.width, self.height);
+        let white = "255 255 255 ";
+        let black = "0 0 0 ";
+
+        for pixel_row in 0..self.height {
+            for pixel_col in 0..self.width {
+                let point = Vector::new(pixel_col as f32, 0.0, (pixel_row as f32) * -1.0) + Vector::new(self.width as f32 / -2.0, 0.0, self.height as f32 / 2.0);
+                let color = match particles.iter().any(|p| p.position.distance(point) <= p.mass / 2.0) {
+                    true => black,
+                    false => white,
+                };
+                s += color;
+            }
+        }
+        s
+    }
+
+    fn save_frame(&self, particles: &Vec<Particle>, time: f32) {
+        let file_name = format!("frame-{:08.3}", time);
+        let path = self.prefix.join(file_name);
+        let mut file = File::create(path).unwrap();
+        let data = self.render_particles(particles);
+        file.write(data.as_bytes()).unwrap();
+    }
+}
+
+
+fn main() {
+    let rotation_matrix = Matrix3::from_axis_angle(Vector::unit_y(), -Rad::full_turn() / 12.0);
+
+    let mut simulation = Simulation {
+        particles: vec![
+            Particle::new(Vector::zero(), 5.0, true),
+            Particle::new(rotation_matrix * Vector::unit_x() * 40.0, 30.0, false),
+            Particle::new(rotation_matrix * Vector::unit_x() * 40.0 + rotation_matrix * rotation_matrix * Vector::unit_x() * 50.0, 15.0, false)
+        ],
+        springs: vec![
+            Spring {
+                particle_ids: [0, 1],
+                rest_length: 40.0,
+                stiffness: 400.0,
+            },
+            Spring {
+                particle_ids: [1, 2],
+                rest_length: 50.0,
+                stiffness: 5.0,
+            },
+        ],
+    };
+
+
+    let ppm = PPM {
+        prefix: PathBuf::from("./out"),
+        width: 150,
+        height: 800,
+    };
+
+    let gravity = Vector::new(0.0, 0.0, -9.8);
+
+    let time_step = 0.01;
+    let mut tick = 1;
+    let mut time = 0.0;
+    while time < 130.0 {
+        {
+            let p = &mut simulation.particles[0];
+            p.apply_force(gravity * -p.mass);
+        }
+        let mut all_particles = Vec::new();
+        for p in &mut simulation.particles {
+            p.apply_force(gravity * p.mass);
+        }
+        simulation.tick(time_step);
+        all_particles.append(&mut simulation.particles.clone());
+
+        if tick % 10 == 0 {
+            ppm.save_frame(&all_particles, time);
+        }
+
+
+        time += time_step;
+        tick +=1;
+    }
+}