rexafs 0.2.14

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

//! Shared inverse-FFT validation and real-spectrum resizing for both array backends.
//!
//! The input stores only the positive-frequency half of a real forward signal.
//! Conjugate symmetry supplies the negative half. The underlying inverse FFT is
//! unnormalized, so `sqrt(pi)/(kstep*nfft)` cancels both its length factor and
//! the public forward multiplier `kstep/sqrt(pi)`. This gives a real result;
//! it is not a complex analytic-signal inverse.
//! A changed nfft preserves available positive bins, pads with zeros or drops
//! high bins, and forces an even-length Nyquist coefficient to be real.
use super::{errors::FFTError, xrayfft::XrayFFTR};
use easyfft::{dyn_size::realfft::DynRealDft, num_complex::Complex, prelude::DynRealIfft};

/// Check the zero-origin uniform R grid and inverse scalar settings.
/// R is in Å. The absolute grid tolerance is `1e-8 * max(abs(R_step), 1)` Å.
/// Explicit kstep must imply the same R spacing through pi/(nfft*kstep).
/// This does not check the separate complex coefficients for finite values.
pub(super) fn validate(r: &[f64], settings: &XrayFFTR) -> Result<(), FFTError> {
    let invalid = |parameter: &str, reason: &str| FFTError::InvalidParameter {
        parameter: parameter.into(),
        reason: reason.into(),
    };
    if r.len() < 2 || r.iter().any(|v| !v.is_finite()) {
        return Err(invalid(
            "r",
            "requires at least two finite, uniformly spaced values starting at zero",
        ));
    }
    let step = r[1] - r[0];
    let tolerance = 1e-8 * step.abs().max(1.0);
    if step <= 0.0
        || r[0].abs() > tolerance
        || r.windows(2)
            .any(|pair| ((pair[1] - pair[0]) - step).abs() > tolerance)
    {
        return Err(invalid(
            "r",
            "must be a uniformly increasing grid starting at zero",
        ));
    }
    if settings.nfft.is_some_and(|n| n < 2) {
        return Err(invalid("nfft", "must be at least 2"));
    }
    for (name, value, positive) in [
        ("kstep", settings.kstep, true),
        ("qmax_out", settings.qmax_out, false),
        ("rweight", settings.rweight, false),
        ("dr", settings.dr, false),
        ("dr2", settings.dr2, false),
        ("rmin", settings.rmin, false),
        ("rmax", settings.rmax, false),
    ] {
        if value.is_some_and(|v| !v.is_finite() || v < 0.0 || (positive && v == 0.0)) {
            return Err(invalid(
                name,
                "must be finite and nonnegative (kstep must be positive)",
            ));
        }
    }
    if settings
        .rmin
        .zip(settings.rmax)
        .is_some_and(|(lo, hi)| lo >= hi)
    {
        return Err(invalid("rmin/rmax", "requires rmin < rmax"));
    }
    if let Some(kstep) = settings.kstep {
        let implied = std::f64::consts::PI / (kstep * settings.nfft.unwrap_or(2048) as f64);
        if (implied - step).abs() > tolerance {
            return Err(invalid(
                "kstep/nfft",
                "must match the input R-grid spacing; leave kstep on Auto when changing nfft",
            ));
        }
    }
    Ok(())
}

/// Allocate q=j*step (Å⁻¹), ending at the largest available index not above maximum.
/// The caller supplies positive finite step and at least one available sample;
/// this helper does not validate those preconditions or stretch the final point.
pub(super) fn q_grid(maximum: f64, step: f64, available: usize) -> Vec<f64> {
    let length = (maximum / step)
        .floor()
        .min(available.saturating_sub(1) as f64) as usize
        + 1;
    (0..length).map(|i| i as f64 * step).collect()
}

/// Rebuild a real DFT for the requested output size. Include DC and every
/// positive-frequency bin; the easyfft frequency-bin accessor excludes endpoints.
/// A new even-length transform's Nyquist bin must be real.
pub(super) fn inverse(chir: &DynRealDft<f64>, nfft: usize, kstep: f64) -> Vec<f64> {
    let mut bins = vec![Complex::new(0.0, 0.0); nfft / 2];
    for (destination, source) in bins.iter_mut().zip(chir.iter().skip(1)) {
        *destination = *source;
    }
    if let Some(nyquist) = bins.last_mut().filter(|_| nfft.is_multiple_of(2)) {
        nyquist.im = 0.0;
    }
    let resized = DynRealDft::new(*chir.get_offset(), &bins, nfft);
    let scale = std::f64::consts::PI.sqrt() / kstep / nfft as f64;
    resized
        .real_ifft()
        .iter()
        .map(|value| value * scale)
        .collect()
}