use crate::atmosphere::ScatteringFactor;
use crate::ext_qtty::Quantity;
use crate::qtty::{Kilometers, Radians};
pub fn van_rhijn_factor(
zenith: Radians,
emission_height: Kilometers,
) -> Quantity<ScatteringFactor> {
van_rhijn_factor_with_radius(
zenith,
emission_height,
crate::bodies::solar_system::EARTH.radius,
)
}
pub fn van_rhijn_factor_with_radius(
zenith: Radians,
emission_height: Kilometers,
body_radius: Kilometers,
) -> Quantity<ScatteringFactor> {
let r = body_radius.value();
let h = emission_height.value();
if !zenith.is_finite() || !h.is_finite() || !r.is_finite() || h <= 0.0 || r <= 0.0 {
return Quantity::<ScatteringFactor>::new(f64::NAN);
}
let ratio = r / (r + h);
let s = zenith.sin();
let inner = 1.0 - (ratio * s) * (ratio * s);
let v = if inner <= 0.0 {
f64::INFINITY
} else {
inner.sqrt().recip()
};
Quantity::<ScatteringFactor>::new(v)
}
#[cfg(test)]
mod tests {
use super::*;
use crate::qtty::{unit::Radian, Degrees};
#[test]
fn van_rhijn_is_unity_at_zenith() {
let f = van_rhijn_factor(Degrees::new(0.0).to::<Radian>(), Kilometers::new(90.0));
assert!((f.value() - 1.0).abs() < 1.0e-12);
}
#[test]
fn van_rhijn_increases_toward_horizon_but_stays_finite() {
let h = Kilometers::new(90.0);
let z45 = van_rhijn_factor(Degrees::new(45.0).to::<Radian>(), h).value();
let z80 = van_rhijn_factor(Degrees::new(80.0).to::<Radian>(), h).value();
let z90 = van_rhijn_factor(Degrees::new(90.0).to::<Radian>(), h).value();
assert!(z45 > 1.0);
assert!(z80 > z45);
assert!(z90 > z80);
assert!(z90.is_finite());
}
}