lindhard/electron/data/
table.rs

1use super::{
2    check_energy_grid, check_in_range, check_nonnegative, checked_text, invalid, lerp_on_grid,
3    parse_error, read_text, ElectronDataError, Result, Subshell, CACHE_FORMAT_VERSION,
4};
5use serde::de::IgnoredAny;
6use serde::{Deserialize, Serialize};
7use std::collections::BTreeMap;
8use std::path::Path;
9
10/// What a [`CrossSectionTable`] samples. Elastic and inelastic tables are
11/// distinct types of data and are never interchangeable.
12#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
13#[serde(rename_all = "snake_case")]
14pub enum SamplingAxis {
15    /// Elastic scattering: the polar deflection angle `θ` in radians,
16    /// `0 <= θ <= π`.
17    ElasticPolarAngle,
18    /// Inelastic scattering: the energy loss `W` in eV, `0 <= W <= E` for an
19    /// incident energy `E`.
20    InelasticEnergyLoss,
21}
22
23impl SamplingAxis {
24    /// Unit of the sampled variable, as written in the cache file.
25    pub fn unit(self) -> &'static str {
26        match self {
27            SamplingAxis::ElasticPolarAngle => "rad",
28            SamplingAxis::InelasticEnergyLoss => "eV",
29        }
30    }
31}
32
33/// The unvalidated contents of a [`CrossSectionTable`], for
34/// [`CrossSectionTable::new`]. See [`CrossSectionTable`] for the meaning and
35/// invariants of each field.
36#[derive(Debug, Clone, PartialEq)]
37pub struct CrossSectionTableParts {
38    /// Identity of the model that produced the table (name and version).
39    pub model: String,
40    /// Material identity.
41    pub material: String,
42    /// Provenance: the inputs of the model and their origin.
43    pub provenance: String,
44    /// The sampled variable.
45    pub axis: SamplingAxis,
46    /// Incident electron energies, eV.
47    pub energy_ev: Vec<f64>,
48    /// Inverse mean free path at each energy, m⁻¹.
49    pub inverse_mfp_per_m: Vec<f64>,
50    /// Cumulative-probability grid of the inverse CDFs.
51    pub probability: Vec<f64>,
52    /// Inverse-CDF ordinates, one row per energy (empty for zero-rate rows).
53    pub quantiles: Vec<Vec<f64>>,
54}
55
56/// A precomputed cross-section table for one material and one process,
57/// cacheable on disk (format version [`CACHE_FORMAT_VERSION`]).
58///
59/// # Contents and dimensions
60///
61/// With `N` incident energies and `M` probability points:
62///
63/// - `energy_ev[N]`: incident electron energies, eV; finite, positive,
64///   strictly increasing, `N >= 2`.
65/// - `inverse_mfp_per_m[N]`: inverse mean free path `λ⁻¹(E) = n σ(E)` for the
66///   process, in m⁻¹; finite and non-negative. Zero means no interaction at
67///   that energy (for example below a threshold).
68/// - `probability[M]`: the cumulative-probability axis shared by every row;
69///   finite, strictly increasing, starting at exactly 0 and ending at exactly
70///   1, `M >= 2`.
71/// - `quantiles[N][..]`: **inverse-CDF ordinates** (the cache stores the
72///   inverse CDF, not the forward CDF). `quantiles[i][j]` is the value `x` of
73///   the sampled variable at which the cumulative distribution at energy
74///   `energy_ev[i]` reaches `probability[j]`. A row is non-decreasing and lies
75///   in the axis's domain (`[0, π]` rad for the elastic polar angle, `[0,
76///   energy_ev[i]]` eV for the inelastic energy loss). A row whose inverse
77///   mean free path is zero must be **empty**: there is no distribution to
78///   store, and [`CrossSectionTable::inverse_cdf`] returns
79///   [`ElectronDataError::ZeroRate`] for it. Every other row has exactly `M`
80///   entries.
81///
82/// The file also records the model identity, the material, the provenance,
83/// the sampling axis and its unit (`ordinate_unit`, which must match the
84/// axis), and the format version.
85///
86/// # Sampling
87///
88/// [`CrossSectionTable::inverse_cdf`] evaluates a stored row at a
89/// probability `u` by linear interpolation between neighbouring probability
90/// points (exact at the stored points). Interpolation between incident
91/// energies is left to the transport code, which owns that choice.
92#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
93#[serde(try_from = "RawCacheIn", into = "RawCacheOut")]
94pub struct CrossSectionTable {
95    parts: CrossSectionTableParts,
96}
97
98/// Deserialization form. Every field is optional so that the format version
99/// is checked first, whatever else the file holds; unknown keys are collected
100/// and rejected after the version check.
101#[derive(Debug, Deserialize)]
102struct RawCacheIn {
103    format_version: Option<u32>,
104    model: Option<String>,
105    material: Option<String>,
106    provenance: Option<String>,
107    axis: Option<SamplingAxis>,
108    ordinate_unit: Option<String>,
109    energy_ev: Option<Vec<f64>>,
110    inverse_mfp_per_m: Option<Vec<f64>>,
111    probability: Option<Vec<f64>>,
112    quantiles: Option<Vec<Vec<f64>>>,
113    #[serde(flatten)]
114    unknown: BTreeMap<String, IgnoredAny>,
115}
116
117#[derive(Debug, Clone, Serialize)]
118struct RawCacheOut {
119    format_version: u32,
120    model: String,
121    material: String,
122    provenance: String,
123    axis: SamplingAxis,
124    ordinate_unit: String,
125    energy_ev: Vec<f64>,
126    inverse_mfp_per_m: Vec<f64>,
127    probability: Vec<f64>,
128    quantiles: Vec<Vec<f64>>,
129}
130
131fn required<T>(field: &'static str, v: Option<T>) -> Result<T> {
132    v.ok_or_else(|| ElectronDataError::Parse(format!("missing field `{field}`")))
133}
134
135impl TryFrom<RawCacheIn> for CrossSectionTable {
136    type Error = ElectronDataError;
137    fn try_from(r: RawCacheIn) -> Result<Self> {
138        if r.format_version != Some(CACHE_FORMAT_VERSION) {
139            return Err(ElectronDataError::UnsupportedVersion {
140                found: r.format_version,
141                supported: CACHE_FORMAT_VERSION,
142            });
143        }
144        if let Some(k) = r.unknown.keys().next() {
145            return Err(ElectronDataError::Parse(format!(
146                "unknown field `{k}` in a version {CACHE_FORMAT_VERSION} cache"
147            )));
148        }
149        let provenance = checked_text("provenance", r.provenance.as_deref())?;
150        let model = checked_text("model", r.model.as_deref())?;
151        let material = checked_text("material", r.material.as_deref())?;
152        let axis = required("axis", r.axis)?;
153        let unit = checked_text("ordinate_unit", r.ordinate_unit.as_deref())?;
154        if unit != axis.unit() {
155            return Err(ElectronDataError::MalformedText {
156                field: "ordinate_unit",
157                reason: format!(
158                    "{unit:?} does not match axis {axis:?} (unit {:?})",
159                    axis.unit()
160                ),
161            });
162        }
163        Self::new(CrossSectionTableParts {
164            model,
165            material,
166            provenance,
167            axis,
168            energy_ev: required("energy_ev", r.energy_ev)?,
169            inverse_mfp_per_m: required("inverse_mfp_per_m", r.inverse_mfp_per_m)?,
170            probability: required("probability", r.probability)?,
171            quantiles: required("quantiles", r.quantiles)?,
172        })
173    }
174}
175
176impl From<CrossSectionTable> for RawCacheOut {
177    fn from(t: CrossSectionTable) -> Self {
178        let p = t.parts;
179        Self {
180            format_version: CACHE_FORMAT_VERSION,
181            ordinate_unit: p.axis.unit().to_string(),
182            model: p.model,
183            material: p.material,
184            provenance: p.provenance,
185            axis: p.axis,
186            energy_ev: p.energy_ev,
187            inverse_mfp_per_m: p.inverse_mfp_per_m,
188            probability: p.probability,
189            quantiles: p.quantiles,
190        }
191    }
192}
193
194impl CrossSectionTable {
195    /// Build a validated table. See the type docs for the invariants.
196    pub fn new(parts: CrossSectionTableParts) -> Result<Self> {
197        let provenance = checked_text("provenance", Some(&parts.provenance))?;
198        let model = checked_text("model", Some(&parts.model))?;
199        let material = checked_text("material", Some(&parts.material))?;
200        let n = parts.energy_ev.len();
201        check_energy_grid("incident energy grid", &parts.energy_ev)?;
202        if parts.inverse_mfp_per_m.len() != n {
203            return Err(invalid(
204                "inverse mean free path",
205                format!(
206                    "need one value per energy: {n} energies, {} values",
207                    parts.inverse_mfp_per_m.len()
208                ),
209            ));
210        }
211        check_nonnegative("inverse mean free path", &parts.inverse_mfp_per_m)?;
212        let p = &parts.probability;
213        if p.len() < 2 {
214            return Err(invalid("probability axis", "needs at least 2 points"));
215        }
216        if p.iter().any(|x| !x.is_finite()) {
217            return Err(invalid("probability axis", "values must be finite"));
218        }
219        if p[0] != 0.0 || p[p.len() - 1] != 1.0 {
220            return Err(invalid(
221                "probability axis",
222                format!(
223                    "must start at 0 and end at 1, got {} .. {}",
224                    p[0],
225                    p[p.len() - 1]
226                ),
227            ));
228        }
229        if let Some(w) = p.windows(2).find(|w| w[1] <= w[0]) {
230            return Err(invalid(
231                "probability axis",
232                format!("must be strictly increasing, got {} then {}", w[0], w[1]),
233            ));
234        }
235        if parts.quantiles.len() != n {
236            return Err(invalid(
237                "inverse-CDF rows",
238                format!(
239                    "need one row per energy: {n} energies, {} rows",
240                    parts.quantiles.len()
241                ),
242            ));
243        }
244        for (i, row) in parts.quantiles.iter().enumerate() {
245            let e = parts.energy_ev[i];
246            if parts.inverse_mfp_per_m[i] == 0.0 {
247                if !row.is_empty() {
248                    return Err(invalid(
249                        "inverse-CDF rows",
250                        format!(
251                            "row {i} ({e} eV) has zero inverse mean free path and must be empty"
252                        ),
253                    ));
254                }
255                continue;
256            }
257            if row.len() != p.len() {
258                return Err(invalid(
259                    "inverse-CDF rows",
260                    format!(
261                        "row {i} ({e} eV) has {} entries; the probability axis has {}",
262                        row.len(),
263                        p.len()
264                    ),
265                ));
266            }
267            let max = match parts.axis {
268                SamplingAxis::ElasticPolarAngle => std::f64::consts::PI,
269                SamplingAxis::InelasticEnergyLoss => e,
270            };
271            if let Some(x) = row
272                .iter()
273                .find(|x| !(x.is_finite() && **x >= 0.0 && **x <= max))
274            {
275                return Err(invalid(
276                    "inverse-CDF rows",
277                    format!(
278                        "row {i} ({e} eV): {x} is outside [0, {max}] {}",
279                        parts.axis.unit()
280                    ),
281                ));
282            }
283            if let Some(w) = row.windows(2).find(|w| w[1] < w[0]) {
284                return Err(invalid(
285                    "inverse-CDF rows",
286                    format!(
287                        "row {i} ({e} eV) must be non-decreasing, got {} then {}",
288                        w[0], w[1]
289                    ),
290                ));
291            }
292        }
293        Ok(Self {
294            parts: CrossSectionTableParts {
295                model,
296                material,
297                provenance,
298                ..parts
299            },
300        })
301    }
302
303    /// Parse and validate the TOML cache form.
304    pub fn from_toml_str(text: &str) -> Result<Self> {
305        // Parse the raw form, then convert, so the caller gets the typed
306        // error (an unsupported version, a missing provenance) rather than
307        // the message toml wraps a conversion error in.
308        let raw: RawCacheIn = toml::from_str(text).map_err(parse_error)?;
309        Self::try_from(raw)
310    }
311
312    /// Serialize to the TOML cache form.
313    pub fn to_toml_string(&self) -> Result<String> {
314        toml::to_string(self).map_err(|e| ElectronDataError::Parse(e.to_string()))
315    }
316
317    /// Load a TOML cache file.
318    pub fn from_toml_file(path: impl AsRef<Path>) -> Result<Self> {
319        Self::from_toml_str(&read_text(path.as_ref())?)
320    }
321
322    /// Write a TOML cache file.
323    pub fn write_toml_file(&self, path: impl AsRef<Path>) -> Result<()> {
324        let path = path.as_ref();
325        std::fs::write(path, self.to_toml_string()?)
326            .map_err(|e| ElectronDataError::Io(format!("{}: {e}", path.display())))
327    }
328
329    /// Format version of this table ([`CACHE_FORMAT_VERSION`]).
330    pub fn format_version(&self) -> u32 {
331        CACHE_FORMAT_VERSION
332    }
333
334    /// Model identity.
335    pub fn model(&self) -> &str {
336        &self.parts.model
337    }
338
339    /// Material identity.
340    pub fn material(&self) -> &str {
341        &self.parts.material
342    }
343
344    /// Provenance string.
345    pub fn provenance(&self) -> &str {
346        &self.parts.provenance
347    }
348
349    /// The sampled variable.
350    pub fn axis(&self) -> SamplingAxis {
351        self.parts.axis
352    }
353
354    /// Incident energies, eV.
355    pub fn energy_ev(&self) -> &[f64] {
356        &self.parts.energy_ev
357    }
358
359    /// Inverse mean free paths, m⁻¹, one per energy.
360    pub fn inverse_mfp_per_m(&self) -> &[f64] {
361        &self.parts.inverse_mfp_per_m
362    }
363
364    /// The cumulative-probability axis.
365    pub fn probability(&self) -> &[f64] {
366        &self.parts.probability
367    }
368
369    /// The inverse-CDF ordinates at grid energy `i`, or `None` for a
370    /// zero-rate row (or an index past the end).
371    pub fn quantiles(&self, i: usize) -> Option<&[f64]> {
372        self.parts
373            .quantiles
374            .get(i)
375            .filter(|r| !r.is_empty())
376            .map(Vec::as_slice)
377    }
378
379    /// The inverse CDF of grid row `i` at cumulative probability `u` in
380    /// `[0, 1]`, in the axis's unit. A zero-rate row is
381    /// [`ElectronDataError::ZeroRate`]: the caller must not sample an
382    /// interaction that cannot happen.
383    pub fn inverse_cdf(&self, i: usize, u: f64) -> Result<f64> {
384        let n = self.parts.energy_ev.len();
385        if i >= n {
386            return Err(ElectronDataError::IndexOutOfRange { index: i, len: n });
387        }
388        check_in_range("cumulative probability", u, 0.0, 1.0)?;
389        let Some(row) = self.quantiles(i) else {
390            return Err(ElectronDataError::ZeroRate {
391                energy_ev: self.parts.energy_ev[i],
392            });
393        };
394        Ok(lerp_on_grid(&self.parts.probability, row, u))
395    }
396
397    /// The validated contents.
398    pub fn parts(&self) -> &CrossSectionTableParts {
399        &self.parts
400    }
401}
402
403/// The format of a [`ShellChannelTable`] file this build reads and writes.
404///
405/// It versions only the wrapper (shell identity, binding energy, binding
406/// provenance); the wrapped [`CrossSectionTable`] keeps its own
407/// [`CACHE_FORMAT_VERSION`], which this type does not change. Bump it whenever
408/// the meaning or layout of the wrapper changes.
409pub const SHELL_CHANNEL_FORMAT_VERSION: u32 = 1;
410
411/// The energy-loss table of one inner-shell ionisation channel, with the
412/// identity of its shell.
413///
414/// The wrapped [`CrossSectionTable`] holds the channel's inverse mean free
415/// path and the inverse CDF of its loss `ω`, built from that shell's own
416/// optical ELF (`crate::electron::inelastic::inner_shell`); the wrapper adds
417/// the shell (`Z`, subshell), its binding energy `B` and the provenance of
418/// `B` (for example [`SubshellBindingTable::provenance`](super::SubshellBindingTable::provenance)). Built by
419/// [`crate::electron::inelastic::table::build_shell_channel_tables`] and
420/// consumed by [`crate::electron::transport::Transport::with_inner_shells`].
421///
422/// # Invariants
423///
424/// - `Z >= 1`; `B` finite and positive; a non-blank one-line binding
425///   provenance.
426/// - The table samples [`SamplingAxis::InelasticEnergyLoss`].
427/// - **Every stored loss is at least `B`**: a shell contributes only above its
428///   edge (de Vera et al. 2022, eq. (2), the step `Θ(E - B)`), so the
429///   secondary energy `ω - B` is never negative. Since a loss is also at most
430///   the row energy, every row with a nonzero rate lies at an energy of at
431///   least `B`.
432///
433/// # File form
434///
435/// TOML: `format_version` ([`SHELL_CHANNEL_FORMAT_VERSION`]), `z`,
436/// `subshell` (the ENDF label, e.g. `"K"`, `"L3"`), `binding_energy_ev`,
437/// `binding_provenance`, and a `[table]` holding the [`CrossSectionTable`]
438/// cache form (with its own `format_version`). The wrapper version is checked
439/// first, unknown keys are rejected, and the same checks as
440/// [`ShellChannelTable::new`] run on read.
441#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
442#[serde(try_from = "RawShellChannelIn", into = "RawShellChannelOut")]
443pub struct ShellChannelTable {
444    z: u8,
445    subshell: Subshell,
446    binding_energy_ev: f64,
447    binding_provenance: String,
448    table: CrossSectionTable,
449}
450
451#[derive(Debug, Deserialize)]
452struct RawShellChannelIn {
453    format_version: Option<u32>,
454    z: Option<u8>,
455    subshell: Option<String>,
456    binding_energy_ev: Option<f64>,
457    binding_provenance: Option<String>,
458    table: Option<RawCacheIn>,
459    #[serde(flatten)]
460    unknown: BTreeMap<String, IgnoredAny>,
461}
462
463#[derive(Debug, Clone, Serialize)]
464struct RawShellChannelOut {
465    format_version: u32,
466    z: u8,
467    subshell: String,
468    binding_energy_ev: f64,
469    binding_provenance: String,
470    table: CrossSectionTable,
471}
472
473impl TryFrom<RawShellChannelIn> for ShellChannelTable {
474    type Error = ElectronDataError;
475    fn try_from(r: RawShellChannelIn) -> Result<Self> {
476        if r.format_version != Some(SHELL_CHANNEL_FORMAT_VERSION) {
477            return Err(ElectronDataError::UnsupportedVersion {
478                found: r.format_version,
479                supported: SHELL_CHANNEL_FORMAT_VERSION,
480            });
481        }
482        if let Some(k) = r.unknown.keys().next() {
483            return Err(ElectronDataError::Parse(format!(
484                "unknown field `{k}` in a version {SHELL_CHANNEL_FORMAT_VERSION} shell-channel table"
485            )));
486        }
487        let label = checked_text("subshell", r.subshell.as_deref())?;
488        let subshell =
489            Subshell::from_label(&label).ok_or_else(|| ElectronDataError::MalformedText {
490                field: "subshell",
491                reason: format!("{label:?} is not an ENDF subshell label"),
492            })?;
493        Self::new(
494            required("z", r.z)?,
495            subshell,
496            required("binding_energy_ev", r.binding_energy_ev)?,
497            r.binding_provenance.as_deref().unwrap_or_default(),
498            CrossSectionTable::try_from(required("table", r.table)?)?,
499        )
500    }
501}
502
503impl From<ShellChannelTable> for RawShellChannelOut {
504    fn from(t: ShellChannelTable) -> Self {
505        Self {
506            format_version: SHELL_CHANNEL_FORMAT_VERSION,
507            z: t.z,
508            subshell: t.subshell.label().to_string(),
509            binding_energy_ev: t.binding_energy_ev,
510            binding_provenance: t.binding_provenance,
511            table: t.table,
512        }
513    }
514}
515
516impl ShellChannelTable {
517    /// A validated channel table: shell `subshell` of element `z` with
518    /// binding energy `binding_energy_ev` (whose origin `binding_provenance`
519    /// states) and its loss `table`. See the type docs for the invariants.
520    pub fn new(
521        z: u8,
522        subshell: Subshell,
523        binding_energy_ev: f64,
524        binding_provenance: &str,
525        table: CrossSectionTable,
526    ) -> Result<Self> {
527        let name = format!("Z = {z} {}", subshell.label());
528        if z == 0 {
529            return Err(invalid("inner-shell channel", "Z must be at least 1"));
530        }
531        if !(binding_energy_ev.is_finite() && binding_energy_ev > 0.0) {
532            return Err(invalid(
533                "inner-shell channel",
534                format!(
535                    "{name}: binding energy must be finite and positive, got {binding_energy_ev} eV"
536                ),
537            ));
538        }
539        let binding_provenance = checked_text("provenance", Some(binding_provenance))?;
540        if table.axis() != SamplingAxis::InelasticEnergyLoss {
541            return Err(invalid(
542                "inner-shell channel",
543                format!(
544                    "{name}: the table samples {:?}, not the energy loss",
545                    table.axis()
546                ),
547            ));
548        }
549        for i in 0..table.energy_ev().len() {
550            if let Some(q) = table.quantiles(i) {
551                if let Some(w) = q.iter().find(|&&w| w < binding_energy_ev) {
552                    return Err(invalid(
553                        "inner-shell channel",
554                        format!(
555                            "{name}: row {i} ({} eV) has the loss {w} eV below the binding \
556                             energy {binding_energy_ev} eV (a shell contributes only above \
557                             its edge)",
558                            table.energy_ev()[i]
559                        ),
560                    ));
561                }
562            }
563        }
564        Ok(Self {
565            z,
566            subshell,
567            binding_energy_ev,
568            binding_provenance,
569            table,
570        })
571    }
572
573    /// Atomic number of the shell's element.
574    pub fn z(&self) -> u8 {
575        self.z
576    }
577
578    /// The subshell.
579    pub fn subshell(&self) -> Subshell {
580        self.subshell
581    }
582
583    /// Binding energy `B`, eV.
584    pub fn binding_energy_ev(&self) -> f64 {
585        self.binding_energy_ev
586    }
587
588    /// Where `B` came from.
589    pub fn binding_provenance(&self) -> &str {
590        &self.binding_provenance
591    }
592
593    /// The channel's loss table.
594    pub fn table(&self) -> &CrossSectionTable {
595        &self.table
596    }
597
598    /// Format version of the wrapper ([`SHELL_CHANNEL_FORMAT_VERSION`]).
599    pub fn format_version(&self) -> u32 {
600        SHELL_CHANNEL_FORMAT_VERSION
601    }
602
603    /// Parse and validate the TOML file form.
604    pub fn from_toml_str(text: &str) -> Result<Self> {
605        let raw: RawShellChannelIn = toml::from_str(text).map_err(parse_error)?;
606        Self::try_from(raw)
607    }
608
609    /// Serialize to the TOML file form.
610    pub fn to_toml_string(&self) -> Result<String> {
611        toml::to_string(self).map_err(|e| ElectronDataError::Parse(e.to_string()))
612    }
613
614    /// Load a TOML file.
615    pub fn from_toml_file(path: impl AsRef<Path>) -> Result<Self> {
616        Self::from_toml_str(&read_text(path.as_ref())?)
617    }
618
619    /// Write a TOML file.
620    pub fn write_toml_file(&self, path: impl AsRef<Path>) -> Result<()> {
621        let path = path.as_ref();
622        std::fs::write(path, self.to_toml_string()?)
623            .map_err(|e| ElectronDataError::Io(format!("{}: {e}", path.display())))
624    }
625}
lindhard/electron/data/table.rs

lindhard/electron/data/
table.rs
