suns 0.0.4

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
153
154
155
156
157
158
159
160
161
162
163
164
165
166

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

use crate::{ELECTRON_CHARGE, ELECTRON_MASS, EPSILON_0, MU_0, PROTON_MASS, SOLAR_RADIUS};

pub const SLOW_WIND_SPEED: f64 = 4.0e5;
pub const FAST_WIND_SPEED: f64 = 8.0e5;
pub const WIND_DENSITY_1AU: f64 = 7e6;
pub const WIND_TEMPERATURE_1AU: f64 = 1e5;
pub const WIND_MAGNETIC_FIELD_1AU: f64 = 5e-9;
pub const PARKER_SPIRAL_OMEGA: f64 = 2.87e-6;
pub const SOLAR_WIND_MASS_LOSS_RATE: f64 = 2e-14;
pub const HELIOSPHERIC_CURRENT_SHEET_TILT_DEG: f64 = 7.25;
pub const TERMINATION_SHOCK_AU: f64 = 85.0;
pub const HELIOPAUSE_AU: f64 = 120.0;

pub struct SolarWindState {
    pub speed_ms: f64,
    pub density_m3: f64,
    pub temperature_k: f64,
    pub magnetic_field_t: f64,
    pub distance_au: f64,
}

impl SolarWindState {
    pub fn slow_wind_1au() -> Self {
        Self {
            speed_ms: SLOW_WIND_SPEED,
            density_m3: WIND_DENSITY_1AU,
            temperature_k: WIND_TEMPERATURE_1AU,
            magnetic_field_t: WIND_MAGNETIC_FIELD_1AU,
            distance_au: 1.0,
        }
    }

    pub fn fast_wind_1au() -> Self {
        Self {
            speed_ms: FAST_WIND_SPEED,
            density_m3: 3e6,
            temperature_k: 2e5,
            magnetic_field_t: WIND_MAGNETIC_FIELD_1AU,
            distance_au: 1.0,
        }
    }

    pub fn at_distance(&self, distance_au: f64) -> Self {
        let r_ratio = self.distance_au / distance_au;
        Self {
            speed_ms: self.speed_ms,
            density_m3: self.density_m3 * r_ratio * r_ratio,
            temperature_k: self.temperature_k * r_ratio.powf(0.7),
            magnetic_field_t: self.magnetic_field_t * r_ratio,
            distance_au,
        }
    }

    pub fn dynamic_pressure(&self) -> f64 {
        0.5 * self.density_m3 * PROTON_MASS * self.speed_ms * self.speed_ms
    }

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

    pub fn thermal_pressure(&self) -> f64 {
        self.density_m3 * K_B * self.temperature_k
    }

    pub fn magnetic_pressure(&self) -> f64 {
        self.magnetic_field_t.powi(2) / (2.0 * MU_0)
    }

    pub fn total_pressure(&self) -> f64 {
        self.dynamic_pressure() + self.thermal_pressure() + self.magnetic_pressure()
    }

    pub fn kinetic_energy_flux(&self) -> f64 {
        0.5 * self.density_m3 * PROTON_MASS * self.speed_ms.powi(3)
    }

    pub fn mass_flux(&self) -> f64 {
        self.density_m3 * PROTON_MASS * self.speed_ms
    }

    pub fn alfven_mach_number(&self) -> f64 {
        let va = self.magnetic_field_t / (MU_0 * self.density_m3 * PROTON_MASS).sqrt();
        self.speed_ms / va
    }

    pub fn sonic_mach_number(&self) -> f64 {
        let cs = (5.0 / 3.0 * K_B * self.temperature_k / PROTON_MASS).sqrt();
        self.speed_ms / cs
    }

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

    pub fn proton_gyroradius(&self) -> f64 {
        let v_th = (2.0 * K_B * self.temperature_k / PROTON_MASS).sqrt();
        PROTON_MASS * v_th / (ELECTRON_CHARGE * self.magnetic_field_t)
    }

    pub fn debye_length(&self) -> f64 {
        (EPSILON_0 * K_B * self.temperature_k
            / (self.density_m3 * ELECTRON_CHARGE * ELECTRON_CHARGE))
            .sqrt()
    }

    pub fn plasma_frequency(&self) -> f64 {
        (self.density_m3 * ELECTRON_CHARGE * ELECTRON_CHARGE / (EPSILON_0 * ELECTRON_MASS)).sqrt()
    }

    pub fn travel_time_to_au(&self, target_au: f64) -> f64 {
        let distance = (target_au - self.distance_au).abs() * AU;
        distance / self.speed_ms
    }
}

pub fn parker_critical_radius() -> f64 {
    let mu = 0.6;
    let cs = (K_B * 1e6 / (mu * PROTON_MASS)).sqrt();
    G * SOLAR_MASS / (2.0 * cs * cs)
}

pub fn parker_solution_speed(r_m: f64) -> f64 {
    let rc = parker_critical_radius();
    let cs = (K_B * 1e6 / (0.6 * PROTON_MASS)).sqrt();
    let ratio = r_m / rc;
    if ratio <= 1.0 {
        cs * (ratio).powf(0.5)
    } else {
        cs * (1.0 + 2.0 * (ratio.ln())).sqrt()
    }
}

pub fn solar_wind_mass_loss_rate() -> f64 {
    let r_1au = AU;
    let area = 4.0 * std::f64::consts::PI * r_1au * r_1au;
    WIND_DENSITY_1AU * PROTON_MASS * SLOW_WIND_SPEED * area
}

pub fn alfven_radius() -> f64 {
    15.0 * SOLAR_RADIUS
}

pub fn heliospheric_magnetic_field_at(distance_au: f64) -> f64 {
    WIND_MAGNETIC_FIELD_1AU / distance_au
}

pub fn bow_shock_standoff(planet_dipole_moment: f64, distance_au: f64) -> f64 {
    let wind = SolarWindState::slow_wind_1au().at_distance(distance_au);
    let p_dyn = wind.dynamic_pressure();
    (planet_dipole_moment / (4.0 * std::f64::consts::PI * p_dyn / MU_0)).powf(1.0 / 6.0)
}

pub fn interplanetary_shock_speed(upstream_speed: f64, compression_ratio: f64) -> f64 {
    upstream_speed * compression_ratio / (compression_ratio - 1.0)
}

pub fn corotation_radius() -> f64 {
    let omega = PARKER_SPIRAL_OMEGA;
    let v = SLOW_WIND_SPEED;
    v / omega
}