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
//! Parity for the Debye-Waller σ² models against larch.
//!
//! Reads `tests/data/ref_sigma2.txt` (generated by `scripts/ref_sigma2.py`):
//! `sigma2_eins` and `rmass` come from larch directly; `sigma2_debye` comes
//! from larch's pure-Python `sigma2_correldebye_py` port of Feff6 `sigms.f`
//! (the C library is x86_64-only and will not load on arm64). The Rust port
//! targets that pure-Python reference.

use std::path::PathBuf;

use feffit::feffdat::{FeffDatFile, sigma2_debye, sigma2_eins};

fn data_dir() -> PathBuf {
    PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("tests/data")
}

fn rel(got: f64, want: f64) -> f64 {
    (got - want).abs() / want.abs().max(1e-300)
}

#[test]
fn sigma2_models_match_larch() {
    let text = std::fs::read_to_string(data_dir().join("ref_sigma2.txt")).unwrap();

    let mut path: Option<FeffDatFile> = None;
    let mut n_eins = 0;
    let mut n_debye = 0;
    let mut n_rmass = 0;

    for line in text.lines() {
        if let Some(fname) = line.strip_prefix("#path ") {
            path = Some(FeffDatFile::from_path(data_dir().join(fname.trim())).unwrap());
        } else if let Some(rest) = line.strip_prefix("#rmass ") {
            let want: f64 = rest.trim().parse().unwrap();
            let got = path.as_ref().unwrap().rmass();
            println!(
                "rmass got={got:.12e} want={want:.12e} rel={:.2e}",
                rel(got, want)
            );
            assert!(rel(got, want) < 1e-12, "rmass: {got} vs {want}");
            n_rmass += 1;
        } else if let Some(rest) = line.strip_prefix("#eins ") {
            let mut it = rest.split_whitespace();
            let t: f64 = it.next().unwrap().parse().unwrap();
            let theta: f64 = it.next().unwrap().parse().unwrap();
            let want: f64 = it.next().unwrap().parse().unwrap();
            let got = sigma2_eins(t, theta, &path.as_ref().unwrap().geom);
            println!(
                "eins t={t} theta={theta} got={got:.12e} want={want:.12e} rel={:.2e}",
                rel(got, want)
            );
            assert!(
                rel(got, want) < 1e-12,
                "sigma2_eins(t={t},theta={theta}): {got} vs {want}"
            );
            n_eins += 1;
        } else if let Some(rest) = line.strip_prefix("#debye ") {
            let mut it = rest.split_whitespace();
            let t: f64 = it.next().unwrap().parse().unwrap();
            let theta: f64 = it.next().unwrap().parse().unwrap();
            let want: f64 = it.next().unwrap().parse().unwrap();
            let p = path.as_ref().unwrap();
            let got = sigma2_debye(t, theta, p.rnorman, &p.geom);
            println!(
                "debye t={t} theta={theta} got={got:.12e} want={want:.12e} rel={:.2e}",
                rel(got, want)
            );
            assert!(
                rel(got, want) < 1e-12,
                "sigma2_debye(t={t},theta={theta}): {got} vs {want}"
            );
            n_debye += 1;
        }
    }

    assert!(
        n_rmass >= 2 && n_eins >= 6 && n_debye >= 6,
        "expected the full reference grid"
    );
}