xraydb 0.1.2

X-ray reference data for the elements
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

/// Linear interpolation (equivalent to numpy.interp).
///
/// Interpolates values from `(xp, fp)` at points `x`.
/// Values outside the range are clamped to the boundary values.
pub fn interp(x: &[f64], xp: &[f64], fp: &[f64]) -> Vec<f64> {
    x.iter().map(|&xi| interp_one(xi, xp, fp)).collect()
}

/// Interpolate a single value.
pub fn interp_one(x: f64, xp: &[f64], fp: &[f64]) -> f64 {
    if x <= xp[0] {
        return fp[0];
    }
    if x >= xp[xp.len() - 1] {
        return fp[fp.len() - 1];
    }

    // Binary search for the bracket
    let idx = xp.partition_point(|&v| v < x);
    if idx == 0 {
        return fp[0];
    }

    // Check for exact match
    if (xp[idx] - x).abs() < f64::EPSILON * xp[idx].abs() {
        return fp[idx];
    }

    let lo = idx - 1;
    let t = (x - xp[lo]) / (xp[idx] - xp[lo]);
    fp[lo] + t * (fp[idx] - fp[lo])
}

/// Log-log linear interpolation.
///
/// Equivalent to `exp(interp(log(x), log(xp), log(fp)))`.
/// Used for Chantler f2 and mu data.
pub fn interp_loglog(x: &[f64], xp: &[f64], fp: &[f64]) -> Vec<f64> {
    let log_xp: Vec<f64> = xp.iter().map(|v| v.ln()).collect();
    let log_fp: Vec<f64> = fp.iter().map(|v| v.ln()).collect();
    x.iter()
        .map(|&xi| interp_one(xi.ln(), &log_xp, &log_fp).exp())
        .collect()
}

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

    #[test]
    fn test_interp_basic() {
        let xp = vec![0.0, 1.0, 2.0];
        let fp = vec![0.0, 10.0, 20.0];

        let result = interp(&[0.5, 1.5], &xp, &fp);
        assert!((result[0] - 5.0).abs() < 1e-10);
        assert!((result[1] - 15.0).abs() < 1e-10);
    }

    #[test]
    fn test_interp_clamping() {
        let xp = vec![1.0, 2.0, 3.0];
        let fp = vec![10.0, 20.0, 30.0];

        let result = interp(&[0.0, 4.0], &xp, &fp);
        assert!((result[0] - 10.0).abs() < 1e-10);
        assert!((result[1] - 30.0).abs() < 1e-10);
    }
}