suns 0.0.3

Sun celestial simulation crate for the MilkyWay SolarSystem workspace
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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152

use sciforge::hub::domain::common::constants::{AU, K_B};

use crate::plasma::solar_wind::{PARKER_SPIRAL_OMEGA, SolarWindState, WIND_MAGNETIC_FIELD_1AU};
use crate::{MU_0, PROTON_MASS};

pub const EARTH_DIPOLE_MOMENT: f64 = 8.22e22;
pub const JUPITER_DIPOLE_MOMENT: f64 = 1.56e27;
pub const SATURN_DIPOLE_MOMENT: f64 = 4.6e25;
pub const MERCURY_DIPOLE_MOMENT: f64 = 4.0e19;
pub const URANUS_DIPOLE_MOMENT: f64 = 3.9e24;
pub const NEPTUNE_DIPOLE_MOMENT: f64 = 2.2e24;
pub const GANYMEDE_DIPOLE_MOMENT: f64 = 1.3e20;

pub const EARTH_MAGNETOPAUSE_R: f64 = 10.0;
pub const JUPITER_MAGNETOPAUSE_R: f64 = 75.0;

pub struct SolarWindInteraction {
    pub planet_dipole_moment: f64,
    pub planet_radius_m: f64,
    pub distance_au: f64,
}

impl SolarWindInteraction {
    pub fn new(planet_dipole_moment: f64, planet_radius_m: f64, distance_au: f64) -> Self {
        Self {
            planet_dipole_moment,
            planet_radius_m,
            distance_au,
        }
    }

    pub fn earth() -> Self {
        Self::new(EARTH_DIPOLE_MOMENT, 6.371e6, 1.0)
    }

    pub fn jupiter() -> Self {
        Self::new(JUPITER_DIPOLE_MOMENT, 7.1492e7, 5.2)
    }

    pub fn mercury() -> Self {
        Self::new(MERCURY_DIPOLE_MOMENT, 2.4397e6, 0.387)
    }

    pub fn wind_state(&self) -> SolarWindState {
        SolarWindState::slow_wind_1au().at_distance(self.distance_au)
    }

    pub fn ram_pressure(&self) -> f64 {
        self.wind_state().dynamic_pressure()
    }

    pub fn magnetopause_distance(&self) -> f64 {
        let p_dyn = self.ram_pressure();
        (self.planet_dipole_moment.powi(2) * MU_0 / (8.0 * std::f64::consts::PI.powi(2) * p_dyn))
            .powf(1.0 / 6.0)
    }

    pub fn magnetopause_distance_in_radii(&self) -> f64 {
        self.magnetopause_distance() / self.planet_radius_m
    }

    pub fn bow_shock_distance(&self) -> f64 {
        1.3 * self.magnetopause_distance()
    }

    pub fn bow_shock_distance_in_radii(&self) -> f64 {
        self.bow_shock_distance() / self.planet_radius_m
    }

    pub fn magnetic_pressure_at_magnetopause(&self) -> f64 {
        let r_mp = self.magnetopause_distance();
        let b = self.planet_dipole_moment * MU_0 / (4.0 * std::f64::consts::PI * r_mp.powi(3));
        b * b / (2.0 * MU_0)
    }

    pub fn magnetosphere_cross_section(&self) -> f64 {
        let r_mp = self.magnetopause_distance();
        std::f64::consts::PI * r_mp * r_mp
    }

    pub fn intercepted_power(&self) -> f64 {
        let wind = self.wind_state();
        let ke_flux = wind.kinetic_energy_flux();
        ke_flux * self.magnetosphere_cross_section()
    }

    pub fn reconnection_rate(&self) -> f64 {
        let wind = self.wind_state();
        let v_a = wind.magnetic_field_t / (MU_0 * wind.density_m3 * PROTON_MASS).sqrt();
        0.1 * v_a
    }

    pub fn polar_cap_potential(&self) -> f64 {
        let wind = self.wind_state();
        let r_mp = self.magnetopause_distance();
        wind.speed_ms * wind.magnetic_field_t * 2.0 * r_mp * 0.1
    }

    pub fn imf_strength(&self) -> f64 {
        WIND_MAGNETIC_FIELD_1AU / self.distance_au
    }

    pub fn parker_spiral_angle_deg(&self) -> f64 {
        let r_m = self.distance_au * AU;
        let wind = self.wind_state();
        (PARKER_SPIRAL_OMEGA * r_m / wind.speed_ms)
            .atan()
            .to_degrees()
    }
}

pub fn magnetopause_standoff(dipole_moment: f64, dynamic_pressure: f64) -> f64 {
    (dipole_moment.powi(2) * MU_0 / (8.0 * std::f64::consts::PI.powi(2) * dynamic_pressure))
        .powf(1.0 / 6.0)
}

pub fn solar_wind_ram_pressure_at(distance_au: f64) -> f64 {
    SolarWindState::slow_wind_1au()
        .at_distance(distance_au)
        .dynamic_pressure()
}

pub fn has_magnetosphere(dipole_moment: f64, planet_radius_m: f64, distance_au: f64) -> bool {
    let p_dyn = solar_wind_ram_pressure_at(distance_au);
    let r_mp = magnetopause_standoff(dipole_moment, p_dyn);
    r_mp > planet_radius_m
}

pub fn atmospheric_ion_escape_rate(
    exosphere_temp_k: f64,
    particle_mass_kg: f64,
    planet_radius_m: f64,
    escape_velocity_ms: f64,
    distance_au: f64,
) -> f64 {
    let v_th = (2.0 * K_B * exosphere_temp_k / particle_mass_kg).sqrt();
    let lambda = (escape_velocity_ms / v_th).powi(2);
    let jeans = (1.0 + lambda) * (-lambda).exp();
    let wind = SolarWindState::slow_wind_1au().at_distance(distance_au);
    let sputtering_factor = wind.density_m3 * wind.speed_ms / 1e10;
    let area = 4.0 * std::f64::consts::PI * planet_radius_m * planet_radius_m;
    let n_exo = 1e12;
    n_exo * v_th * area * (jeans + sputtering_factor)
}

pub fn energy_input_to_magnetosphere(dipole_moment: f64, distance_au: f64) -> f64 {
    let wind = SolarWindState::slow_wind_1au().at_distance(distance_au);
    let r_mp = magnetopause_standoff(dipole_moment, wind.dynamic_pressure());
    let area = std::f64::consts::PI * r_mp * r_mp;
    let epsilon = 0.1;
    epsilon * wind.kinetic_energy_flux() * area
}