feffit 0.1.0

//! `diffkk` vs `larch` (z=29, Cu K edge, scalar KK kernel).
//!
//! Reference generated by `scripts/ref_diffkk.py`. The Chantler `f1`/`f2` larch
//! computes internally (via `xraydb`) are dumped at full resolution and fed
//! identically into the Rust port, isolating the diffKK algorithm from the
//! atomic-data table source so the comparison is bit-exact.

use std::collections::HashMap;
use std::path::PathBuf;

use feffit::xasproc::{MbackParams, diffkk};

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

fn load_xmu() -> (Vec<f64>, Vec<f64>) {
    let text = std::fs::read_to_string(data("cu.xmu")).unwrap();
    let mut e = Vec::new();
    let mut m = Vec::new();
    for line in text.lines() {
        let s = line.trim();
        if s.is_empty() || s.starts_with('#') {
            continue;
        }
        let mut it = s.split_whitespace();
        e.push(it.next().unwrap().parse().unwrap());
        m.push(it.next().unwrap().parse().unwrap());
    }
    (e, m)
}

struct Ref {
    f1: Vec<f64>,
    f2: Vec<f64>,
    scalars: HashMap<String, f64>,
    arrays: HashMap<String, Vec<(usize, f64)>>,
}

fn load_ref(npts: usize) -> Ref {
    let text = std::fs::read_to_string(data("ref_diffkk.txt")).unwrap();
    let mut f1 = vec![0.0; npts];
    let mut f2 = vec![0.0; npts];
    let mut scalars = HashMap::new();
    let mut arrays: HashMap<String, Vec<(usize, f64)>> = HashMap::new();
    for line in text.lines() {
        let s = line.trim();
        if s.is_empty() || s.starts_with('#') {
            continue;
        }
        let parts: Vec<&str> = s.split_whitespace().collect();
        match parts.as_slice() {
            [key, val] => {
                if *key != "npts" && *key != "grid_npts" {
                    scalars.insert(key.to_string(), val.parse().unwrap());
                }
            }
            [key, i, val] => {
                let i: usize = i.parse().unwrap();
                let v: f64 = val.parse().unwrap();
                match *key {
                    "f1" => f1[i] = v,
                    "f2" => f2[i] = v,
                    _ => arrays.entry(key.to_string()).or_default().push((i, v)),
                }
            }
            _ => panic!("bad line: {s}"),
        }
    }
    Ref {
        f1,
        f2,
        scalars,
        arrays,
    }
}

fn check_arr(name: &str, got: &[f64], want: &[(usize, f64)], tol: f64) -> f64 {
    let mut worst = 0.0_f64;
    for &(i, w) in want {
        let g = got[i];
        let d = (g - w).abs();
        worst = worst.max(d);
        assert!(d < tol + tol * w.abs(), "{name}[{i}] {g} vs {w} (abs {d})");
    }
    worst
}

#[test]
fn diffkk_matches_larch() {
    let (energy, mu) = load_xmu();
    let r = load_ref(energy.len());

    let dkk = diffkk(&energy, &mu, &r.f1, &r.f2, &MbackParams::default());

    let d_e0 = (dkk.e0 - r.scalars["e0"]).abs();
    assert!(d_e0 < 1e-12, "e0 {} vs {}", dkk.e0, r.scalars["e0"]);

    let wfpp = check_arr("fpp", &dkk.fpp, &r.arrays["fpp"], 1e-10);
    let wfp = check_arr("fp", &dkk.fp, &r.arrays["fp"], 1e-10);
    eprintln!("diffkk worst abs delta: fpp={wfpp:.2e} fp={wfp:.2e}");
}