1use crate::constants::{
2 EARTH_MASS, EARTH_ORBITAL_VELOCITY, EARTH_RADIUS, G, JUPITER_AXIAL_TILT, JUPITER_BOND_ALBEDO,
3 JUPITER_ECCENTRICITY, JUPITER_ESCAPE_VELOCITY, JUPITER_FLATTENING, JUPITER_INCLINATION,
4 JUPITER_MASS, JUPITER_MEAN_DENSITY, JUPITER_ORBITAL_PERIOD, JUPITER_ORBITAL_VELOCITY,
5 JUPITER_RADIUS, JUPITER_SEMI_MAJOR_AXIS, JUPITER_SIDEREAL_DAY, JUPITER_SURFACE_GRAVITY, K_B,
6 MARS_AXIAL_TILT, MARS_BOND_ALBEDO, MARS_ECCENTRICITY, MARS_ESCAPE_VELOCITY, MARS_FLATTENING,
7 MARS_INCLINATION, MARS_MASS, MARS_MEAN_DENSITY, MARS_ORBITAL_PERIOD, MARS_ORBITAL_VELOCITY,
8 MARS_RADIUS, MARS_SEMI_MAJOR_AXIS, MARS_SIDEREAL_DAY, MARS_SURFACE_GRAVITY, MERCURY_AXIAL_TILT,
9 MERCURY_BOND_ALBEDO, MERCURY_ECCENTRICITY, MERCURY_ESCAPE_VELOCITY, MERCURY_FLATTENING,
10 MERCURY_INCLINATION, MERCURY_MASS, MERCURY_MEAN_DENSITY, MERCURY_ORBITAL_PERIOD,
11 MERCURY_ORBITAL_VELOCITY, MERCURY_RADIUS, MERCURY_SEMI_MAJOR_AXIS, MERCURY_SIDEREAL_DAY,
12 MERCURY_SURFACE_GRAVITY, MU_0, NEPTUNE_AXIAL_TILT, NEPTUNE_BOND_ALBEDO, NEPTUNE_ECCENTRICITY,
13 NEPTUNE_ESCAPE_VELOCITY, NEPTUNE_FLATTENING, NEPTUNE_INCLINATION, NEPTUNE_MASS,
14 NEPTUNE_MEAN_DENSITY, NEPTUNE_ORBITAL_PERIOD, NEPTUNE_ORBITAL_VELOCITY, NEPTUNE_RADIUS,
15 NEPTUNE_SEMI_MAJOR_AXIS, NEPTUNE_SIDEREAL_DAY, NEPTUNE_SURFACE_GRAVITY, SATURN_AXIAL_TILT,
16 SATURN_BOND_ALBEDO, SATURN_ECCENTRICITY, SATURN_ESCAPE_VELOCITY, SATURN_FLATTENING,
17 SATURN_INCLINATION, SATURN_MASS, SATURN_MEAN_DENSITY, SATURN_ORBITAL_PERIOD,
18 SATURN_ORBITAL_VELOCITY, SATURN_RADIUS, SATURN_SEMI_MAJOR_AXIS, SATURN_SIDEREAL_DAY,
19 SATURN_SURFACE_GRAVITY, SIGMA_SB, URANUS_AXIAL_TILT, URANUS_BOND_ALBEDO, URANUS_ECCENTRICITY,
20 URANUS_ESCAPE_VELOCITY, URANUS_FLATTENING, URANUS_INCLINATION, URANUS_MASS,
21 URANUS_MEAN_DENSITY, URANUS_ORBITAL_PERIOD, URANUS_ORBITAL_VELOCITY, URANUS_RADIUS,
22 URANUS_SEMI_MAJOR_AXIS, URANUS_SIDEREAL_DAY, URANUS_SURFACE_GRAVITY, VENUS_AXIAL_TILT,
23 VENUS_BOND_ALBEDO, VENUS_ECCENTRICITY, VENUS_ESCAPE_VELOCITY, VENUS_FLATTENING,
24 VENUS_INCLINATION, VENUS_MASS, VENUS_MEAN_DENSITY, VENUS_ORBITAL_PERIOD,
25 VENUS_ORBITAL_VELOCITY, VENUS_RADIUS, VENUS_SEMI_MAJOR_AXIS, VENUS_SIDEREAL_DAY,
26 VENUS_SURFACE_GRAVITY,
27};
28
29pub fn hydrostatic_pressure(density: f64, g_surface: f64, depth: f64) -> f64 {
30 density * g_surface * depth
31}
32
33pub fn central_pressure(density: f64, g_surface: f64, radius: f64) -> f64 {
34 3.0 * g_surface * density * radius / 2.0
35}
36
37pub fn adiabatic_temperature_gradient(g_local: f64, specific_heat: f64) -> f64 {
38 g_local / specific_heat
39}
40
41pub fn planetary_moment_of_inertia_factor(
42 core_density: f64,
43 mantle_density: f64,
44 core_radius: f64,
45 total_radius: f64,
46) -> f64 {
47 let rc5 = core_radius.powi(5);
48 let rt5 = total_radius.powi(5);
49 let rc3 = core_radius.powi(3);
50 let rt3 = total_radius.powi(3);
51 let numerator = core_density * rc5 + mantle_density * (rt5 - rc5);
52 let denominator = core_density * rc3 + mantle_density * (rt3 - rc3);
53 0.4 * numerator / (denominator * total_radius.powi(2) / rt3)
54}
55
56pub fn love_number_k2(rigidity: f64, density: f64, radius: f64) -> f64 {
57 1.5 / (1.0 + 19.0 * rigidity / (2.0 * G * density.powi(2) * radius.powi(2) / 3.0))
58}
59
60pub fn tidal_heating(
61 radius: f64,
62 eccentricity: f64,
63 mean_motion: f64,
64 k2: f64,
65 tidal_q: f64,
66 perturber_mass: f64,
67 semi_major_axis: f64,
68) -> f64 {
69 21.0 / 2.0 * k2 / tidal_q * G * perturber_mass.powi(2) * radius.powi(5) * mean_motion
70 / semi_major_axis.powi(6)
71 * eccentricity.powi(2)
72}
73
74pub fn tidal_locking_timescale(
75 mass: f64,
76 radius: f64,
77 semi_major_axis: f64,
78 perturber_mass: f64,
79 tidal_q: f64,
80 rigidity: f64,
81 initial_spin: f64,
82) -> f64 {
83 let k2 = love_number_k2(
84 rigidity,
85 mass / (4.0 / 3.0 * std::f64::consts::PI * radius.powi(3)),
86 radius,
87 );
88 let i = 0.4 * mass * radius.powi(2);
89 i * tidal_q * semi_major_axis.powi(6) * initial_spin
90 / (3.0 * k2 * G * perturber_mass.powi(2) * radius.powi(3))
91}
92
93pub fn roche_limit_fluid(primary_radius: f64, primary_density: f64, secondary_density: f64) -> f64 {
94 2.456 * primary_radius * (primary_density / secondary_density).powf(1.0 / 3.0)
95}
96
97pub fn roche_limit_rigid(primary_radius: f64, primary_density: f64, secondary_density: f64) -> f64 {
98 1.26 * primary_radius * (primary_density / secondary_density).powf(1.0 / 3.0)
99}
100
101pub fn equilibrium_temperature(stellar_luminosity: f64, semi_major_axis: f64, albedo: f64) -> f64 {
102 (stellar_luminosity * (1.0 - albedo)
103 / (16.0 * std::f64::consts::PI * SIGMA_SB * semi_major_axis.powi(2)))
104 .powf(0.25)
105}
106
107pub fn greenhouse_surface_temperature(equilibrium_temp: f64, optical_depth_ir: f64) -> f64 {
108 equilibrium_temp * (1.0 + 0.75 * optical_depth_ir).powf(0.25)
109}
110
111pub fn scale_height(temperature: f64, mean_molecular_mass: f64, g_surface: f64) -> f64 {
112 K_B * temperature / (mean_molecular_mass * g_surface)
113}
114
115pub fn atmospheric_escape_jeans(
116 temperature: f64,
117 molecular_mass: f64,
118 escape_velocity: f64,
119) -> f64 {
120 let v_thermal = (2.0 * K_B * temperature / molecular_mass).sqrt();
121 let lambda = (escape_velocity / v_thermal).powi(2);
122 (1.0 + lambda) * (-lambda).exp()
123}
124
125pub fn magnetopause_standoff(dipole_moment: f64, solar_wind_pressure: f64) -> f64 {
126 (MU_0 * dipole_moment.powi(2) / (8.0 * std::f64::consts::PI.powi(2) * solar_wind_pressure))
127 .powf(1.0 / 6.0)
128}
129
130pub fn ring_roche_inner(planet_mass: f64, planet_radius: f64, particle_density: f64) -> f64 {
131 let planet_density = planet_mass / (4.0 / 3.0 * std::f64::consts::PI * planet_radius.powi(3));
132 roche_limit_fluid(planet_radius, planet_density, particle_density)
133}
134
135pub fn ring_optical_depth(
136 surface_density: f64,
137 particle_radius: f64,
138 particle_density: f64,
139) -> f64 {
140 3.0 * surface_density / (4.0 * particle_density * particle_radius)
141}
142
143pub fn synchronous_orbit_radius(mass: f64, rotation_period: f64) -> f64 {
144 let omega = 2.0 * std::f64::consts::PI / rotation_period;
145 (G * mass / omega.powi(2)).powf(1.0 / 3.0)
146}
147
148pub fn oblateness(rotation_rate: f64, equatorial_radius: f64, mass: f64) -> f64 {
149 let q = rotation_rate.powi(2) * equatorial_radius.powi(3) / (G * mass);
150 q / 2.0
151}
152
153pub fn precession_rate(
154 oblateness_j2: f64,
155 planet_radius: f64,
156 semi_major_axis: f64,
157 mean_motion: f64,
158) -> f64 {
159 -1.5 * mean_motion * oblateness_j2 * (planet_radius / semi_major_axis).powi(2)
160}
161
162pub fn bond_albedo_from_geometric(geometric_albedo: f64, phase_integral: f64) -> f64 {
163 geometric_albedo * phase_integral
164}
165
166pub fn thermal_inertia(thermal_conductivity: f64, density: f64, specific_heat: f64) -> f64 {
167 (thermal_conductivity * density * specific_heat).sqrt()
168}
169
170pub fn diurnal_skin_depth(thermal_diffusivity: f64, rotation_period: f64) -> f64 {
171 (thermal_diffusivity * rotation_period / std::f64::consts::PI).sqrt()
172}
173
174pub fn subsolar_temperature(stellar_flux: f64, albedo: f64, emissivity: f64) -> f64 {
175 (stellar_flux * (1.0 - albedo) / (emissivity * SIGMA_SB)).powf(0.25)
176}
177
178pub fn nightside_temperature(
179 thermal_inertia_val: f64,
180 subsolar_temp: f64,
181 rotation_period: f64,
182) -> f64 {
183 let theta = thermal_inertia_val
184 / (SIGMA_SB * subsolar_temp.powi(3) * (rotation_period / std::f64::consts::PI).sqrt());
185 subsolar_temp * theta.powf(0.5).min(1.0)
186}
187
188pub fn sputtering_loss_rate(
189 stellar_wind_flux: f64,
190 sputtering_yield: f64,
191 cross_section: f64,
192) -> f64 {
193 stellar_wind_flux * sputtering_yield * cross_section
194}
195
196pub fn hill_sphere_atmospheric(planet_mass: f64, stellar_mass: f64, semi_major_axis: f64) -> f64 {
197 semi_major_axis * (planet_mass / (3.0 * stellar_mass)).powf(1.0 / 3.0)
198}
199
200pub struct PlanetData {
201 pub mass: f64,
202 pub radius: f64,
203 pub flattening: f64,
204 pub orbital_period: f64,
205 pub semi_major_axis: f64,
206 pub eccentricity: f64,
207 pub inclination: f64,
208 pub axial_tilt: f64,
209 pub sidereal_day: f64,
210 pub surface_gravity: f64,
211 pub escape_velocity: f64,
212 pub mean_density: f64,
213 pub bond_albedo: f64,
214 pub orbital_velocity: f64,
215}
216
217pub fn planet_data(name: &str) -> Option<PlanetData> {
218 match name.to_ascii_lowercase().as_str() {
219 "mercury" => Some(PlanetData {
220 mass: MERCURY_MASS,
221 radius: MERCURY_RADIUS,
222 flattening: MERCURY_FLATTENING,
223 orbital_period: MERCURY_ORBITAL_PERIOD,
224 semi_major_axis: MERCURY_SEMI_MAJOR_AXIS,
225 eccentricity: MERCURY_ECCENTRICITY,
226 inclination: MERCURY_INCLINATION,
227 axial_tilt: MERCURY_AXIAL_TILT,
228 sidereal_day: MERCURY_SIDEREAL_DAY,
229 surface_gravity: MERCURY_SURFACE_GRAVITY,
230 escape_velocity: MERCURY_ESCAPE_VELOCITY,
231 mean_density: MERCURY_MEAN_DENSITY,
232 bond_albedo: MERCURY_BOND_ALBEDO,
233 orbital_velocity: MERCURY_ORBITAL_VELOCITY,
234 }),
235 "venus" => Some(PlanetData {
236 mass: VENUS_MASS,
237 radius: VENUS_RADIUS,
238 flattening: VENUS_FLATTENING,
239 orbital_period: VENUS_ORBITAL_PERIOD,
240 semi_major_axis: VENUS_SEMI_MAJOR_AXIS,
241 eccentricity: VENUS_ECCENTRICITY,
242 inclination: VENUS_INCLINATION,
243 axial_tilt: VENUS_AXIAL_TILT,
244 sidereal_day: VENUS_SIDEREAL_DAY,
245 surface_gravity: VENUS_SURFACE_GRAVITY,
246 escape_velocity: VENUS_ESCAPE_VELOCITY,
247 mean_density: VENUS_MEAN_DENSITY,
248 bond_albedo: VENUS_BOND_ALBEDO,
249 orbital_velocity: VENUS_ORBITAL_VELOCITY,
250 }),
251 "earth" => Some(PlanetData {
252 mass: EARTH_MASS,
253 radius: EARTH_RADIUS,
254 flattening: 0.003_353_36,
255 orbital_period: 365.256_363_004,
256 semi_major_axis: crate::constants::AU,
257 eccentricity: 0.016_708_6,
258 inclination: 0.0,
259 axial_tilt: 23.439_3,
260 sidereal_day: 86164.1,
261 surface_gravity: crate::constants::EARTH_GRAVITY,
262 escape_velocity: 11186.0,
263 mean_density: 5514.0,
264 bond_albedo: 0.306,
265 orbital_velocity: EARTH_ORBITAL_VELOCITY,
266 }),
267 "mars" => Some(PlanetData {
268 mass: MARS_MASS,
269 radius: MARS_RADIUS,
270 flattening: MARS_FLATTENING,
271 orbital_period: MARS_ORBITAL_PERIOD,
272 semi_major_axis: MARS_SEMI_MAJOR_AXIS,
273 eccentricity: MARS_ECCENTRICITY,
274 inclination: MARS_INCLINATION,
275 axial_tilt: MARS_AXIAL_TILT,
276 sidereal_day: MARS_SIDEREAL_DAY,
277 surface_gravity: MARS_SURFACE_GRAVITY,
278 escape_velocity: MARS_ESCAPE_VELOCITY,
279 mean_density: MARS_MEAN_DENSITY,
280 bond_albedo: MARS_BOND_ALBEDO,
281 orbital_velocity: MARS_ORBITAL_VELOCITY,
282 }),
283 "jupiter" => Some(PlanetData {
284 mass: JUPITER_MASS,
285 radius: JUPITER_RADIUS,
286 flattening: JUPITER_FLATTENING,
287 orbital_period: JUPITER_ORBITAL_PERIOD,
288 semi_major_axis: JUPITER_SEMI_MAJOR_AXIS,
289 eccentricity: JUPITER_ECCENTRICITY,
290 inclination: JUPITER_INCLINATION,
291 axial_tilt: JUPITER_AXIAL_TILT,
292 sidereal_day: JUPITER_SIDEREAL_DAY,
293 surface_gravity: JUPITER_SURFACE_GRAVITY,
294 escape_velocity: JUPITER_ESCAPE_VELOCITY,
295 mean_density: JUPITER_MEAN_DENSITY,
296 bond_albedo: JUPITER_BOND_ALBEDO,
297 orbital_velocity: JUPITER_ORBITAL_VELOCITY,
298 }),
299 "saturn" => Some(PlanetData {
300 mass: SATURN_MASS,
301 radius: SATURN_RADIUS,
302 flattening: SATURN_FLATTENING,
303 orbital_period: SATURN_ORBITAL_PERIOD,
304 semi_major_axis: SATURN_SEMI_MAJOR_AXIS,
305 eccentricity: SATURN_ECCENTRICITY,
306 inclination: SATURN_INCLINATION,
307 axial_tilt: SATURN_AXIAL_TILT,
308 sidereal_day: SATURN_SIDEREAL_DAY,
309 surface_gravity: SATURN_SURFACE_GRAVITY,
310 escape_velocity: SATURN_ESCAPE_VELOCITY,
311 mean_density: SATURN_MEAN_DENSITY,
312 bond_albedo: SATURN_BOND_ALBEDO,
313 orbital_velocity: SATURN_ORBITAL_VELOCITY,
314 }),
315 "uranus" => Some(PlanetData {
316 mass: URANUS_MASS,
317 radius: URANUS_RADIUS,
318 flattening: URANUS_FLATTENING,
319 orbital_period: URANUS_ORBITAL_PERIOD,
320 semi_major_axis: URANUS_SEMI_MAJOR_AXIS,
321 eccentricity: URANUS_ECCENTRICITY,
322 inclination: URANUS_INCLINATION,
323 axial_tilt: URANUS_AXIAL_TILT,
324 sidereal_day: URANUS_SIDEREAL_DAY,
325 surface_gravity: URANUS_SURFACE_GRAVITY,
326 escape_velocity: URANUS_ESCAPE_VELOCITY,
327 mean_density: URANUS_MEAN_DENSITY,
328 bond_albedo: URANUS_BOND_ALBEDO,
329 orbital_velocity: URANUS_ORBITAL_VELOCITY,
330 }),
331 "neptune" => Some(PlanetData {
332 mass: NEPTUNE_MASS,
333 radius: NEPTUNE_RADIUS,
334 flattening: NEPTUNE_FLATTENING,
335 orbital_period: NEPTUNE_ORBITAL_PERIOD,
336 semi_major_axis: NEPTUNE_SEMI_MAJOR_AXIS,
337 eccentricity: NEPTUNE_ECCENTRICITY,
338 inclination: NEPTUNE_INCLINATION,
339 axial_tilt: NEPTUNE_AXIAL_TILT,
340 sidereal_day: NEPTUNE_SIDEREAL_DAY,
341 surface_gravity: NEPTUNE_SURFACE_GRAVITY,
342 escape_velocity: NEPTUNE_ESCAPE_VELOCITY,
343 mean_density: NEPTUNE_MEAN_DENSITY,
344 bond_albedo: NEPTUNE_BOND_ALBEDO,
345 orbital_velocity: NEPTUNE_ORBITAL_VELOCITY,
346 }),
347 _ => None,
348 }
349}
350
351pub fn planet_mass(name: &str) -> Option<f64> {
352 planet_data(name).map(|p| p.mass)
353}
354
355pub fn planet_radius(name: &str) -> Option<f64> {
356 planet_data(name).map(|p| p.radius)
357}
358
359pub fn planet_flattening(name: &str) -> Option<f64> {
360 planet_data(name).map(|p| p.flattening)
361}
362
363pub fn planet_orbital_period(name: &str) -> Option<f64> {
364 planet_data(name).map(|p| p.orbital_period)
365}
366
367pub fn planet_semi_major_axis(name: &str) -> Option<f64> {
368 planet_data(name).map(|p| p.semi_major_axis)
369}
370
371pub fn planet_eccentricity(name: &str) -> Option<f64> {
372 planet_data(name).map(|p| p.eccentricity)
373}
374
375pub fn planet_inclination(name: &str) -> Option<f64> {
376 planet_data(name).map(|p| p.inclination)
377}
378
379pub fn planet_axial_tilt(name: &str) -> Option<f64> {
380 planet_data(name).map(|p| p.axial_tilt)
381}
382
383pub fn planet_sidereal_day(name: &str) -> Option<f64> {
384 planet_data(name).map(|p| p.sidereal_day)
385}
386
387pub fn planet_surface_gravity(name: &str) -> Option<f64> {
388 planet_data(name).map(|p| p.surface_gravity)
389}
390
391pub fn planet_escape_velocity(name: &str) -> Option<f64> {
392 planet_data(name).map(|p| p.escape_velocity)
393}
394
395pub fn planet_mean_density(name: &str) -> Option<f64> {
396 planet_data(name).map(|p| p.mean_density)
397}
398
399pub fn planet_bond_albedo(name: &str) -> Option<f64> {
400 planet_data(name).map(|p| p.bond_albedo)
401}
402
403pub fn planet_orbital_velocity(name: &str) -> Option<f64> {
404 planet_data(name).map(|p| p.orbital_velocity)
405}
406
407pub fn planet_volume(name: &str) -> Option<f64> {
408 planet_data(name).map(|p| 4.0 / 3.0 * std::f64::consts::PI * p.radius.powi(3))
409}
410
411pub fn planet_circumference(name: &str) -> Option<f64> {
412 planet_data(name).map(|p| 2.0 * std::f64::consts::PI * p.radius)
413}
414
415pub fn planet_surface_area(name: &str) -> Option<f64> {
416 planet_data(name).map(|p| 4.0 * std::f64::consts::PI * p.radius.powi(2))
417}
418
419pub fn planet_gravitational_parameter(name: &str) -> Option<f64> {
420 planet_data(name).map(|p| G * p.mass)
421}
422
423pub fn planet_hill_sphere(name: &str, stellar_mass: f64) -> Option<f64> {
424 planet_data(name).map(|p| {
425 p.semi_major_axis * (1.0 - p.eccentricity) * (p.mass / (3.0 * stellar_mass)).powf(1.0 / 3.0)
426 })
427}
428
429pub fn planet_roche_limit(name: &str, secondary_density: f64) -> Option<f64> {
430 planet_data(name)
431 .map(|p| 2.456 * p.radius * (p.mean_density / secondary_density).powf(1.0 / 3.0))
432}
433
434pub fn planet_synchronous_orbit(name: &str) -> Option<f64> {
435 planet_data(name).map(|p| {
436 let omega = 2.0 * std::f64::consts::PI / p.sidereal_day;
437 (G * p.mass / omega.powi(2)).powf(1.0 / 3.0)
438 })
439}