use sciforge::hub::domain::common::constants::EARTH_RADIUS;
use sciforge::hub::prelude::constants::N_A;
use sciforge::hub::prelude::constants::elements::atomic_mass;
fn rayleigh_beta_rgb() -> [f64; 3] {
let m_n2 = 2.0 * atomic_mass(7) * 1e-3;
let m_o2 = 2.0 * atomic_mass(8) * 1e-3;
let m_air = 0.78084 * m_n2 + 0.20946 * m_o2 + 0.00934 * atomic_mass(18) * 1e-3;
let n_density = N_A * *crate::SEA_LEVEL_AIR_DENSITY / m_air;
let n_minus_1 = 2.9e-4;
let n_sq_minus_1_sq = (2.0 * n_minus_1) * (2.0 * n_minus_1);
let coeff = 8.0 * std::f64::consts::PI.powi(3) * n_sq_minus_1_sq / (3.0 * n_density);
let lambdas = [680e-9_f64, 550e-9, 440e-9];
[
coeff / lambdas[0].powi(4),
coeff / lambdas[1].powi(4),
coeff / lambdas[2].powi(4),
]
}
pub struct AtmosphereParams {
pub planet_radius_m: f64,
pub atmosphere_height_m: f64,
pub rayleigh_coefficients: [f64; 3],
pub rayleigh_scale_height_m: f64,
pub mie_coefficient: f64,
pub mie_scale_height_m: f64,
pub mie_asymmetry: f64,
pub sun_intensity: f64,
pub num_samples: u32,
pub num_light_samples: u32,
}
impl Default for AtmosphereParams {
fn default() -> Self {
let beta = rayleigh_beta_rgb();
Self {
planet_radius_m: EARTH_RADIUS,
atmosphere_height_m: 100_000.0,
rayleigh_coefficients: beta,
rayleigh_scale_height_m: 8_500.0,
mie_coefficient: 21e-6,
mie_scale_height_m: 1_200.0,
mie_asymmetry: 0.76,
sun_intensity: 22.0,
num_samples: 16,
num_light_samples: 8,
}
}
}
impl AtmosphereParams {
pub fn rayleigh_density(&self, altitude_m: f64) -> f64 {
(-altitude_m / self.rayleigh_scale_height_m).exp()
}
pub fn mie_density(&self, altitude_m: f64) -> f64 {
(-altitude_m / self.mie_scale_height_m).exp()
}
pub fn scatter_rayleigh(&self, cos_theta: f64) -> f64 {
3.0 / (16.0 * std::f64::consts::PI) * (1.0 + cos_theta * cos_theta)
}
pub fn scatter_mie(&self, cos_theta: f64) -> f64 {
let g = self.mie_asymmetry;
let g2 = g * g;
let num = 3.0 * (1.0 - g2) * (1.0 + cos_theta * cos_theta);
let den =
8.0 * std::f64::consts::PI * (2.0 + g2) * (1.0 + g2 - 2.0 * g * cos_theta).powf(1.5);
num / den
}
pub fn sky_color(
&self,
sun_dir: [f64; 3],
view_dir: [f64; 3],
camera_altitude_m: f64,
) -> [f64; 3] {
let cos_theta =
sun_dir[0] * view_dir[0] + sun_dir[1] * view_dir[1] + sun_dir[2] * view_dir[2];
let rayleigh_phase = self.scatter_rayleigh(cos_theta);
let mie_phase = self.scatter_mie(cos_theta);
const LAMBDA_R: f64 = 680.0;
const LAMBDA_G: f64 = 550.0;
const LAMBDA_B: f64 = 440.0;
let inv4_r = (LAMBDA_G / LAMBDA_R).powi(4);
let inv4_g = 1.0;
let inv4_b = (LAMBDA_G / LAMBDA_B).powi(4);
let n_steps = self.num_samples as usize;
let atmo_top = self.planet_radius_m + self.atmosphere_height_m;
let ray_len = atmo_top - (self.planet_radius_m + camera_altitude_m);
let step_len = ray_len / n_steps as f64;
let mut color = [0.0_f64; 3];
let mut optical_r = 0.0;
let mut optical_g = 0.0;
let mut optical_b = 0.0;
for i in 0..n_steps {
let h = camera_altitude_m + step_len * (i as f64 + 0.5);
let rho_ray = self.rayleigh_density(h);
let rho_mie = self.mie_density(h);
optical_r += (self.rayleigh_coefficients[0] * inv4_r * rho_ray
+ self.mie_coefficient * rho_mie)
* step_len;
optical_g += (self.rayleigh_coefficients[1] * inv4_g * rho_ray
+ self.mie_coefficient * rho_mie)
* step_len;
optical_b += (self.rayleigh_coefficients[2] * inv4_b * rho_ray
+ self.mie_coefficient * rho_mie)
* step_len;
let n_light = self.num_light_samples as usize;
let light_step = (atmo_top - (self.planet_radius_m + h)) / n_light as f64;
let mut sun_r = 0.0;
let mut sun_g = 0.0;
let mut sun_b = 0.0;
for j in 0..n_light {
let lh = h + light_step * (j as f64 + 0.5);
let lr = self.rayleigh_density(lh);
let lm = self.mie_density(lh);
sun_r += (self.rayleigh_coefficients[0] * inv4_r * lr + self.mie_coefficient * lm)
* light_step;
sun_g += (self.rayleigh_coefficients[1] * inv4_g * lr + self.mie_coefficient * lm)
* light_step;
sun_b += (self.rayleigh_coefficients[2] * inv4_b * lr + self.mie_coefficient * lm)
* light_step;
}
let atten_r = (-sun_r - optical_r).exp();
let atten_g = (-sun_g - optical_g).exp();
let atten_b = (-sun_b - optical_b).exp();
color[0] += (self.rayleigh_coefficients[0] * inv4_r * rho_ray * rayleigh_phase
+ self.mie_coefficient * rho_mie * mie_phase)
* atten_r
* step_len;
color[1] += (self.rayleigh_coefficients[1] * inv4_g * rho_ray * rayleigh_phase
+ self.mie_coefficient * rho_mie * mie_phase)
* atten_g
* step_len;
color[2] += (self.rayleigh_coefficients[2] * inv4_b * rho_ray * rayleigh_phase
+ self.mie_coefficient * rho_mie * mie_phase)
* atten_b
* step_len;
}
[
self.sun_intensity * color[0],
self.sun_intensity * color[1],
self.sun_intensity * color[2],
]
}
}