spectrafit_models/
moffat.rs

1use crate::Model;
2
3/// Moffat peak: `A / (((x − c)/σ)² + 1)^β`.
4///
5/// Parameters (in order): `[amplitude, center, sigma, beta]`. `amplitude` is the peak
6/// height at `x == center`; `beta` controls the wings (β→∞ Gaussian-like, small β heavier
7/// tails). numpy oracle identical: `A / (((x-c)/σ)**2 + 1)**β`.
8pub struct Moffat;
9
10impl Model for Moffat {
11    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
12        let (a, c, sigma, beta) = (params[0], params[1], params[2], params[3]);
13        let z = (x[0] - c) / sigma;
14        a / (z * z + 1.0).powf(beta)
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![
34            "amplitude".into(),
35            "center".into(),
36            "sigma".into(),
37            "beta".into(),
38        ]
39    }
40}
41
42#[cfg(test)]
43mod tests {
44    use super::*;
45    use approx::assert_relative_eq;
46
47    #[test]
48    fn eval_at_center_equals_amplitude() {
49        assert_relative_eq!(
50            Moffat.eval(&[1.5], &[3.0, 1.5, 0.8, 2.0]),
51            3.0,
52            epsilon = 1e-12
53        );
54    }
55
56    #[test]
57    fn param_names_and_jacobian() {
58        assert_eq!(
59            Moffat
60                .param_names()
61                .iter()
62                .map(|c| c.as_ref())
63                .collect::<Vec<_>>(),
64            &["amplitude", "center", "sigma", "beta"]
65        );
66        let j = Moffat.jacobian(&[1.5], &[3.0, 1.5, 0.8, 2.0]);
67        assert_eq!(j.len(), 4);
68        assert_relative_eq!(j[0], 1.0, epsilon = 1e-6); // ∂/∂A at center = 1
69    }
70}
spectrafit_models/moffat.rs

spectrafit_models/
moffat.rs
