lindhard/electron/inelastic/table/
sampler.rs

1//! The momentum-transfer sampler of the inelastic table (parent module
2//! docs).
3
4use super::super::penn::{hartree_ev, SinglePolePenn};
5use super::density::{adaptive_density, LinearDensity};
6use super::{check_no_exchange, invalid, InelasticTableError, DEFAULT_Q_TOLERANCE};
7use crate::constants::BOHR_RADIUS;
8
9/// Samples the momentum transfer `q` of an inelastic collision with incident
10/// kinetic energy `E` and loss `W` (module docs). Built once per `(E, W)`.
11#[derive(Debug, Clone, PartialEq)]
12pub struct MomentumTransferSampler {
13    density: LinearDensity,
14    /// Hartree-unit wavenumbers `k` and `k'`, and `T' = E + E_F`.
15    k: f64,
16    kp: f64,
17    tp: f64,
18}
19
20impl MomentumTransferSampler {
21    /// Build the sampler with the default tolerance [`DEFAULT_Q_TOLERANCE`].
22    pub fn new(
23        penn: &SinglePolePenn,
24        energy_ev: f64,
25        loss_ev: f64,
26    ) -> Result<Self, InelasticTableError> {
27        Self::with_tolerance(penn, energy_ev, loss_ev, DEFAULT_Q_TOLERANCE)
28    }
29
30    /// Build the sampler with density-model tolerance `tol` in
31    /// `[1e-9, 1e-2]`. Fails if `0 < W <= E` does not hold or the slice has
32    /// no support (`W` below the ELF table, or no `q` allowed).
33    pub fn with_tolerance(
34        penn: &SinglePolePenn,
35        energy_ev: f64,
36        loss_ev: f64,
37        tol: f64,
38    ) -> Result<Self, InelasticTableError> {
39        check_no_exchange(penn)?;
40        if !(energy_ev.is_finite() && energy_ev > 0.0) {
41            return Err(invalid("electron energy", format!("{energy_ev} eV")));
42        }
43        if !(loss_ev.is_finite() && loss_ev > 0.0 && loss_ev <= energy_ev) {
44            return Err(invalid(
45                "energy loss",
46                format!("{loss_ev} eV (need 0 < W <= E = {energy_ev} eV)"),
47            ));
48        }
49        if !(tol.is_finite() && (1.0e-9..=1.0e-2).contains(&tol)) {
50            return Err(invalid(
51                "q tolerance",
52                format!("{tol} (need 1e-9 ..= 1e-2)"),
53            ));
54        }
55        let h = hartree_ev();
56        let (t, w) = (energy_ev / h, loss_ev / h);
57        let breaks = penn.q_slice_breaks_au(t, w).ok_or_else(|| {
58            invalid(
59                "momentum-transfer slice",
60                format!("no allowed q at E = {energy_ev} eV, W = {loss_ev} eV"),
61            )
62        })?;
63        let f = |u: f64| penn.q_integrand_au(w, u);
64        // The node bound is an error for the loss rows of a table (#345); a
65        // momentum-transfer slice keeps the density it reached, as before.
66        let (density, _) = adaptive_density(&f, &breaks, 8, tol, false, &[], 0.0);
67        if !(density.total() > 0.0 && density.total().is_finite()) {
68            return Err(invalid(
69                "momentum-transfer slice",
70                format!("zero loss function at E = {energy_ev} eV, W = {loss_ev} eV"),
71            ));
72        }
73        let tp = t + penn.fermi_energy_ev() / h;
74        Ok(Self {
75            density,
76            k: (2.0 * tp).sqrt(),
77            kp: (2.0 * (tp - w)).max(0.0).sqrt(),
78            tp,
79        })
80    }
81
82    /// The momentum transfer (as a wavenumber, m⁻¹) at cumulative
83    /// probability `u` in `[0, 1]` (clamped).
84    pub fn sample(&self, u: f64) -> f64 {
85        self.density.quantile(u.clamp(0.0, 1.0)).exp() / BOHR_RADIUS
86    }
87
88    /// The support `(q_min, q_max)` of the sampler, m⁻¹: the kinematic
89    /// limits `q±` cut to the tabulated ELF range.
90    pub fn q_range_per_m(&self) -> (f64, f64) {
91        let n = self.density.x.len();
92        (
93            self.density.x[0].exp() / BOHR_RADIUS,
94            self.density.x[n - 1].exp() / BOHR_RADIUS,
95        )
96    }
97
98    /// The kinematic limits `q- = k - k'` and `q+ = k + k'`, m⁻¹.
99    pub fn kinematic_limits_per_m(&self) -> (f64, f64) {
100        (
101            (self.k - self.kp) / BOHR_RADIUS,
102            (self.k + self.kp) / BOHR_RADIUS,
103        )
104    }
105
106    /// The cosine of the primary's deflection angle for momentum transfer
107    /// `q` (m⁻¹), `(k² + k'² - q²)/(2 k k')`, clamped to `[-1, 1]`.
108    pub fn cos_theta(&self, q_per_m: f64) -> f64 {
109        let q = q_per_m * BOHR_RADIUS;
110        let c = (self.k * self.k + self.kp * self.kp - q * q) / (2.0 * self.k * self.kp);
111        c.clamp(-1.0, 1.0)
112    }
113
114    /// The DIIMFP `p(E, W)` (m⁻¹ eV⁻¹) of the sampler's own density model:
115    /// the integral of the density over `ln q`, divided by `π T'`. It agrees
116    /// with [`SinglePolePenn::diimfp_per_m_ev`] to the model tolerance.
117    pub fn diimfp_per_m_ev(&self) -> f64 {
118        let h = hartree_ev();
119        self.density.total() / (std::f64::consts::PI * self.tp) / (BOHR_RADIUS * h)
120    }
121}
lindhard/electron/inelastic/table/sampler.rs

lindhard/electron/inelastic/table/
sampler.rs
