satellitesfactory 0.0.3

Satellite factory — classify, build and catalogue natural satellites for any planetary system: Solar System moons (Moon, Galileans, Titan, Triton…) or custom configurations.
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

use crate::config::parameters::*;

pub struct SatelliteEvolution {
    pub mass: f64,
    pub radius: f64,
    pub semi_major_axis: f64,
    pub eccentricity: f64,
    pub parent_mass: f64,
    pub parent_radius: f64,
}

impl SatelliteEvolution {
    pub fn new(
        mass: f64,
        radius: f64,
        semi_major_axis: f64,
        eccentricity: f64,
        parent_mass: f64,
        parent_radius: f64,
    ) -> Self {
        Self {
            mass,
            radius,
            semi_major_axis,
            eccentricity,
            parent_mass,
            parent_radius,
        }
    }

    pub fn tidal_migration_rate(&self, q_planet: f64) -> f64 {
        let k2p = 0.38;
        let n = (G * self.parent_mass / self.semi_major_axis.powi(3)).sqrt();
        3.0 * k2p * self.mass * n * self.parent_radius.powi(5)
            / (q_planet * self.parent_mass * self.semi_major_axis.powi(4))
    }

    pub fn orbital_decay_timescale(&self, q_planet: f64) -> f64 {
        let da_dt = self.tidal_migration_rate(q_planet);
        if da_dt.abs() < 1e-30 {
            return f64::INFINITY;
        }
        self.semi_major_axis / da_dt.abs()
    }

    pub fn circularization_timescale(&self, q_sat: f64, k2_sat: f64) -> f64 {
        let n = (G * self.parent_mass / self.semi_major_axis.powi(3)).sqrt();
        if n < 1e-30 {
            return f64::INFINITY;
        }
        (2.0 * q_sat * self.mass * self.semi_major_axis.powi(5) * n)
            / (63.0 * k2_sat * G * self.parent_mass * self.parent_mass * self.radius.powi(5))
    }

    pub fn mean_motion_resonance_period_ratio(&self, other_a: f64) -> f64 {
        (self.semi_major_axis / other_a).powf(1.5)
    }

    pub fn capture_into_resonance_probability(mass_ratio: f64, migration_rate: f64) -> f64 {
        let mu = mass_ratio;
        let tau = migration_rate.abs();
        if tau < 1e-30 {
            return 1.0;
        }
        (mu.powf(0.3) / (tau * 1e15 + 1.0)).min(1.0)
    }

    pub fn sputtering_rate(&self, ion_flux: f64) -> f64 {
        let cross_section = std::f64::consts::PI * self.radius * self.radius;
        let yield_factor = 0.1;
        ion_flux * cross_section * yield_factor * 18.0 * 1.66e-27
    }

    pub fn sublimation_mass_loss(&self, temperature: f64) -> f64 {
        if temperature < 100.0 {
            return 0.0;
        }
        let p_vap = 6.11e2 * (-5370.0 / temperature + 19.8).exp();
        let area = 4.0 * std::f64::consts::PI * self.radius * self.radius;
        let m_h2o = 18.0 * 1.66e-27;
        area * p_vap * (m_h2o / (2.0 * std::f64::consts::PI * K_B * temperature)).sqrt()
    }

    pub fn impact_gardening_depth(flux: f64, age: f64) -> f64 {
        0.1 * (flux * age).powf(0.5)
    }
}

pub fn secular_acceleration(a: f64, da_dt: f64) -> f64 {
    if a < 1e-3 {
        return 0.0;
    }
    1.5 * da_dt / a
}

pub fn kozai_lidov_max_eccentricity(i_mutual: f64) -> f64 {
    let cos_i = i_mutual.cos();
    (1.0 - (5.0 / 3.0) * cos_i * cos_i).sqrt().max(0.0)
}

pub fn yarkovsky_drift(radius: f64, density: f64, luminosity: f64, distance: f64) -> f64 {
    let area = std::f64::consts::PI * radius * radius;
    let flux = luminosity / (4.0 * std::f64::consts::PI * distance * distance);
    let force = area * flux / SPEED_OF_LIGHT;
    let mass = 4.0 / 3.0 * std::f64::consts::PI * radius.powi(3) * density;
    force / (mass + 1e-30)
}