qualia_core_db/specialized_libs/physics_simulation/
nbody.rs

1use super::*;
2
3impl PhysicsSimulationLibrary {
4    /// Astrophysics — Newtonian N-body gravitation in 2D by direct force summation.
5    ///
6    /// `positions` and `velocities` are flat `[x0,y0,x1,y1,…]` (length `2·N`), `masses`
7    /// length `N`. Accelerations `aᵢ = Σⱼ G·mⱼ·(rⱼ−rᵢ)/(|rⱼ−rᵢ|²+ε²)^{3/2}` are assembled
8    /// here; the time integration is `integrate_dopri5`. Total energy and angular momentum
9    /// are reported for conservation checks.
10    pub fn run_nbody_gravitation(
11        &self,
12        masses: Vec<f64>,
13        positions: Vec<f64>,
14        velocities: Vec<f64>,
15        g: f64,
16        softening: f64,
17        total_time: f64,
18        num_samples: usize,
19    ) -> Result<NBodyResult, PhysicsError> {
20        let n = masses.len();
21        if n == 0 || positions.len() != 2 * n || velocities.len() != 2 * n {
22            return Err(PhysicsError::InvalidConfiguration(
23                "masses length N, positions/velocities length 2N required".to_string(),
24            ));
25        }
26        let eps2 = softening * softening;
27        let masses_e = masses.clone();
28        // Layout: [pos(2N), vel(2N)].
29        let mut state = Vec::with_capacity(4 * n);
30        state.extend_from_slice(&positions);
31        state.extend_from_slice(&velocities);
32
33        let energy = |st: &[f64]| -> (f64, f64) {
34            let mut ke = 0.0;
35            let mut pe = 0.0;
36            let mut angmom = 0.0;
37            for i in 0..n {
38                let (vx, vy) = (st[2 * n + 2 * i], st[2 * n + 2 * i + 1]);
39                ke += 0.5 * masses[i] * (vx * vx + vy * vy);
40                let (x, y) = (st[2 * i], st[2 * i + 1]);
41                angmom += masses[i] * (x * vy - y * vx);
42                for j in (i + 1)..n {
43                    let dx = st[2 * j] - x;
44                    let dy = st[2 * j + 1] - y;
45                    let r = (dx * dx + dy * dy + eps2).sqrt();
46                    pe -= g * masses[i] * masses[j] / r;
47                }
48            }
49            (ke + pe, angmom)
50        };
51        let (energy_initial, angmom_initial) = energy(&state);
52
53        let deriv = move |_t: f64, y: &[f64], dy: &mut [f64]| -> Result<(), OdeError> {
54            // Positions' derivative = velocities.
55            for i in 0..(2 * n) {
56                dy[i] = y[2 * n + i];
57            }
58            // Velocities' derivative = accelerations (direct sum).
59            for i in 0..n {
60                let (xi, yi) = (y[2 * i], y[2 * i + 1]);
61                let mut ax = 0.0;
62                let mut ay = 0.0;
63                for j in 0..n {
64                    if j == i {
65                        continue;
66                    }
67                    let dx = y[2 * j] - xi;
68                    let dyj = y[2 * j + 1] - yi;
69                    let r2 = dx * dx + dyj * dyj + eps2;
70                    let inv_r3 = 1.0 / (r2 * r2.sqrt());
71                    ax += g * masses_e[j] * dx * inv_r3;
72                    ay += g * masses_e[j] * dyj * inv_r3;
73                }
74                dy[2 * n + 2 * i] = ax;
75                dy[2 * n + 2 * i + 1] = ay;
76            }
77            Ok(())
78        };
79        let (final_state, snapshots, accepted, rejected) =
80            self.integrate_ode_samples(state, total_time, num_samples, deriv)?;
81        let (energy_final, angmom_final) = energy(&final_state);
82        let position_snapshots: Vec<Vec<f64>> =
83            snapshots.iter().map(|s| s[..2 * n].to_vec()).collect();
84        let n_pts = snapshots.len();
85        let times: Vec<f64> = (0..n_pts)
86            .map(|k| total_time * k as f64 / (n_pts - 1).max(1) as f64)
87            .collect();
88        let energy_drift_rel = if energy_initial.abs() > f64::MIN_POSITIVE {
89            (energy_final - energy_initial).abs() / energy_initial.abs()
90        } else {
91            (energy_final - energy_initial).abs()
92        };
93        Ok(NBodyResult {
94            num_bodies: n,
95            times,
96            position_snapshots,
97            final_positions: final_state[..2 * n].to_vec(),
98            final_velocities: final_state[2 * n..].to_vec(),
99            energy_initial,
100            energy_final,
101            energy_drift_rel,
102            angular_momentum_initial: angmom_initial,
103            angular_momentum_final: angmom_final,
104            steps_accepted: accepted,
105            steps_rejected: rejected,
106        })
107    }
108}
qualia_core_db/specialized_libs/physics_simulation/nbody.rs

qualia_core_db/specialized_libs/physics_simulation/
nbody.rs
