use rand::RngCore;
use crate::consts::*;
use crate::structs::planetesimal::*;
use crate::utils::*;
pub fn orbital_zone(luminosity: &f64, distance_to_primary_star: f64) -> i32 {
if distance_to_primary_star < (4.0 * luminosity.sqrt()) {
return 1;
} else if distance_to_primary_star < (15.0 * luminosity.sqrt()) {
return 2;
}
3
}
pub fn volume_radius(mass: &f64, density: &f64) -> f64 {
let volume = mass * SOLAR_MASS_IN_GRAMS / density;
((3.0 * volume) / (4.0 * PI)).powf(0.33) / CM_PER_KM
}
pub fn kothari_radius(mass: &f64, giant: &bool, zone: &i32) -> f64 {
let mut atomic_weight = 0.0;
let mut atomic_num = 0.0;
match (zone, giant) {
(1, true) => {
atomic_weight = 9.5;
atomic_num = 4.5;
}
(1, false) => {
atomic_weight = 15.0;
atomic_num = 8.0;
}
(2, true) => {
atomic_weight = 2.47;
atomic_num = 2.0;
}
(2, false) => {
atomic_weight = 10.0;
atomic_num = 5.0;
}
(3, true) => {
atomic_weight = 7.0;
atomic_num = 4.0;
}
(3, false) => {
atomic_weight = 10.0;
atomic_num = 5.0;
}
(_, _) => (),
}
let mut temp: f64 = atomic_weight * atomic_num;
temp = 2.0 * BETA_20 * SOLAR_MASS_IN_GRAMS.powf(1.0 / 3.0) / (A1_20 * temp.powf(1.0 / 3.0));
let mut temp2 = A2_20 * atomic_weight.powf(4.0 / 3.0) * SOLAR_MASS_IN_GRAMS.powf(2.0 / 3.0);
temp2 *= mass.powf(2.0 / 3.0);
temp2 /= A1_20 * atomic_num.powf(2.0);
temp2 += 1.0;
temp /= temp2;
temp = (temp * mass.powf(1.0 / 3.0)) / CM_PER_KM;
temp /= JIMS_FUDGE;
temp
}
pub fn empirical_density(
mass: &f64,
distance_to_primary_star: &f64,
ecosphere_radius: &f64,
is_gas_giant: &bool,
) -> f64 {
let mut density = (mass * EARTH_MASSES_PER_SOLAR_MASS).powf(1.0 / 8.0);
density *= (ecosphere_radius / distance_to_primary_star).powf(0.25);
match is_gas_giant {
true => density * 1.2,
false => density * 5.5,
}
}
pub fn volume_density(mass: &f64, equat_radius: &f64) -> f64 {
let equat_radius = equat_radius * CM_PER_KM;
let volume = (4.0 * PI * equat_radius.powf(3.0)) / 3.0;
mass * SOLAR_MASS_IN_GRAMS / volume
}
pub fn period(separation: &f64, small_mass: &f64, large_mass: &f64) -> f64 {
let period_in_years = (separation.powf(3.0) / (small_mass + large_mass)).sqrt();
period_in_years * DAYS_IN_A_YEAR
}
pub fn day_length(planet: &mut Planetesimal, stellar_mass: &f64, main_sequence_age: &f64) -> f64 {
let planet_mass_in_grams = planet.mass * SOLAR_MASS_IN_GRAMS;
let equatorial_radius_in_cm = planet.radius * CM_PER_KM;
let year_in_hours = planet.orbital_period_days;
let k2 = match planet.is_gas_giant {
true => 0.24,
false => 0.33,
};
let base_angular_velocity =
(2.0 * J * planet_mass_in_grams).sqrt() / (k2 * equatorial_radius_in_cm.powf(2.0));
let change_in_angular_velocity = CHANGE_IN_EARTH_ANG_VEL
* (planet.density / EARTH_DENSITY)
* (equatorial_radius_in_cm / EARTH_RADIUS_IN_CM)
* (EARTH_MASS_IN_GRAMS / planet_mass_in_grams)
* stellar_mass.powf(2.0)
* (1.0 / planet.a.powf(6.0));
let ang_velocity = base_angular_velocity + (change_in_angular_velocity * main_sequence_age);
let (stopped, day_in_hours) = match ang_velocity <= 0.0 {
true => (true, INCREDIBLY_LARGE_NUMBER),
false => (
false,
RADIANS_PER_ROTATION / (SECONDS_PER_HOUR * ang_velocity),
),
};
if day_in_hours >= year_in_hours || stopped {
if planet.e > 0.1 {
let spin_resonance_factor = (1.0 - planet.e) / (1.0 + planet.e);
planet.resonant_period = true;
return spin_resonance_factor * year_in_hours;
} else {
return year_in_hours;
}
}
day_in_hours
}
pub fn inclination(orbital_radius: &f64, rng: &mut dyn RngCore) -> f64 {
let inclination = orbital_radius.powf(0.2) * about(EARTH_AXIAL_TILT, 0.4, rng);
inclination % 360.0
}
pub fn escape_vel(mass: &f64, radius: &f64) -> f64 {
let mass_in_grams = mass * SOLAR_MASS_IN_GRAMS;
let radius_in_cm = radius * CM_PER_KM;
(2.0 * GRAV_CONSTANT * mass_in_grams / radius_in_cm).sqrt()
}
pub fn rms_vel(molecular_weight: &f64, orbital_radius: &f64) -> f64 {
let exospheric_temp = EARTH_EXOSPHERE_TEMP / orbital_radius.powf(2.0);
((3.0 * MOLAR_GAS_CONST * exospheric_temp) / molecular_weight).sqrt() * CM_PER_METER
}
pub fn molecule_limit(mass: &f64, equatorial_radius: &f64) -> f64 {
let escape_velocity = escape_vel(mass, equatorial_radius);
3.0 * (GAS_RETENTION_THRESHOLD * CM_PER_METER).powf(2.0)
* MOLAR_GAS_CONST
* EARTH_EXOSPHERE_TEMP
/ escape_velocity.powf(2.0)
}
pub fn acceleration(mass: &f64, radius: &f64) -> f64 {
GRAV_CONSTANT * mass * SOLAR_MASS_IN_GRAMS / (radius * CM_PER_KM).powf(2.0)
}
pub fn gravity(acceleration: &f64) -> f64 {
acceleration / EARTH_ACCELERATION
}
pub fn greenhouse(
distance_to_primary_star: &f64,
orbit_zone: &i32,
surface_pressure_bar: &f64,
ecosphere_radius: &f64,
) -> bool {
let greenhouse_radius = ecosphere_radius * GREENHOUSE_EFFECT_CONST;
*distance_to_primary_star < greenhouse_radius && *orbit_zone == 1 && *surface_pressure_bar > 0.0
}
pub fn vol_inventory(
mass: &f64,
escape_vel: &f64,
rms_vel: &f64,
stellar_mass: &f64,
zone: &i32,
greenhouse_effect: &bool,
rng: &mut dyn RngCore,
) -> f64 {
let velocity_ratio = escape_vel / rms_vel;
if velocity_ratio < GAS_RETENTION_THRESHOLD {
return 0.0;
}
let proportion_const = match zone {
1 => 100000.0,
2 => 75000.0,
3 => 250.0,
_ => 10.0,
};
let mass_in_earth_units = mass * EARTH_MASSES_PER_SOLAR_MASS;
let temp1 = proportion_const * mass_in_earth_units / stellar_mass;
let temp2 = about(temp1, 0.2, rng);
if *greenhouse_effect {
return temp2;
}
temp2 / 100.0
}
pub fn pressure(volatile_gas_inventory: &f64, equatorial_radius: &f64, gravity: &f64) -> f64 {
let equatorial_radius = EARTH_RADIUS_IN_KM / equatorial_radius;
(volatile_gas_inventory * gravity / equatorial_radius.powf(2.0)) / EARTH_SURF_PRES_IN_MILLIBARS
}
pub fn boiling_point_kelvin(surface_pressure_bar: &f64) -> f64 {
1.0 / (surface_pressure_bar.log(std::f64::consts::E) / -5050.5 + 1.0 / 373.0)
}
pub fn hydrosphere_fraction(volatile_gas_inventory: &f64, planetary_radius: &f64) -> f64 {
let hydrosphere_fraction =
0.75 * volatile_gas_inventory / 1000.0 * (EARTH_RADIUS_IN_KM / planetary_radius).powf(2.0);
match hydrosphere_fraction >= 1.0 {
true => 1.0,
false => hydrosphere_fraction,
}
}
pub fn cloud_fraction(
surface_temp_kelvin: f64,
smallest_mw_retained: f64,
equatorial_radius: f64,
hydrosphere_fraction: f64,
) -> f64 {
if smallest_mw_retained > WATER_VAPOR {
return 0.0;
}
let surface_area = 4.0 * PI * equatorial_radius.powf(2.0);
let hydrosphere_mass = hydrosphere_fraction * surface_area * EARTH_WATER_MASS_PER_AREA;
let water_vapor_in_kg =
(0.00000001 * hydrosphere_mass) * (Q2_36 * (surface_temp_kelvin - 288.0)).exp();
let mut fraction = CLOUD_COVERAGE_FACTOR * water_vapor_in_kg / surface_area;
if fraction >= 1.0 {
fraction = 1.0;
}
fraction
}
pub fn ice_fraction(hydrosphere_fraction: &f64, surface_temp_kelvin: &f64) -> f64 {
let surface_temp_kelvin = match *surface_temp_kelvin > 328.0 {
true => 328.0,
false => *surface_temp_kelvin,
};
let mut temp = ((328.0 - surface_temp_kelvin) / 70.0).powf(5.0);
if temp > 1.5 * hydrosphere_fraction {
temp = 1.5 * hydrosphere_fraction;
}
if temp >= 1.0 {
return 1.0;
}
temp
}
pub fn eff_temp(ecosphere_radius: &f64, orbital_radius: &f64, albedo: &f64) -> f64 {
(ecosphere_radius / orbital_radius).sqrt()
* ((1.0 - albedo) / 0.7).powf(0.25)
* EARTH_EFFECTIVE_TEMP
}
pub fn green_rise(optical_depth: f64, effective_temp: f64, surface_pressure_bar: f64) -> f64 {
let convection_factor = EARTH_CONVECTION_FACTOR * surface_pressure_bar.powf(0.25);
((1.0 + 0.75 * optical_depth).powf(0.25) - 1.0) * effective_temp * convection_factor
}
pub fn planet_albedo(
water_fraction: &f64,
cloud_fraction: &f64,
ice_fraction: &f64,
surface_pressure_bar: &f64,
rng: &mut dyn RngCore,
) -> f64 {
let mut rock_fraction = 1.0 - *water_fraction - *ice_fraction;
let mut components: f64 = 0.0;
if *water_fraction > 0.0 {
components += 1.0;
}
if *ice_fraction > 0.0 {
components += 1.0;
}
if rock_fraction > 0.0 {
components += 1.0;
}
let cloud_adjustment = *cloud_fraction / components;
if rock_fraction >= cloud_adjustment {
rock_fraction -= cloud_adjustment;
} else {
rock_fraction = 0.0;
}
let water_fraction = match *water_fraction > cloud_adjustment {
true => *water_fraction - cloud_adjustment,
false => 0.0,
};
let ice_fraction = match *ice_fraction > cloud_adjustment {
true => *ice_fraction - cloud_adjustment,
false => 0.0,
};
let cloud_contribution = *cloud_fraction * about(CLOUD_ALBEDO, 0.2, rng);
let water_contribution = water_fraction * about(WATER_ALBEDO, 0.2, rng);
let (rock_contribution, ice_contribution) = match *surface_pressure_bar == 0.0 {
true => (
rock_fraction * about(AIRLESS_ROCKY_ALBEDO, 0.3, rng),
ice_fraction * about(AIRLESS_ICE_ALBEDO, 0.4, rng),
),
false => (
rock_fraction * about(ROCKY_ALBEDO, 0.1, rng),
ice_fraction * about(ICE_ALBEDO, 0.1, rng),
),
};
cloud_contribution + rock_contribution + water_contribution + ice_contribution
}
pub fn opacity(molecular_weight: f64, surface_pressure_bar: f64) -> f64 {
let mut optical_depth = 0.0;
if (0.0..10.0).contains(&molecular_weight) {
optical_depth += 3.0;
}
if (10.0..20.0).contains(&molecular_weight) {
optical_depth += 2.34;
}
if (20.0..30.0).contains(&molecular_weight) {
optical_depth += 1.0;
}
if (30.0..45.0).contains(&molecular_weight) {
optical_depth += 0.15;
}
if (45.0..100.0).contains(&molecular_weight) {
optical_depth += 0.05;
}
if surface_pressure_bar >= 0.07 {
optical_depth *= 8.333;
} else if surface_pressure_bar >= 0.05 {
optical_depth *= 6.666;
} else if surface_pressure_bar >= 0.03 {
optical_depth *= 3.333;
} else if surface_pressure_bar >= 0.01 {
optical_depth *= 2.0;
} else if surface_pressure_bar >= 0.005 {
optical_depth *= 1.5;
}
optical_depth
}
pub fn get_earth_mass(mass: f64) -> f64 {
mass * EARTH_MASSES_PER_SOLAR_MASS
}
pub fn iterate_surface_temp(
planet: &mut Planetesimal,
ecosphere_radius: &f64,
rng: &mut dyn RngCore,
) -> f64 {
let mut albedo = 0.0;
let mut water = 0.0;
let mut clouds = 0.0;
let mut ice = 0.0;
let mut optical_depth = opacity(planet.molecule_weight, planet.surface_pressure_bar);
let mut effective_temp = eff_temp(ecosphere_radius, &planet.a, &EARTH_ALBEDO);
let mut greenhouse_rise =
green_rise(optical_depth, effective_temp, planet.surface_pressure_bar);
let mut surface_temp_kelvin = effective_temp + greenhouse_rise;
let mut previous_temp = surface_temp_kelvin - 5.0;
while (surface_temp_kelvin - previous_temp).abs() > 1.0 {
previous_temp = surface_temp_kelvin;
water = hydrosphere_fraction(&planet.volatile_gas_inventory, &planet.radius);
clouds = cloud_fraction(
surface_temp_kelvin,
planet.molecule_weight,
planet.radius,
water,
);
ice = ice_fraction(&water, &surface_temp_kelvin);
if surface_temp_kelvin >= planet.boiling_point_kelvin
|| surface_temp_kelvin <= FREEZING_POINT_OF_WATER
{
water = 0.0;
}
albedo = planet_albedo(&water, &clouds, &ice, &planet.surface_pressure_bar, rng);
optical_depth = opacity(planet.molecule_weight, planet.surface_pressure_bar);
effective_temp = eff_temp(ecosphere_radius, &planet.a, &albedo);
greenhouse_rise = green_rise(optical_depth, effective_temp, planet.surface_pressure_bar);
surface_temp_kelvin = effective_temp + greenhouse_rise;
}
planet.hydrosphere = water;
planet.cloud_cover = clouds;
planet.ice_cover = ice;
planet.albedo = albedo;
planet.surface_temp_kelvin = surface_temp_kelvin;
surface_temp_kelvin
}
pub fn check_tidal_lock(day_length: f64, orbital_period: f64) -> bool {
let error_margin = f64::EPSILON;
(day_length - orbital_period * day_length).abs() < error_margin
}
fn lim(x: f64) -> f64 {
x / (1.0 + x.powf(4.0)).sqrt().sqrt()
}
fn soft(v: &f64, max: &f64, min: &f64) -> f64 {
let dv = v - min;
let dm = max - min;
(lim(2.0 * dv / dm - 1.0) + 1.0) / 2.0 * dm + min
}
pub fn get_day_night_temp_kelvin(planet: &mut Planetesimal) {
let Planetesimal {
surface_pressure_bar,
surface_temp_kelvin,
axial_tilt,
e,
day_hours,
..
} = planet;
let pressmod: f64 = 1.0 / (1.0 + 20.0 * *surface_pressure_bar).sqrt();
let ppmod = 1.0 / (10.0 + 5.0 * *surface_pressure_bar).sqrt();
let tiltmod = ((*axial_tilt * PI / 180.0).cos() * (1.0 + *e).powf(2.0)).abs();
let daymod = 1.0 / (200.0 / *day_hours + 1.0);
let mh = (1.0 + daymod).powf(pressmod);
let ml = (1.0 - daymod).powf(pressmod);
let hi = mh * *surface_temp_kelvin;
let mut lo = ml * *surface_temp_kelvin;
let sh = hi + ((100.0 + hi) * tiltmod).powf((ppmod).sqrt());
let mut wl = lo - ((150.0 + lo) * tiltmod).powf((ppmod).sqrt());
let max = *surface_temp_kelvin + surface_temp_kelvin.sqrt() * 10.0;
let min = *surface_temp_kelvin / (*day_hours + 24.0).sqrt();
if lo < min {
lo = min;
}
if wl < 0.0 {
wl = 0.0;
}
let high_temp = soft(&hi, &max, &min);
let low_temp = soft(&lo, &max, &min);
let max_temp = soft(&sh, &max, &min);
let min_temp = soft(&wl, &max, &min);
if (high_temp - max_temp).abs() > 10.0 || (low_temp - min_temp).abs() > 10.0 {
planet.day_temp_kelvin = high_temp;
planet.night_temp_kelvin = low_temp;
}
planet.day_temp_kelvin = max_temp;
planet.night_temp_kelvin = min_temp;
}