solarsystems 0.0.1

N-body solar system engine — gravitational dynamics, orbital mechanics, perturbations, event detection, and full celestial orchestration
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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119

use crate::gravity::interactions::{j2_acceleration, relativistic_correction, total_acceleration};
use crate::{CelestialBody, Vec3};

/// Mass threshold below which a body is treated as a test particle.
/// Test particles receive gravity from massive bodies but do NOT contribute
/// gravity to others. This reduces the complexity from O(N²) to O(N_massive × N_total)
/// while remaining physically accurate (asteroids don't perturb planets).
const TEST_PARTICLE_MASS: f64 = 1e20;

fn dm_halo_acceleration(body: &CelestialBody) -> Vec3 {
    let r_vec = body.position;
    let r = r_vec.magnitude();
    if r < 1.0 {
        return Vec3::ZERO;
    }
    let rho_local = crate::darkmatter::halo::local_dm_density();
    let m_enc = crate::darkmatter::halo::dm_enclosed_mass_uniform(rho_local, r);
    let factor = darkmatter::constants::physical::G * m_enc / (r * r * r);
    r_vec * (-factor)
}

pub fn compute_accelerations(bodies: &mut [CelestialBody]) {
    let snapshot: Vec<CelestialBody> = bodies.to_vec();

    // Pre-collect indices of massive bodies for the test-particle optimisation.
    let massive_indices: Vec<usize> = snapshot
        .iter()
        .enumerate()
        .filter(|(_, b)| b.mass >= TEST_PARTICLE_MASS)
        .map(|(i, _)| i)
        .collect();

    for body in bodies.iter_mut() {
        let mut acc = Vec3::ZERO;

        if body.mass >= TEST_PARTICLE_MASS {
            // Massive body: only sum gravity from other massive bodies.
            for &mi in &massive_indices {
                if snapshot[mi].id == body.id {
                    continue;
                }
                acc +=
                    crate::gravity::interactions::gravitational_acceleration(body, &snapshot[mi]);
            }
        } else {
            // Test particle: receive gravity from massive bodies only.
            for &mi in &massive_indices {
                acc +=
                    crate::gravity::interactions::gravitational_acceleration(body, &snapshot[mi]);
            }
        }

        if let Some(parent_id) = body.id.parent()
            && let Some(parent) = snapshot.iter().find(|b| b.id == parent_id)
        {
            acc += j2_acceleration(body, parent);
        }
        acc += dm_halo_acceleration(body);
        body.acceleration = acc;
    }
}

pub fn compute_accelerations_relativistic(bodies: &mut [CelestialBody]) {
    let snapshot: Vec<CelestialBody> = bodies.to_vec();

    let massive_indices: Vec<usize> = snapshot
        .iter()
        .enumerate()
        .filter(|(_, b)| b.mass >= TEST_PARTICLE_MASS)
        .map(|(i, _)| i)
        .collect();

    for body in bodies.iter_mut() {
        let mut acc = Vec3::ZERO;

        if body.mass >= TEST_PARTICLE_MASS {
            for &mi in &massive_indices {
                if snapshot[mi].id == body.id {
                    continue;
                }
                acc +=
                    crate::gravity::interactions::gravitational_acceleration(body, &snapshot[mi]);
            }
        } else {
            for &mi in &massive_indices {
                acc +=
                    crate::gravity::interactions::gravitational_acceleration(body, &snapshot[mi]);
            }
        }

        if let Some(parent_id) = body.id.parent()
            && let Some(parent) = snapshot.iter().find(|b| b.id == parent_id)
        {
            acc += j2_acceleration(body, parent);
            acc += relativistic_correction(body, parent);
        }
        acc += dm_halo_acceleration(body);
        body.acceleration = acc;
    }
}

pub fn compute_pure_gravitational(bodies: &mut [CelestialBody]) {
    let snapshot: Vec<CelestialBody> = bodies.to_vec();
    for (i, body) in bodies.iter_mut().enumerate() {
        body.acceleration = total_acceleration(&snapshot[i], &snapshot);
    }
}

pub fn total_system_energy(bodies: &[CelestialBody]) -> f64 {
    crate::config::parameters::total_energy(bodies)
}

pub fn total_system_angular_momentum(bodies: &[CelestialBody]) -> Vec3 {
    crate::config::parameters::total_angular_momentum(bodies)
}

pub fn barycenter(bodies: &[CelestialBody]) -> (Vec3, Vec3) {
    crate::config::parameters::system_barycenter(bodies)
}