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
167
168
169
170
171
172
173
174
175
176

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

use crate::{
    ELECTRON_CHARGE, ELECTRON_MASS, MU_0, PROTON_MASS, SECONDS_PER_DAY, SOLAR_CYCLE_YEARS,
    SOLAR_RADIUS, SOLAR_ROTATION_EQUATORIAL_DAYS, SOLAR_ROTATION_POLAR_DAYS,
};

pub const QUIET_SUN_FIELD_T: f64 = 1e-4;
pub const ACTIVE_REGION_FIELD_T: f64 = 0.1;
pub const SUNSPOT_FIELD_T: f64 = 0.3;
pub const CORE_FIELD_ESTIMATE_T: f64 = 30.0;
pub const TACHOCLINE_FIELD_T: f64 = 10.0;
pub const MAGNETIC_REYNOLDS_NUMBER: f64 = 1e10;
pub const PHOTOSPHERIC_DIFFUSIVITY: f64 = 1e3;

pub struct MagneticFieldRegion {
    pub name: &'static str,
    pub field_strength_t: f64,
    pub temperature_k: f64,
    pub density_kgm3: f64,
    pub scale_length_m: f64,
}

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

    pub fn thermal_pressure(&self) -> f64 {
        let n = self.density_kgm3 / PROTON_MASS;
        n * K_B * self.temperature_k
    }

    pub fn plasma_beta(&self) -> f64 {
        self.thermal_pressure() / self.magnetic_pressure().max(1e-30)
    }

    pub fn alfven_speed(&self) -> f64 {
        self.field_strength_t / (MU_0 * self.density_kgm3).sqrt()
    }

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

    pub fn magnetic_diffusion_time(&self) -> f64 {
        MU_0 * self.scale_length_m.powi(2) / PHOTOSPHERIC_DIFFUSIVITY
    }

    pub fn magnetic_reynolds_at_scale(&self, velocity: f64) -> f64 {
        velocity * self.scale_length_m / PHOTOSPHERIC_DIFFUSIVITY
    }
}

pub fn quiet_sun() -> MagneticFieldRegion {
    MagneticFieldRegion {
        name: "Quiet Sun",
        field_strength_t: QUIET_SUN_FIELD_T,
        temperature_k: 5778.0,
        density_kgm3: 2e-4,
        scale_length_m: 1e6,
    }
}

pub fn active_region() -> MagneticFieldRegion {
    MagneticFieldRegion {
        name: "Active Region",
        field_strength_t: ACTIVE_REGION_FIELD_T,
        temperature_k: 5000.0,
        density_kgm3: 3e-4,
        scale_length_m: 3e7,
    }
}

pub fn sunspot_umbra() -> MagneticFieldRegion {
    MagneticFieldRegion {
        name: "Sunspot Umbra",
        field_strength_t: SUNSPOT_FIELD_T,
        temperature_k: 3700.0,
        density_kgm3: 5e-4,
        scale_length_m: 1e7,
    }
}

pub fn coronal_field() -> MagneticFieldRegion {
    MagneticFieldRegion {
        name: "Corona",
        field_strength_t: 1e-4,
        temperature_k: 1e6,
        density_kgm3: 1e-12,
        scale_length_m: 1e8,
    }
}

pub fn differential_rotation_rate(latitude_deg: f64) -> f64 {
    let a = 14.713;
    let b = -2.396;
    let c = -1.787;
    let sin_lat = latitude_deg.to_radians().sin();
    let sin2 = sin_lat * sin_lat;
    let sin4 = sin2 * sin2;
    (a + b * sin2 + c * sin4) * std::f64::consts::PI / 180.0 / SECONDS_PER_DAY
}

pub fn differential_rotation_period_days(latitude_deg: f64) -> f64 {
    let omega = differential_rotation_rate(latitude_deg);
    2.0 * std::f64::consts::PI / (omega * SECONDS_PER_DAY)
}

pub fn omega_effect_shear() -> f64 {
    let omega_eq = 2.0 * std::f64::consts::PI / (SOLAR_ROTATION_EQUATORIAL_DAYS * SECONDS_PER_DAY);
    let omega_pole = 2.0 * std::f64::consts::PI / (SOLAR_ROTATION_POLAR_DAYS * SECONDS_PER_DAY);
    (omega_eq - omega_pole).abs()
}

pub fn dynamo_number(alpha: f64, shear: f64, diffusivity: f64, length: f64) -> f64 {
    alpha * shear * length.powi(3) / (diffusivity * diffusivity)
}

pub fn solar_dynamo_number_estimate() -> f64 {
    let alpha = 1.0;
    let shear = omega_effect_shear();
    let eta = PHOTOSPHERIC_DIFFUSIVITY;
    let l = 0.3 * SOLAR_RADIUS;
    dynamo_number(alpha, shear, eta, l)
}

pub fn magnetic_flux(field_t: f64, area_m2: f64) -> f64 {
    field_t * area_m2
}

pub fn magnetic_energy(field_t: f64, volume_m3: f64) -> f64 {
    field_t.powi(2) * volume_m3 / (2.0 * MU_0)
}

pub fn larmor_radius(velocity: f64, mass: f64, charge: f64, field_t: f64) -> f64 {
    mass * velocity / (charge.abs() * field_t)
}

pub fn proton_larmor_radius(temperature_k: f64, field_t: f64) -> f64 {
    let v_thermal = (2.0 * K_B * temperature_k / PROTON_MASS).sqrt();
    larmor_radius(v_thermal, PROTON_MASS, ELECTRON_CHARGE, field_t)
}

pub fn electron_larmor_radius(temperature_k: f64, field_t: f64) -> f64 {
    let v_thermal = (2.0 * K_B * temperature_k / ELECTRON_MASS).sqrt();
    larmor_radius(v_thermal, ELECTRON_MASS, ELECTRON_CHARGE, field_t)
}

pub fn cyclotron_frequency(charge: f64, mass: f64, field_t: f64) -> f64 {
    charge.abs() * field_t / mass
}

pub fn proton_cyclotron_frequency(field_t: f64) -> f64 {
    cyclotron_frequency(ELECTRON_CHARGE, PROTON_MASS, field_t)
}

pub fn electron_cyclotron_frequency(field_t: f64) -> f64 {
    cyclotron_frequency(ELECTRON_CHARGE, ELECTRON_MASS, field_t)
}

pub fn sweet_parker_reconnection_rate(alfven_speed: f64, lundquist_number: f64) -> f64 {
    alfven_speed / lundquist_number.sqrt()
}

pub fn magnetic_helicity(field_t: f64, length_m: f64) -> f64 {
    field_t * field_t * length_m.powi(3)
}

pub fn force_free_parameter(current_density: f64, field_t: f64) -> f64 {
    MU_0 * current_density / field_t.max(1e-30)
}

pub fn hale_cycle_period() -> f64 {
    2.0 * SOLAR_CYCLE_YEARS
}