satellitesfactory 0.0.2

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

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

pub struct SatelliteInterior {
    pub mass: f64,
    pub radius: f64,
    pub ice_fraction: f64,
    pub rock_fraction: f64,
    pub core_density: f64,
    pub mantle_density: f64,
}

impl SatelliteInterior {
    pub fn new(mass: f64, radius: f64, ice_fraction: f64, rock_fraction: f64) -> Self {
        let metal_fraction = (1.0 - ice_fraction - rock_fraction).max(0.0);
        let core_density = IRON_DENSITY * metal_fraction + ROCK_DENSITY * (1.0 - metal_fraction);
        let mantle_density = if ice_fraction > 0.01 {
            ice_fraction * WATER_ICE_DENSITY + (1.0 - ice_fraction) * ROCK_DENSITY
        } else {
            ROCK_DENSITY
        };
        Self {
            mass,
            radius,
            ice_fraction,
            rock_fraction,
            core_density,
            mantle_density,
        }
    }

    pub fn mean_density(&self) -> f64 {
        mean_density(self.mass, self.radius)
    }

    pub fn core_radius(&self) -> f64 {
        let metal_fraction = (1.0 - self.ice_fraction - self.rock_fraction).max(0.0);
        self.radius * metal_fraction.cbrt()
    }

    pub fn mantle_thickness(&self) -> f64 {
        self.radius - self.core_radius()
    }

    pub fn central_pressure(&self) -> f64 {
        let rho = self.mean_density();
        2.0 * std::f64::consts::PI / 3.0 * G * rho * rho * self.radius * self.radius
    }

    pub fn moment_of_inertia_factor(&self) -> f64 {
        let r_c = self.core_radius();
        let r = self.radius;
        if r < 1e-3 {
            return 0.4;
        }
        let x = r_c / r;
        let rho_c = self.core_density;
        let rho_m = self.mantle_density;
        let top = rho_c * x.powi(5) + rho_m * (1.0 - x.powi(5));
        let bot = rho_c * x.powi(3) + rho_m * (1.0 - x.powi(3));
        0.4 * top / bot
    }

    pub fn gravitational_binding_energy(&self) -> f64 {
        3.0 * G * self.mass * self.mass / (5.0 * self.radius)
    }

    pub fn ocean_depth_estimate(&self) -> f64 {
        if self.ice_fraction < 0.01 {
            return 0.0;
        }
        let ice_shell = 0.1 * self.radius * self.ice_fraction;
        let ocean = self.radius * self.ice_fraction * 0.3;
        if ocean > ice_shell { ocean } else { 0.0 }
    }
}

pub fn differentiation_energy(mass: f64, radius: f64, core_fraction: f64) -> f64 {
    let e_bind = 3.0 * G * mass * mass / (5.0 * radius);
    e_bind * 0.1 * core_fraction
}

pub fn ice_melting_depth(surface_temp: f64, tidal_flux: f64) -> f64 {
    if tidal_flux < 1e-10 {
        return f64::INFINITY;
    }
    let k_ice = 2.2;
    let t_melt = 273.0;
    if surface_temp >= t_melt {
        return 0.0;
    }
    k_ice * (t_melt - surface_temp) / tidal_flux
}

pub fn convective_vigor(tidal_flux: f64, viscosity: f64, alpha: f64, depth: f64) -> f64 {
    let g_local = 1.0;
    let rho = WATER_ICE_DENSITY;
    let kappa = 1e-6;
    let delta_t = tidal_flux * depth / 2.2;
    rho * g_local * alpha * delta_t * depth.powi(3) / (viscosity * kappa)
}