use crate::constants::physical::G;
use std::f64::consts::PI;
#[derive(Clone, Debug, PartialEq)]
pub struct Particle {
pub position: [f64; 3],
pub velocity: [f64; 3],
pub mass: f64,
}
pub fn gravitational_acceleration_direct(
particles: &[Particle],
index: usize,
softening: f64,
) -> [f64; 3] {
let mut acc = [0.0; 3];
let pi = &particles[index];
for (j, pj) in particles.iter().enumerate() {
if j == index {
continue;
}
let dx = pj.position[0] - pi.position[0];
let dy = pj.position[1] - pi.position[1];
let dz = pj.position[2] - pi.position[2];
let r2 = dx * dx + dy * dy + dz * dz + softening * softening;
let r = r2.sqrt();
let r3 = r * r2;
let f = G * pj.mass / r3;
acc[0] += f * dx;
acc[1] += f * dy;
acc[2] += f * dz;
}
acc
}
pub fn all_accelerations(particles: &[Particle], softening: f64) -> Vec<[f64; 3]> {
let n = particles.len();
let mut accs = vec![[0.0; 3]; n];
for i in 0..n {
for j in (i + 1)..n {
let dx = particles[j].position[0] - particles[i].position[0];
let dy = particles[j].position[1] - particles[i].position[1];
let dz = particles[j].position[2] - particles[i].position[2];
let r2 = dx * dx + dy * dy + dz * dz + softening * softening;
let r3 = r2.sqrt() * r2;
let fi = G * particles[j].mass / r3;
let fj = G * particles[i].mass / r3;
accs[i][0] += fi * dx;
accs[i][1] += fi * dy;
accs[i][2] += fi * dz;
accs[j][0] -= fj * dx;
accs[j][1] -= fj * dy;
accs[j][2] -= fj * dz;
}
}
accs
}
pub fn leapfrog_step_nbody(particles: &mut [Particle], dt: f64, softening: f64) {
let accs = all_accelerations(particles, softening);
for (p, a) in particles.iter_mut().zip(accs.iter()) {
p.velocity[0] += 0.5 * dt * a[0];
p.velocity[1] += 0.5 * dt * a[1];
p.velocity[2] += 0.5 * dt * a[2];
}
for p in particles.iter_mut() {
p.position[0] += dt * p.velocity[0];
p.position[1] += dt * p.velocity[1];
p.position[2] += dt * p.velocity[2];
}
let accs_new = all_accelerations(particles, softening);
for (p, a) in particles.iter_mut().zip(accs_new.iter()) {
p.velocity[0] += 0.5 * dt * a[0];
p.velocity[1] += 0.5 * dt * a[1];
p.velocity[2] += 0.5 * dt * a[2];
}
}
pub fn nbody_simulate(
particles: &mut [Particle],
dt: f64,
n_steps: usize,
softening: f64,
) -> Vec<Vec<[f64; 3]>> {
let mut snapshots = Vec::with_capacity(n_steps + 1);
snapshots.push(particles.iter().map(|p| p.position).collect());
for _ in 0..n_steps {
leapfrog_step_nbody(particles, dt, softening);
snapshots.push(particles.iter().map(|p| p.position).collect());
}
snapshots
}
pub fn total_kinetic_energy(particles: &[Particle]) -> f64 {
particles
.iter()
.map(|p| {
let v2 = p.velocity[0].powi(2) + p.velocity[1].powi(2) + p.velocity[2].powi(2);
0.5 * p.mass * v2
})
.sum()
}
pub fn total_potential_energy(particles: &[Particle], softening: f64) -> f64 {
let n = particles.len();
let mut pe = 0.0;
for i in 0..n {
for j in (i + 1)..n {
let dx = particles[j].position[0] - particles[i].position[0];
let dy = particles[j].position[1] - particles[i].position[1];
let dz = particles[j].position[2] - particles[i].position[2];
let r = (dx * dx + dy * dy + dz * dz + softening * softening).sqrt();
pe -= G * particles[i].mass * particles[j].mass / r;
}
}
pe
}
pub fn center_of_mass(particles: &[Particle]) -> [f64; 3] {
let mut com = [0.0; 3];
let mut total_mass = 0.0;
for p in particles {
com[0] += p.mass * p.position[0];
com[1] += p.mass * p.position[1];
com[2] += p.mass * p.position[2];
total_mass += p.mass;
}
if total_mass > 0.0 {
com[0] /= total_mass;
com[1] /= total_mass;
com[2] /= total_mass;
}
com
}
pub fn velocity_dispersion(particles: &[Particle]) -> f64 {
let n = particles.len() as f64;
if n < 1.0 {
return 0.0;
}
let mut mean_v = [0.0; 3];
for p in particles {
mean_v[0] += p.velocity[0];
mean_v[1] += p.velocity[1];
mean_v[2] += p.velocity[2];
}
mean_v[0] /= n;
mean_v[1] /= n;
mean_v[2] /= n;
let mut sigma2 = 0.0;
for p in particles {
sigma2 += (p.velocity[0] - mean_v[0]).powi(2)
+ (p.velocity[1] - mean_v[1]).powi(2)
+ (p.velocity[2] - mean_v[2]).powi(2);
}
(sigma2 / (3.0 * n)).sqrt()
}
pub fn optimal_softening(n_particles: usize, r_system: f64) -> f64 {
r_system / (n_particles as f64).powf(1.0 / 3.0)
}
pub fn plummer_softened_potential(m: f64, r: f64, epsilon: f64) -> f64 {
-G * m / (r * r + epsilon * epsilon).sqrt()
}
pub fn adaptive_timestep(acc_max: f64, softening: f64, eta: f64) -> f64 {
if acc_max < 1e-50 {
return 1e10;
}
eta * (softening / acc_max).sqrt()
}
pub fn virial_ratio(particles: &[Particle], softening: f64) -> f64 {
let ke = total_kinetic_energy(particles);
let pe = total_potential_energy(particles, softening);
if pe.abs() < 1e-100 {
return 0.0;
}
2.0 * ke / pe.abs()
}
pub fn half_mass_radius(particles: &[Particle]) -> f64 {
let com = center_of_mass(particles);
let total_mass: f64 = particles.iter().map(|p| p.mass).sum();
let mut radii: Vec<(f64, f64)> = particles
.iter()
.map(|p| {
let dx = p.position[0] - com[0];
let dy = p.position[1] - com[1];
let dz = p.position[2] - com[2];
((dx * dx + dy * dy + dz * dz).sqrt(), p.mass)
})
.collect();
radii.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
let half = total_mass * 0.5;
let mut cumulative = 0.0;
for (r, m) in &radii {
cumulative += m;
if cumulative >= half {
return *r;
}
}
radii.last().map_or(0.0, |(r, _)| *r)
}
pub fn crossing_time(particles: &[Particle]) -> f64 {
let r_h = half_mass_radius(particles);
let sigma = velocity_dispersion(particles);
if sigma < 1e-100 {
return f64::INFINITY;
}
2.0 * PI * r_h / sigma
}