feffit 0.1.0

Pure-Rust EXAFS toolkit — data reduction (pre-edge/normalize/AUTOBK), Fourier transforms, FEFF path fitting (feffit), and feff.inp build/run; a port of larch.xafs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

//! Modified Bessel function of the first kind, order zero, `I0(x)`.
//!
//! Direct port of the Cephes `i0` routine (public domain) — the same code
//! `scipy.special.i0` calls — so the Kaiser-Bessel FT window matches larch to
//! machine precision.

// The coefficient tables below are transcribed verbatim from Cephes; the digit
// groups are kept un-separated and at full upstream precision so they can be
// diffed line-by-line against the public-domain source. The excess digits past
// f64's ~17 sig-figs round to the same bits the compiler would pick anyway, so
// both readability lints are intentionally disabled for this reference data.
#![allow(clippy::unreadable_literal, clippy::excessive_precision)]

// Chebyshev coefficients for exp(-x) I0(x) on [0, 8].
#[rustfmt::skip]
const A: [f64; 30] = [
    -4.41534164647933937950E-18, 3.33079451882223809783E-17,
    -2.43127984654795469359E-16, 1.71539128555513303061E-15,
    -1.16853328779934516808E-14, 7.67618549860493561688E-14,
    -4.85644678311192946090E-13, 2.95505266312963983461E-12,
    -1.72682629144155570723E-11, 9.67580903537323691224E-11,
    -5.18979560163526290666E-10, 2.65982372468238665035E-9,
    -1.30002500998624804212E-8,  6.04699502254191894932E-8,
    -2.67079385394061173391E-7,  1.11738753912010371815E-6,
    -4.41673835845875056359E-6,  1.64484480707288970893E-5,
    -5.75419501008210370398E-5,  1.88502885095841655729E-4,
    -5.76375574538582365885E-4,  1.63947561694133579842E-3,
    -4.32430999505057594430E-3,  1.05464603945949983183E-2,
    -2.37374148058994688156E-2,  4.93052842396707084878E-2,
    -9.49010970480476444210E-2,  1.71620901522208775349E-1,
    -3.04682672343198398683E-1,  6.76795274409476084995E-1,
];

// Chebyshev coefficients for exp(-x) sqrt(x) I0(x) on [8, inf].
#[rustfmt::skip]
const B: [f64; 25] = [
    -7.23318048787475395456E-18, -4.83050448594418207126E-18,
    4.46562142029675999901E-17,  3.46122286769746109310E-17,
    -2.82762398051658348494E-16, -3.42548561967721913462E-16,
    1.77256013305652638360E-15,  3.81168066935262242075E-15,
    -9.55484669882830764870E-15, -4.15056934728722208663E-14,
    1.54008621752140982691E-14,  3.85277838274214270114E-13,
    7.18012445138366623367E-13,  -1.79417853150680611778E-12,
    -1.32158118404477131188E-11, -3.14991652796324136454E-11,
    1.18891471078464383424E-11,  4.94060238822496958910E-10,
    3.39623202570838634515E-9,   2.26666899049817806459E-8,
    2.04891858946906374183E-7,   2.89137052083475648297E-6,
    6.88975834691682398426E-5,   3.36911647825569408990E-3,
    8.04490411014108831608E-1,
];

/// Evaluate a Chebyshev series (Cephes `chbevl`).
fn chbevl(x: f64, coeffs: &[f64]) -> f64 {
    let mut b0 = coeffs[0];
    let mut b1 = 0.0_f64;
    let mut b2 = 0.0_f64;
    for &c in &coeffs[1..] {
        b2 = b1;
        b1 = b0;
        b0 = x * b1 - b2 + c;
    }
    0.5 * (b0 - b2)
}

/// Modified Bessel function `I0(x)` (Cephes parity with `scipy.special.i0`).
pub fn i0(x: f64) -> f64 {
    let x = x.abs();
    if x <= 8.0 {
        let y = (x / 2.0) - 2.0;
        x.exp() * chbevl(y, &A)
    } else {
        x.exp() * chbevl(32.0 / x - 2.0, &B) / x.sqrt()
    }
}

#[cfg(test)]
mod tests {
    use super::i0;

    #[test]
    fn known_values() {
        // reference values from scipy.special.i0
        assert!((i0(0.0) - 1.0).abs() < 1e-15);
        assert!((i0(1.0) - 1.2660658777520084).abs() < 1e-13);
        assert!((i0(4.0) - 11.301921952136330).abs() < 1e-11);
        assert!((i0(10.0) - 2815.7166284662544).abs() < 1e-8);
    }
}