xraydb 0.4.0

X-ray reference data for the elements: absorption edges, emission lines, and Elam/Chantler cross-sections
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

//! Internal cubic-spline evaluation for the Elam tables.

/// Evaluate a cubic spline at one point using pre-computed second derivatives.
///
/// This is the core interpolation behind every Elam photoabsorption and scattering
/// cross-section lookup; it is a direct port of upstream XrayDB's `elam_spline`.
///
/// * `xin` — knot abscissae, strictly increasing
/// * `yin` — knot ordinates
/// * `yspl` — pre-computed second derivatives at the knots
///
/// Returns `NaN` for a degenerate (empty or single-knot) table rather than panicking.
pub(crate) fn elam_spline_one(xin: &[f64], yin: &[f64], yspl: &[f64], x: f64) -> f64 {
    let n = xin.len().min(yin.len()).min(yspl.len());
    if n == 0 {
        return f64::NAN;
    }
    if n == 1 {
        return yin[0];
    }

    // Bracket [lo, hi] with 1 <= hi <= n-1.
    let hi = xin[..n].partition_point(|&v| v < x).clamp(1, n - 1);
    let lo = hi - 1;

    let diff = xin[hi] - xin[lo];
    if diff.is_nan() || diff <= 0.0 {
        // Non-increasing knots: fall back to the nearest ordinate instead of
        // dividing by zero.
        return yin[hi];
    }

    let a = (xin[hi] - x) / diff;
    let b = (x - xin[lo]) / diff;

    a * yin[lo]
        + b * yin[hi]
        + (diff * diff / 6.0) * ((a * a - 1.0) * a * yspl[lo] + (b * b - 1.0) * b * yspl[hi])
}

/// Evaluate a cubic spline at each point of `xout`.
#[cfg(test)]
pub(crate) fn elam_spline(xin: &[f64], yin: &[f64], yspl: &[f64], xout: &[f64]) -> Vec<f64> {
    xout.iter()
        .map(|&x| elam_spline_one(xin, yin, yspl, x))
        .collect()
}

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

    #[test]
    fn test_spline_at_knot_points() {
        let x = [1.0, 2.0, 3.0, 4.0, 5.0];
        let y = [1.0, 4.0, 9.0, 16.0, 25.0];
        let spl = [0.0; 5]; // zero second derivatives = linear

        for (&xi, &yi) in x.iter().zip(y.iter()) {
            let result = elam_spline(&x, &y, &spl, &[xi]);
            assert!(
                (result[0] - yi).abs() < 1e-10,
                "at x={xi}: got {} expected {yi}",
                result[0]
            );
        }
    }

    #[test]
    fn test_spline_linear_interpolation() {
        // With zero spline coefficients, should do linear interpolation
        let x = [0.0, 1.0, 2.0];
        let y = [0.0, 1.0, 2.0];
        let spl = [0.0; 3];

        let result = elam_spline(&x, &y, &spl, &[0.5, 1.5]);
        assert!((result[0] - 0.5).abs() < 1e-10);
        assert!((result[1] - 1.5).abs() < 1e-10);
    }

    #[test]
    fn empty_table_yields_nan_not_panic() {
        let out = elam_spline(&[], &[], &[], &[1.0]);
        assert!(out[0].is_nan());
    }

    #[test]
    fn single_knot_returns_that_knot() {
        assert_eq!(elam_spline_one(&[2.0], &[9.0], &[0.0], 100.0), 9.0);
    }

    #[test]
    fn mismatched_lengths_do_not_read_out_of_bounds() {
        // yspl shorter than xin: n collapses to 2, no panic.
        let v = elam_spline_one(&[0.0, 1.0, 2.0], &[0.0, 1.0, 2.0], &[0.0, 0.0], 0.5);
        assert!((v - 0.5).abs() < 1e-12, "got {v}");
    }

    #[test]
    fn extrapolates_beyond_the_table_without_panicking() {
        let x = [0.0, 1.0, 2.0];
        let y = [0.0, 1.0, 2.0];
        let spl = [0.0; 3];
        assert!(elam_spline_one(&x, &y, &spl, -10.0).is_finite());
        assert!(elam_spline_one(&x, &y, &spl, 10.0).is_finite());
    }
}