darkmatter 0.0.2

Dark matter simulation engine — gravitational fields, particle dynamics, halo stability, and cosmological constants
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

use darkmatter::dynamics::barnes_hut::{
    LeapfrogParams, all_accelerations_barnes_hut, build_octree, leapfrog_step_barnes_hut,
    optimal_half_width, tree_potential_energy,
};

fn main() {
    let mut positions = [
        [1.0, 0.0, 0.0],
        [-1.0, 0.0, 0.0],
        [0.0, 1.0, 0.0],
        [0.0, -1.0, 0.0],
    ];
    let mut velocities = [
        [0.0, 0.01, 0.0],
        [0.0, -0.01, 0.0],
        [-0.01, 0.0, 0.0],
        [0.01, 0.0, 0.0],
    ];
    let masses = [1e30, 1e30, 1e30, 1e30];

    let hw = optimal_half_width(&positions);
    let center = [0.0, 0.0, 0.0];
    let theta = 0.5;
    let softening = 0.01;

    let tree = build_octree(&positions, &masses, center, hw);
    let total_mass = tree.total_mass;

    let accs = all_accelerations_barnes_hut(&positions, &masses, center, hw, theta, softening);
    let acc_magnitude: f64 = accs
        .iter()
        .map(|a| (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt())
        .sum();

    let pe = tree_potential_energy(&positions, &masses, center, hw, theta, softening);

    let params = LeapfrogParams {
        dt: 1e4,
        center,
        half_width: hw,
        theta,
        softening,
    };

    for step in 0..10 {
        leapfrog_step_barnes_hut(&mut positions, &mut velocities, &masses, &params);
        let ke: f64 = velocities
            .iter()
            .zip(masses.iter())
            .map(|(v, m)| 0.5 * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]))
            .sum();
        let check = ke + step as f64;
        assert!(check.is_finite());
    }

    let check = total_mass + acc_magnitude + pe;
    assert!(check.is_finite());
}