cometsfactory 0.0.1

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

use crate::config::parameters;

/// Interstellar comets on hyperbolic orbits (eccentricity > 1.0).
/// Examples: 1I/'Oumuamua, 2I/Borisov.
#[derive(Debug)]
pub struct Interstellar {
    pub nucleus_radius: f64,
    pub density: f64,
    pub albedo: f64,
    pub eccentricity: f64,
    pub excess_velocity: f64,
    pub perihelion_distance: f64,
    pub inclination: f64,
    pub dust_to_gas_ratio: f64,
    pub active_fraction: f64,
}

impl Interstellar {
    pub fn new(nucleus_radius: f64, excess_velocity: f64, perihelion_distance: f64) -> Self {
        Self {
            nucleus_radius: nucleus_radius.max(1.0),
            density: parameters::COMET_NUCLEUS_DENSITY,
            albedo: parameters::COMET_ALBEDO,
            eccentricity: 1.0
                + (excess_velocity.max(0.1)).powi(2) * parameters::AU
                    / (parameters::G * parameters::SOLAR_MASS
                        / (perihelion_distance.max(0.01) * parameters::AU)),
            excess_velocity: excess_velocity.max(0.1),
            perihelion_distance: perihelion_distance.max(0.01),
            inclination: 90.0,
            dust_to_gas_ratio: 0.5,
            active_fraction: 0.05,
        }
    }

    pub fn with_density(mut self, density: f64) -> Self {
        self.density = density.clamp(100.0, 2000.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_dust_to_gas_ratio(mut self, ratio: f64) -> Self {
        self.dust_to_gas_ratio = ratio.clamp(0.01, 100.0);
        self
    }

    pub fn with_active_fraction(mut self, fraction: f64) -> Self {
        self.active_fraction = fraction.clamp(0.001, 1.0);
        self
    }

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

    /// Velocity at infinity (km/s)
    pub fn velocity_at_infinity_kms(&self) -> f64 {
        self.excess_velocity
    }

    /// Velocity at perihelion (km/s) via vis-viva for hyperbolic orbit
    pub fn perihelion_velocity(&self) -> f64 {
        let r = self.perihelion_distance * parameters::AU;
        let v_inf = self.excess_velocity * 1000.0;
        (v_inf * v_inf + 2.0 * parameters::G * parameters::SOLAR_MASS / r).sqrt()
    }

    /// Is the object on a truly hyperbolic (unbound) orbit?
    pub fn is_unbound(&self) -> bool {
        self.eccentricity > 1.0
    }

    pub fn perihelion_temperature(&self) -> f64 {
        let r_m = self.perihelion_distance * parameters::AU;
        parameters::equilibrium_temperature(r_m, self.albedo)
    }

    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)
    }
}