rexafs 0.2.4

Rust-powered X-ray absorption spectroscopy analysis and EXAFS fitting
Documentation 
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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124

const A: &[f64] = &[
    -4.415_341_646_479_339_5E-18,
    3.330_794_518_822_238_4E-17,
    -2.431_279_846_547_955E-16,
    1.715_391_285_555_133E-15,
    -1.168_533_287_799_345_1E-14,
    7.676_185_498_604_936E-14,
    -4.856_446_783_111_929E-13,
    2.955_052_663_129_64E-12,
    -1.726_826_291_441_556E-11,
    9.675_809_035_373_237E-11,
    -5.189_795_601_635_263E-10,
    2.659_823_724_682_386_6E-9,
    -1.300_025_009_986_248E-8,
    6.046_995_022_541_919E-8,
    -2.670_793_853_940_612E-7,
    1.117_387_539_120_103_7E-6,
    -4.416_738_358_458_750_5E-6,
    1.644_844_807_072_889_6E-5,
    -5.754_195_010_082_104E-5,
    1.885_028_850_958_416_5E-4,
    -5.763_755_745_385_824E-4,
    1.639_475_616_941_335_7E-3,
    -4.324_309_995_050_576E-3,
    1.054_646_039_459_499_8E-2,
    -2.373_741_480_589_947E-2,
    4.930_528_423_967_071E-2,
    -9.490_109_704_804_764E-2,
    1.716_209_015_222_087_7E-1,
    -3.046_826_723_431_984E-1,
    6.767_952_744_094_761E-1,
];

const B: &[f64] = &[
    -7.233_180_487_874_754E-18,
    -4.830_504_485_944_182E-18,
    4.465_621_420_296_76E-17,
    3.461_222_867_697_461E-17,
    -2.827_623_980_516_583_6E-16,
    -3.425_485_619_677_219E-16,
    1.772_560_133_056_526_3E-15,
    3.811_680_669_352_622_4E-15,
    -9.554_846_698_828_307E-15,
    -4.150_569_347_287_222E-14,
    1.540_086_217_521_41E-14,
    3.852_778_382_742_142_6E-13,
    7.180_124_451_383_666E-13,
    -1.794_178_531_506_806_2E-12,
    -1.321_581_184_044_771_3E-11,
    -3.149_916_527_963_241_6E-11,
    1.188_914_710_784_643_9E-11,
    4.940_602_388_224_97E-10,
    3.396_232_025_708_386_5E-9,
    2.266_668_990_498_178E-8,
    2.048_918_589_469_063_8E-7,
    2.891_370_520_834_756_7E-6,
    6.889_758_346_916_825E-5,
    3.369_116_478_255_694_3E-3,
    8.044_904_110_141_088E-1,
];

pub fn chvevl(x: f64, array: &[f64]) -> f64 {
    let mut b0: f64 = array[0];
    let mut b1: f64 = 0.0;
    let mut b2: f64 = 0.0;

    for p0 in &array[1..] {
        b2 = b1;
        b1 = b0;
        b0 = x * b1 - b2 + *p0;
    }

    0.5 * (b0 - b2)
}

pub fn bessel_i0(mut x: f64) -> f64 {
    let y: f64;

    if x < 0.0 {
        x = -x;
    }

    if x <= 8.0 {
        y = (&x / 2.0) - 2.0;
        return x.exp() * chvevl(y, A);
    }

    x.exp() * chvevl(32.0 / x - 2.0, B) / x.sqrt()
}

#[cfg(test)]
mod tests {
    use itertools::Itertools;
    use std::vec;

    use super::*;
    use crate::xafs::tests::TEST_TOL;
    use approx::assert_abs_diff_eq;

    #[test]
    fn test_bessel_i0() {
        let x: Vec<f64> = (0..=10).map(|x| x as f64).collect();
        let y = x.iter().map(|x| bessel_i0(*x)).collect_vec();

        // Expected values from scipy.special.i0
        let expected = [
            1.0,
            1.2660658777520082,
            2.279585302336067,
            4.880792585865024,
            11.30192195213633,
            27.239871823604442,
            67.23440697647797,
            168.59390851028968,
            427.56411572180474,
            1093.5883545113745,
            2815.716628466254,
        ];

        y.iter()
            .zip(expected.iter())
            .for_each(|(y, expected)| assert_abs_diff_eq!(y, expected, epsilon = TEST_TOL));
    }
}