lindhard/geometry/
stack.rs

1//! The 1D layered [`Stack`]: ordered homogeneous layers, laterally infinite.
2//! See the [module docs](super) for the coordinates and the `Geometry`
3//! contract.
4
5use super::{Exit, ExitOutcome, Face, Geometry, GeometryError};
6use crate::material::Material;
7
8/// One homogeneous layer.
9#[derive(Debug, Clone, PartialEq)]
10pub struct Layer {
11    material: Material,
12    /// Depth of the front face, m.
13    front_m: f64,
14    /// Depth of the back face, m (`INFINITY` for a substrate).
15    back_m: f64,
16}
17
18impl Layer {
19    /// The layer material.
20    pub fn material(&self) -> &Material {
21        &self.material
22    }
23
24    /// Depth of the front face, m.
25    pub fn front_m(&self) -> f64 {
26        self.front_m
27    }
28
29    /// Depth of the back face, m (`f64::INFINITY` for a semi-infinite substrate).
30    pub fn back_m(&self) -> f64 {
31        self.back_m
32    }
33
34    /// Thickness, m (`f64::INFINITY` for a semi-infinite substrate).
35    pub fn thickness_m(&self) -> f64 {
36        self.back_m - self.front_m
37    }
38}
39
40/// A 1D layered target: finite layers front to back, optionally ending in a
41/// semi-infinite substrate.
42#[derive(Debug, Clone, PartialEq)]
43pub struct Stack {
44    layers: Vec<Layer>,
45}
46
47impl Stack {
48    /// Build from `(material, thickness in m)` pairs, front first, and an
49    /// optional semi-infinite substrate behind them.
50    pub fn new(
51        layers: Vec<(Material, f64)>,
52        substrate: Option<Material>,
53    ) -> Result<Self, GeometryError> {
54        if layers.is_empty() && substrate.is_none() {
55            return Err(GeometryError::Empty);
56        }
57        let mut out = Vec::with_capacity(layers.len() + 1);
58        let mut x = 0.0f64;
59        for (index, (material, thickness_m)) in layers.into_iter().enumerate() {
60            if !(thickness_m.is_finite() && thickness_m > 0.0) {
61                return Err(GeometryError::InvalidThickness { index, thickness_m });
62            }
63            let back = x + thickness_m;
64            out.push(Layer {
65                material,
66                front_m: x,
67                back_m: back,
68            });
69            x = back;
70        }
71        if let Some(material) = substrate {
72            out.push(Layer {
73                material,
74                front_m: x,
75                back_m: f64::INFINITY,
76            });
77        }
78        Ok(Self { layers: out })
79    }
80
81    /// A single semi-infinite material.
82    pub fn semi_infinite(material: Material) -> Self {
83        Self {
84            layers: vec![Layer {
85                material,
86                front_m: 0.0,
87                back_m: f64::INFINITY,
88            }],
89        }
90    }
91
92    /// Layers front to back.
93    pub fn layers(&self) -> &[Layer] {
94        &self.layers
95    }
96
97    /// Depth of the back face, m, or `None` if the last layer is semi-infinite.
98    pub fn back_face_m(&self) -> Option<f64> {
99        let b = self.layers.last().expect("non-empty").back_m;
100        b.is_finite().then_some(b)
101    }
102
103    /// Index of the layer containing depth `x`, or `None` outside the target.
104    /// A depth exactly on an interface belongs to the deeper layer.
105    pub fn layer_at(&self, x: f64) -> Option<usize> {
106        if x.is_nan() || x < 0.0 {
107            return None;
108        }
109        self.layers.iter().position(|l| x < l.back_m)
110    }
111}
112
113impl Geometry for Stack {
114    fn as_stack(&self) -> Option<&Stack> {
115        Some(self)
116    }
117
118    fn n_materials(&self) -> usize {
119        self.layers.len()
120    }
121
122    fn material(&self, index: usize) -> &Material {
123        &self.layers[index].material
124    }
125
126    fn n_regions(&self) -> usize {
127        self.layers.len()
128    }
129
130    /// Each layer is its own material entry, so layer data is cached per
131    /// layer.
132    fn material_index(&self, region: usize) -> usize {
133        region
134    }
135
136    fn locate(&self, pos: [f64; 3]) -> Option<usize> {
137        self.layer_at(pos[0])
138    }
139
140    fn entry_point(&self) -> [f64; 3] {
141        [0.0; 3]
142    }
143
144    fn exit(&self, region: usize, pos: [f64; 3], dir: [f64; 3], limit: f64) -> Option<Exit> {
145        let layer = &self.layers[region];
146        let (d, boundary) = if dir[0] > 0.0 {
147            ((layer.back_m - pos[0]) / dir[0], layer.back_m)
148        } else if dir[0] < 0.0 {
149            ((layer.front_m - pos[0]) / dir[0], layer.front_m)
150        } else {
151            return None;
152        };
153        if d.is_nan() || d > limit {
154            return None;
155        }
156        let at = [boundary, pos[1] + d * dir[1], pos[2] + d * dir[2]];
157        let outward = dir[0] > 0.0;
158        let outcome = if !outward && region == 0 {
159            ExitOutcome::Escape {
160                face: Face::Front,
161                normal: [-1.0, 0.0, 0.0],
162            }
163        } else if outward && region + 1 == self.layers.len() {
164            ExitOutcome::Escape {
165                face: Face::Back,
166                normal: [1.0, 0.0, 0.0],
167            }
168        } else {
169            ExitOutcome::Enter {
170                region: if outward { region + 1 } else { region - 1 },
171                pos: at,
172            }
173        };
174        Some(Exit {
175            distance: d,
176            at,
177            outcome,
178        })
179    }
180
181    fn in_vacuum(&self, pos: [f64; 3]) -> bool {
182        pos[0] < 0.0
183    }
184}
185
186#[cfg(test)]
187mod tests {
188    use super::super::test_util::{c, enter, si, unit};
189    use super::*;
190
191    #[test]
192    fn boundaries_accumulate() {
193        let s = Stack::new(vec![(si(), 1e-8), (si(), 2e-8)], Some(si())).unwrap();
194        let b: Vec<_> = s
195            .layers()
196            .iter()
197            .map(|l| (l.front_m(), l.back_m()))
198            .collect();
199        assert_eq!(b[0], (0.0, 1e-8));
200        assert_eq!(b[1], (1e-8, 1e-8 + 2e-8));
201        assert_eq!(b[2].1, f64::INFINITY);
202        assert_eq!(s.back_face_m(), None);
203        assert_eq!(s.layer_at(0.0), Some(0));
204        assert_eq!(s.layer_at(1e-8), Some(1));
205        assert_eq!(s.layer_at(1.0), Some(2));
206        assert_eq!(s.layer_at(-1e-12), None);
207        let f = Stack::new(vec![(si(), 1e-8)], None).unwrap();
208        assert_eq!(f.back_face_m(), Some(1e-8));
209        assert_eq!(f.layer_at(1e-8), None);
210    }
211
212    #[test]
213    fn rejects_bad_input() {
214        assert_eq!(Stack::new(vec![], None), Err(GeometryError::Empty));
215        for t in [0.0, -1.0, f64::NAN, f64::INFINITY] {
216            assert!(matches!(
217                Stack::new(vec![(si(), t)], None),
218                Err(GeometryError::InvalidThickness { index: 0, .. })
219            ));
220        }
221    }
222
223    #[test]
224    fn stack_exit_matches_the_layer_contract() {
225        let s = Stack::new(vec![(si(), 1.0), (c(), 2.0)], None).unwrap();
226        assert_eq!(s.n_materials(), 2);
227        assert_eq!(s.locate([0.0; 3]), Some(0));
228        assert_eq!(s.entry_point(), [0.0; 3]);
229        assert!(s.in_vacuum([-1e-9, 0.0, 0.0]) && !s.in_vacuum([4.0, 0.0, 0.0]));
230        let d = unit([1.0, 1.0, 0.0]);
231        let e = s.exit(0, [0.25, 0.0, 0.0], d, 10.0).unwrap();
232        assert_eq!(e.at[0], 1.0);
233        assert_eq!(enter(&e).0, 1);
234        let e = s.exit(1, [1.0, 0.0, 0.0], [1.0, 0.0, 0.0], 10.0).unwrap();
235        assert_eq!(e.distance, 2.0);
236        assert!(matches!(
237            e.outcome,
238            ExitOutcome::Escape {
239                face: Face::Back,
240                ..
241            }
242        ));
243        let e = s.exit(0, [0.5, 0.0, 0.0], [-1.0, 0.0, 0.0], 10.0).unwrap();
244        assert!(matches!(
245            e.outcome,
246            ExitOutcome::Escape {
247                face: Face::Front,
248                normal: [-1.0, 0.0, 0.0]
249            }
250        ));
251        // Parallel to the interfaces, or too short: no event.
252        assert!(s.exit(0, [0.5; 3], [0.0, 1.0, 0.0], 1e9).is_none());
253        assert!(s.exit(0, [0.5; 3], [1.0, 0.0, 0.0], 0.1).is_none());
254        // Substrate: no back face.
255        let sub = Stack::semi_infinite(si());
256        assert!(sub.exit(0, [0.5; 3], [1.0, 0.0, 0.0], 1e9).is_none());
257    }
258}
lindhard/geometry/stack.rs

lindhard/geometry/
stack.rs
