refeff-io 0.2.0

//! FEFF `xmu.dat` spectrum text codec.
//!
//! FEFF writes `xmu.dat` as a comment-rich header followed by six numeric
//! columns: photon energy, edge-relative energy, photoelectron wave number,
//! normalized total absorption, normalized atomic background, and chi. The
//! full-spectrum reader in `FULLSPECTRUM/rdxmu.f90` uses the `xsedge+ 50`
//! header scalar to convert normalized `mu` and `mu0` to absolute cross
//! sections.

use std::fmt::Write as _;
use std::path::Path;

use ndarray::{Array1, ArrayView1};
use refeff_core::{
    FEFF_ALPHA_INV, FEFF_BOHR_ANGSTROM, FEFF_HARTREE_EV, FullSpectrumBackgroundSegmentInput,
    FullSpectrumFineStructureSegmentInput, FullSpectrumValenceInput,
    full_spectrum_valence_epsilon2,
};

use crate::error::{IoError, Result};
use crate::format::{FortranField, write_fortran_row};

const XMU_DAT_ROW_WIDTH: usize = 6;
const COMPACT_FIXED_PRECISION: i32 = 3;
const COLUMN_EQUALITY_TOLERANCE: f64 = 1.0e-12;

/// FEFF's compact `xmu.dat` row layout: `omega`, edge-relative energy, and
/// wave number as fixed-point columns, followed by three `E13.5` value
/// columns, all written back-to-back with no separator.
const XMU_COMPACT_OMEGA: FortranField = FortranField::F {
    width: 12,
    precision: 3,
};
const XMU_COMPACT_EDGE: FortranField = FortranField::F {
    width: 11,
    precision: 3,
};
const XMU_COMPACT_K: FortranField = FortranField::F {
    width: 8,
    precision: 3,
};
const XMU_COMPACT_VALUE: FortranField = FortranField::E {
    width: 13,
    precision: 5,
};

/// FEFF's FPRIME `xmu.dat` row layout: `omega` and edge-relative energy as
/// fixed-point columns, followed by four zero-scaled `E13.5` value columns.
const XMU_FPRIME_VALUE: FortranField = FortranField::ZeroScaledE {
    width: 13,
    precision: 5,
    exp_width: 2,
};

/// FEFF's wide `xmu.dat` row layout, used when the compact fixed-point
/// columns would lose precision.
const XMU_WIDE_OMEGA: FortranField = FortranField::F {
    width: 21,
    precision: 10,
};
const XMU_WIDE_EDGE: FortranField = FortranField::F {
    width: 20,
    precision: 10,
};
const XMU_WIDE_K: FortranField = FortranField::F {
    width: 20,
    precision: 10,
};
const XMU_WIDE_VALUE: FortranField = FortranField::E {
    width: 20,
    precision: 10,
};

/// Parsed FEFF `xmu.dat` contents.
#[derive(Debug, Clone, PartialEq)]
pub struct XmuDatData {
    /// Header lines before and around the numeric data table.
    pub header_lines: Vec<String>,
    /// `xsedge+ 50` normalization scalar when the header provides it.
    pub normalization: Option<f64>,
    /// Photon energy in eV.
    pub photon_energy_ev: Array1<f64>,
    /// Photoelectron energy relative to the edge in eV.
    pub relative_energy_ev: Array1<f64>,
    /// Photoelectron wave number in inverse Angstrom.
    pub wave_number: Array1<f64>,
    /// Normalized total absorption coefficient.
    pub mu: Array1<f64>,
    /// Normalized atomic background absorption coefficient.
    pub mu0: Array1<f64>,
    /// Fine structure, `chi = mu - mu0` in normalized units.
    pub chi: Array1<f64>,
}

/// Units carried by a FEFF FULLSPECTRUM view of `xmu.dat`.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum FullSpectrumXmuUnits {
    /// Normalized `mu` and `mu0`, matching `FULLSPECTRUM/rdxmunorm.f90`.
    Normalized,
    /// Absolute cross sections in square Angstrom, matching
    /// `FULLSPECTRUM/rdxmu.f90`.
    SquareAngstrom,
}

/// FEFF FULLSPECTRUM view of an `xmu.dat` spectrum.
#[derive(Debug, Clone, PartialEq)]
pub struct FullSpectrumXmuData {
    /// Photon energy from column one, in eV.
    pub photon_energy_ev: Array1<f64>,
    /// Edge-relative photoelectron energy from column two, in eV.
    pub relative_energy_ev: Array1<f64>,
    /// Photoelectron wave number from column three, in inverse Angstrom.
    pub wave_number_inverse_angstrom: Array1<f64>,
    /// Total absorption in [`Self::units`].
    pub mu: Array1<f64>,
    /// Atomic background absorption in [`Self::units`].
    pub mu0: Array1<f64>,
    /// Units of `mu` and `mu0`.
    pub units: FullSpectrumXmuUnits,
}

/// FEFF FULLSPECTRUM FPRIME background segment read from one `xmu.dat`.
#[derive(Debug, Clone, PartialEq)]
pub struct FullSpectrumBackgroundSegmentData {
    /// Photon energy from column one, in eV.
    pub photon_energy_ev: Array1<f64>,
    /// Real anomalous scattering factor from column four.
    pub f_prime: Array1<f64>,
    /// Imaginary anomalous scattering factor from column five.
    pub f_double_prime: Array1<f64>,
}

/// FEFF FULLSPECTRUM fine-structure segment read from one `xmu.dat`.
#[derive(Debug, Clone, PartialEq)]
pub struct FullSpectrumFineStructureSegmentData {
    /// Photon energy from column one, in eV.
    pub photon_energy_ev: Array1<f64>,
    /// Photoelectron wave number from column three, in inverse Angstrom.
    pub wave_number_inverse_angstrom: Array1<f64>,
    /// Real or imaginary scattering-factor component.
    pub scattering_factor: Array1<f64>,
    /// Atomic-background component matching [`Self::scattering_factor`].
    pub background: Array1<f64>,
}

impl XmuDatData {
    /// Number of spectrum rows.
    #[must_use]
    pub fn point_count(&self) -> usize {
        self.photon_energy_ev.len()
    }

    /// Total absorption converted to absolute units when normalization exists.
    #[must_use]
    pub fn absolute_mu(&self) -> Option<Array1<f64>> {
        self.normalization
            .map(|normalization| self.mu.mapv(|value| value * normalization))
    }

    /// Atomic background converted to absolute units when normalization exists.
    #[must_use]
    pub fn absolute_mu0(&self) -> Option<Array1<f64>> {
        self.normalization
            .map(|normalization| self.mu0.mapv(|value| value * normalization))
    }
}

impl FullSpectrumXmuData {
    /// Number of spectrum rows.
    #[must_use]
    pub fn point_count(&self) -> usize {
        self.photon_energy_ev.len()
    }
}

impl FullSpectrumBackgroundSegmentData {
    /// Number of spectrum rows.
    #[must_use]
    pub fn point_count(&self) -> usize {
        self.photon_energy_ev.len()
    }

    /// Borrow this segment as a core FULLSPECTRUM `rdbkg` input.
    #[must_use]
    pub fn as_core_input(&self) -> FullSpectrumBackgroundSegmentInput<'_> {
        FullSpectrumBackgroundSegmentInput {
            photon_energy_ev: self.photon_energy_ev.view(),
            f_prime: self.f_prime.view(),
            f_double_prime: self.f_double_prime.view(),
        }
    }
}

impl FullSpectrumFineStructureSegmentData {
    /// Number of spectrum rows.
    #[must_use]
    pub fn point_count(&self) -> usize {
        self.photon_energy_ev.len()
    }

    /// Borrow this segment as a core FULLSPECTRUM `rdst` input.
    #[must_use]
    pub fn as_core_input(&self) -> FullSpectrumFineStructureSegmentInput<'_> {
        FullSpectrumFineStructureSegmentInput {
            photon_energy_ev: self.photon_energy_ev.view(),
            wave_number_inverse_angstrom: self.wave_number_inverse_angstrom.view(),
            scattering_factor: self.scattering_factor.view(),
            background: self.background.view(),
        }
    }
}

/// Render FEFF-compatible `xmu.dat` text.
pub fn xmu_dat_string(data: &XmuDatData) -> Result<String> {
    validate_xmu_dat(data)?;

    let mut out = String::new();
    for line in &data.header_lines {
        writeln!(out, "{line}")?;
    }
    let format = xmu_render_format(data);
    for (((((omega, edge), k), mu), mu0), chi) in data
        .photon_energy_ev
        .iter()
        .zip(data.relative_energy_ev.iter())
        .zip(data.wave_number.iter())
        .zip(data.mu.iter())
        .zip(data.mu0.iter())
        .zip(data.chi.iter())
    {
        match format {
            XmuRenderFormat::Compact => {
                write_fortran_row(
                    &mut out,
                    "",
                    [
                        (XMU_COMPACT_OMEGA, *omega),
                        (XMU_COMPACT_EDGE, *edge),
                        (XMU_COMPACT_K, *k),
                        (XMU_COMPACT_VALUE, *mu),
                        (XMU_COMPACT_VALUE, *mu0),
                        (XMU_COMPACT_VALUE, *chi),
                    ],
                )?;
            }
            XmuRenderFormat::FPrime => {
                write_fortran_row(
                    &mut out,
                    "",
                    [
                        (XMU_COMPACT_OMEGA, *omega),
                        (XMU_COMPACT_EDGE, *edge),
                        (XMU_FPRIME_VALUE, *k),
                        (XMU_FPRIME_VALUE, *mu),
                        (XMU_FPRIME_VALUE, *mu0),
                        (XMU_FPRIME_VALUE, *chi),
                    ],
                )?;
            }
            XmuRenderFormat::Wide => {
                write_fortran_row(
                    &mut out,
                    "",
                    [
                        (XMU_WIDE_OMEGA, *omega),
                        (XMU_WIDE_EDGE, *edge),
                        (XMU_WIDE_K, *k),
                        (XMU_WIDE_VALUE, *mu),
                        (XMU_WIDE_VALUE, *mu0),
                        (XMU_WIDE_VALUE, *chi),
                    ],
                )?;
            }
        }
        out.push('\n');
    }
    Ok(out)
}

/// Parse FEFF `xmu.dat` text.
pub fn parse_xmu_dat(text: &str) -> Result<XmuDatData> {
    let mut header_lines = Vec::new();
    let mut normalization = None;
    let mut photon_energy_ev = Vec::new();
    let mut relative_energy_ev = Vec::new();
    let mut wave_number = Vec::new();
    let mut mu = Vec::new();
    let mut mu0 = Vec::new();
    let mut chi = Vec::new();

    for (index, raw) in text.lines().enumerate() {
        let line_number = index + 1;
        let line = raw.trim_end();
        let tokens = line.split_whitespace().collect::<Vec<_>>();
        if tokens.first().is_some_and(|token| is_numeric_token(token)) {
            if tokens.len() != XMU_DAT_ROW_WIDTH {
                return Err(IoError::XmuDatRowWidth {
                    line: line_number,
                    actual: tokens.len(),
                    expected: XMU_DAT_ROW_WIDTH,
                });
            }
            photon_energy_ev.push(parse_f64(line_number, "omega", tokens[0])?);
            relative_energy_ev.push(parse_f64(line_number, "edge-relative energy", tokens[1])?);
            wave_number.push(parse_f64(line_number, "wave number", tokens[2])?);
            mu.push(parse_f64(line_number, "mu", tokens[3])?);
            mu0.push(parse_f64(line_number, "mu0", tokens[4])?);
            chi.push(parse_f64(line_number, "chi", tokens[5])?);
        } else {
            if let Some(value) = parse_normalization(line, line_number)? {
                normalization = Some(value);
            }
            header_lines.push(raw.to_string());
        }
    }

    let data = XmuDatData {
        header_lines,
        normalization,
        photon_energy_ev: Array1::from_vec(photon_energy_ev),
        relative_energy_ev: Array1::from_vec(relative_energy_ev),
        wave_number: Array1::from_vec(wave_number),
        mu: Array1::from_vec(mu),
        mu0: Array1::from_vec(mu0),
        chi: Array1::from_vec(chi),
    };
    validate_xmu_dat(&data)?;
    Ok(data)
}

/// Write FEFF `xmu.dat` text to a file.
pub fn write_xmu_dat(path: impl AsRef<Path>, data: &XmuDatData) -> Result<()> {
    let path = path.as_ref();
    std::fs::write(path, xmu_dat_string(data)?).map_err(|source| IoError::io(path, source))
}

/// Read FEFF `xmu.dat` text from a file.
pub fn read_xmu_dat(path: impl AsRef<Path>) -> Result<XmuDatData> {
    let path = path.as_ref();
    let text = std::fs::read_to_string(path).map_err(|source| IoError::io(path, source))?;
    parse_xmu_dat(&text)
}

/// Convert parsed `xmu.dat` to FEFF `FULLSPECTRUM/rdxmu.f90` absolute units.
///
/// FEFF stores normalized `mu` and `mu0` in `xmu.dat` and multiplies them by
/// the `xsedge+ 50` header scalar before the full-spectrum assembly consumes
/// the cross sections.
pub fn fullspectrum_absolute_xmu_from_xmu_dat(data: &XmuDatData) -> Result<FullSpectrumXmuData> {
    validate_xmu_dat(data)?;
    let normalization = data
        .normalization
        .ok_or_else(|| invalid_xmu_dat("normalization", "missing xsedge normalization"))?;
    let mu = data.mu.mapv(|value| value * normalization);
    let mu0 = data.mu0.mapv(|value| value * normalization);
    validate_fullspectrum_xmu_values("mu", &mu)?;
    validate_fullspectrum_xmu_values("mu0", &mu0)?;
    Ok(fullspectrum_xmu_data(
        data,
        mu,
        mu0,
        FullSpectrumXmuUnits::SquareAngstrom,
    ))
}

/// Convert parsed `xmu.dat` to FEFF `FULLSPECTRUM/rdxmunorm.f90` columns.
///
/// This keeps `mu` and `mu0` in FEFF's normalized units while preserving the
/// photon energy, edge-relative energy, and wave-number columns.
pub fn fullspectrum_normalized_xmu_from_xmu_dat(data: &XmuDatData) -> Result<FullSpectrumXmuData> {
    validate_xmu_dat(data)?;
    Ok(fullspectrum_xmu_data(
        data,
        data.mu.clone(),
        data.mu0.clone(),
        FullSpectrumXmuUnits::Normalized,
    ))
}

/// Convert a FPRIME `xmu.dat` table to a FEFF `FULLSPECTRUM/rdbkg.f90` segment.
///
/// FEFF reads FPRIME `xmu.dat` through `rdxmunorm.f90` and consumes columns
/// four and five directly as `f'` and `f''`, without `xsedge` normalization.
pub fn fullspectrum_background_segment_from_fprime_xmu_dat(
    data: &XmuDatData,
) -> Result<FullSpectrumBackgroundSegmentData> {
    let normalized = fullspectrum_normalized_xmu_from_xmu_dat(data)?;
    Ok(FullSpectrumBackgroundSegmentData {
        photon_energy_ev: normalized.photon_energy_ev,
        f_prime: normalized.mu,
        f_double_prime: normalized.mu0,
    })
}

/// Convert parsed `xmu.dat` to a FEFF `FULLSPECTRUM/rdst.f90` real segment.
///
/// FEFF uses `rdxmunorm.f90` for the real `fms_re`/`path_re` branches, so
/// column four is consumed as `f'` and column five as its atomic background
/// without applying the `xsedge` cross-section normalization.
pub fn fullspectrum_real_fine_structure_segment_from_xmu_dat(
    data: &XmuDatData,
) -> Result<FullSpectrumFineStructureSegmentData> {
    let normalized = fullspectrum_normalized_xmu_from_xmu_dat(data)?;
    Ok(FullSpectrumFineStructureSegmentData {
        photon_energy_ev: normalized.photon_energy_ev,
        wave_number_inverse_angstrom: normalized.wave_number_inverse_angstrom,
        scattering_factor: normalized.mu,
        background: normalized.mu0,
    })
}

/// Convert parsed `xmu.dat` to a FEFF `FULLSPECTRUM/rdst.f90` imaginary segment.
///
/// FEFF uses `rdxmu.f90` for the imaginary `fms_im`/`path_im` branches and
/// converts the absolute square-Angstrom absorption columns to `f''` with
/// `mu * alpha_inv * omega_hartree * bohr^2`.
pub fn fullspectrum_imaginary_fine_structure_segment_from_xmu_dat(
    data: &XmuDatData,
) -> Result<FullSpectrumFineStructureSegmentData> {
    let absolute = fullspectrum_absolute_xmu_from_xmu_dat(data)?;
    let bohr_squared = FEFF_BOHR_ANGSTROM.powi(2);
    let scattering_factor = Array1::from_iter(
        absolute
            .photon_energy_ev
            .iter()
            .copied()
            .zip(absolute.mu.iter().copied())
            .map(|(energy_ev, mu)| {
                mu * FEFF_ALPHA_INV * energy_ev / FEFF_HARTREE_EV * bohr_squared
            }),
    );
    let background = Array1::from_iter(
        absolute
            .photon_energy_ev
            .iter()
            .copied()
            .zip(absolute.mu0.iter().copied())
            .map(|(energy_ev, mu0)| {
                mu0 * FEFF_ALPHA_INV * energy_ev / FEFF_HARTREE_EV * bohr_squared
            }),
    );
    validate_fullspectrum_xmu_values("fine_structure scattering_factor", &scattering_factor)?;
    validate_fullspectrum_xmu_values("fine_structure background", &background)?;
    Ok(FullSpectrumFineStructureSegmentData {
        photon_energy_ev: absolute.photon_energy_ev,
        wave_number_inverse_angstrom: absolute.wave_number_inverse_angstrom,
        scattering_factor,
        background,
    })
}

/// Compute FEFF `FULLSPECTRUM/rdval.f90` valence eps2 from parsed `xmu.dat`.
///
/// The parsed `xmu.dat` must include FEFF's `xsedge+50` normalization scalar so
/// normalized `mu` can be converted to the absolute square-Angstrom cross
/// section consumed by the FULLSPECTRUM valence projection.
pub fn valence_epsilon2_from_xmu_dat(
    number_density: f64,
    omega: ArrayView1<'_, f64>,
    data: &XmuDatData,
) -> Result<Array1<f64>> {
    let absolute_mu = data
        .absolute_mu()
        .ok_or_else(|| invalid_xmu_dat("normalization", "missing xsedge normalization"))?;
    full_spectrum_valence_epsilon2(FullSpectrumValenceInput {
        number_density,
        omega,
        source_energy_ev: data.photon_energy_ev.view(),
        source_absorption_angstrom2: absolute_mu.view(),
    })
    .map_err(|source| invalid_xmu_dat("valence_epsilon2", source.to_string()))
}

fn fullspectrum_xmu_data(
    data: &XmuDatData,
    mu: Array1<f64>,
    mu0: Array1<f64>,
    units: FullSpectrumXmuUnits,
) -> FullSpectrumXmuData {
    FullSpectrumXmuData {
        photon_energy_ev: data.photon_energy_ev.clone(),
        relative_energy_ev: data.relative_energy_ev.clone(),
        wave_number_inverse_angstrom: data.wave_number.clone(),
        mu,
        mu0,
        units,
    }
}

fn parse_normalization(line: &str, line_number: usize) -> Result<Option<f64>> {
    let lower = line.to_ascii_lowercase();
    if !(lower.contains("xsedge") && lower.contains("normalize")) {
        return Ok(None);
    }
    let Some(token) = line.split_whitespace().last() else {
        return Ok(None);
    };
    Ok(Some(parse_f64(line_number, "xsedge normalization", token)?))
}

pub(crate) fn validate_xmu_dat(data: &XmuDatData) -> Result<()> {
    let point_count = data.point_count();
    if point_count == 0 {
        return Err(invalid_xmu_dat(
            "rows",
            "at least one spectrum row is required",
        ));
    }
    validate_len(
        "relative_energy_ev",
        data.relative_energy_ev.len(),
        point_count,
    )?;
    validate_len("wave_number", data.wave_number.len(), point_count)?;
    validate_len("mu", data.mu.len(), point_count)?;
    validate_len("mu0", data.mu0.len(), point_count)?;
    validate_len("chi", data.chi.len(), point_count)?;

    if let Some(normalization) = data.normalization {
        validate_finite("xsedge normalization", normalization)?;
    }
    for (row, (((((omega, edge), k), mu), mu0), chi)) in data
        .photon_energy_ev
        .iter()
        .zip(data.relative_energy_ev.iter())
        .zip(data.wave_number.iter())
        .zip(data.mu.iter())
        .zip(data.mu0.iter())
        .zip(data.chi.iter())
        .enumerate()
    {
        let row = row + 1;
        validate_finite_row("omega", *omega, row)?;
        validate_finite_row("edge-relative energy", *edge, row)?;
        validate_finite_row("wave number", *k, row)?;
        validate_finite_row("mu", *mu, row)?;
        validate_finite_row("mu0", *mu0, row)?;
        validate_finite_row("chi", *chi, row)?;
    }
    Ok(())
}

fn validate_len(field: &'static str, actual: usize, expected: usize) -> Result<()> {
    if actual == expected {
        Ok(())
    } else {
        Err(IoError::XmuDatShape {
            field,
            actual,
            expected,
        })
    }
}

fn validate_fullspectrum_xmu_values(field: &'static str, values: &Array1<f64>) -> Result<()> {
    for (row, value) in values.iter().copied().enumerate() {
        validate_finite_row(field, value, row + 1)?;
    }
    Ok(())
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum XmuRenderFormat {
    Compact,
    FPrime,
    Wide,
}

fn xmu_render_format(data: &XmuDatData) -> XmuRenderFormat {
    if looks_like_fprime_xmu(data) || header_looks_like_fprime_xmu(data) {
        XmuRenderFormat::FPrime
    } else if needs_wide_xmu_format(data) {
        XmuRenderFormat::Wide
    } else {
        XmuRenderFormat::Compact
    }
}

fn looks_like_fprime_xmu(data: &XmuDatData) -> bool {
    data.wave_number
        .iter()
        .zip(data.mu.iter())
        .all(|(wave_number, mu)| (*wave_number - *mu).abs() <= COLUMN_EQUALITY_TOLERANCE)
        && data
            .mu0
            .iter()
            .zip(data.chi.iter())
            .all(|(mu0, chi)| (*mu0 - *chi).abs() <= COLUMN_EQUALITY_TOLERANCE)
}

fn header_looks_like_fprime_xmu(data: &XmuDatData) -> bool {
    data.header_lines.iter().any(|line| {
        let line = line.to_ascii_lowercase();
        line.contains("f'") && line.contains("f''")
    })
}

fn needs_wide_xmu_format(data: &XmuDatData) -> bool {
    data.photon_energy_ev
        .iter()
        .chain(data.relative_energy_ev.iter())
        .chain(data.wave_number.iter())
        .any(|value| has_more_decimal_precision(*value, COMPACT_FIXED_PRECISION))
}

fn has_more_decimal_precision(value: f64, precision: i32) -> bool {
    let scale = 10.0_f64.powi(precision);
    let rounded = (value * scale).round() / scale;
    (value - rounded).abs() > 1.0e-9
}

fn parse_f64(line: usize, field: &'static str, token: &str) -> Result<f64> {
    token
        .replace(['D', 'd'], "E")
        .parse::<f64>()
        .map_err(|_| IoError::XmuDatParse {
            field,
            line,
            token: token.to_string(),
        })
}

fn validate_finite(field: &'static str, value: f64) -> Result<()> {
    if value.is_finite() {
        Ok(())
    } else {
        Err(invalid_xmu_dat(field, "value must be finite"))
    }
}

fn validate_finite_row(field: &'static str, value: f64, row: usize) -> Result<()> {
    if value.is_finite() {
        Ok(())
    } else {
        Err(IoError::InvalidXmuDat {
            field,
            message: format!("row {row} value must be finite"),
        })
    }
}

fn invalid_xmu_dat(field: &'static str, message: impl Into<String>) -> IoError {
    IoError::InvalidXmuDat {
        field,
        message: message.into(),
    }
}

fn is_numeric_token(token: &str) -> bool {
    token.replace(['D', 'd'], "E").parse::<f64>().is_ok()
}

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

    #[test]
    fn parses_feff_xmu_reference_shape() -> Result<()> {
        let data = parse_xmu_dat(XMU_DAT)?;
        assert_eq!(data.point_count(), 3);
        assert_eq!(data.normalization, Some(1.2667e-4));
        assert_eq!(data.photon_energy_ev[0], 11076.317);
        assert_eq!(data.relative_energy_ev[1], -39.429);
        assert_eq!(data.wave_number[2], -2.965);
        assert_eq!(data.mu[0], 9.93209e-3);
        assert_eq!(data.mu0[1], 8.38540e-3);
        assert_eq!(data.chi[2], 3.54700e-4);
        let absolute = data
            .absolute_mu()
            .ok_or_else(|| invalid_xmu_dat("mu", "missing norm"))?;
        assert!((absolute[0] - 9.93209e-3 * 1.2667e-4).abs() < 1.0e-14);
        Ok(())
    }

    #[test]
    fn converts_xmu_to_feff_fullspectrum_rdxmu_views() -> Result<()> {
        let data = parse_xmu_dat(XMU_DAT)?;

        let absolute = fullspectrum_absolute_xmu_from_xmu_dat(&data)?;
        assert_eq!(absolute.point_count(), 3);
        assert_eq!(absolute.units, FullSpectrumXmuUnits::SquareAngstrom);
        assert_eq!(absolute.photon_energy_ev[0], data.photon_energy_ev[0]);
        assert_eq!(absolute.relative_energy_ev[1], data.relative_energy_ev[1]);
        assert_eq!(
            absolute.wave_number_inverse_angstrom[2],
            data.wave_number[2]
        );
        assert_close(absolute.mu[0], 9.93209e-3 * 1.2667e-4);
        assert_close(absolute.mu0[1], 8.38540e-3 * 1.2667e-4);

        let normalized = fullspectrum_normalized_xmu_from_xmu_dat(&data)?;
        assert_eq!(normalized.units, FullSpectrumXmuUnits::Normalized);
        assert_eq!(normalized.mu[0], data.mu[0]);
        assert_eq!(normalized.mu0[1], data.mu0[1]);
        Ok(())
    }

    #[test]
    fn converts_xmu_to_feff_fullspectrum_rdst_segments() -> Result<()> {
        let data = parse_xmu_dat(XMU_DAT)?;

        let real = fullspectrum_real_fine_structure_segment_from_xmu_dat(&data)?;
        assert_eq!(real.point_count(), 3);
        assert_eq!(real.photon_energy_ev[0], data.photon_energy_ev[0]);
        assert_eq!(real.wave_number_inverse_angstrom[2], data.wave_number[2]);
        assert_close(real.scattering_factor[0], data.mu[0]);
        assert_close(real.background[1], data.mu0[1]);
        let real_input = real.as_core_input();
        assert_eq!(real_input.scattering_factor[0], data.mu[0]);

        let imaginary = fullspectrum_imaginary_fine_structure_segment_from_xmu_dat(&data)?;
        let normalization = data
            .normalization
            .ok_or_else(|| invalid_xmu_dat("normalization", "missing norm"))?;
        let bohr_squared = FEFF_BOHR_ANGSTROM.powi(2);
        let expected_fpp = data.mu[0] * normalization * FEFF_ALPHA_INV * data.photon_energy_ev[0]
            / FEFF_HARTREE_EV
            * bohr_squared;
        let expected_background =
            data.mu0[1] * normalization * FEFF_ALPHA_INV * data.photon_energy_ev[1]
                / FEFF_HARTREE_EV
                * bohr_squared;
        assert_close(imaginary.scattering_factor[0], expected_fpp);
        assert_close(imaginary.background[1], expected_background);
        let imaginary_input = imaginary.as_core_input();
        assert_eq!(imaginary_input.background[1], imaginary.background[1]);
        Ok(())
    }

    #[test]
    fn converts_fprime_xmu_to_feff_fullspectrum_rdbkg_segments() -> Result<()> {
        let data = parse_xmu_dat(FPRIME_XMU_DAT)?;

        let segment = fullspectrum_background_segment_from_fprime_xmu_dat(&data)?;
        assert_eq!(segment.point_count(), 3);
        assert_eq!(segment.photon_energy_ev[1], data.photon_energy_ev[1]);
        assert_close(segment.f_prime[0], data.mu[0]);
        assert_close(segment.f_double_prime[2], data.mu0[2]);
        let core_input = segment.as_core_input();
        assert_eq!(core_input.f_prime[0], segment.f_prime[0]);
        assert_eq!(core_input.f_double_prime[2], segment.f_double_prime[2]);

        let rendered = xmu_dat_string(&data)?;
        assert_eq!(parse_xmu_dat(&rendered)?, data);
        Ok(())
    }

    #[test]
    fn renders_non_duplicate_fprime_columns_with_fprime_format() -> Result<()> {
        let mut data = parse_xmu_dat(FPRIME_XMU_DAT)?;
        data.wave_number[1] = -1.875;
        data.mu[1] = -1.625;
        data.mu0[1] = 0.625;
        data.chi[1] = 0.375;

        let rendered = xmu_dat_string(&data)?;
        let row = rendered
            .lines()
            .find(|line| line.contains("100.500"))
            .ok_or_else(|| invalid_xmu_dat("fprime", "missing rendered middle row"))?;
        let tokens = row.split_whitespace().collect::<Vec<_>>();

        for token in &tokens[2..=5] {
            assert!(token.contains('E'));
        }
        assert_close(tokens[2].parse::<f64>().unwrap(), -1.875);
        assert_close(tokens[3].parse::<f64>().unwrap(), -1.625);
        assert_close(tokens[4].parse::<f64>().unwrap(), 0.625);
        assert_close(tokens[5].parse::<f64>().unwrap(), 0.375);
        assert_eq!(parse_xmu_dat(&rendered)?, data);
        Ok(())
    }

    #[test]
    fn derives_valence_epsilon2_from_xmu_dat() -> Result<()> {
        let data = parse_xmu_dat(VALENCE_XMU_DAT)?;
        let omega = Array1::from_vec(vec![
            5.0 / 27.211_396,
            10.0 / 27.211_396,
            15.0 / 27.211_396,
            25.0 / 27.211_396,
            40.0 / 27.211_396,
        ]);

        let epsilon2 = valence_epsilon2_from_xmu_dat(0.075, omega.view(), &data)?;

        assert_eq!(epsilon2.len(), omega.len());
        assert_eq!(epsilon2[0], 0.0);
        assert_eq!(epsilon2[1], 0.0);
        assert!((epsilon2[2] - 131.219_281_455_964_96).abs() < 1.0e-12);
        assert!((epsilon2[3] - 157.463_137_747_157_93).abs() < 1.0e-12);
        assert_eq!(epsilon2[4], 0.0);
        Ok(())
    }

    #[test]
    fn roundtrips_xmu_text() -> Result<()> {
        let data = parse_xmu_dat(XMU_DAT)?;
        let rendered = xmu_dat_string(&data)?;
        assert_eq!(rendered, XMU_DAT);
        assert_eq!(parse_xmu_dat(&rendered)?, data);
        Ok(())
    }

    #[test]
    fn rejects_bad_xmu_inputs() -> Result<()> {
        assert!(parse_xmu_dat("# no data\n").is_err());
        assert!(parse_xmu_dat("1 2 3\n").is_err());
        assert!(parse_xmu_dat("1 2 3 NaN 5 6\n").is_err());
        assert!(parse_xmu_dat("# xsedge+ 50, used to normalize mu nope\n1 2 3 4 5 6\n").is_err());
        let missing_normalization = parse_xmu_dat("# omega e k mu mu0 chi\n1 2 3 4 5 6\n")?;
        assert!(fullspectrum_absolute_xmu_from_xmu_dat(&missing_normalization).is_err());
        assert!(
            fullspectrum_imaginary_fine_structure_segment_from_xmu_dat(&missing_normalization)
                .is_err()
        );
        assert!(
            fullspectrum_real_fine_structure_segment_from_xmu_dat(&missing_normalization).is_ok()
        );
        assert!(
            fullspectrum_background_segment_from_fprime_xmu_dat(&missing_normalization).is_ok()
        );
        Ok(())
    }

    fn assert_close(actual: f64, expected: f64) {
        assert!(
            (actual - expected).abs() <= 1.0e-14,
            "actual {actual} expected {expected}"
        );
    }

    const XMU_DAT: &str = r#"# # Cu                                                           FEFF 10.0.0
#  S02=1.000  Temp=   0.00  Debye_temp=   0.00  Global_sig2= 0.00000
#     0/   0 paths used
#  xsedge+ 50, used to normalize mu           1.2667E-04
#  -----------------------------------------------------------------------
#  omega    e    k    mu    mu0     chi     @#
   11076.317    -40.000  -3.016  9.93209E-03  9.60242E-03  3.29662E-04
   11076.888    -39.429  -2.991  8.72601E-03  8.38540E-03  3.40613E-04
   11077.459    -38.858  -2.965  7.66539E-03  7.31069E-03  3.54700E-04
"#;

    const VALENCE_XMU_DAT: &str = r#"# FEFF valence xmu.dat
#     0/   0 paths used
#  xsedge+ 50, used to normalize mu           2.0000E+00
#  omega    e    k    mu    mu0     chi     @#
      10.000      0.000   0.000  5.00000E-01  1.00000E-01  0.00000E+00
      20.000     10.000   0.000  1.50000E+00  2.00000E-01  0.00000E+00
      40.000     30.000   0.000  3.50000E+00  3.00000E-01  0.00000E+00
"#;

    const FPRIME_XMU_DAT: &str = r#"# FEFF FPRIME xmu.dat
#  omega    e    f'    f'    f''    f''     @#
       0.000      0.000 -2.50000E+00 -2.50000E+00  0.00000E+00  0.00000E+00
     100.500    100.500 -1.75000E+00 -1.75000E+00  5.00000E-01  5.00000E-01
     250.000    250.000 -5.00000E-01 -5.00000E-01  1.25000E+00  1.25000E+00
"#;

    /// Round-trip property coverage (F7): generators snap floats to the
    /// exact decimals the Compact `xmu.dat` layout renders (`F12.3`/`F11.3`
    /// fixed columns, `E13.5` value columns) and steer clear of the
    /// `looks_like_fprime_xmu` column-equality heuristic, so
    /// `parse_xmu_dat(xmu_dat_string(data)) == data` holds for the Compact
    /// render path without expected precision loss from an unsnapped `f64`
    /// or an accidental switch to the FPrime/Wide layouts (already covered
    /// by the example-based tests above).
    mod proptests {
        use super::*;
        use crate::format::fortran_exp;
        use proptest::prelude::*;

        fn snap_fixed(value: f64, precision: usize) -> f64 {
            format!("{value:.precision$}")
                .parse::<f64>()
                .unwrap_or(value)
        }

        fn snap_exp(value: f64) -> f64 {
            fortran_exp(value, 13, 5)
                .trim()
                .parse::<f64>()
                .unwrap_or(value)
        }

        fn omega_strategy() -> impl Strategy<Value = f64> {
            (-9_999_999_i64..9_999_999).prop_map(|n| snap_fixed(n as f64 / 1000.0, 3))
        }

        fn edge_strategy() -> impl Strategy<Value = f64> {
            (-999_999_i64..999_999).prop_map(|n| snap_fixed(n as f64 / 1000.0, 3))
        }

        fn k_strategy() -> impl Strategy<Value = f64> {
            (-99_999_i64..99_999).prop_map(|n| snap_fixed(n as f64 / 1000.0, 3))
        }

        fn value_strategy() -> impl Strategy<Value = f64> {
            (-999_999_i64..999_999).prop_map(|n| snap_exp(n as f64 / 100.0))
        }

        fn compact_row_strategy() -> impl Strategy<Value = (f64, f64, f64, f64, f64, f64)> {
            (
                omega_strategy(),
                edge_strategy(),
                k_strategy(),
                value_strategy(),
                value_strategy(),
                value_strategy(),
            )
                .prop_map(|(omega, edge, k, mu, mu0, chi)| {
                    // Nudge away exact wave_number/mu or mu0/chi coincidences
                    // so the render format detector picks Compact rather than
                    // FPrime, matching what this property targets.
                    let k = if (k - mu).abs() <= COLUMN_EQUALITY_TOLERANCE {
                        snap_fixed(k + 1.0, 3)
                    } else {
                        k
                    };
                    let mu0 = if (mu0 - chi).abs() <= COLUMN_EQUALITY_TOLERANCE {
                        snap_exp(mu0 + 1.0)
                    } else {
                        mu0
                    };
                    (omega, edge, k, mu, mu0, chi)
                })
        }

        proptest! {
            #[test]
            fn roundtrips_compact_rows(
                rows in prop::collection::vec(compact_row_strategy(), 1..6),
            ) {
                let data = XmuDatData {
                    header_lines: vec!["# proptest compact header".to_string()],
                    normalization: None,
                    photon_energy_ev: Array1::from_iter(rows.iter().map(|row| row.0)),
                    relative_energy_ev: Array1::from_iter(rows.iter().map(|row| row.1)),
                    wave_number: Array1::from_iter(rows.iter().map(|row| row.2)),
                    mu: Array1::from_iter(rows.iter().map(|row| row.3)),
                    mu0: Array1::from_iter(rows.iter().map(|row| row.4)),
                    chi: Array1::from_iter(rows.iter().map(|row| row.5)),
                };
                prop_assume!(!looks_like_fprime_xmu(&data));
                let rendered = xmu_dat_string(&data)?;
                let reparsed = parse_xmu_dat(&rendered)?;
                prop_assert_eq!(reparsed, data);
            }
        }
    }
}

/// Checked, unit-bearing access to one absorption-spectrum row.
#[derive(Debug, Clone, Copy)]
pub struct XmuPoint {
    pub photon_energy: refeff_core::units::ElectronVolts,
    pub relative_energy: refeff_core::units::ElectronVolts,
    pub wave_number: refeff_core::units::InverseAngstroms,
    pub mu: f64,
    pub mu0: f64,
    pub chi: f64,
}
impl XmuDatData {
    /// None for out-of-range, inconsistent-length or non-finite rows.
    pub fn point(&self, index: usize) -> Option<XmuPoint> {
        use refeff_core::units::{ElectronVolts, InverseAngstroms};
        let point = XmuPoint {
            photon_energy: ElectronVolts::new(*self.photon_energy_ev.get(index)?).ok()?,
            relative_energy: ElectronVolts::new(*self.relative_energy_ev.get(index)?).ok()?,
            wave_number: InverseAngstroms::new(*self.wave_number.get(index)?).ok()?,
            mu: *self.mu.get(index)?,
            mu0: *self.mu0.get(index)?,
            chi: *self.chi.get(index)?,
        };
        [point.mu, point.mu0, point.chi]
            .iter()
            .all(|value| value.is_finite())
            .then_some(point)
    }
}