use super::*;
use crate::units::Nanolamberts;
pub(crate) struct KrisciunasSchaefer1991 {
location: Geodetic<ECEF>,
k_ext: MagnitudesPerAirmass,
}
impl KrisciunasSchaefer1991 {
pub(crate) fn published_reference(location: Geodetic<ECEF>) -> Self {
Self {
location,
k_ext: DEFAULT_K_EXT,
}
}
pub(crate) fn compute(
&self,
time: Time<UTC>,
target: SphericalDirection<EquatorialMeanJ2000>,
) -> Result<MoonOutputs> {
self.compute_with_sky(time, target, &MoonlightSky::at(time, self.location))
}
pub(crate) fn compute_with_sky(
&self,
time: Time<UTC>,
target: SphericalDirection<EquatorialMeanJ2000>,
sky: &MoonlightSky,
) -> Result<MoonOutputs> {
let geometry = sky.geometry(time.to::<TT>().to::<JD>(), self.location, target);
compute_krisciunas_schaefer_1991(&geometry, self.k_ext)
}
}
pub(super) fn compute_krisciunas_schaefer_1991(
inp: &MoonlightGeometry,
k_ext: MagnitudesPerAirmass,
) -> Result<MoonOutputs> {
inp.validate()?;
if !k_ext.is_finite() {
return Err(crate::NsbError::OutOfRange(
"non-finite Moonlight extinction".into(),
));
}
if inp.moon_zenith >= Degrees::new(90.0)
|| inp.source_zenith >= Degrees::new(90.0)
|| inp.separation <= Degrees::new(0.0)
{
return Ok(zero_outputs());
}
let b_nl = scattered_brightness_nanolamberts(
inp.phase.phase_angle,
inp.separation,
inp.moon_zenith,
inp.source_zenith,
k_ext,
);
if !b_nl.is_finite() {
return Err(crate::NsbError::OutOfRange(
"non-finite K&S lunar brightness".into(),
));
}
if b_nl <= Nanolamberts::new(0.0) {
return Ok(zero_outputs());
}
let v_mag_arcsec2 = v_mag_per_arcsec2_from_nl(b_nl);
let v_s10 = 10f64.powf(0.4 * (NSB_S10_ZP.value() - v_mag_arcsec2.value()));
let integrated = S10_V_TO_INTEGRATED_PH * v_s10;
Ok(MoonOutputs {
integrated,
b_flux_s10: radiometry::S10s::new(v_s10),
v_flux_s10: radiometry::S10s::new(v_s10),
})
}
fn lunar_illuminance_outside_atmosphere(alpha: Radians) -> f64 {
let a = alpha.abs().to::<Degree>().value();
let exponent = -0.4 * (3.84 + 0.026 * a + 4.0e-9 * a.powi(4));
10f64.powf(exponent)
}
fn scattering_function(rho: Degrees) -> f64 {
let cos_rho = rho.cos();
let rayleigh = 10f64.powf(5.36) * (1.06 + cos_rho * cos_rho);
let aureole = 10f64.powf(6.15 - rho.value() / 40.0);
rayleigh + aureole
}
fn v_mag_per_arcsec2_from_nl(b_nl: Nanolamberts) -> crate::units::photometry::SurfaceBrightness {
crate::units::photometry::SurfaceBrightness::new(
(20.7233 - (b_nl.value() / 34.08).ln()) / 0.92104,
)
}
fn scattered_brightness_nanolamberts(
alpha: Radians,
rho: Degrees,
z_moon: Degrees,
z_src: Degrees,
k_ext: MagnitudesPerAirmass,
) -> Nanolamberts {
let k_ext = k_ext.value();
let i_star = lunar_illuminance_outside_atmosphere(alpha);
let f_rho = scattering_function(rho);
let am_moon = airmass::<KrisciunasSchaeferAirmass>(z_moon.to::<Radian>());
let am_src = airmass::<KrisciunasSchaeferAirmass>(z_src.to::<Radian>());
let trans_moon = 10f64.powf(-0.4 * k_ext * am_moon.value());
let absorb_path = 1.0 - 10f64.powf(-0.4 * k_ext * am_src.value());
Nanolamberts::new(f_rho * i_star * trans_moon * absorb_path)
}
#[cfg(test)]
mod tests {
use super::*;
use crate::units::angular::Radians;
use crate::units::photometry::s10_to_surface_brightness;
use crate::units::radiometry::S10s;
use siderust::qtty::IlluminationFractions;
fn make_phase(alpha_deg: f64) -> MoonPhaseGeometry {
MoonPhaseGeometry {
phase_angle: Radians::new(alpha_deg.to_radians()),
illuminated_fraction: IlluminationFractions::new(
0.5 * (1.0 + alpha_deg.to_radians().cos()),
),
elongation: Radians::new(0.0),
waxing: true,
}
}
fn inputs(alpha_deg: f64, rho_deg: f64, z_moon: f64, z_src: f64) -> MoonlightGeometry {
MoonlightGeometry {
separation: Degrees::new(rho_deg),
moon_zenith: Degrees::new(z_moon),
phase: make_phase(alpha_deg),
source_zenith: Degrees::new(z_src),
moon_distance: siderust::event::lunar::photometry::MEAN_MOON_DISTANCE,
}
}
fn compute(inp: &MoonlightGeometry) -> Result<MoonOutputs> {
compute_krisciunas_schaefer_1991(inp, DEFAULT_K_EXT)
}
fn v_mag_arcsec2(out: &MoonOutputs) -> f64 {
s10_to_surface_brightness(out.v_flux_s10, NSB_S10_ZP.value()).value()
}
#[test]
fn full_moon_reference_geometry_matches_published_brightness() {
let fixture = include_str!("../../../tests/data/external_reference_cases.csv");
let row = fixture.lines().nth(1).expect("external reference row");
let columns: Vec<_> = row.split(',').collect();
assert_eq!(columns.len(), 10, "external fixture schema changed");
assert_eq!(columns[0], "ks91_full_moon_zenith_90deg");
assert!(columns[1].contains("Krisciunas and Schaefer 1991"));
assert_eq!(columns[4], "mag_per_arcsec2");
assert_eq!(columns[5], "Johnson-V-approximation");
assert_eq!(columns[9], "model-choice");
let expected: f64 = columns[6].parse().expect("expected value");
let tolerance: f64 = columns[7].parse().expect("absolute tolerance");
let out = compute(&inputs(0.0, 90.0, 45.0, 45.0)).unwrap();
assert!(out.v_flux_s10 > S10s::zero());
let v_mag = v_mag_arcsec2(&out);
assert!(
(v_mag - expected).abs() < tolerance,
"V_sky = {v_mag:.2} mag/arcsec² not within {tolerance} of published {expected}"
);
}
#[test]
fn new_moon_contribution_is_negligible_relative_to_full() {
let out_new = compute(&inputs(180.0, 90.0, 45.0, 45.0)).unwrap();
let out_full = compute(&inputs(0.0, 90.0, 45.0, 45.0)).unwrap();
assert!(out_full.v_flux_s10 > S10s::zero());
assert!(
out_new.v_flux_s10 < out_full.v_flux_s10 * 1e-3,
"new-moon V S10 ({}) should be << full-moon ({})",
out_new.v_flux_s10.value(),
out_full.v_flux_s10.value()
);
}
#[test]
fn scattering_function_exhibits_expected_behavior() {
let inputs_at_rho = |rho| compute(&inputs(45.0, rho, 30.0, 60.0)).unwrap();
let b5 = inputs_at_rho(5.0).v_flux_s10;
let b90 = inputs_at_rho(90.0).v_flux_s10;
let b120 = inputs_at_rho(120.0).v_flux_s10;
let b175 = inputs_at_rho(175.0).v_flux_s10;
assert!(
b5 > S10s::zero(),
"brightness at rho=5 deg must be positive"
);
assert!(
b90 > S10s::zero(),
"brightness at rho=90 deg must be positive"
);
assert!(
b90 < b5,
"brightness at rho=90 deg should be less than forward scattering"
);
assert!(
b120 > b90,
"brightness should increase from rho=90 deg to rho=120 deg"
);
assert!(b175 > b120, "brightness should increase toward rho=180 deg");
assert!(
b175 < b5,
"backscattering should be weaker than forward scattering"
);
}
#[test]
fn moon_below_horizon_returns_zero() {
let out = compute(&inputs(0.0, 30.0, 95.0, 30.0)).unwrap();
assert_eq!(out.v_flux_s10, S10s::zero());
assert_eq!(out.integrated.value(), 0.0);
}
#[test]
fn airmass_at_zenith_is_unity() {
let am = airmass::<KrisciunasSchaeferAirmass>(Degrees::new(0.0).to::<Radian>());
assert!((am.value() - 1e0).abs() < 1e-12, "X(0) = {:?}", am);
}
}