siderust 0.7.0

// SPDX-License-Identifier: AGPL-3.0-or-later
// Copyright (C) 2026 Vallés Puig, Ramon

//! # Mie / aerosol optical depth — power-law parameterisation
//!
//! ## Scientific scope
//!
//! Aerosol particles in the lower atmosphere scatter and absorb optical
//! photons through Mie theory; for routine sky-brightness work the
//! wavelength dependence of the resulting column-integrated optical depth
//! is well-approximated by the Patat (2011) power law
//!
//! ```text
//! τ_M(λ) = τ₀ · (λ / λ_ref)^α
//! ```
//!
//! where `τ₀` is the optical depth at a reference wavelength `λ_ref`
//! (typically 550 nm) and `α` is a (typically negative) wavelength
//! exponent. This module exposes [`MieParams`] for the three numbers and
//! [`mie_optical_depth`] for the evaluation, plus four observatory
//! presets covering Cerro Paranal, La Palma (ORM), Mauna Kea, and La
//! Silla.
//!
//! ## Technical scope
//!
//! - `τ₀` is stored as a typed [`OpticalDepths`].
//! - `λ_ref` is stored as a typed [`Nanometers`].
//! - `α` is left as a raw `f64` because it is a true dimensionless
//!   exponent.
//! - The evaluator returns [`OpticalDepths`].
//!
//! ## References
//!
//! - Patat, F., et al. (2011). "Optical atmospheric extinction over Cerro
//!   Paranal: 6 years of night-sky observations". *A&A* 527, A91.
//! - Lombardi, G., et al. (2008). *PASP* 120, 212.
//! - García-Gil, A., et al. (2010). *PASP* 122, 1109.
//! - Krisciunas, K. (1990). *PASP* 102, 1235.
//! - Burki, G., et al. (1995). *A&AS* 112, 383.

use crate::ext_qtty::length::Nanometers;
use crate::provenance::{DataSource, Provenance};
use crate::qtty::OpticalDepths;

/// Aerosol optical-depth model parameters.
///
/// The optical depth is `τ_M(λ) = τ₀ · (λ / λ_ref)^α`. Negative `α`
/// encodes the usual aerosol behaviour (more extinction at shorter
/// wavelengths).
///
/// # Reference
///
/// Patat, F. (2011), "The dancing sky: 6 years of night-sky observations
/// at Cerro Paranal", *A&A* 527, A91.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct MieParams {
    /// Optical depth at the reference wavelength.
    pub tau0: OpticalDepths,
    /// Wavelength exponent `α` (typically negative). True dimensionless
    /// exponent — kept as `f64` because raising a typed quantity to a
    /// non-integer power has no `qtty` representation.
    pub alpha: f64,
    /// Reference wavelength `λ_ref`.
    pub lambda_ref: Nanometers,
}

impl MieParams {
    /// Cerro Paranal default used by `darknsb` / NSB
    /// (`τ₀ = 0.05`, `α = -1.38`, `λ_ref = 550 nm`).
    ///
    /// See [`crate::observatories::EL_PARANAL`].
    pub const PARANAL: MieParams = MieParams {
        tau0: OpticalDepths::new(0.05),
        alpha: -1.38,
        lambda_ref: Nanometers::new(550.0),
    };

    /// Roque de los Muchachos Observatory, La Palma (ORM, Canary Islands,
    /// Spain).
    ///
    /// `τ₀ = 0.02` at `λ_ref = 550 nm`, `α = −1.38`, giving
    /// AOD(500 nm) ≈ 0.022. Lombardi et al. 2008 + García-Gil et al. 2010
    /// median; not a single-night calibration.
    ///
    /// See [`crate::observatories::ROQUE_DE_LOS_MUCHACHOS`].
    pub const LA_PALMA: MieParams = MieParams {
        tau0: OpticalDepths::new(0.02),
        alpha: -1.38,
        lambda_ref: Nanometers::new(550.0),
    };

    /// Mauna Kea Observatories (Hawaiʻi, USA) — CFHT / Subaru site.
    ///
    /// `τ₀ = 0.03` at `λ_ref = 550 nm`, `α = −1.38`, giving
    /// AOD(500 nm) ≈ 0.034. Krisciunas 1990 representative median.
    ///
    /// See [`crate::observatories::MAUNA_KEA`].
    pub const MAUNA_KEA: MieParams = MieParams {
        tau0: OpticalDepths::new(0.03),
        alpha: -1.38,
        lambda_ref: Nanometers::new(550.0),
    };

    /// La Silla Observatory (ESO, Chile).
    ///
    /// `τ₀ = 0.04` at `λ_ref = 550 nm`, `α = −1.38`, giving
    /// AOD(500 nm) ≈ 0.045. Burki et al. 1995 representative median.
    ///
    /// See [`crate::observatories::LA_SILLA_OBSERVATORY`].
    pub const LA_SILLA: MieParams = MieParams {
        tau0: OpticalDepths::new(0.04),
        alpha: -1.38,
        lambda_ref: Nanometers::new(550.0),
    };

    /// Provenance record for [`MieParams::PARANAL`].
    pub fn paranal_provenance() -> Provenance {
        Provenance {
            source: Some(DataSource::LiteratureCitation {
                bibkey: "patat2011".to_string(),
                doi: Some("10.1051/0004-6361/201015537".to_string()),
            }),
            notes: Some(
                "Cerro Paranal aerosol τ₀ = 0.05 at 550 nm (α = −1.38). \
                 Patat 2011 (A&A 527, A91) Table 1 / darknsb model. \
                 Representative median of photometric nights over 6 years."
                    .to_string(),
            ),
            ..Provenance::default()
        }
    }

    /// Provenance record for [`MieParams::LA_PALMA`].
    pub fn la_palma_provenance() -> Provenance {
        Provenance {
            source: Some(DataSource::LiteratureCitation {
                bibkey: "lombardi2008".to_string(),
                doi: Some("10.1086/526338".to_string()),
            }),
            notes: Some(
                "La Palma ORM aerosol τ₀ = 0.02 at 550 nm (α = −1.38). \
                 Lombardi et al. 2008 (PASP 120, 212) median k_V ≈ 0.11 mag/airmass; \
                 Mie component ≈ 0.02 after Rayleigh subtraction. \
                 Confirmed by García-Gil et al. 2010 (PASP 122, 1109; DOI 10.1086/656169). \
                 Representative median, not a single-night calibration."
                    .to_string(),
            ),
            ..Provenance::default()
        }
    }

    /// Provenance record for [`MieParams::MAUNA_KEA`].
    pub fn mauna_kea_provenance() -> Provenance {
        Provenance {
            source: Some(DataSource::LiteratureCitation {
                bibkey: "krisciunas1990".to_string(),
                doi: Some("10.1086/132757".to_string()),
            }),
            notes: Some(
                "Mauna Kea aerosol τ₀ = 0.03 at 550 nm (α = −1.38). \
                 Krisciunas 1990 (PASP 102, 1235) V-band extinction at Mauna Kea; \
                 Mie component ≈ 0.02–0.04 after Rayleigh subtraction at 615 hPa. \
                 Consistent with CFHT sky-quality monitoring (Cuillandre et al.). \
                 Representative median, not a single-night calibration."
                    .to_string(),
            ),
            ..Provenance::default()
        }
    }

    /// Provenance record for [`MieParams::LA_SILLA`].
    pub fn la_silla_provenance() -> Provenance {
        Provenance {
            source: Some(DataSource::LiteratureCitation {
                bibkey: "burki1995".to_string(),
                doi: None,
            }),
            notes: Some(
                "La Silla aerosol τ₀ = 0.04 at 550 nm (α = −1.38). \
                 Burki et al. 1995 (A&AS 112, 383; bibcode 1995A&AS..112..383B) \
                 broad-band extinction at La Silla; Mie component ≈ 0.03–0.06 \
                 after Rayleigh subtraction at 760 hPa. See also Rufener 1986 \
                 (A&A 165, 275; bibcode 1986A&A...165..275R). \
                 Representative median, not a single-night calibration."
                    .to_string(),
            ),
            ..Provenance::default()
        }
    }
}

/// Mie / aerosol optical depth at the given wavelength using the
/// power-law parameterisation stored in `params`.
///
/// Returns a typed [`OpticalDepths`].
#[inline]
pub fn mie_optical_depth(params: &MieParams, wavelength: Nanometers) -> OpticalDepths {
    let ratio: f64 = wavelength / params.lambda_ref;
    params.tau0 * ratio.powf(params.alpha)
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn at_lambda_ref_returns_tau0() {
        let p = MieParams::PARANAL;
        assert_eq!(
            mie_optical_depth(&p, Nanometers::new(550.0)).value(),
            p.tau0.value()
        );
    }

    #[test]
    fn shorter_wavelength_higher_tau() {
        let p = MieParams::PARANAL;
        let blue = mie_optical_depth(&p, Nanometers::new(400.0)).value();
        let red = mie_optical_depth(&p, Nanometers::new(700.0)).value();
        assert!(blue > p.tau0.value());
        assert!(red < p.tau0.value());
    }

    #[test]
    fn site_presets_differ_from_paranal() {
        assert_ne!(
            MieParams::LA_PALMA.tau0.value(),
            MieParams::PARANAL.tau0.value()
        );
        assert_ne!(
            MieParams::MAUNA_KEA.tau0.value(),
            MieParams::PARANAL.tau0.value()
        );
        assert_ne!(
            MieParams::LA_SILLA.tau0.value(),
            MieParams::PARANAL.tau0.value()
        );
    }

    #[test]
    fn site_presets_optical_depth_ordering() {
        let lambda = Nanometers::new(550.0);
        let tau_lp = mie_optical_depth(&MieParams::LA_PALMA, lambda).value();
        let tau_mk = mie_optical_depth(&MieParams::MAUNA_KEA, lambda).value();
        let tau_ls = mie_optical_depth(&MieParams::LA_SILLA, lambda).value();
        let tau_pn = mie_optical_depth(&MieParams::PARANAL, lambda).value();

        for tau in [tau_lp, tau_mk, tau_ls, tau_pn] {
            assert!(
                tau > 0.0 && tau.is_finite(),
                "tau must be positive finite: {tau}"
            );
        }

        assert!(tau_lp < tau_mk, "La Palma should be cleaner than Mauna Kea");
        assert!(tau_mk < tau_ls, "Mauna Kea should be cleaner than La Silla");
        assert!(tau_ls < tau_pn, "La Silla should be cleaner than Paranal");
    }

    #[test]
    fn provenance_non_empty() {
        assert!(MieParams::paranal_provenance().source.is_some());
        assert!(MieParams::la_palma_provenance().source.is_some());
        assert!(MieParams::mauna_kea_provenance().source.is_some());
        assert!(MieParams::la_silla_provenance().source.is_some());
    }

    #[test]
    fn aod_500nm_within_published_ranges() {
        let lambda = Nanometers::new(500.0);
        let aod_lp = mie_optical_depth(&MieParams::LA_PALMA, lambda).value();
        let aod_mk = mie_optical_depth(&MieParams::MAUNA_KEA, lambda).value();
        let aod_ls = mie_optical_depth(&MieParams::LA_SILLA, lambda).value();

        assert!(
            (0.01..=0.05).contains(&aod_lp),
            "La Palma AOD(500nm) = {aod_lp:.4} outside 0.01–0.05"
        );
        assert!(
            (0.02..=0.05).contains(&aod_mk),
            "Mauna Kea AOD(500nm) = {aod_mk:.4} outside 0.02–0.05"
        );
        assert!(
            (0.03..=0.07).contains(&aod_ls),
            "La Silla AOD(500nm) = {aod_ls:.4} outside 0.03–0.07"
        );
    }
}