planetsfactory 0.0.3

Planet factory — classify, build and catalogue planets for any star system: Solar System, TRAPPIST-1, Kepler-90, Proxima Centauri, or fully custom worlds.
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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

use crate::config::parameters::*;
use std::f64::consts::PI;

#[derive(Clone, Copy)]
pub enum FormationChannel {
    CoreAccretion,
    DiskInstability,
    PebbleAccretion,
}

pub struct FormationModel {
    pub channel: FormationChannel,
    pub disk_mass: f64,
    pub stellar_mass: f64,
    pub metallicity: f64,
    pub snow_line_au: f64,
}

impl FormationModel {
    pub fn new(
        channel: FormationChannel,
        disk_mass: f64,
        stellar_mass: f64,
        metallicity: f64,
        snow_line_au: f64,
    ) -> Self {
        Self {
            channel,
            disk_mass,
            stellar_mass,
            metallicity,
            snow_line_au,
        }
    }

    pub fn solar_like() -> Self {
        Self {
            channel: FormationChannel::CoreAccretion,
            disk_mass: 0.01 * SOLAR_MASS,
            stellar_mass: SOLAR_MASS,
            metallicity: 0.02,
            snow_line_au: 2.7,
        }
    }

    pub fn snow_line_radius(&self) -> f64 {
        self.snow_line_au * AU
    }

    pub fn disk_surface_density(&self, r_au: f64) -> f64 {
        let sigma_0 = self.disk_mass / (2.0 * PI * (AU * AU) * 100.0);
        sigma_0 * r_au.powf(-1.5)
    }

    pub fn isolation_mass(&self, r_au: f64) -> f64 {
        let sigma = self.disk_surface_density(r_au);
        let a = r_au * AU;
        let feeding_zone = 10.0 * hill_radius(a, EARTH_MASS, self.stellar_mass);
        2.0 * PI * a * feeding_zone * sigma
    }

    pub fn critical_core_mass(&self) -> f64 {
        10.0 * EARTH_MASS
    }

    pub fn runaway_gas_accretion_rate(&self, core_mass: f64) -> f64 {
        let m_crit = self.critical_core_mass();
        if core_mass < m_crit {
            return 0.0;
        }
        let excess = (core_mass - m_crit) / m_crit;
        1e-3 * EARTH_MASS / YEAR * excess.powi(2)
    }

    pub fn pebble_accretion_rate(&self, planet_mass: f64, r_au: f64) -> f64 {
        let a = r_au * AU;
        let v_k = orbital_velocity(self.stellar_mass, a);
        let r_h = hill_radius(a, planet_mass, self.stellar_mass);
        let pebble_flux = 1e-4 * EARTH_MASS / YEAR;
        let tau_s = 0.1;
        let r_acc = r_h * (tau_s / 0.1_f64).sqrt();
        pebble_flux * (2.0 * r_acc / a).min(1.0) * (v_k / (v_k + 1.0e-30))
    }

    pub fn disk_instability_mass(&self, r_au: f64) -> f64 {
        let sigma = self.disk_surface_density(r_au);
        let a = r_au * AU;
        let cs = 1000.0;
        let q = cs * (G * self.stellar_mass / (a * a * a)).sqrt() / (PI * G * sigma);
        if q < 1.0 {
            PI * sigma * a * a / 4.0
        } else {
            0.0
        }
    }

    pub fn formation_timescale(&self, r_au: f64) -> f64 {
        match self.channel {
            FormationChannel::CoreAccretion => 1e5 * YEAR * r_au.powf(3.0),
            FormationChannel::DiskInstability => 1e3 * YEAR,
            FormationChannel::PebbleAccretion => 1e4 * YEAR * r_au.powf(1.5),
        }
    }
}

pub fn minimum_mass_solar_nebula_density(r_au: f64) -> f64 {
    1700.0 * r_au.powf(-1.5)
}

pub fn oligarchic_growth_time(planet_mass: f64, sigma: f64, r_au: f64) -> f64 {
    let a = r_au * AU;
    let rho_s = 3000.0;
    let r_p = (3.0 * planet_mass / (4.0 * PI * rho_s)).cbrt();
    let v_esc = (2.0 * G * planet_mass / r_p).sqrt();
    let omega = (G * SOLAR_MASS / (a * a * a)).sqrt();
    planet_mass / (PI * r_p * r_p * sigma * omega * (v_esc / (a * omega)).powi(2))
}