cometsfactory 0.0.3

Comet factory — classify, build and catalogue comets of any type: short-period, long-period, Halley-type, sungrazer, interstellar, main-belt comet, centaur-transition, and extinct.
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

use crate::config::parameters;

/// Extinct (or dormant) comets: nuclei that have exhausted their volatiles
/// or developed a thick mantle suppressing outgassing.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum DormancyState {
    Dormant,
    Extinct,
}

#[derive(Debug)]
pub struct Extinct {
    pub nucleus_radius: f64,
    pub density: f64,
    pub albedo: f64,
    pub semi_major_axis: f64,
    pub eccentricity: f64,
    pub inclination: f64,
    pub mantle_thickness: f64,
    pub residual_ice_fraction: f64,
    pub state: DormancyState,
}

impl Extinct {
    pub fn new(nucleus_radius: f64, semi_major_axis: f64, eccentricity: f64) -> Self {
        Self {
            nucleus_radius: nucleus_radius.max(10.0),
            density: 1200.0,
            albedo: 0.04,
            semi_major_axis: semi_major_axis.clamp(1.0, 50.0),
            eccentricity: eccentricity.clamp(0.0, 0.99),
            inclination: 15.0,
            mantle_thickness: 1.0,
            residual_ice_fraction: 0.0,
            state: DormancyState::Extinct,
        }
    }

    pub fn with_density(mut self, density: f64) -> Self {
        self.density = density.clamp(500.0, 3500.0);
        self
    }

    pub fn with_albedo(mut self, albedo: f64) -> Self {
        self.albedo = albedo.clamp(0.01, 0.5);
        self
    }

    pub fn with_inclination(mut self, inclination: f64) -> Self {
        self.inclination = inclination.clamp(0.0, 180.0);
        self
    }

    pub fn with_mantle_thickness(mut self, thickness: f64) -> Self {
        self.mantle_thickness = thickness.clamp(0.01, 100.0);
        self
    }

    pub fn with_residual_ice(mut self, fraction: f64) -> Self {
        self.residual_ice_fraction = fraction.clamp(0.0, 0.3);
        if self.residual_ice_fraction > 0.001 {
            self.state = DormancyState::Dormant;
        }
        self
    }

    pub fn mass(&self) -> f64 {
        parameters::sphere_mass(self.nucleus_radius, self.density)
    }

    pub fn perihelion(&self) -> f64 {
        parameters::perihelion_au(self.semi_major_axis, self.eccentricity)
    }

    pub fn aphelion(&self) -> f64 {
        parameters::aphelion_au(self.semi_major_axis, self.eccentricity)
    }

    pub fn orbital_period_years(&self) -> f64 {
        parameters::orbital_period(self.semi_major_axis, parameters::SOLAR_MASS) / parameters::YEAR
    }

    /// Tisserand parameter w.r.t. Jupiter (a_J = 5.2 AU)
    pub fn tisserand_jupiter(&self) -> f64 {
        let a_j = 5.2;
        let ratio = a_j / self.semi_major_axis;
        let inc_rad = self.inclination.to_radians();
        ratio
            + 2.0
                * ((1.0 - self.eccentricity * self.eccentricity) * (self.semi_major_axis / a_j))
                    .sqrt()
                * inc_rad.cos()
    }

    /// Could this object be a "stealth" comet on a near-Earth orbit?
    pub fn is_neo_like(&self) -> bool {
        self.perihelion() < 1.3
    }

    /// Could it reactivate if the mantle is breached?
    pub fn can_reactivate(&self) -> bool {
        self.residual_ice_fraction > 0.001
    }

    pub fn surface_gravity(&self) -> f64 {
        parameters::surface_gravity(self.mass(), self.nucleus_radius)
    }

    pub fn escape_velocity(&self) -> f64 {
        parameters::escape_velocity(self.mass(), self.nucleus_radius)
    }
}