treams_core/coeffs/mie_cyl/
saved.rs

1//! Fixed numerical state for a cylindrical coefficient record, without recomputation.
2//!
3//! treams-rs extension.
4
5use super::{Interface, Material, Matrix2, Matrix4, Matrix42, MieCylResidual, RadialJet, Side};
6use crate::{
7    Error, MAX_DEGREE, Result,
8    saved::{Reader, SavedState, Writer},
9};
10
11fn invalid() -> Error {
12    Error::InvalidInput("invalid cylindrical Mie derivative state".into())
13}
14
15impl MieCylResidual {
16    /// Byte count of the saved numerical state for `boundaries` layer boundaries.
17    pub fn state_size(boundaries: usize) -> Result<usize> {
18        if boundaries == 0 {
19            return Err(invalid());
20        }
21        // Two integer headers, n+2 real numbers, and 81*n+5 complex numbers.
22        boundaries
23            .checked_mul(1304)
24            .and_then(|size| size.checked_add(112))
25            .ok_or_else(invalid)
26    }
27
28    #[allow(clippy::cast_sign_loss)] // The recorded order is in -MAX_DEGREE..=MAX_DEGREE.
29    pub(crate) fn write_state(&self, writer: &mut Writer) {
30        writer.usize(self.radii.len());
31        writer.usize((self.order + MAX_DEGREE) as usize);
32        writer.f64(self.kz);
33        writer.f64(self.k0);
34        for &radius in &self.radii {
35            writer.f64(radius);
36        }
37        for material in &self.materials {
38            writer.complex(material.epsilon);
39            writer.complex(material.mu);
40            writer.complex(material.kappa);
41        }
42        for interface in &self.interfaces {
43            interface.inside.write_state(writer);
44            interface.outside.write_state(writer);
45            for &value in interface
46                .inverse_outer
47                .iter()
48                .chain(interface.transfer.iter())
49                .chain(interface.before.iter())
50            {
51                writer.complex(value);
52            }
53        }
54        for &value in self.inverse.iter().chain(self.value.iter()) {
55            writer.complex(value);
56        }
57    }
58
59    pub(crate) fn read_state(reader: &mut Reader<'_>) -> Result<Self> {
60        let boundaries = reader.usize()?;
61        let size = Self::state_size(boundaries)?;
62        if reader.remaining_len() < size - 8 {
63            return Err(invalid());
64        }
65        let encoded = i32::try_from(reader.usize()?).map_err(|_| invalid())?;
66        if encoded > 2 * MAX_DEGREE {
67            return Err(invalid());
68        }
69        let order = encoded - MAX_DEGREE;
70        let kz = reader.f64()?;
71        let k0 = reader.f64()?;
72        let radii = (0..boundaries)
73            .map(|_| reader.f64())
74            .collect::<Result<Vec<_>>>()?;
75        let materials = (0..=boundaries)
76            .map(|_| {
77                Ok(Material {
78                    epsilon: reader.complex()?,
79                    mu: reader.complex()?,
80                    kappa: reader.complex()?,
81                })
82            })
83            .collect::<Result<Vec<_>>>()?;
84        let mut interfaces = Vec::with_capacity(boundaries);
85        for i in 0..boundaries {
86            let inside = Side::read_state(reader, i == 0)?;
87            let outside = Side::read_state(reader, false)?;
88            let mut inverse_outer = Matrix4::zeros();
89            let mut transfer = Matrix4::zeros();
90            let mut before = Matrix42::zeros();
91            for value in inverse_outer
92                .iter_mut()
93                .chain(transfer.iter_mut())
94                .chain(before.iter_mut())
95            {
96                *value = reader.complex()?;
97            }
98            interfaces.push(Interface {
99                inside,
100                outside,
101                inverse_outer,
102                transfer,
103                before,
104            });
105        }
106        let mut inverse = Matrix2::zeros();
107        let mut value = Matrix2::zeros();
108        for entry in inverse.iter_mut().chain(value.iter_mut()) {
109            *entry = reader.complex()?;
110        }
111        Ok(Self {
112            order,
113            kz,
114            k0,
115            radii,
116            materials,
117            interfaces,
118            inverse,
119            value,
120        })
121    }
122}
123
124impl Side {
125    fn write_state(&self, writer: &mut Writer) {
126        for value in self
127            .ks
128            .into_iter()
129            .chain([self.impedance])
130            .chain(self.kr)
131            .chain(self.x)
132        {
133            writer.complex(value);
134        }
135        // The first inside side has only regular waves. Its absent singular
136        // entries follow from the boundary topology and need no tags.
137        for radial in self.radial.iter().flatten().flatten() {
138            writer.complex(radial.value);
139            writer.complex(radial.first);
140            writer.complex(radial.second);
141        }
142    }
143
144    fn read_state(reader: &mut Reader<'_>, regular_only: bool) -> Result<Self> {
145        let ks = [reader.complex()?, reader.complex()?];
146        let impedance = reader.complex()?;
147        let kr = [reader.complex()?, reader.complex()?];
148        let x = [reader.complex()?, reader.complex()?];
149        let mut radial = [[None; 2]; 2];
150        for polarization in &mut radial {
151            for entry in polarization
152                .iter_mut()
153                .take(if regular_only { 1 } else { 2 })
154            {
155                *entry = Some(RadialJet {
156                    value: reader.complex()?,
157                    first: reader.complex()?,
158                    second: reader.complex()?,
159                });
160            }
161        }
162        Ok(Self {
163            ks,
164            impedance,
165            kr,
166            x,
167            radial,
168        })
169    }
170}
171
172impl SavedState for MieCylResidual {
173    fn save_state(&self) -> Result<Vec<u8>> {
174        let mut writer = Writer::new(Self::state_size(self.boundaries())?);
175        self.write_state(&mut writer);
176        Ok(writer.finish())
177    }
178
179    fn from_state(bytes: &[u8]) -> Result<Self> {
180        let mut reader = Reader::new(bytes);
181        let result = Self::read_state(&mut reader)?;
182        reader.finish()?;
183        Ok(result)
184    }
185}
186
187#[cfg(test)]
188mod tests {
189    use super::*;
190    use crate::{Complex, coeffs::mie_cyl};
191
192    fn record(order: i32, boundaries: usize) -> MieCylResidual {
193        let radii: Vec<_> = (0..boundaries)
194            .map(|i| 0.2 + 0.3 * f64::from(u32::try_from(i).unwrap()))
195            .collect();
196        let materials: Vec<_> = (0..=boundaries)
197            .map(|i| Material {
198                epsilon: Complex::new(2.0 + 0.4 * f64::from(u32::try_from(i).unwrap()), 0.1),
199                mu: Complex::new(1.2, 0.03),
200                kappa: Complex::new(0.02, 0.004),
201            })
202            .collect();
203        mie_cyl(0.3, order, 1.2, &radii, &materials).unwrap()
204    }
205
206    #[test]
207    fn saved_state_preserves_values_and_both_derivatives() {
208        for boundaries in [1, 2, 3] {
209            for order in [-3, 0, 2] {
210                let original = record(order, boundaries);
211                let bytes = original.save_state().unwrap();
212                assert_eq!(bytes.len(), MieCylResidual::state_size(boundaries).unwrap());
213                let restored = MieCylResidual::from_state(&bytes).unwrap();
214                assert_eq!(restored.save_state().unwrap(), bytes);
215                assert_eq!(restored.value(), original.value());
216                let radii = vec![0.07; boundaries];
217                let materials = vec![
218                    Material {
219                        epsilon: Complex::new(0.1, 0.03),
220                        mu: Complex::new(-0.02, 0.01),
221                        kappa: Complex::new(0.04, -0.02),
222                    };
223                    boundaries + 1
224                ];
225                assert_eq!(
226                    restored
227                        .pushforward(0.03, -0.04, &radii, &materials)
228                        .unwrap(),
229                    original
230                        .pushforward(0.03, -0.04, &radii, &materials)
231                        .unwrap()
232                );
233                let cotangent = Matrix2::from_fn(|i, j| {
234                    Complex::new(
235                        0.3 + f64::from(u32::try_from(i).unwrap()) * 0.1,
236                        -0.2 + f64::from(u32::try_from(j).unwrap()) * 0.15,
237                    )
238                });
239                let actual = restored.pullback(&cotangent).unwrap();
240                let expected = original.pullback(&cotangent).unwrap();
241                assert_eq!(actual.kz.to_bits(), expected.kz.to_bits());
242                assert_eq!(actual.k0.to_bits(), expected.k0.to_bits());
243                assert_eq!(actual.layers.radii, expected.layers.radii);
244                assert_eq!(actual.layers.epsilon, expected.layers.epsilon);
245                assert_eq!(actual.layers.mu, expected.layers.mu);
246                assert_eq!(actual.layers.kappa, expected.layers.kappa);
247            }
248        }
249    }
250
251    #[test]
252    fn saved_state_rejects_invalid_counts_orders_and_lengths() {
253        let bytes = record(1, 2).save_state().unwrap();
254        for length in [0, 7, 8, 15, bytes.len() - 1] {
255            assert!(MieCylResidual::from_state(&bytes[..length]).is_err());
256        }
257        let mut trailing = bytes.clone();
258        trailing.push(0);
259        assert!(MieCylResidual::from_state(&trailing).is_err());
260        for boundaries in [0_u64, 1, u64::MAX] {
261            let mut invalid = bytes.clone();
262            invalid[..8].copy_from_slice(&boundaries.to_le_bytes());
263            assert!(MieCylResidual::from_state(&invalid).is_err());
264        }
265        for boundaries in [0_u64, 3, u64::MAX] {
266            let mut invalid = bytes.clone();
267            invalid[..8].copy_from_slice(&boundaries.to_le_bytes());
268            assert!(MieCylResidual::read_state(&mut Reader::new(&invalid)).is_err());
269        }
270        for order in [257_u64, u64::MAX] {
271            let mut invalid = bytes.clone();
272            invalid[8..16].copy_from_slice(&order.to_le_bytes());
273            assert!(MieCylResidual::from_state(&invalid).is_err());
274        }
275        assert!(MieCylResidual::state_size(usize::MAX).is_err());
276    }
277}
treams_core/coeffs/mie_cyl/saved.rs

treams_core/coeffs/mie_cyl/
saved.rs
