suns 0.0.1

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

use sciforge::hub::domain::common::constants::AU;

use crate::{
    CME_MASS_RANGE_HIGH, CME_MASS_RANGE_LOW, CME_SPEED_FAST, CME_SPEED_SLOW, MU_0,
};

pub const CME_RATE_SOLAR_MIN_PER_DAY: f64 = 0.5;
pub const CME_RATE_SOLAR_MAX_PER_DAY: f64 = 4.0;
pub const CME_ANGULAR_WIDTH_DEG: f64 = 60.0;
pub const HALO_CME_FRACTION: f64 = 0.04;
pub const CME_MAGNETIC_CLOUD_FIELD_T: f64 = 2e-8;
pub const CME_SHEATH_COMPRESSION: f64 = 4.0;
pub const ICME_DURATION_HOURS: f64 = 24.0;
pub const DST_STORM_THRESHOLD_NT: f64 = -50.0;

pub struct CoronalMassEjection {
    pub speed_ms: f64,
    pub mass_kg: f64,
    pub angular_width_deg: f64,
    pub magnetic_field_t: f64,
    pub source_latitude_deg: f64,
}

impl CoronalMassEjection {
    pub fn slow() -> Self {
        Self {
            speed_ms: CME_SPEED_SLOW,
            mass_kg: CME_MASS_RANGE_LOW,
            angular_width_deg: 40.0,
            magnetic_field_t: 1e-8,
            source_latitude_deg: 30.0,
        }
    }

    pub fn fast() -> Self {
        Self {
            speed_ms: CME_SPEED_FAST,
            mass_kg: CME_MASS_RANGE_HIGH,
            angular_width_deg: 90.0,
            magnetic_field_t: 5e-8,
            source_latitude_deg: 15.0,
        }
    }

    pub fn halo() -> Self {
        Self {
            speed_ms: 2.0e6,
            mass_kg: 5e12,
            angular_width_deg: 360.0,
            magnetic_field_t: 3e-8,
            source_latitude_deg: 10.0,
        }
    }

    pub fn kinetic_energy(&self) -> f64 {
        0.5 * self.mass_kg * self.speed_ms * self.speed_ms
    }

    pub fn magnetic_energy(&self) -> f64 {
        let volume = self.volume_at_1au();
        self.magnetic_field_t.powi(2) * volume / (2.0 * MU_0)
    }

    pub fn total_energy(&self) -> f64 {
        self.kinetic_energy() + self.magnetic_energy()
    }

    pub fn travel_time_to_earth(&self) -> f64 {
        AU / self.speed_ms
    }

    pub fn travel_time_to_au(&self, distance_au: f64) -> f64 {
        distance_au * AU / self.speed_ms
    }

    pub fn momentum(&self) -> f64 {
        self.mass_kg * self.speed_ms
    }

    pub fn solid_angle_sr(&self) -> f64 {
        2.0 * std::f64::consts::PI
            * (1.0 - (self.angular_width_deg / 2.0).to_radians().cos())
    }

    pub fn volume_at_1au(&self) -> f64 {
        let r = AU;
        self.solid_angle_sr() * r.powi(3) / 3.0
    }

    pub fn density_at_1au(&self) -> f64 {
        self.mass_kg / self.volume_at_1au()
    }

    pub fn dynamic_pressure_at_1au(&self) -> f64 {
        0.5 * self.density_at_1au() * self.speed_ms * self.speed_ms
    }

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

    pub fn is_supercritical(&self) -> bool {
        self.alfven_mach_number() > 1.0
    }

    pub fn sheath_thickness_m(&self) -> f64 {
        AU * 0.01
    }

    pub fn magnetic_cloud_radius_m(&self) -> f64 {
        0.1 * AU
    }

    pub fn dst_index_estimate(&self) -> f64 {
        let v = self.speed_ms / 1e3;
        let bz = self.magnetic_field_t * 1e9;
        -0.01 * v * bz
    }

    pub fn is_geoeffective(&self) -> bool {
        self.dst_index_estimate() < DST_STORM_THRESHOLD_NT
    }

    pub fn forbush_decrease_percent(&self) -> f64 {
        let e = self.kinetic_energy();
        (e / 1e24).min(20.0)
    }
}

pub fn cme_rate_at_cycle_phase(phase: f64) -> f64 {
    let min = CME_RATE_SOLAR_MIN_PER_DAY;
    let max = CME_RATE_SOLAR_MAX_PER_DAY;
    min + (max - min) * (std::f64::consts::PI * phase).sin().powi(2)
}

pub fn snow_plow_speed(cme_mass: f64, cme_speed: f64, ambient_density: f64, distance: f64) -> f64 {
    let swept_mass = ambient_density * 4.0 * std::f64::consts::PI * distance.powi(3) / 3.0;
    cme_mass * cme_speed / (cme_mass + swept_mass)
}

pub fn drag_based_cme_speed(
    initial_speed: f64,
    solar_wind_speed: f64,
    drag_coeff: f64,
    time_s: f64,
) -> f64 {
    let dv = initial_speed - solar_wind_speed;
    solar_wind_speed + dv * (-drag_coeff * time_s).exp()
}

pub fn cme_arrival_probability(angular_width_deg: f64) -> f64 {
    (angular_width_deg / 360.0).min(1.0)
}

pub fn interplanetary_shock_compression(mach: f64) -> f64 {
    let gamma = 5.0 / 3.0;
    (gamma + 1.0) * mach * mach / ((gamma - 1.0) * mach * mach + 2.0)
}