use sciforge::hub::prelude::constants::N_A;
use sciforge::hub::prelude::constants::elements::atomic_mass;
pub struct MolecularSpecies {
pub name: &'static str,
pub symbol: &'static str,
pub molar_mass_kg_mol: f64,
pub volume_fraction: f64,
pub refractive_index_stp: f64,
pub depolarization_factor: f64,
}
pub struct AtmosphereEndpoint {
pub planet_radius_m: f64,
pub atmosphere_height_m: f64,
pub sea_level_pressure_pa: f64,
pub sea_level_temperature_k: f64,
pub sea_level_number_density_m3: f64,
pub mean_molar_mass_kg_mol: f64,
pub rayleigh_scale_height_m: f64,
pub mie_scale_height_m: f64,
pub mie_asymmetry_g: f64,
pub mie_coefficient: f64,
pub ammonia_absorption: [f64; 3],
pub rayleigh_coefficients_rgb: [f64; 3],
pub species: Vec<MolecularSpecies>,
pub sun_irradiance_w_m2: f64,
}
fn h2_molar() -> f64 {
2.0 * atomic_mass(1) * 1e-3
}
fn he_molar() -> f64 {
atomic_mass(2) * 1e-3
}
fn ch4_molar() -> f64 {
(atomic_mass(6) + 4.0 * atomic_mass(1)) * 1e-3
}
fn nh3_molar() -> f64 {
(atomic_mass(7) + 3.0 * atomic_mass(1)) * 1e-3
}
fn h2o_molar() -> f64 {
(2.0 * atomic_mass(1) + atomic_mass(8)) * 1e-3
}
fn h2s_molar() -> f64 {
(2.0 * atomic_mass(1) + atomic_mass(16)) * 1e-3
}
fn ph3_molar() -> f64 {
(atomic_mass(15) + 3.0 * atomic_mass(1)) * 1e-3
}
fn mean_molar_mass() -> f64 {
0.898 * h2_molar()
+ 0.102 * he_molar()
+ 3.0e-3 * ch4_molar()
+ 2.6e-4 * nh3_molar()
+ 6.0e-6 * h2o_molar()
+ 7.7e-5 * h2s_molar()
+ 6.0e-7 * ph3_molar()
}
fn sea_level_number_density() -> f64 {
let m_air = mean_molar_mass();
let rho = *crate::ONEBARDENSITY;
N_A * rho / m_air
}
fn rayleigh_beta_rgb(n_density: f64) -> [f64; 3] {
let n_minus_1 = 1.32e-4;
let ns2 = (2.0 * n_minus_1) * (2.0 * n_minus_1);
let coeff = 8.0 * std::f64::consts::PI.powi(3) * ns2 / (3.0 * n_density);
[
coeff / (680e-9_f64).powi(4),
coeff / (550e-9_f64).powi(4),
coeff / (440e-9_f64).powi(4),
]
}
impl AtmosphereEndpoint {
pub fn jupiter() -> Self {
let n_density = sea_level_number_density();
let beta = rayleigh_beta_rgb(n_density);
let species = vec![
MolecularSpecies {
name: "Dihydrogen",
symbol: "H2",
molar_mass_kg_mol: h2_molar(),
volume_fraction: 0.898,
refractive_index_stp: 1.000_132_0,
depolarization_factor: 0.020,
},
MolecularSpecies {
name: "Helium",
symbol: "He",
molar_mass_kg_mol: he_molar(),
volume_fraction: 0.102,
refractive_index_stp: 1.000_035_0,
depolarization_factor: 0.0,
},
MolecularSpecies {
name: "Methane",
symbol: "CH4",
molar_mass_kg_mol: ch4_molar(),
volume_fraction: 3.0e-3,
refractive_index_stp: 1.000_444_0,
depolarization_factor: 0.0,
},
MolecularSpecies {
name: "Ammonia",
symbol: "NH3",
molar_mass_kg_mol: nh3_molar(),
volume_fraction: 2.6e-4,
refractive_index_stp: 1.000_376_0,
depolarization_factor: 0.0,
},
MolecularSpecies {
name: "Water vapor",
symbol: "H2O",
molar_mass_kg_mol: h2o_molar(),
volume_fraction: 6.0e-6,
refractive_index_stp: 1.000_256_0,
depolarization_factor: 0.17,
},
MolecularSpecies {
name: "Hydrogen sulfide",
symbol: "H2S",
molar_mass_kg_mol: h2s_molar(),
volume_fraction: 7.7e-5,
refractive_index_stp: 1.000_644_0,
depolarization_factor: 0.0,
},
MolecularSpecies {
name: "Phosphine",
symbol: "PH3",
molar_mass_kg_mol: ph3_molar(),
volume_fraction: 6.0e-7,
refractive_index_stp: 1.000_500_0,
depolarization_factor: 0.0,
},
];
Self {
planet_radius_m: crate::JUPITEREQUATORIALRADIUS,
atmosphere_height_m: 1_000_000.0,
sea_level_pressure_pa: 100_000.0,
sea_level_temperature_k: 165.0,
sea_level_number_density_m3: n_density,
mean_molar_mass_kg_mol: mean_molar_mass(),
rayleigh_scale_height_m: *crate::SCALEHEIGHT,
mie_scale_height_m: 20_000.0,
mie_asymmetry_g: 0.78,
mie_coefficient: 8e-6,
ammonia_absorption: [0.001, 0.0005, 0.0002],
rayleigh_coefficients_rgb: beta,
species,
sun_irradiance_w_m2: 1_361.0 / (5.2044 * 5.2044),
}
}
pub fn sky_luminance(&self, observer_alt: f64, cos_zenith: f64) -> f64 {
let cos_z = cos_zenith.max(0.01);
let molecular_rayleigh: f64 = self
.species
.iter()
.map(|s| {
let n_s = self.sea_level_number_density_m3 * s.volume_fraction;
let delta = s.depolarization_factor;
let king = (6.0 + 3.0 * delta) / (6.0 - 7.0 * delta);
let ns = s.refractive_index_stp - 1.0;
8.0 * std::f64::consts::PI.powi(3) * (2.0 * ns).powi(2) * king / (3.0 * n_s)
})
.sum();
let altitude_factor = (-observer_alt / self.rayleigh_scale_height_m).exp();
let mie_factor = (-observer_alt / self.mie_scale_height_m).exp() * self.mie_coefficient;
let g = self.mie_asymmetry_g;
let hg_phase = (1.0 - g * g)
/ (4.0 * std::f64::consts::PI * (1.0 + g * g - 2.0 * g * cos_z).powf(1.5));
let ammonia_abs = self.ammonia_absorption[0] * 0.3
+ self.ammonia_absorption[1] * 0.59
+ self.ammonia_absorption[2] * 0.11;
let scatter_550 = molecular_rayleigh / (550e-9_f64).powi(4) * altitude_factor;
let sky_r = self.rayleigh_coefficients_rgb[0] * cos_z * altitude_factor * 1e4;
let sky_g = self.rayleigh_coefficients_rgb[1] * cos_z * altitude_factor * 1e4;
let sky_b = self.rayleigh_coefficients_rgb[2] * cos_z * altitude_factor * 1e4;
sky_r * 0.3
+ sky_g * 0.59
+ sky_b * 0.11
+ scatter_550 * 1e-20
+ mie_factor * 1e-6
+ ammonia_abs * 1e-8
+ hg_phase * mie_factor * 1e-4
}
pub fn barometric_density(&self) -> f64 {
self.sea_level_pressure_pa
/ (self.sea_level_temperature_k * 8.314_462_618 / self.mean_molar_mass_kg_mol)
}
pub fn shell_volume(&self) -> f64 {
4.0 / 3.0
* std::f64::consts::PI
* ((self.planet_radius_m + self.atmosphere_height_m).powi(3)
- self.planet_radius_m.powi(3))
}
}