Skip to main content

phyz_gravity/
particle.rs

1//! Gravity particle representation.
2
3use phyz_math::Vec3;
4
5/// A massive particle subject to gravitational forces.
6#[derive(Clone, Debug)]
7pub struct GravityParticle {
8    /// Position (m).
9    pub x: Vec3,
10    /// Velocity (m/s).
11    pub v: Vec3,
12    /// Force accumulator (N).
13    pub f: Vec3,
14    /// Mass (kg).
15    pub m: f64,
16}
17
18impl GravityParticle {
19    /// Create a new gravity particle.
20    pub fn new(x: Vec3, v: Vec3, m: f64) -> Self {
21        Self {
22            x,
23            v,
24            f: Vec3::zeros(),
25            m,
26        }
27    }
28
29    /// Reset force accumulator.
30    pub fn reset_force(&mut self) {
31        self.f = Vec3::zeros();
32    }
33
34    /// Add force to accumulator.
35    pub fn add_force(&mut self, f: Vec3) {
36        self.f += f;
37    }
38
39    /// Integrate using velocity Verlet (kick-drift-kick).
40    pub fn velocity_verlet_step(&mut self, dt: f64) {
41        // Half-step velocity update
42        self.v += self.f / self.m * (dt / 2.0);
43
44        // Full-step position update
45        self.x += self.v * dt;
46    }
47
48    /// Complete velocity Verlet after force recomputation.
49    pub fn velocity_verlet_complete(&mut self, dt: f64) {
50        // Second half-step velocity update
51        self.v += self.f / self.m * (dt / 2.0);
52    }
53
54    /// Kinetic energy: 0.5 * m * v².
55    pub fn kinetic_energy(&self) -> f64 {
56        0.5 * self.m * self.v.norm_squared()
57    }
58}
59
60#[cfg(test)]
61mod tests {
62    use super::*;
63
64    #[test]
65    fn test_particle_creation() {
66        let p = GravityParticle::new(Vec3::new(1.0, 2.0, 3.0), Vec3::zeros(), 1.0);
67        assert_eq!(p.x, Vec3::new(1.0, 2.0, 3.0));
68        assert_eq!(p.m, 1.0);
69    }
70
71    #[test]
72    fn test_kinetic_energy() {
73        let p = GravityParticle::new(Vec3::zeros(), Vec3::new(1.0, 0.0, 0.0), 2.0);
74        assert!((p.kinetic_energy() - 1.0).abs() < 1e-10);
75    }
76
77    #[test]
78    fn test_velocity_verlet() {
79        let mut p = GravityParticle::new(Vec3::zeros(), Vec3::new(1.0, 0.0, 0.0), 1.0);
80        p.f = Vec3::new(1.0, 0.0, 0.0); // 1 N force
81
82        let dt = 0.1;
83        p.velocity_verlet_step(dt);
84
85        // x = v0*t + 0.5*a*t^2 = 1.0*0.1 + 0.5*1.0*0.01 = 0.105
86        assert!((p.x.x - 0.105).abs() < 1e-10);
87    }
88}