spectrafit_models/
tauc.rs

1use crate::Model;
2
3/// Tauc optical band-gap edge: `A · ((x − e_gap) · H(x − e_gap))^p`.
4///
5/// Parameters (in order): `[amplitude, e_gap, exponent]`
6///
7/// - `amplitude` (`A`) scales the absorption above the gap.
8/// - `e_gap` is the optical band-gap energy: the edge onset. Below it the model is
9///   exactly `0` (Heaviside `H(x − e_gap)`); at and above it the absorption rises as a
10///   power law in the excess energy `(x − e_gap)`.
11/// - `exponent` (`p`) is the Tauc power (`p = 2` for an allowed indirect transition,
12///   `p = 1/2` for an allowed direct one). The classic Tauc plot fits `(α·hν)^{1/p}`
13///   linearly in `hν`; here we expose the forward edge `A·(hν − E_g)^p` directly.
14///
15/// The kernel is `0` for `x ≤ e_gap` (where the base `(x − e_gap)` is non-positive and
16/// a fractional `exponent` would otherwise be NaN). The numpy benchmark formula is
17/// identical — `np.where(x > e_gap, A·(x − e_gap)^p, 0)` — so numpy↔Rust parity is exact.
18pub struct Tauc;
19
20impl Model for Tauc {
21    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
22        let (a, e_gap, exponent) = (params[0], params[1], params[2]);
23        let excess = x[0] - e_gap;
24        if excess > 0.0 {
25            a * excess.powf(exponent)
26        } else {
27            0.0
28        }
29    }
30
31    /// Central finite-difference Jacobian. The Heaviside cut-off makes the closed-form
32    /// derivative piecewise (and the power law's `∂/∂e_gap` diverges as `x → e_gap⁺` for
33    /// `exponent < 1`), so a numerical Jacobian is used — matching `log_normal`.
34    fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64> {
35        let mut p = params.to_vec();
36        (0..params.len())
37            .map(|i| {
38                let h = 1e-7_f64 * params[i].abs().max(1e-7);
39                p[i] = params[i] + h;
40                let f_plus = self.eval(x, &p);
41                p[i] = params[i] - h;
42                let f_minus = self.eval(x, &p);
43                p[i] = params[i];
44                (f_plus - f_minus) / (2.0 * h)
45            })
46            .collect()
47    }
48
49    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
50        vec!["amplitude".into(), "e_gap".into(), "exponent".into()]
51    }
52}
53
54#[cfg(test)]
55mod tests {
56    use super::*;
57    use approx::assert_relative_eq;
58
59    #[test]
60    fn below_gap_is_zero() {
61        let m = Tauc;
62        assert_eq!(m.eval(&[1.0], &[2.0, 1.5, 2.0]), 0.0);
63        assert_eq!(m.eval(&[1.5], &[2.0, 1.5, 2.0]), 0.0); // exactly at the gap
64    }
65
66    #[test]
67    fn above_gap_power_law() {
68        // (x − e_gap) = 2.0, p = 2 ⇒ A·4.0.
69        let m = Tauc;
70        assert_relative_eq!(m.eval(&[3.5], &[2.0, 1.5, 2.0]), 2.0 * 4.0, epsilon = 1e-12);
71    }
72
73    #[test]
74    fn param_names_are_canonical() {
75        assert_eq!(
76            Tauc.param_names()
77                .iter()
78                .map(|c| c.as_ref())
79                .collect::<Vec<_>>(),
80            &["amplitude", "e_gap", "exponent"]
81        );
82    }
83
84    #[test]
85    fn jacobian_shape_and_amplitude() {
86        let m = Tauc;
87        let j = m.jacobian(&[3.5], &[2.0, 1.5, 2.0]);
88        assert_eq!(j.len(), 3);
89        // ∂/∂amplitude above the gap is (x − e_gap)^p = 4.0.
90        assert_relative_eq!(j[0], 4.0, epsilon = 1e-6);
91    }
92}
spectrafit_models/tauc.rs

spectrafit_models/
tauc.rs
