asteroidsfactory 0.0.1

Asteroid factory — classify, build and catalogue asteroids of any type: near-Earth, main belt, trojan, centaur, binary, rubble pile, metallic, and potentially hazardous.
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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133

use crate::config::parameters::*;

pub struct PotentiallyHazardousAsteroid {
    pub radius: f64,
    pub density: f64,
    pub mass: f64,
    pub semi_major_axis: f64,
    pub eccentricity: f64,
    pub inclination: f64,
    pub moid: f64,
    pub albedo: f64,
    pub absolute_mag: f64,
    pub velocity_infinity: f64,
}

impl PotentiallyHazardousAsteroid {
    pub fn new(
        radius_m: f64,
        density: f64,
        semi_major_au: f64,
        eccentricity: f64,
        moid_au: f64,
    ) -> Self {
        let radius = radius_m.clamp(70.0, MAX_ASTEROID_RADIUS);
        let density = density.clamp(500.0, 8000.0);
        let mass = sphere_mass(radius, density);
        let ecc = eccentricity.clamp(0.0, 0.999);
        let moid = moid_au.clamp(0.0, 0.05);
        let albedo: f64 = 0.15;
        let diameter_km: f64 = 2.0 * radius / 1000.0;
        let absolute_mag = 15.618 - 5.0 * diameter_km.log10() - 2.5 * albedo.log10();

        let v_circ = (G * SOLAR_MASS / (1.0 * AU)).sqrt();
        let v_ast = vis_viva(SOLAR_MASS, 1.0 * AU, semi_major_au * AU);
        let velocity_infinity = (v_ast - v_circ).abs();

        Self {
            radius,
            density,
            mass,
            semi_major_axis: semi_major_au,
            eccentricity: ecc,
            inclination: 0.0,
            moid,
            albedo,
            absolute_mag,
            velocity_infinity,
        }
    }

    pub fn with_inclination(mut self, inc_deg: f64) -> Self {
        self.inclination = inc_deg.to_radians();
        self
    }

    pub fn with_albedo(mut self, a: f64) -> Self {
        self.albedo = a.clamp(0.01, 0.99);
        let diameter_km = 2.0 * self.radius / 1000.0;
        self.absolute_mag = 15.618 - 5.0 * diameter_km.log10() - 2.5 * self.albedo.log10();
        self
    }

    pub fn is_pha(&self) -> bool {
        self.moid < 0.05 && self.absolute_mag < 22.0
    }

    pub fn impact_energy_joules(&self) -> f64 {
        0.5 * self.mass * self.velocity_infinity * self.velocity_infinity
    }

    pub fn impact_energy_megatons(&self) -> f64 {
        self.impact_energy_joules() / 4.184e15
    }

    pub fn torino_scale_estimate(&self) -> u8 {
        let e_mt = self.impact_energy_megatons();
        if self.moid > 0.05 {
            0
        } else if e_mt < 1.0 {
            1
        } else if e_mt < 100.0 {
            3
        } else if e_mt < 1e5 {
            5
        } else {
            8
        }
    }

    pub fn crater_diameter_estimate(&self) -> f64 {
        let d = 2.0 * self.radius;
        let v = self.velocity_infinity;
        1.161 * d.powf(0.78) * v.powf(0.44) * self.density.powf(0.33) / 2500.0_f64.powf(0.33)
    }

    pub fn surface_gravity(&self) -> f64 {
        surface_gravity(self.mass, self.radius)
    }

    pub fn escape_velocity(&self) -> f64 {
        escape_velocity(self.mass, self.radius)
    }

    pub fn orbital_period(&self) -> f64 {
        orbital_period(self.semi_major_axis, SOLAR_MASS)
    }

    pub fn orbital_period_years(&self) -> f64 {
        self.orbital_period() / YEAR
    }

    pub fn perihelion_au(&self) -> f64 {
        self.semi_major_axis * (1.0 - self.eccentricity)
    }

    pub fn aphelion_au(&self) -> f64 {
        self.semi_major_axis * (1.0 + self.eccentricity)
    }

    pub fn warning_time_estimate(&self, detection_distance_au: f64) -> f64 {
        let v_approach = self.velocity_infinity;
        (detection_distance_au * AU) / v_approach
    }

    pub fn deflection_delta_v(&self, lead_time_years: f64) -> f64 {
        let t = lead_time_years * YEAR;
        if t < 1.0 {
            return f64::INFINITY;
        }
        let miss = self.moid * AU;
        miss / t
    }
}