Skip to main content

phyz_gravity/
constant.rs

1//! Constant uniform gravity solver (Layer 1).
2
3use crate::{GravityParticle, GravitySolver};
4use phyz_math::Vec3;
5
6/// Constant uniform gravity field.
7#[derive(Debug, Clone)]
8pub struct ConstantGravity {
9    /// Gravitational acceleration vector (m/s²).
10    pub g: Vec3,
11}
12
13impl ConstantGravity {
14    /// Create constant gravity pointing down (Earth-like).
15    pub fn new(g: Vec3) -> Self {
16        Self { g }
17    }
18
19    /// Standard Earth gravity (-9.81 m/s² in z).
20    pub fn earth() -> Self {
21        Self {
22            g: Vec3::new(0.0, 0.0, -crate::G_EARTH),
23        }
24    }
25}
26
27impl GravitySolver for ConstantGravity {
28    fn compute_forces(&mut self, particles: &mut [GravityParticle]) {
29        for p in particles.iter_mut() {
30            p.reset_force();
31            p.add_force(self.g * p.m);
32        }
33    }
34
35    fn potential_energy(&self, particles: &[GravityParticle]) -> f64 {
36        // U = -m * g · x (assuming g points down, so U = m*|g|*h)
37        particles.iter().map(|p| -p.m * self.g.dot(&p.x)).sum()
38    }
39}
40
41#[cfg(test)]
42mod tests {
43    use super::*;
44
45    #[test]
46    fn test_constant_gravity() {
47        let mut solver = ConstantGravity::earth();
48        let mut particles = vec![GravityParticle::new(
49            Vec3::new(0.0, 0.0, 10.0),
50            Vec3::zeros(),
51            1.0,
52        )];
53
54        solver.compute_forces(&mut particles);
55
56        // Force should be -9.81 N in z
57        assert!((particles[0].f.z + crate::G_EARTH).abs() < 1e-10);
58        assert!(particles[0].f.x.abs() < 1e-10);
59        assert!(particles[0].f.y.abs() < 1e-10);
60    }
61
62    #[test]
63    fn test_potential_energy() {
64        let solver = ConstantGravity::earth();
65        let particles = vec![GravityParticle::new(
66            Vec3::new(0.0, 0.0, 10.0),
67            Vec3::zeros(),
68            2.0,
69        )];
70
71        let u = solver.potential_energy(&particles);
72        // U = m*g*h = 2.0 * 9.80665 * 10.0 = 196.133
73        assert!((u - 196.133).abs() < 1e-3);
74    }
75}