xraydb 0.1.1

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
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
108
109
110
111
112
113
114
115
116
117
118
119
120
121

use crate::chemparser::chemparse;
use crate::constants::{AVOGADRO, R_ELECTRON_CM};
use crate::db::XrayDb;
use crate::elam::CrossSectionKind;
use crate::error::{Result, XrayDbError};

impl XrayDb {
    /// Returns mass attenuation coefficient for a material in 1/cm.
    ///
    /// # Arguments
    /// * `formula` - Chemical formula (e.g., "H2O", "SiO2")
    /// * `density` - Density in g/cm³
    /// * `energies` - X-ray energies in eV
    /// * `kind` - Cross-section kind (total, photo, coherent, incoherent)
    pub fn material_mu(
        &self,
        formula: &str,
        density: f64,
        energies: &[f64],
        kind: CrossSectionKind,
    ) -> Result<Vec<f64>> {
        let composition = chemparse(formula)?;
        let mut species = Vec::with_capacity(composition.len());
        for (sym, &count) in &composition {
            let molar_mass = self.molar_mass(sym)?;
            species.push((sym.as_str(), count, molar_mass));
        }

        // Calculate total formula weight
        let total_weight: f64 = species
            .iter()
            .map(|(_, count, molar_mass)| count * molar_mass)
            .sum();

        if total_weight <= 0.0 {
            return Err(XrayDbError::InvalidFormula(format!(
                "zero weight formula: {formula}"
            )));
        }

        // Sum weighted cross-sections
        let mut mu = vec![0.0_f64; energies.len()];
        for &(sym, count, molar_mass) in &species {
            let frac = count * molar_mass / total_weight;
            let elem_mu = self.mu_elam(sym, energies, kind)?;
            for (i, &val) in elem_mu.iter().enumerate() {
                mu[i] += frac * val;
            }
        }

        // Convert from cm²/g to 1/cm by multiplying by density
        for val in &mut mu {
            *val *= density;
        }

        Ok(mu)
    }

    /// Returns X-ray refractive index components (delta, beta, attenuation_length_cm).
    ///
    /// The complex index of refraction is: n = 1 - delta - i*beta
    ///
    /// # Arguments
    /// * `formula` - Chemical formula
    /// * `density` - Density in g/cm³
    /// * `energy` - X-ray energy in eV
    pub fn xray_delta_beta(
        &self,
        formula: &str,
        density: f64,
        energy: f64,
    ) -> Result<(f64, f64, f64)> {
        let composition = chemparse(formula)?;
        let mut species = Vec::with_capacity(composition.len());
        for (sym, &count) in &composition {
            let molar_mass = self.molar_mass(sym)?;
            species.push((sym.as_str(), count, molar_mass));
        }

        // wavelength in cm
        let wavelength = 1.0e-7 * crate::constants::PLANCK_HC / energy;

        let total_weight: f64 = species
            .iter()
            .map(|(_, count, molar_mass)| count * molar_mass)
            .sum();

        if total_weight <= 0.0 {
            return Err(XrayDbError::InvalidFormula(format!(
                "zero weight formula: {formula}"
            )));
        }

        // Sum f1 and f2 contributions weighted by composition
        let mut sum_f1 = 0.0_f64;
        let mut sum_f2 = 0.0_f64;

        for &(sym, count, _) in &species {
            let z = self.atomic_number(sym)? as f64;
            let f1_vals = self.f1_chantler(sym, &[energy])?;
            let f2_vals = self.f2_chantler(sym, &[energy])?;
            // f1 from Chantler is f' (anomalous correction), so full f1 = Z + f'
            sum_f1 += count * (z + f1_vals[0]);
            sum_f2 += count * f2_vals[0];
        }

        // Number density factor
        let prefactor = R_ELECTRON_CM * wavelength * wavelength * density * AVOGADRO
            / (2.0 * std::f64::consts::PI * total_weight);

        let delta = prefactor * sum_f1;
        let beta = prefactor * sum_f2;
        let atlen = if beta > 0.0 {
            wavelength / (4.0 * std::f64::consts::PI * beta)
        } else {
            f64::INFINITY
        };

        Ok((delta, beta, atlen))
    }
}