spectrafit_models/
fano.rs

1use crate::Model;
2
3/// Fano resonance lineshape used in XPS asymmetric peaks.
4///
5/// Formula: `A · (q + ε)² / (1 + ε²)`,  ε = (x − x₀) / Γ
6///
7/// Parameters (in order): `[amplitude, center, gamma, q]`
8///
9/// - `amplitude` — overall scale factor `A`
10/// - `center`    — resonance centre `x₀`
11/// - `gamma`     — half-width at half-maximum `Γ` (> 0)
12/// - `q`         — Fano asymmetry parameter (dimensionless)
13///
14/// Jacobians are fully analytical.
15pub struct Fano;
16
17impl Model for Fano {
18    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
19        let (a, x0, gamma, q) = (params[0], params[1], params[2], params[3]);
20        let eps = (x[0] - x0) / gamma;
21        let num = (q + eps) * (q + eps);
22        let den = 1.0 + eps * eps;
23        a * num / den
24    }
25
26    fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64> {
27        let (a, x0, gamma, q) = (params[0], params[1], params[2], params[3]);
28        let dx = x[0] - x0;
29        let eps = dx / gamma;
30        let qe = q + eps;
31        let e2_1 = 1.0 + eps * eps;
32        let fano = qe * qe / e2_1; // normalised Fano value
33
34        // ∂/∂amplitude
35        let da = fano;
36
37        // d(fano)/d(eps): use quotient rule
38        // f(ε) = (q+ε)²/(1+ε²)
39        // f'(ε) = [2(q+ε)(1+ε²) − (q+ε)²·2ε] / (1+ε²)²
40        //       = 2(q+ε)[(1+ε²) − (q+ε)ε] / (1+ε²)²
41        let dfano_deps = 2.0 * qe * (e2_1 - qe * eps) / (e2_1 * e2_1);
42
43        // ε = (x − x₀)/γ  →  ∂ε/∂x₀ = −1/γ,  ∂ε/∂γ = −ε/γ
44        let dc = a * dfano_deps * (-1.0 / gamma);
45        let dg = a * dfano_deps * (-eps / gamma);
46
47        // ∂/∂q:  d/dq [(q+ε)²/(1+ε²)] = 2(q+ε)/(1+ε²)
48        let dq = a * 2.0 * qe / e2_1;
49
50        vec![da, dc, dg, dq]
51    }
52
53    #[inline]
54    fn jacobian_into(&self, x: &[f64], params: &[f64], out: &mut [f64]) {
55        let (a, x0, gamma, q) = (params[0], params[1], params[2], params[3]);
56        let dx = x[0] - x0;
57        let eps = dx / gamma;
58        let qe = q + eps;
59        let e2_1 = 1.0 + eps * eps;
60        let fano = qe * qe / e2_1;
61        let dfano_deps = 2.0 * qe * (e2_1 - qe * eps) / (e2_1 * e2_1);
62        out[0] = fano;
63        out[1] = a * dfano_deps * (-1.0 / gamma);
64        out[2] = a * dfano_deps * (-eps / gamma);
65        out[3] = a * 2.0 * qe / e2_1;
66    }
67
68    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
69        vec![
70            "amplitude".into(),
71            "center".into(),
72            "gamma".into(),
73            "q".into(),
74        ]
75    }
76}
77
78#[cfg(test)]
79mod tests {
80    use super::*;
81    use approx::assert_relative_eq;
82
83    #[test]
84    fn at_center_eps_zero() {
85        // ε=0 → (q+0)²/(1+0) = q², so f = A·q²
86        let v = Fano.eval(&[0.0], &[1.0, 0.0, 1.0, 2.0]);
87        assert_relative_eq!(v, 4.0, epsilon = 1e-12);
88    }
89
90    #[test]
91    fn zero_at_eps_minus_q() {
92        // ε = −q → numerator = 0 → f = 0
93        let v = Fano.eval(&[-2.0], &[1.0, 0.0, 1.0, 2.0]);
94        assert_relative_eq!(v, 0.0, epsilon = 1e-12);
95    }
96
97    #[test]
98    fn jacobian_shape() {
99        let j = Fano.jacobian(&[0.5], &[1.0, 0.0, 1.0, 1.5]);
100        assert_eq!(j.len(), 4);
101    }
102
103    #[test]
104    fn jacobian_finite_diff_check() {
105        // Upgraded from a single-point forward difference (h=1e-5,
106        // epsilon=1e-4, O(h) truncation error the size of its own tolerance)
107        // to a central difference (O(h^2)) swept over three parameter
108        // regimes and several x-points relative to the resonance centre.
109        let param_sets = [
110            [2.0, 0.2, 0.8, 1.5],   // nominal
111            [1e-3, 0.0, 0.05, 0.5], // small amplitude/width
112            [5.0, -2.0, 3.0, 10.0], // offset centre, wide, large asymmetry
113        ];
114        for p in param_sets {
115            let (center, gamma) = (p[1], p[2]);
116            for &mult in &[-3.0_f64, -1.0, -0.1, 0.0, 0.1, 1.0, 3.0] {
117                let x = [center + mult * gamma];
118                let j_anal = Fano.jacobian(&x, &p);
119                for i in 0..p.len() {
120                    let h = 1e-6 * p[i].abs().max(1.0);
121                    let (mut a, mut b) = (p, p);
122                    a[i] += h;
123                    b[i] -= h;
124                    let fd = (Fano.eval(&x, &a) - Fano.eval(&x, &b)) / (2.0 * h);
125                    assert_relative_eq!(j_anal[i], fd, epsilon = 1e-7, max_relative = 1e-6);
126                }
127            }
128        }
129    }
130
131    #[test]
132    fn da_equals_fano_value_over_amplitude() {
133        // ∂/∂amplitude is just the normalised lineshape
134        let p = [3.0, 0.0, 1.0, 1.0];
135        let x = [0.5];
136        let j = Fano.jacobian(&x, &p);
137        let expected_da = Fano.eval(&x, &[1.0, p[1], p[2], p[3]]);
138        assert_relative_eq!(j[0], expected_da, epsilon = 1e-12);
139    }
140
141    // ----- Limiting-case asymptotic (ground-truth verification) -----
142    //
143    // For very large q, the Fano lineshape collapses to a Lorentzian scaled by q²:
144    //
145    //     A·(q+ε)²/(1+ε²)  →  A·q²/(1+ε²)         as q → ∞
146    //
147    // So Fano(A=1, q=Q)/Q² should converge to a Lorentzian of unit amplitude as
148    // Q grows. The remaining error is the 2qε + ε² subleading term, of relative
149    // order 1/q. Test at q = 1e6 → relative error ≤ ~1e-5.
150
151    #[test]
152    fn fano_q_large_normalised_matches_lorentzian() {
153        use crate::lorentzian::Lorentzian;
154        let x = [0.4_f64];
155        let q = 1.0e6;
156        let p_fano = [1.0, 0.0, 1.0, q]; // A=1, c=0, γ=1, q huge
157        let p_lor = [1.0, 0.0, 1.0]; // A=1, c=0, σ=γ=1
158        let fano_normalised = Fano.eval(&x, &p_fano) / (q * q);
159        let lor = Lorentzian.eval(&x, &p_lor);
160        // Subleading 1/q term gives ~1e-6 relative error at q=1e6.
161        assert_relative_eq!(fano_normalised, lor, epsilon = 1e-5);
162    }
163}
spectrafit_models/fano.rs

spectrafit_models/
fano.rs
