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

use std::f64::consts::PI;

use asteroidsfactory::types::main_belt::{MainBeltAsteroid, SpectralType};

use crate::config::parameters::MU_SUN;
use crate::orbits::keplerian::OrbitalElements;
use crate::{BodyId, CelestialBody, Vec3};

const AU: f64 = 1.495978707e11;

/// Simple deterministic PRNG (xorshift64) for reproducible belt generation.
struct Rng(u64);

impl Rng {
    fn next_u64(&mut self) -> u64 {
        self.0 ^= self.0 << 13;
        self.0 ^= self.0 >> 7;
        self.0 ^= self.0 << 17;
        self.0
    }

    fn f64(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 / ((1u64 << 53) as f64)
    }

    fn range(&mut self, lo: f64, hi: f64) -> f64 {
        lo + self.f64() * (hi - lo)
    }
}

/// Spectral type distribution: C ≈40%, S ≈30%, M ≈10%, B ≈5%, D ≈5%, V ≈5%, X ≈5%
fn random_spectral_type(rng: &mut Rng) -> SpectralType {
    let r = rng.f64();
    if r < 0.40 {
        SpectralType::C
    } else if r < 0.70 {
        SpectralType::S
    } else if r < 0.80 {
        SpectralType::M
    } else if r < 0.85 {
        SpectralType::B
    } else if r < 0.90 {
        SpectralType::D
    } else if r < 0.95 {
        SpectralType::V
    } else {
        SpectralType::X
    }
}

/// Density by spectral type (kg/m³).
fn density_for_type(st: SpectralType, rng: &mut Rng) -> f64 {
    match st {
        SpectralType::C | SpectralType::B | SpectralType::D => rng.range(1200.0, 2500.0),
        SpectralType::S | SpectralType::V => rng.range(2500.0, 3800.0),
        SpectralType::M | SpectralType::X => rng.range(3500.0, 5500.0),
    }
}

/// Generate `count` main-belt asteroids using `asteroidsfactory::MainBeltAsteroid`.
///
/// Semi-major axes: 2.1–3.3 AU (rejecting Kirkwood gaps via the factory method).
/// Eccentricities: 0–0.3.
/// Inclinations: 0–20°.
/// Positions & velocities computed from Keplerian elements.
pub fn generate_main_belt(count: u16) -> Vec<CelestialBody> {
    let mut rng = Rng(0xCE1E_5714_A573_B317);
    let mut bodies = Vec::with_capacity(count as usize);

    for idx in 0..count {
        let spectral = random_spectral_type(&mut rng);
        let density = density_for_type(spectral, &mut rng);
        let radius = 10.0_f64.powf(rng.range(3.0, 5.5)); // 1 km – 316 km

        // Semi-major axis: reject Kirkwood gaps using the factory
        let a_au = loop {
            let candidate = rng.range(2.1, 3.3);
            let test = MainBeltAsteroid::new(radius, density, candidate);
            if !test.is_in_kirkwood_gap() {
                break candidate;
            }
        };

        let ecc = rng.range(0.0, 0.3);
        let inc_deg = rng.range(0.0, 20.0);

        let asteroid = MainBeltAsteroid::new(radius, density, a_au)
            .with_eccentricity(ecc)
            .with_inclination(inc_deg)
            .with_spectral_type(spectral)
            .with_rotation_period(rng.range(2.0, 24.0));

        // Convert orbital elements to cartesian state vectors
        let raan = rng.range(0.0, 2.0 * PI);
        let arg_periapsis = rng.range(0.0, 2.0 * PI);
        let true_anomaly = rng.range(0.0, 2.0 * PI);

        let elements = OrbitalElements {
            semi_major_axis: asteroid.semi_major_axis * AU,
            eccentricity: asteroid.eccentricity,
            inclination: asteroid.inclination, // already in radians from with_inclination
            longitude_ascending_node: raan,
            argument_periapsis: arg_periapsis,
            true_anomaly,
        };

        let (position, velocity) = elements.to_state_vectors(MU_SUN);

        bodies.push(CelestialBody {
            id: BodyId::BeltAsteroid(idx),
            name: "Belt Asteroid",
            mass: asteroid.mass,
            radius: asteroid.radius,
            position,
            velocity,
            acceleration: Vec3::ZERO,
            j2: 0.0,
        });
    }

    bodies
}