dwarfplanetsfactory 0.0.3

Dwarf planet factory — classify, build and catalogue dwarf planets of any type: Kuiper belt, scattered disk, plutino, cold classical, detached, binary, Ceres-type, and sednoid.
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

use crate::config::parameters;

/// Equilibrium temperature at distance r (AU) for a given albedo
pub fn equilibrium_temp(r_au: f64, albedo: f64) -> f64 {
    parameters::equilibrium_temperature(r_au * parameters::AU, albedo)
}

/// Equilibrium temperature at perihelion
pub fn temp_at_perihelion(semi_major_au: f64, eccentricity: f64, albedo: f64) -> f64 {
    let q = parameters::perihelion_au(semi_major_au, eccentricity);
    equilibrium_temp(q, albedo)
}

/// Equilibrium temperature at aphelion
pub fn temp_at_aphelion(semi_major_au: f64, eccentricity: f64, albedo: f64) -> f64 {
    let big_q = parameters::aphelion_au(semi_major_au, eccentricity);
    equilibrium_temp(big_q, albedo)
}

/// Thermal inertia indicator: ratio of orbital period to thermal
/// relaxation timescale (dimensionless). High values → isothermal surface.
pub fn thermal_parameter(
    radius: f64,
    thermal_inertia: f64,
    albedo: f64,
    semi_major_au: f64,
) -> f64 {
    let t_eq = equilibrium_temp(semi_major_au, albedo);
    if t_eq < 1.0 || thermal_inertia < 1.0e-10 {
        return 0.0;
    }
    let period = parameters::orbital_period(semi_major_au, parameters::SOLAR_MASS);
    let omega = 2.0 * parameters::PI / period;
    let skin_depth = radius * 1.0e-6;
    thermal_inertia * (omega).sqrt() / (parameters::STEFAN_BOLTZMANN * t_eq.powi(3))
        + skin_depth * 0.0
}

/// Sublimation rate for exposed ices at temperature T (kg/m²/s)
/// Using Clausius-Clapeyron approximation for N₂ ice
pub fn n2_sublimation_rate(temp_k: f64) -> f64 {
    if temp_k < 10.0 {
        return 0.0;
    }
    let latent_heat: f64 = 2.56e5; // J/kg for N2 ice
    let mol_mass: f64 = 0.028; // kg/mol
    let p_ref: f64 = 1.0e5; // Pa reference
    let t_ref: f64 = 77.36; // K boiling point at 1 atm
    let r_gas: f64 = 8.314;
    let p_vap = p_ref * (-(latent_heat * mol_mass / r_gas) * (1.0 / temp_k - 1.0 / t_ref)).exp();
    p_vap * (mol_mass / (2.0 * parameters::PI * r_gas * temp_k)).sqrt()
}

/// CO ice sublimation rate (kg/m²/s)
pub fn co_sublimation_rate(temp_k: f64) -> f64 {
    if temp_k < 10.0 {
        return 0.0;
    }
    let latent_heat: f64 = 2.27e5;
    let mol_mass: f64 = 0.028;
    let p_ref: f64 = 1.0e5;
    let t_ref: f64 = 81.6;
    let r_gas: f64 = 8.314;
    let p_vap = p_ref * (-(latent_heat * mol_mass / r_gas) * (1.0 / temp_k - 1.0 / t_ref)).exp();
    p_vap * (mol_mass / (2.0 * parameters::PI * r_gas * temp_k)).sqrt()
}

/// Surface recession rate (m/s) from sublimation
pub fn surface_recession(sublimation_rate: f64, ice_density: f64) -> f64 {
    if ice_density < 1.0 {
        return 0.0;
    }
    sublimation_rate / ice_density
}