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

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

pub struct PlanetEvolution {
    pub mass: f64,
    pub radius: f64,
    pub age: f64,
    pub semi_major: f64,
    pub stellar_luminosity: f64,
}

impl PlanetEvolution {
    pub fn new(mass: f64, radius: f64, semi_major: f64, stellar_luminosity: f64) -> Self {
        Self {
            mass,
            radius,
            age: 0.0,
            semi_major,
            stellar_luminosity,
        }
    }

    pub fn cooling_luminosity(&self) -> f64 {
        let t_eff = 150.0
            * (JUPITER_MASS / self.mass).powf(0.3)
            * (1e9 * YEAR / (self.age + 1e6 * YEAR)).powf(0.3);
        let area = 4.0 * PI * self.radius * self.radius;
        SIGMA_SB * t_eff.powi(4) * area
    }

    pub fn radius_at_age(&self, age: f64) -> f64 {
        if self.mass > 50.0 * EARTH_MASS {
            self.radius * (1.0 + 0.1 * (1e9 * YEAR / (age + 1e6 * YEAR)).powf(0.1))
        } else {
            self.radius
        }
    }

    pub fn xuv_mass_loss_rate(&self, stellar_xuv_flux: f64) -> f64 {
        let eta = 0.15;
        let r_xuv = self.radius * 1.5;
        let area = PI * r_xuv * r_xuv;
        eta * stellar_xuv_flux * area * self.radius / (G * self.mass)
    }

    pub fn atmospheric_loss_timescale(&self, stellar_xuv_flux: f64) -> f64 {
        let m_atm = 0.01 * self.mass;
        let mdot = self.xuv_mass_loss_rate(stellar_xuv_flux);
        if mdot > 0.0 {
            m_atm / mdot
        } else {
            f64::INFINITY
        }
    }

    pub fn tidal_migration_rate(&self, star_mass: f64, q_star: f64) -> f64 {
        let a = self.semi_major;
        let n = (G * star_mass / (a * a * a)).sqrt();
        -9.0 / 2.0
            * n
            * (self.mass / star_mass)
            * (star_mass / self.mass)
            * (self.radius / a).powi(5)
            / q_star
    }

    pub fn orbital_circularization_time(&self, star_mass: f64, q_planet: f64) -> f64 {
        let a = self.semi_major;
        let n = (G * star_mass / (a * a * a)).sqrt();
        2.0 / 63.0 * q_planet / n * (self.mass / star_mass) * (a / self.radius).powi(5)
    }

    pub fn impact_erosion_fraction(impact_velocity: f64, escape_velocity: f64) -> f64 {
        let ratio = impact_velocity / escape_velocity;
        if ratio > 1.0 {
            0.5 * (1.0 - 1.0 / ratio)
        } else {
            0.0
        }
    }
}

pub fn runaway_greenhouse_distance(stellar_luminosity: f64) -> f64 {
    (stellar_luminosity / (4.0 * PI * 340.0 * SIGMA_SB * 340.0_f64.powi(3))).sqrt()
}

pub fn maximum_greenhouse_distance(stellar_luminosity: f64) -> f64 {
    (stellar_luminosity / (4.0 * PI * 175.0 * SIGMA_SB * 175.0_f64.powi(3))).sqrt()
}

pub fn heavy_bombardment_rate(age: f64) -> f64 {
    let t0 = 4.0e9 * YEAR;
    let peak = 3.9e9 * YEAR;
    if age < t0 {
        1e20 * (-(age - peak).abs() / (0.1e9 * YEAR)).exp()
    } else {
        1e14 * (-age / (1e9 * YEAR)).exp()
    }
}