lindhard/ion/stopping/
table.rs

1//! User-supplied `S_e(E)` tables, loaded from TOML, with mandatory provenance.
2//!
3//! ```toml
4//! provenance = "Author, Journal vol, page (year), Table N"   # required
5//! ion_z = 14
6//! ion_mass_amu = 28.0855            # optional; default: standard weight
7//! target_z = 14
8//! energy_ev = [1.0e3, 1.0e4, 1.0e5]
9//! stopping_ev_1e15_cm2 = [10.0, 30.0, 60.0]   # eV 1e-15 cm^2 per atom
10//! ```
11//!
12//! A table without a non-empty `provenance` is rejected at load time
13//! (see `docs/data-provenance.md`: data without an origin is not admitted).
14//! Do not load SRIM- or ICRU-derived tables for anything committed to the
15//! repository.
16//!
17//! The table is tied to the projectile mass it was declared for (the standard
18//! atomic weight when `ion_mass_amu` is omitted): a query for an ion whose mass
19//! differs by more than [`MASS_RELATIVE_TOLERANCE`] is an error
20//! ([`StoppingError::TableMassMismatch`]), because the same energy at a
21//! different mass is a different ion speed. No energy-axis conversion is
22//! attempted.
23//!
24//! Interpolation is piecewise linear in `ln S` versus `ln E`, which preserves
25//! the monotonicity of each segment of the data. **Out of range** queries are
26//! an error ([`StoppingError::OutOfTableRange`]); nothing is extrapolated or
27//! clamped silently.
28
29use super::{check_mass, from_ev_1e15_cm2, ElectronicStopping, Ion, StoppingError, ValidityRange};
30use serde::Deserialize;
31use std::path::Path;
32
33/// Relative tolerance when comparing a table's declared projectile mass with
34/// the queried ion's mass. Tight enough to separate neighbouring isotopes,
35/// loose enough to absorb rounding in a quoted mass.
36pub const MASS_RELATIVE_TOLERANCE: f64 = 1.0e-4;
37
38#[derive(Deserialize)]
39struct RawTable {
40    provenance: Option<String>,
41    ion_z: u8,
42    ion_mass_amu: Option<f64>,
43    target_z: u8,
44    energy_ev: Vec<f64>,
45    stopping_ev_1e15_cm2: Vec<f64>,
46}
47
48/// A validated user stopping table for one ion/target pair.
49#[derive(Debug, Clone, PartialEq)]
50pub struct StoppingTable {
51    provenance: String,
52    ion_z: u8,
53    ion_mass_amu: f64,
54    target_z: u8,
55    /// ln(E / eV)
56    ln_e: Vec<f64>,
57    /// ln(S / (J m^2))
58    ln_s: Vec<f64>,
59    energy_ev: Vec<f64>,
60}
61
62impl StoppingTable {
63    /// Parse and validate a TOML table.
64    pub fn from_toml_str(text: &str) -> Result<Self, StoppingError> {
65        let raw: RawTable =
66            toml::from_str(text).map_err(|e| StoppingError::InvalidTable(e.to_string()))?;
67        let provenance = raw.provenance.unwrap_or_default().trim().to_string();
68        if provenance.is_empty() {
69            return Err(StoppingError::MissingProvenance);
70        }
71        let ion_element = super::target(raw.ion_z)?;
72        let ion_mass_amu = match raw.ion_mass_amu {
73            Some(m) => {
74                check_mass(m)?;
75                m
76            }
77            None => ion_element.atomic_weight,
78        };
79        super::target(raw.target_z)?;
80        let (e, s) = (raw.energy_ev, raw.stopping_ev_1e15_cm2);
81        if e.len() != s.len() || e.len() < 2 {
82            return Err(StoppingError::InvalidTable(
83                "energy and stopping need the same length, at least 2 points".into(),
84            ));
85        }
86        if e.iter().any(|x| !(x.is_finite() && *x > 0.0))
87            || s.iter().any(|x| !(x.is_finite() && *x > 0.0))
88        {
89            return Err(StoppingError::InvalidTable(
90                "energies and stopping values must be finite and positive".into(),
91            ));
92        }
93        if e.windows(2).any(|w| w[1] <= w[0]) {
94            return Err(StoppingError::InvalidTable(
95                "energies must be strictly increasing".into(),
96            ));
97        }
98        Ok(Self {
99            provenance,
100            ion_z: raw.ion_z,
101            ion_mass_amu,
102            target_z: raw.target_z,
103            ln_e: e.iter().map(|x| x.ln()).collect(),
104            ln_s: s.iter().map(|x| from_ev_1e15_cm2(*x).ln()).collect(),
105            energy_ev: e,
106        })
107    }
108
109    /// Load from a TOML file.
110    pub fn from_toml_file(path: impl AsRef<Path>) -> Result<Self, StoppingError> {
111        let text = std::fs::read_to_string(path.as_ref())
112            .map_err(|e| StoppingError::InvalidTable(e.to_string()))?;
113        Self::from_toml_str(&text)
114    }
115
116    /// The provenance string (to be copied into run metadata).
117    pub fn provenance(&self) -> &str {
118        &self.provenance
119    }
120
121    /// Projectile atomic number the table is for.
122    pub fn ion_z(&self) -> u8 {
123        self.ion_z
124    }
125
126    /// Projectile mass the table is for, u.
127    pub fn ion_mass_amu(&self) -> f64 {
128        self.ion_mass_amu
129    }
130
131    /// Target atomic number the table is for.
132    pub fn target_z(&self) -> u8 {
133        self.target_z
134    }
135
136    /// The tabulated energy range `(min, max)`, eV.
137    pub fn energy_range_ev(&self) -> (f64, f64) {
138        (self.energy_ev[0], self.energy_ev[self.energy_ev.len() - 1])
139    }
140
141    /// Whether the table is for this projectile and target element (the
142    /// projectile mass is checked at query time).
143    pub fn covers_pair(&self, ion_z: u8, target_z: u8) -> bool {
144        self.ion_z == ion_z && self.target_z == target_z
145    }
146}
147
148/// User tables layered over a base model: a query for a pair some table
149/// declares goes to that table (so its range and mass checks apply, and it is
150/// never silently replaced by the base model); every other pair goes to the
151/// base model.
152///
153/// `new` does not check for duplicate pairs; the first matching table wins.
154/// The CLI input layer rejects duplicates with a field-naming error.
155pub struct TableOverride {
156    tables: Vec<StoppingTable>,
157    base: Box<dyn ElectronicStopping + Send + Sync>,
158}
159
160// Manual: the base model is a trait object without `Debug`.
161impl std::fmt::Debug for TableOverride {
162    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
163        f.debug_struct("TableOverride")
164            .field("tables", &self.tables)
165            .finish_non_exhaustive()
166    }
167}
168
169impl TableOverride {
170    /// `tables` take precedence over `base` for the pairs they declare.
171    pub fn new(
172        tables: Vec<StoppingTable>,
173        base: Box<dyn ElectronicStopping + Send + Sync>,
174    ) -> Self {
175        Self { tables, base }
176    }
177}
178
179impl ElectronicStopping for TableOverride {
180    fn name(&self) -> &'static str {
181        "user-table-override"
182    }
183
184    fn stopping(&self, ion: &Ion, target_z: u8, energy_ev: f64) -> Result<f64, StoppingError> {
185        match self
186            .tables
187            .iter()
188            .find(|t| t.covers_pair(ion.z(), target_z))
189        {
190            Some(t) => t.stopping(ion, target_z, energy_ev),
191            None => self.base.stopping(ion, target_z, energy_ev),
192        }
193    }
194
195    /// The base model's advisory range; a table's own range is enforced as an
196    /// error at query time instead.
197    fn validity(&self, ion: &Ion) -> ValidityRange {
198        self.base.validity(ion)
199    }
200}
201
202impl ElectronicStopping for StoppingTable {
203    fn name(&self) -> &'static str {
204        "user-table"
205    }
206
207    fn stopping(&self, ion: &Ion, target_z: u8, energy_ev: f64) -> Result<f64, StoppingError> {
208        super::check_energy(energy_ev)?;
209        if ion.z() != self.ion_z || target_z != self.target_z {
210            return Err(StoppingError::TableMismatch {
211                table_ion: self.ion_z,
212                table_target: self.target_z,
213                ion: ion.z(),
214                target: target_z,
215            });
216        }
217        if (ion.mass_amu() - self.ion_mass_amu).abs() > MASS_RELATIVE_TOLERANCE * self.ion_mass_amu
218        {
219            return Err(StoppingError::TableMassMismatch {
220                table_mass_amu: self.ion_mass_amu,
221                mass_amu: ion.mass_amu(),
222            });
223        }
224        let (lo, hi) = (self.energy_ev[0], self.energy_ev[self.energy_ev.len() - 1]);
225        if energy_ev < lo || energy_ev > hi {
226            return Err(StoppingError::OutOfTableRange {
227                energy_ev,
228                min_ev: lo,
229                max_ev: hi,
230            });
231        }
232        let x = energy_ev.ln();
233        let i = self
234            .ln_e
235            .partition_point(|&v| v <= x)
236            .clamp(1, self.ln_e.len() - 1);
237        let t = (x - self.ln_e[i - 1]) / (self.ln_e[i] - self.ln_e[i - 1]);
238        Ok((self.ln_s[i - 1] + t * (self.ln_s[i] - self.ln_s[i - 1])).exp())
239    }
240
241    fn validity(&self, _ion: &Ion) -> ValidityRange {
242        ValidityRange {
243            min_energy_ev: self.energy_ev[0],
244            max_energy_ev: self.energy_ev[self.energy_ev.len() - 1],
245        }
246    }
247}
248
249#[cfg(test)]
250mod tests {
251    use super::*;
252
253    const GOOD: &str = r#"
254provenance = "synthetic test data, made up for unit tests"
255ion_z = 14
256target_z = 14
257energy_ev = [1.0e3, 1.0e4, 1.0e5]
258stopping_ev_1e15_cm2 = [10.0, 40.0, 90.0]
259"#;
260
261    #[test]
262    fn override_routes_declared_pairs_to_the_table_and_others_to_the_base() {
263        use crate::ion::stopping::lindhard_scharff::LindhardScharff;
264        let t = StoppingTable::from_toml_str(GOOD).unwrap();
265        let o = TableOverride::new(vec![t.clone()], Box::new(LindhardScharff::new()));
266        let si = Ion::new(14).unwrap();
267        assert_eq!(
268            o.stopping(&si, 14, 1.0e4).unwrap(),
269            t.stopping(&si, 14, 1.0e4).unwrap()
270        );
271        // The declared pair is never replaced by the base model: out of range
272        // stays an error.
273        assert!(matches!(
274            o.stopping(&si, 14, 10.0),
275            Err(StoppingError::OutOfTableRange { .. })
276        ));
277        // An undeclared pair uses the base model.
278        let ls = LindhardScharff::new();
279        assert_eq!(
280            o.stopping(&si, 8, 1.0e4).unwrap(),
281            ls.stopping(&si, 8, 1.0e4).unwrap()
282        );
283    }
284
285    #[test]
286    fn rejects_missing_or_blank_provenance() {
287        let no = GOOD.replace(
288            "provenance = \"synthetic test data, made up for unit tests\"\n",
289            "",
290        );
291        let msg = StoppingTable::from_toml_str(&no).unwrap_err().to_string();
292        assert!(msg.contains("provenance"), "{msg}");
293        let blank = GOOD.replace("synthetic test data, made up for unit tests", "  ");
294        assert_eq!(
295            StoppingTable::from_toml_str(&blank).unwrap_err(),
296            StoppingError::MissingProvenance
297        );
298    }
299
300    #[test]
301    fn interpolates_exactly_at_nodes_and_log_log_between() {
302        let t = StoppingTable::from_toml_str(GOOD).unwrap();
303        let ion = Ion::new(14).unwrap();
304        let at = |e| super::super::to_ev_1e15_cm2(t.stopping(&ion, 14, e).unwrap());
305        assert!((at(1.0e4) / 40.0 - 1.0).abs() < 1e-12);
306        assert!((at(1.0e5) / 90.0 - 1.0).abs() < 1e-12);
307        let mid = at(10f64.powf(3.5));
308        assert!((mid / (10.0f64 * 40.0).sqrt() - 1.0).abs() < 1e-12);
309        let mut prev = 0.0;
310        for i in 0..200 {
311            let v = at(1.0e3 * 10f64.powf(f64::from(i) / 100.0));
312            assert!(v > prev);
313            prev = v;
314        }
315    }
316
317    #[test]
318    fn out_of_range_and_mismatch_are_errors() {
319        let t = StoppingTable::from_toml_str(GOOD).unwrap();
320        let ion = Ion::new(14).unwrap();
321        assert!(matches!(
322            t.stopping(&ion, 14, 10.0),
323            Err(StoppingError::OutOfTableRange { .. })
324        ));
325        assert!(matches!(
326            t.stopping(&ion, 8, 1.0e4),
327            Err(StoppingError::TableMismatch { .. })
328        ));
329    }
330
331    #[test]
332    fn declared_mass_must_match_projectile_isotope() {
333        let h1 = GOOD.replace("ion_z = 14", "ion_z = 1\nion_mass_amu = 1.00794");
334        let t = StoppingTable::from_toml_str(&h1).unwrap();
335        // matching isotope (same mass) is accepted
336        let ion = Ion::with_mass(1, 1.00794).unwrap();
337        assert!(t.stopping(&ion, 14, 1.0e4).is_ok());
338        // a deuteron at the same energy is a different ion speed: rejected
339        let d = Ion::with_mass(1, 2.014).unwrap();
340        assert!(matches!(
341            t.stopping(&d, 14, 1.0e4),
342            Err(StoppingError::TableMassMismatch { .. })
343        ));
344    }
345
346    #[test]
347    fn omitted_mass_defaults_to_standard_weight() {
348        let t = StoppingTable::from_toml_str(GOOD).unwrap();
349        assert!(t.stopping(&Ion::new(14).unwrap(), 14, 1.0e4).is_ok());
350        let si28 = Ion::with_mass(14, 27.9769).unwrap();
351        assert!(matches!(
352            t.stopping(&si28, 14, 1.0e4),
353            Err(StoppingError::TableMassMismatch { .. })
354        ));
355    }
356
357    #[test]
358    fn rejects_invalid_declared_mass() {
359        for m in ["0.0", "-1.0", "nan", "inf"] {
360            let bad = GOOD.replace("ion_z = 14", &format!("ion_z = 14\nion_mass_amu = {m}"));
361            assert!(
362                matches!(
363                    StoppingTable::from_toml_str(&bad),
364                    Err(StoppingError::InvalidMass(_))
365                ),
366                "mass {m}"
367            );
368        }
369    }
370
371    #[test]
372    fn rejects_unsorted() {
373        let bad = GOOD.replace("[1.0e3, 1.0e4, 1.0e5]", "[1.0e4, 1.0e3, 1.0e5]");
374        assert!(StoppingTable::from_toml_str(&bad).is_err());
375    }
376}
lindhard/ion/stopping/table.rs

lindhard/ion/stopping/
table.rs
