lindhard/ion/stopping/
mod.rs

1//! Electronic stopping of ions in matter, from published formulas only.
2//!
3//! # Conventions
4//!
5//! * The stopping quantity returned by every model is the **stopping cross
6//!   section per target atom**, `S_e = (1/N) (-dE/dx)`, in **J m² per atom**
7//!   (SI). [`to_ev_1e15_cm2`] converts to the customary eV·10⁻¹⁵ cm²/atom.
8//! * Energies are the ion's kinetic energy in the laboratory frame, in eV
9//!   (`_ev` suffix at the API boundary, as in [`crate::units`]).
10//! * Reduced quantities (`ε`, `ρ`) are the Lindhard-Scharff ones, defined in
11//!   [`lindhard_scharff`].
12//! * Every model is a pure function of `(ion, target Z, energy)`: no RNG, no
13//!   state, so results are bitwise reproducible across threads.
14//!
15//! Validity ranges per model are in `docs/stopping-models.md` and are exposed
16//! at run time by [`ElectronicStopping::validity`].
17//!
18//! No tabulated SRIM or ICRU data enters this module. Terms whose published
19//! coefficients are only available from such tables are omitted, and the
20//! omission is listed in `docs/stopping-models.md`.
21
22pub mod bethe;
23pub mod bragg;
24pub mod correction;
25pub mod dataset;
26pub mod lindhard_scharff;
27pub mod mix;
28pub mod none;
29pub mod oen_robinson;
30pub mod straggling;
31pub mod table;
32
33use crate::constants::{ATOMIC_MASS_UNIT, FINE_STRUCTURE, SPEED_OF_LIGHT};
34use crate::elements::element;
35use crate::units::J_PER_EV;
36
37/// Errors from stopping models, tables and Bragg sums.
38#[derive(Debug, Clone, thiserror::Error, PartialEq)]
39pub enum StoppingError {
40    /// Energy was not finite and positive.
41    #[error("energy must be finite and positive, got {0} eV")]
42    InvalidEnergy(f64),
43    /// Atomic number outside the element table (1..=92).
44    #[error("unknown element Z = {0}")]
45    UnknownElement(u8),
46    /// The model's formula is not meaningful at this energy (for example the
47    /// Bethe logarithm bracket is not positive far below its validity range).
48    #[error("{model} is not applicable at {energy_ev} eV for this ion")]
49    NotApplicable {
50        /// Model name.
51        model: &'static str,
52        /// Energy, eV.
53        energy_ev: f64,
54    },
55    /// A model parameter (for example a mean excitation energy) was not finite,
56    /// or not positive where it must be.
57    #[error("invalid model parameter {name} = {value}")]
58    InvalidParameter {
59        /// Parameter name.
60        name: &'static str,
61        /// Offending value.
62        value: f64,
63    },
64    /// A user table has no (or an empty) provenance field.
65    #[error(
66        "stopping table has no provenance; a citation or description of its origin is required"
67    )]
68    MissingProvenance,
69    /// A user table is malformed.
70    #[error("invalid stopping table: {0}")]
71    InvalidTable(String),
72    /// A line of a stopping dataset ([`dataset`]) is malformed or fails
73    /// validation.
74    #[error("stopping dataset, line {line}: {reason}")]
75    InvalidDataset {
76        /// 1-based line number in the source text.
77        line: usize,
78        /// What is wrong with the line.
79        reason: String,
80    },
81    /// A correction fit ([`correction`]) could not be made, with the reason.
82    #[error("invalid correction fit: {0}")]
83    InvalidFit(String),
84    /// A user table was queried outside its energy range.
85    #[error("energy {energy_ev} eV outside table range [{min_ev}, {max_ev}] eV")]
86    OutOfTableRange {
87        /// Queried energy, eV.
88        energy_ev: f64,
89        /// Table minimum, eV.
90        min_ev: f64,
91        /// Table maximum, eV.
92        max_ev: f64,
93    },
94    /// A user table was queried for a different ion or target than it holds.
95    #[error("table is for Z1={table_ion}, Z2={table_target}; asked for Z1={ion}, Z2={target}")]
96    TableMismatch {
97        /// Table ion Z.
98        table_ion: u8,
99        /// Table target Z.
100        table_target: u8,
101        /// Requested ion Z.
102        ion: u8,
103        /// Requested target Z.
104        target: u8,
105    },
106    /// Projectile mass was not finite and positive.
107    #[error("projectile mass must be finite and positive, got {0} u")]
108    InvalidMass(f64),
109    /// A user table was queried for a projectile mass different from the one
110    /// it was declared for (a different isotope).
111    #[error("table is for an ion of mass {table_mass_amu} u; asked for {mass_amu} u")]
112    TableMassMismatch {
113        /// Mass the table was declared for, u.
114        table_mass_amu: f64,
115        /// Requested projectile mass, u.
116        mass_amu: f64,
117    },
118    /// A material error while summing a compound.
119    #[error(transparent)]
120    Material(#[from] crate::material::MaterialError),
121}
122
123/// A projectile: atomic number and mass.
124#[derive(Debug, Clone, Copy, PartialEq)]
125///
126/// Fields are private so that every `Ion` satisfies its invariants: a known
127/// element and a finite, positive mass.
128pub struct Ion {
129    z: u8,
130    mass_amu: f64,
131}
132
133impl Ion {
134    /// An ion of element `z` with the standard atomic weight as its mass.
135    pub fn new(z: u8) -> Result<Self, StoppingError> {
136        let e = element(z).ok_or(StoppingError::UnknownElement(z))?;
137        Ok(Self {
138            z,
139            mass_amu: e.atomic_weight,
140        })
141    }
142
143    /// An ion with an explicit (isotopic) mass in amu.
144    ///
145    /// The mass must be finite and positive.
146    pub fn with_mass(z: u8, mass_amu: f64) -> Result<Self, StoppingError> {
147        element(z).ok_or(StoppingError::UnknownElement(z))?;
148        check_mass(mass_amu)?;
149        Ok(Self { z, mass_amu })
150    }
151
152    /// Atomic number `Z1`.
153    pub fn z(&self) -> u8 {
154        self.z
155    }
156
157    /// Mass in unified atomic mass units.
158    pub fn mass_amu(&self) -> f64 {
159        self.mass_amu
160    }
161
162    /// A proton (CODATA proton mass).
163    pub fn proton() -> Self {
164        Self {
165            z: 1,
166            mass_amu: crate::constants::PROTON_MASS / ATOMIC_MASS_UNIT,
167        }
168    }
169
170    /// Mass in kg.
171    pub fn mass_kg(&self) -> f64 {
172        self.mass_amu * ATOMIC_MASS_UNIT
173    }
174}
175
176/// Advisory energy range over which a model is meant to be used, eV.
177#[derive(Debug, Clone, Copy, PartialEq)]
178pub struct ValidityRange {
179    /// Lower bound, eV (0 if none).
180    pub min_energy_ev: f64,
181    /// Upper bound, eV (`INFINITY` if none).
182    pub max_energy_ev: f64,
183}
184
185impl ValidityRange {
186    /// True if `energy_ev` is inside the range.
187    pub fn contains(&self, energy_ev: f64) -> bool {
188        energy_ev >= self.min_energy_ev && energy_ev <= self.max_energy_ev
189    }
190}
191
192/// An electronic stopping model.
193pub trait ElectronicStopping {
194    /// Short stable model name (used in output metadata).
195    fn name(&self) -> &'static str;
196
197    /// Stopping cross section per atom of element `target_z`, J m².
198    fn stopping(&self, ion: &Ion, target_z: u8, energy_ev: f64) -> Result<f64, StoppingError>;
199
200    /// Advisory validity range for this ion.
201    fn validity(&self, ion: &Ion) -> ValidityRange;
202
203    /// If this model is exactly `S(E) = c sqrt(E)` for this ion and target
204    /// (the velocity-proportional Lindhard-Scharff form), return `c` in
205    /// J m² eV^(-1/2); `None` (the default) otherwise.
206    ///
207    /// A model that returns `Some(c)` promises that
208    /// [`stopping`](Self::stopping) succeeds for every finite positive energy
209    /// and equals `c sqrt(E)` to rounding. A caller that evaluates the model
210    /// many times for one ion and target (the collision loop) can then
211    /// compute `c`, which involves several fractional powers, once. Return
212    /// `None` whenever [`stopping`](Self::stopping) could fail, so the error
213    /// is still reported at evaluation.
214    fn sqrt_energy_coefficient(&self, _ion: &Ion, _target_z: u8) -> Option<f64> {
215        None
216    }
217}
218
219/// Convert J m² per atom to eV·10⁻¹⁵ cm² per atom (1e-15 cm² = 1e-19 m²).
220pub fn to_ev_1e15_cm2(s_j_m2: f64) -> f64 {
221    s_j_m2 / J_PER_EV / 1.0e-19
222}
223
224/// Convert eV·10⁻¹⁵ cm² per atom to J m² per atom.
225pub fn from_ev_1e15_cm2(s: f64) -> f64 {
226    s * 1.0e-19 * J_PER_EV
227}
228
229/// Bohr velocity `v0 = α c`, m/s.
230pub fn bohr_velocity() -> f64 {
231    FINE_STRUCTURE * SPEED_OF_LIGHT
232}
233
234/// Ion speed from kinetic energy, non-relativistic `sqrt(2E/M)`, m/s.
235pub fn speed_nonrel(ion: &Ion, energy_ev: f64) -> f64 {
236    (2.0 * energy_ev * J_PER_EV / ion.mass_kg()).sqrt()
237}
238
239pub(crate) fn check_mass(mass_amu: f64) -> Result<(), StoppingError> {
240    if mass_amu.is_finite() && mass_amu > 0.0 {
241        Ok(())
242    } else {
243        Err(StoppingError::InvalidMass(mass_amu))
244    }
245}
246
247pub(crate) fn check_energy(energy_ev: f64) -> Result<(), StoppingError> {
248    if energy_ev.is_finite() && energy_ev > 0.0 {
249        Ok(())
250    } else {
251        Err(StoppingError::InvalidEnergy(energy_ev))
252    }
253}
254
255pub(crate) fn target(z: u8) -> Result<&'static crate::elements::Element, StoppingError> {
256    element(z).ok_or(StoppingError::UnknownElement(z))
257}
258
259/// Energy per amu (eV/u) at which the ion speed equals `n_v0 * v0 * Z1^(2/3)`:
260/// the Lindhard-Scharff velocity scale that separates the low-energy and
261/// Bethe regimes.
262pub(crate) fn velocity_scale_energy_per_amu_ev(z1: u8, n_v0: f64) -> f64 {
263    let v = n_v0 * bohr_velocity() * f64::from(z1).powf(2.0 / 3.0);
264    0.5 * ATOMIC_MASS_UNIT * v * v / J_PER_EV
265}
266
267#[cfg(test)]
268mod tests {
269    use super::*;
270
271    #[test]
272    fn with_mass_rejects_invalid_mass() {
273        for m in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
274            assert!(
275                matches!(Ion::with_mass(1, m), Err(StoppingError::InvalidMass(_))),
276                "mass {m}"
277            );
278        }
279        assert!(Ion::with_mass(1, 2.014).is_ok());
280    }
281
282    #[test]
283    fn unit_round_trip() {
284        let s = 12.5;
285        assert!((to_ev_1e15_cm2(from_ev_1e15_cm2(s)) - s).abs() < 1e-12);
286    }
287
288    #[test]
289    fn models_are_bitwise_deterministic_across_threads() {
290        let ion = Ion::new(15).unwrap();
291        let m = bethe::BetheBloch::new();
292        let r0 = m.stopping(&ion, 14, 1.2e8).unwrap();
293        let hs: Vec<_> = (0..4)
294            .map(|_| {
295                std::thread::spawn(move || {
296                    let ion = Ion::new(15).unwrap();
297                    bethe::BetheBloch::new().stopping(&ion, 14, 1.2e8).unwrap()
298                })
299            })
300            .collect();
301        for h in hs {
302            assert_eq!(h.join().unwrap().to_bits(), r0.to_bits());
303        }
304        assert_eq!(m.stopping(&ion, 14, 1.2e8).unwrap().to_bits(), r0.to_bits());
305    }
306}
lindhard/ion/stopping/mod.rs

lindhard/ion/stopping/
mod.rs
