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

use crate::config::parameters::PI;

/// Absolute total magnitude H_T of a comet using the standard law:
/// m = H_T + 5 log10(Δ) + 2.5 n log10(r)
/// Returns H_T given an apparent magnitude, geocentric distance Δ (AU),
/// heliocentric distance r (AU), and activity index n.
pub fn absolute_total_magnitude(
    apparent_mag: f64,
    delta_au: f64,
    r_au: f64,
    activity_index: f64,
) -> f64 {
    apparent_mag
        - 5.0 * delta_au.max(0.001).log10()
        - 2.5 * activity_index * r_au.max(0.001).log10()
}

/// Apparent total magnitude from absolute magnitude H_T.
pub fn apparent_total_magnitude(h_t: f64, delta_au: f64, r_au: f64, activity_index: f64) -> f64 {
    h_t + 5.0 * delta_au.max(0.001).log10() + 2.5 * activity_index * r_au.max(0.001).log10()
}

/// Nuclear magnitude (point source) assuming geometric albedo p and
/// nucleus radius R (m).
/// H_nuc = M_sun - 2.5 log10(p (R/1AU_m)²)
/// We use M_sun_V ≈ -26.74 and 1 AU = 1.496e11 m.
pub fn nuclear_magnitude(nucleus_radius: f64, albedo: f64) -> f64 {
    let au_m = 1.496e11;
    let m_sun = -26.74;
    let cross_section = PI * nucleus_radius * nucleus_radius;
    let reflected = albedo * cross_section / (PI * au_m * au_m);
    if reflected < 1.0e-50 {
        return 99.0;
    }
    m_sun - 2.5 * reflected.log10()
}

/// Phase angle correction (mag) using the HG system β function.
/// Very simplified: Δm ≈ β * phase_angle_deg / 100
pub fn phase_correction(phase_angle_deg: f64, beta_coeff: f64) -> f64 {
    beta_coeff * phase_angle_deg / 100.0
}

/// Total visual cross-section of the coma (m²) from Afρ (cm) and aperture ρ (m).
/// σ_coma = π ρ² × (Afρ / ρ_cm)
pub fn coma_cross_section(afrho_cm: f64, aperture_m: f64) -> f64 {
    let aperture_cm = aperture_m * 100.0;
    if aperture_cm < 0.01 {
        return 0.0;
    }
    PI * aperture_m * aperture_m * (afrho_cm / aperture_cm)
}

/// Colour index (B-V) estimate from dust-to-gas ratio.
/// Dusty comets tend to be redder; gas-rich comets bluer.
pub fn color_index_bv(dust_to_gas: f64) -> f64 {
    let solar_bv = 0.65;
    let reddening = 0.1 * (dust_to_gas / (1.0 + dust_to_gas));
    solar_bv + reddening
}