Skip to main content

sciforge_lib/maths/vector/
sim.rs

1use super::fields::{gravitational_force, lorentz_force};
2use super::types::{Particle, Vec3};
3
4pub struct VectorFieldSim {
5    pub particles: Vec<Particle>,
6    pub dt: f64,
7}
8
9impl VectorFieldSim {
10    pub fn new(dt: f64) -> Self {
11        Self {
12            particles: Vec::new(),
13            dt,
14        }
15    }
16
17    pub fn add_particle(&mut self, p: Particle) {
18        self.particles.push(p);
19    }
20
21    pub fn step_gravity(&mut self) {
22        let n = self.particles.len();
23        let mut forces = vec![Vec3::zero(); n];
24        for i in 0..n {
25            for j in (i + 1)..n {
26                let f = gravitational_force(&self.particles[i], &self.particles[j]);
27                forces[i] = &forces[i] + &f;
28                forces[j] = &forces[j] + &(-&f);
29            }
30        }
31        for (fi, pi) in forces.iter().zip(self.particles.iter_mut()) {
32            let accel = fi.scale(1.0 / pi.mass);
33            pi.velocity = &pi.velocity + &accel.scale(self.dt);
34            pi.position = &pi.position + &pi.velocity.scale(self.dt);
35        }
36    }
37
38    pub fn step_em(&mut self, e_field: &Vec3, b_field: &Vec3) {
39        for p in &mut self.particles {
40            let f = lorentz_force(p.charge, &p.velocity, e_field, b_field);
41            let accel = f.scale(1.0 / p.mass);
42            p.velocity = &p.velocity + &accel.scale(self.dt);
43            p.position = &p.position + &p.velocity.scale(self.dt);
44        }
45    }
46
47    pub fn total_kinetic_energy(&self) -> f64 {
48        self.particles.iter().map(|p| p.kinetic_energy()).sum()
49    }
50
51    pub fn total_momentum(&self) -> Vec3 {
52        self.particles
53            .iter()
54            .fold(Vec3::zero(), |acc, p| &acc + &p.momentum())
55    }
56
57    pub fn center_of_mass(&self) -> Vec3 {
58        let total_mass: f64 = self.particles.iter().map(|p| p.mass).sum();
59        if total_mass < 1e-30 {
60            return Vec3::zero();
61        }
62        let weighted = self
63            .particles
64            .iter()
65            .fold(Vec3::zero(), |acc, p| &acc + &p.position.scale(p.mass));
66        weighted.scale(1.0 / total_mass)
67    }
68}