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

use crate::config::parameters::{
    AU, K_B, LATENT_HEAT_WATER, PI, SOLAR_LUMINOSITY, STEFAN_BOLTZMANN, WATER_MOLECULE_MASS,
};

/// Water-ice sublimation rate per unit area (kg/m²/s) at heliocentric
/// distance r_au, assuming energy-limited sublimation.
/// Z ≈ F_absorbed / L   with F = L_sun (1-A) / (16 π r²)
pub fn sublimation_rate(r_au: f64, albedo: f64) -> f64 {
    let r = r_au * AU;
    if r < 1.0 {
        return 0.0;
    }
    let flux = SOLAR_LUMINOSITY * (1.0 - albedo) / (16.0 * PI * r * r);
    flux / LATENT_HEAT_WATER
}

/// Total gas production rate Q (molecules/s) from a spherical nucleus.
/// Q = Z * active_area / m_H2O
pub fn gas_production_rate(
    nucleus_radius: f64,
    active_fraction: f64,
    r_au: f64,
    albedo: f64,
) -> f64 {
    let z = sublimation_rate(r_au, albedo);
    let area = 4.0 * PI * nucleus_radius * nucleus_radius * active_fraction;
    z * area / WATER_MOLECULE_MASS
}

/// Mass loss rate (kg/s) from sublimation.
pub fn mass_loss_rate(nucleus_radius: f64, active_fraction: f64, r_au: f64, albedo: f64) -> f64 {
    let z = sublimation_rate(r_au, albedo);
    let area = 4.0 * PI * nucleus_radius * nucleus_radius * active_fraction;
    z * area
}

/// Non-gravitational acceleration magnitude (m/s²) on a nucleus
/// of mass m due to asymmetric outgassing, rough Marsden model:
/// a_ng ≈ v_gas * dm/dt / m
pub fn nongrav_acceleration(
    nucleus_radius: f64,
    density: f64,
    active_fraction: f64,
    r_au: f64,
    albedo: f64,
) -> f64 {
    let mass = (4.0 / 3.0) * PI * nucleus_radius.powi(3) * density;
    if mass < 1.0 {
        return 0.0;
    }
    let dm = mass_loss_rate(nucleus_radius, active_fraction, r_au, albedo);
    let v_gas = thermal_gas_velocity(r_au, albedo);
    v_gas * dm / mass
}

/// Thermal velocity of sublimated water molecules (m/s)
/// v = sqrt(8 k_B T / (π m_H2O))
pub fn thermal_gas_velocity(r_au: f64, albedo: f64) -> f64 {
    let r = r_au * AU;
    if r < 1.0 {
        return 0.0;
    }
    let t_eq =
        (SOLAR_LUMINOSITY * (1.0 - albedo) / (16.0 * PI * STEFAN_BOLTZMANN * r * r)).powf(0.25);
    (8.0 * K_B * t_eq / (PI * WATER_MOLECULE_MASS)).sqrt()
}

/// Sublimation-driven lifetime estimate (s).
/// τ ≈ M / dm_dt   evaluated at characteristic distance r_au.
pub fn sublimation_lifetime(
    nucleus_radius: f64,
    density: f64,
    active_fraction: f64,
    r_au: f64,
    albedo: f64,
) -> f64 {
    let mass = (4.0 / 3.0) * PI * nucleus_radius.powi(3) * density;
    let dm = mass_loss_rate(nucleus_radius, active_fraction, r_au, albedo);
    if dm < 1.0e-30 {
        return f64::INFINITY;
    }
    mass / dm
}

/// Equilibrium surface temperature (K) at heliocentric distance r_au.
pub fn surface_temperature(r_au: f64, albedo: f64) -> f64 {
    let r = r_au * AU;
    if r < 1.0 {
        return 0.0;
    }
    (SOLAR_LUMINOSITY * (1.0 - albedo) / (16.0 * PI * STEFAN_BOLTZMANN * r * r)).powf(0.25)
}