dwarfplanetsfactory 0.0.3

Dwarf planet factory — classify, build and catalogue dwarf planets of any type: Kuiper belt, scattered disk, plutino, cold classical, detached, binary, Ceres-type, and sednoid.
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

use crate::config::parameters;

/// Mean motion (rad/s) from Kepler's third law
pub fn mean_motion(semi_major_au: f64) -> f64 {
    let period = parameters::orbital_period(semi_major_au, parameters::SOLAR_MASS);
    if period < 1.0 {
        return 0.0;
    }
    2.0 * parameters::PI / period
}

/// Orbital period in Julian years
pub fn period_years(semi_major_au: f64) -> f64 {
    parameters::orbital_period(semi_major_au, parameters::SOLAR_MASS) / parameters::YEAR
}

/// Tisserand parameter relative to Neptune
pub fn tisserand_neptune(semi_major_au: f64, eccentricity: f64, inclination_deg: f64) -> f64 {
    let a_n = parameters::NEPTUNE_A;
    let a_ratio = a_n / semi_major_au;
    let inc = inclination_deg.to_radians();
    a_ratio + 2.0 * ((semi_major_au / a_n) * (1.0 - eccentricity * eccentricity)).sqrt() * inc.cos()
}

/// Tisserand parameter relative to Jupiter (a_j ≈ 5.2 AU)
pub fn tisserand_jupiter(semi_major_au: f64, eccentricity: f64, inclination_deg: f64) -> f64 {
    let a_j: f64 = 5.2;
    let a_ratio = a_j / semi_major_au;
    let inc = inclination_deg.to_radians();
    a_ratio + 2.0 * ((semi_major_au / a_j) * (1.0 - eccentricity * eccentricity)).sqrt() * inc.cos()
}

/// Solve Kepler's equation M = E - e sin(E) by Newton iteration
pub fn kepler_solve(mean_anomaly: f64, eccentricity: f64) -> f64 {
    let mut ecc_anom = mean_anomaly;
    for _ in 0..50 {
        let de = (ecc_anom - eccentricity * ecc_anom.sin() - mean_anomaly)
            / (1.0 - eccentricity * ecc_anom.cos());
        ecc_anom -= de;
        if de.abs() < 1.0e-12 {
            break;
        }
    }
    ecc_anom
}

/// True anomaly from eccentric anomaly
pub fn true_anomaly(ecc_anomaly: f64, eccentricity: f64) -> f64 {
    2.0 * ((1.0 + eccentricity).sqrt() * (ecc_anomaly / 2.0).sin())
        .atan2((1.0 - eccentricity).sqrt() * (ecc_anomaly / 2.0).cos())
}

/// Heliocentric distance (AU) at a given true anomaly
pub fn heliocentric_distance(semi_major_au: f64, eccentricity: f64, true_anom: f64) -> f64 {
    semi_major_au * (1.0 - eccentricity * eccentricity) / (1.0 + eccentricity * true_anom.cos())
}

/// Synodic period (s) between the body and Neptune
pub fn synodic_period_neptune(semi_major_au: f64) -> f64 {
    let n_body = mean_motion(semi_major_au);
    let n_neptune = mean_motion(parameters::NEPTUNE_A);
    let dn = (n_body - n_neptune).abs();
    if dn < 1.0e-20 {
        return f64::INFINITY;
    }
    2.0 * parameters::PI / dn
}

/// Kozai-Lidov timescale (s) for a binary perturbed by the Sun
pub fn kozai_timescale(a_binary_m: f64, a_helio_au: f64, m_primary: f64) -> f64 {
    let a_sun = a_helio_au * parameters::AU;
    if a_binary_m < 1.0 || m_primary < 1.0 {
        return 0.0;
    }
    let p_bin = 2.0 * parameters::PI * (a_binary_m.powi(3) / (parameters::G * m_primary)).sqrt();
    let p_sun =
        2.0 * parameters::PI * (a_sun.powi(3) / (parameters::G * parameters::SOLAR_MASS)).sqrt();
    if p_bin < 1.0 {
        return 0.0;
    }
    (p_sun * p_sun / p_bin)
        * (parameters::SOLAR_MASS / m_primary)
        * (1.0 - (0.3_f64).powi(2)).powf(1.5) // assume e_outer ~ 0.3
}