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

use crate::config::parameters::{PI, YEAR};

/// Mantle growth rate (m/s): how fast a refractory lag deposit builds.
/// Very rough model: dm/dt ∝ sublimation rate × dust-to-gas ratio.
/// Returns thickness growth per orbit passage (m).
pub fn mantle_growth_per_orbit(
    nucleus_radius: f64,
    active_fraction: f64,
    dust_to_gas: f64,
    period_years: f64,
) -> f64 {
    let area = 4.0 * PI * nucleus_radius * nucleus_radius * active_fraction;
    let erosion_rate_m_per_s = 1.0e-9;
    let dust_retention = dust_to_gas / (1.0 + dust_to_gas);
    let growth = erosion_rate_m_per_s * dust_retention * period_years * YEAR;
    growth * area / (4.0 * PI * nucleus_radius * nucleus_radius)
}

/// Radius shrinkage per orbit (m) from sublimation.
/// ΔR ≈ dm_dt * P / (4π R² ρ)   where dm_dt is evaluated at perihelion.
pub fn radius_loss_per_orbit(
    dm_dt_perihelion: f64,
    period_years: f64,
    nucleus_radius: f64,
    density: f64,
) -> f64 {
    let area = 4.0 * PI * nucleus_radius * nucleus_radius;
    if area < 1.0 || density < 1.0 {
        return 0.0;
    }
    dm_dt_perihelion * period_years * YEAR / (area * density)
}

/// Number of orbits until the nucleus is fully eroded.
pub fn orbits_to_exhaustion(nucleus_radius: f64, radius_loss_per_pass: f64) -> f64 {
    if radius_loss_per_pass < 1.0e-30 {
        return f64::INFINITY;
    }
    nucleus_radius / radius_loss_per_pass
}

/// Spin-up timescale (s) from asymmetric outgassing torques.
/// τ ≈ I ω / Torque,   very rough: Torque ~ F_jet * R * moment arm fraction.
pub fn spinup_timescale(
    nucleus_radius: f64,
    density: f64,
    rotation_period_s: f64,
    dm_dt: f64,
    gas_velocity: f64,
) -> f64 {
    let mass = (4.0 / 3.0) * PI * nucleus_radius.powi(3) * density;
    let inertia = 0.4 * mass * nucleus_radius * nucleus_radius;
    let omega = 2.0 * PI / rotation_period_s;
    let moment_arm = 0.01;
    let torque = dm_dt * gas_velocity * nucleus_radius * moment_arm;
    if torque < 1.0e-30 {
        return f64::INFINITY;
    }
    inertia * omega / torque
}

/// Splitting probability heuristic: comets with perihelion < q_crit AU
/// and rotation period < P_crit s are at risk.
pub fn splitting_risk(perihelion_au: f64, rotation_period_s: f64) -> f64 {
    let q_factor = if perihelion_au < 1.0 {
        (1.0 / perihelion_au.max(0.01)).sqrt()
    } else {
        1.0 / perihelion_au
    };
    let p_factor = if rotation_period_s < 7200.0 {
        7200.0 / rotation_period_s.max(100.0)
    } else {
        1.0
    };
    (q_factor * p_factor).min(1.0)
}

/// Fading parameter: fraction of brightness retained after N perihelion passages.
/// Simple exponential fading model: f = exp(-N / N_half)
pub fn brightness_fading(passages: f64, half_life_passages: f64) -> f64 {
    if half_life_passages < 1.0e-10 {
        return 0.0;
    }
    (-passages * (2.0_f64).ln() / half_life_passages).exp()
}