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
120
121
122
123
124

use crate::CelestialBody;
use crate::gravity::nbody::compute_accelerations;

pub enum Integrator {
    LeapfrogKDK,
    VelocityVerlet,
    Yoshida4,
}

impl Integrator {
    pub fn step(&self, bodies: &mut [CelestialBody], dt: f64) {
        match self {
            Integrator::LeapfrogKDK => leapfrog_kdk(bodies, dt),
            Integrator::VelocityVerlet => velocity_verlet(bodies, dt),
            Integrator::Yoshida4 => yoshida4(bodies, dt),
        }
    }
}

pub fn leapfrog_kdk(bodies: &mut [CelestialBody], dt: f64) {
    // Kick (half)
    for b in bodies.iter_mut() {
        b.velocity += b.acceleration * (0.5 * dt);
    }
    // Drift
    for b in bodies.iter_mut() {
        b.position += b.velocity * dt;
    }
    // Update accelerations
    compute_accelerations(bodies);
    // Kick (half)
    for b in bodies.iter_mut() {
        b.velocity += b.acceleration * (0.5 * dt);
    }
}

pub fn velocity_verlet(bodies: &mut [CelestialBody], dt: f64) {
    // Position update
    for b in bodies.iter_mut() {
        b.position += b.velocity * dt + b.acceleration * (0.5 * dt * dt);
    }
    // Store old accelerations
    let old_acc: Vec<_> = bodies.iter().map(|b| b.acceleration).collect();
    // Compute new accelerations
    compute_accelerations(bodies);
    // Velocity update
    for (b, a_old) in bodies.iter_mut().zip(old_acc.iter()) {
        b.velocity += (*a_old + b.acceleration) * (0.5 * dt);
    }
}

pub fn yoshida4(bodies: &mut [CelestialBody], dt: f64) {
    const W0: f64 = -1.702414383919315;
    const W1: f64 = 1.351207191959657;
    let c1 = W1 / 2.0;
    let c2 = (W0 + W1) / 2.0;
    let c3 = c2;
    let c4 = c1;
    let d1 = W1;
    let d2 = W0;
    let d3 = W1;

    // x1
    for b in bodies.iter_mut() {
        b.position += b.velocity * (c1 * dt);
    }
    compute_accelerations(bodies);
    for b in bodies.iter_mut() {
        b.velocity += b.acceleration * (d1 * dt);
    }
    // x2
    for b in bodies.iter_mut() {
        b.position += b.velocity * (c2 * dt);
    }
    compute_accelerations(bodies);
    for b in bodies.iter_mut() {
        b.velocity += b.acceleration * (d2 * dt);
    }
    // x3
    for b in bodies.iter_mut() {
        b.position += b.velocity * (c3 * dt);
    }
    compute_accelerations(bodies);
    for b in bodies.iter_mut() {
        b.velocity += b.acceleration * (d3 * dt);
    }
    // x4
    for b in bodies.iter_mut() {
        b.position += b.velocity * (c4 * dt);
    }
}

pub fn adaptive_step(
    bodies: &mut [CelestialBody],
    dt: f64,
    tolerance: f64,
    integrator: &Integrator,
) -> f64 {
    let snapshot: Vec<CelestialBody> = bodies.to_vec();

    // Full step
    integrator.step(bodies, dt);
    let full: Vec<_> = bodies.iter().map(|b| b.position).collect();

    // Two half steps
    bodies.clone_from_slice(&snapshot);
    integrator.step(bodies, dt / 2.0);
    integrator.step(bodies, dt / 2.0);

    let max_err = bodies
        .iter()
        .zip(full.iter())
        .map(|(b, f)| b.position.distance(f) / b.position.magnitude().max(1.0))
        .fold(0.0_f64, f64::max);

    if max_err > tolerance && dt > 1.0 {
        bodies.clone_from_slice(&snapshot);
        adaptive_step(bodies, dt / 2.0, tolerance, integrator)
    } else if max_err < tolerance * 0.01 {
        dt * 2.0
    } else {
        dt
    }
}