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

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

/// Coma radius estimate (m) from gas production rate Q (molecules/s).
/// R_coma ≈ v_gas * τ_photo,  where τ_photo ≈ 8.2e4 s for water at 1 AU,
/// scaled as r².
pub fn coma_radius(gas_production_q: f64, r_au: f64, albedo: f64) -> f64 {
    let v_gas = gas_thermal_velocity(r_au, albedo);
    let tau_photo_1au = 8.2e4;
    let tau_photo = tau_photo_1au * r_au * r_au;
    let radius = v_gas * tau_photo;
    let _ = gas_production_q;
    radius.max(0.0)
}

/// Dust tail length (m) — very rough Finson-Probstein approximation.
/// L ∝ (1 - β) × acceleration_time² , simplified to geometric scaling.
pub fn dust_tail_length(beta: f64, r_au: f64, orbital_velocity_ms: f64) -> f64 {
    let r = r_au * AU;
    if r < 1.0 || orbital_velocity_ms < 1.0 {
        return 0.0;
    }
    let a_rad = beta * 5.93e-3 / (r_au * r_au);
    let time = r / orbital_velocity_ms;
    0.5 * a_rad * time * time
}

/// Ion tail length (m) — proportional to solar wind interaction region.
/// Rough estimate: L_ion ~ 1 AU × (Q / Q_ref)^0.5  at 1 AU, scaled by r.
pub fn ion_tail_length(gas_production_q: f64, r_au: f64) -> f64 {
    let q_ref = 1.0e28;
    let l_ref = AU;
    l_ref * (gas_production_q / q_ref).sqrt() / r_au.max(0.01)
}

/// Afρ parameter (cm) — proxy for dust production.
/// Afρ = (2 Δ r / ρ_aperture)² × 10^(0.4*(m_sun - m_comet))
/// Simplified: proportional to Q_dust / r²
pub fn afrho(dust_production_rate: f64, r_au: f64) -> f64 {
    if r_au < 0.01 {
        return 0.0;
    }
    dust_production_rate / (r_au * r_au)
}

/// Gas thermal velocity (m/s) for water at equilibrium temperature.
fn gas_thermal_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()
}

/// Optical depth of the coma at distance d from nucleus center.
/// τ ≈ Q σ / (4 π v_gas d)  where σ is the scattering cross-section.
pub fn coma_optical_depth(gas_production_q: f64, distance_m: f64, r_au: f64, albedo: f64) -> f64 {
    let v = gas_thermal_velocity(r_au, albedo);
    if v < 1.0 || distance_m < 1.0 {
        return 0.0;
    }
    let sigma = 1.0e-20;
    gas_production_q * sigma / (4.0 * PI * v * distance_m)
}

/// Hydrogen Lyman-alpha coma brightness (photons/s) — proportional to Q(H₂O).
pub fn lyman_alpha_production(gas_production_q: f64) -> f64 {
    gas_production_q * 0.85
}