spectrafit_models/
harmonic_ir.rs

1use crate::Model;
2
3/// Driven damped harmonic-oscillator IR absorption: `A / ((c² − x²)² + (σ·x)²)`.
4///
5/// Parameters (in order): `[amplitude, center, sigma]` — `center` is the resonance frequency,
6/// `sigma` the damping; `amplitude` is a scale (peak ≠ A). Reuses the canonical
7/// amplitude/center/sigma names. numpy oracle identical: `A / ((c**2 - x**2)**2 + (σ*x)**2)`.
8pub struct HarmonicIr;
9
10impl Model for HarmonicIr {
11    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
12        let (a, c, sigma) = (params[0], params[1], params[2]);
13        let d = c * c - x[0] * x[0];
14        a / (d * d + (sigma * x[0]) * (sigma * x[0]))
15    }
16
17    fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64> {
18        let mut p = params.to_vec();
19        (0..params.len())
20            .map(|i| {
21                let h = 1e-7_f64 * params[i].abs().max(1e-7);
22                p[i] = params[i] + h;
23                let f_plus = self.eval(x, &p);
24                p[i] = params[i] - h;
25                let f_minus = self.eval(x, &p);
26                p[i] = params[i];
27                (f_plus - f_minus) / (2.0 * h)
28            })
29            .collect()
30    }
31
32    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
33        vec!["amplitude".into(), "center".into(), "sigma".into()]
34    }
35}
36
37#[cfg(test)]
38mod tests {
39    use super::*;
40    use approx::assert_relative_eq;
41
42    #[test]
43    fn eval_at_resonance() {
44        // At x == center: A / ((0)² + (σ·c)²) = A / (σ·c)². With A=1,c=2,σ=0.5 → 1/1 = 1.
45        assert_relative_eq!(
46            HarmonicIr.eval(&[2.0], &[1.0, 2.0, 0.5]),
47            1.0,
48            epsilon = 1e-12
49        );
50    }
51
52    #[test]
53    fn param_names_and_jacobian() {
54        assert_eq!(
55            HarmonicIr
56                .param_names()
57                .iter()
58                .map(|c| c.as_ref())
59                .collect::<Vec<_>>(),
60            &["amplitude", "center", "sigma"]
61        );
62        assert_eq!(HarmonicIr.jacobian(&[1.0], &[1.0, 2.0, 0.5]).len(), 3);
63    }
64
65    // ----- Limiting-case asymptotic (ground-truth verification) -----
66    //
67    // The undamped limit (σ = 0) reduces to a closed-form rational:
68    //
69    //     HarmonicIr(A, c, σ=0)  =  A / (c² − x²)²   for x ≠ ±c
70    //
71    // The denominator vanishes at x = ±c (the resonance), so test only at
72    // off-resonance points. Catches: wrong squaring of the damping term,
73    // missing parenthesization that would couple σ into the undamped form.
74
75    #[test]
76    fn sigma_zero_undamped_matches_closed_form() {
77        let m = HarmonicIr;
78        let a = 2.0_f64;
79        let c = 1.5_f64;
80        let p = [a, c, 0.0]; // undamped
81        for &xi in &[0.5_f64, 1.0, 2.0, 3.0, 5.0] {
82            let d = c * c - xi * xi;
83            let expected = a / (d * d);
84            assert_relative_eq!(m.eval(&[xi], &p), expected, epsilon = 1e-12);
85        }
86    }
87}
spectrafit_models/harmonic_ir.rs

spectrafit_models/
harmonic_ir.rs
