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

use sciforge::hub::domain::common::constants::G;

use crate::config::parameters::GRAVITATIONAL_SOFTENING;
use crate::{CelestialBody, Vec3};

pub fn gravitational_acceleration(body: &CelestialBody, source: &CelestialBody) -> Vec3 {
    let r = source.position - body.position;
    let dist_sq = r.magnitude_squared() + GRAVITATIONAL_SOFTENING * GRAVITATIONAL_SOFTENING;
    let dist = dist_sq.sqrt();
    r * (G * source.mass / (dist_sq * dist))
}

pub fn total_acceleration(body: &CelestialBody, bodies: &[CelestialBody]) -> Vec3 {
    let mut acc = Vec3::ZERO;
    for other in bodies {
        if std::ptr::eq(body, other) {
            continue;
        }
        acc += gravitational_acceleration(body, other);
    }
    acc
}

pub fn tidal_acceleration(
    body_pos: Vec3,
    perturber_mass: f64,
    perturber_pos: Vec3,
    central_pos: Vec3,
) -> Vec3 {
    let r_body = body_pos - perturber_pos;
    let r_central = central_pos - perturber_pos;
    let rb3 = r_body.magnitude().powi(3);
    let rc3 = r_central.magnitude().powi(3);
    (r_central / rc3 - r_body / rb3) * (G * perturber_mass)
}

pub fn mutual_gravitational_force(a: &CelestialBody, b: &CelestialBody) -> f64 {
    let r = a.position.distance(&b.position);
    G * a.mass * b.mass / (r * r)
}

pub fn tidal_force_magnitude(perturber_mass: f64, distance: f64, body_radius: f64) -> f64 {
    2.0 * G * perturber_mass * body_radius / distance.powi(3)
}

pub fn j2_acceleration(body: &CelestialBody, central: &CelestialBody) -> Vec3 {
    if central.j2 == 0.0 {
        return Vec3::ZERO;
    }
    let r_vec = body.position - central.position;
    let r = r_vec.magnitude();
    let mu = G * central.mass;
    let re = central.radius;
    let j2 = central.j2;
    let z = r_vec.z;
    let r2 = r * r;
    let factor = 1.5 * j2 * mu * re * re / (r2 * r2 * r);
    let zr2 = (z * z) / r2;
    Vec3::new(
        r_vec.x * factor * (5.0 * zr2 - 1.0),
        r_vec.y * factor * (5.0 * zr2 - 1.0),
        r_vec.z * factor * (5.0 * zr2 - 3.0),
    )
}

pub fn relativistic_correction(body: &CelestialBody, central: &CelestialBody) -> Vec3 {
    let c = 2.99792458e8;
    let r_vec = body.position - central.position;
    let r = r_vec.magnitude();
    let v = body.velocity;
    let mu = G * central.mass;
    let v2 = v.magnitude_squared();
    let r_hat = r_vec.normalize();
    let rdot = r_vec.dot(&v) / r;
    let factor = mu / (r * r * c * c);
    let term1 = r_hat * ((4.0 * mu / r - v2) * factor);
    let term2 = v * (4.0 * rdot * factor);
    term1 + term2
}