Skip to main content

spectrafit_models/
breit_wigner.rs

1use crate::Model;
2
3/// Breit-Wigner-Fano (BWF) resonance: `A·(q·g + (x − c))² / (g² + (x − c)²)`, `g = σ/2`.
4///
5/// Parameters (in order): `[amplitude, center, sigma, q]`. `amplitude` is a scale factor
6/// (peak ≠ A — the asymmetric Fano family, like `fano`); `q` is the asymmetry. numpy oracle
7/// identical: `g = σ/2; A*(q*g + (x-c))**2 / (g*g + (x-c)**2)`.
8pub struct BreitWigner;
9
10impl Model for BreitWigner {
11    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
12        let (a, c, sigma, q) = (params[0], params[1], params[2], params[3]);
13        let g = sigma / 2.0;
14        let dx = x[0] - c;
15        a * (q * g + dx) * (q * g + dx) / (g * g + dx * dx)
16    }
17
18    fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64> {
19        let mut p = params.to_vec();
20        (0..params.len())
21            .map(|i| {
22                let h = 1e-7_f64 * params[i].abs().max(1e-7);
23                p[i] = params[i] + h;
24                let f_plus = self.eval(x, &p);
25                p[i] = params[i] - h;
26                let f_minus = self.eval(x, &p);
27                p[i] = params[i];
28                (f_plus - f_minus) / (2.0 * h)
29            })
30            .collect()
31    }
32
33    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
34        vec![
35            "amplitude".into(),
36            "center".into(),
37            "sigma".into(),
38            "q".into(),
39        ]
40    }
41}
42
43#[cfg(test)]
44mod tests {
45    use super::*;
46    use approx::assert_relative_eq;
47
48    #[test]
49    fn eval_at_center_is_amp_times_q_squared() {
50        // At x == center, dx = 0 ⇒ A·(q·g)² / g² = A·q².
51        let v = BreitWigner.eval(&[1.5], &[2.0, 1.5, 1.0, 1.7]);
52        assert_relative_eq!(v, 2.0 * 1.7 * 1.7, epsilon = 1e-12);
53    }
54
55    #[test]
56    fn param_names_and_jacobian() {
57        assert_eq!(
58            BreitWigner
59                .param_names()
60                .iter()
61                .map(|c| c.as_ref())
62                .collect::<Vec<_>>(),
63            &["amplitude", "center", "sigma", "q"]
64        );
65        assert_eq!(BreitWigner.jacobian(&[1.0], &[2.0, 1.5, 1.0, 1.7]).len(), 4);
66    }
67}