treams_core/coeffs/
material.rs

1//! Isotropic, optionally chiral materials and the validation of concentric layers.
2//!
3//! Upstream: `treams.Material` and `treams.misc.refractive_index`.
4#![allow(clippy::indexing_slicing)] // Adjacent boundaries from `windows(2)`.
5
6use std::ops::AddAssign;
7
8use crate::{Complex, Error, Result, numerics::finite};
9
10/// Isotropic reciprocal material, optionally chiral.
11///
12/// Upstream: `treams.Material`. Differences: [`index`](Self::index) keeps the
13/// principal square root, as `treams.coeffs.mie` does, where `Material.n` flips its
14/// sign for a negative imaginary part.
15#[derive(Clone, Copy, Debug)]
16pub struct Material {
17    /// Relative permittivity.
18    pub epsilon: Complex,
19    /// Relative permeability.
20    pub mu: Complex,
21    /// Chirality parameter.
22    pub kappa: Complex,
23}
24
25impl Default for Material {
26    fn default() -> Self {
27        Self {
28            epsilon: Complex::new(1.0, 0.0),
29            mu: Complex::new(1.0, 0.0),
30            kappa: Complex::default(),
31        }
32    }
33}
34
35impl Material {
36    /// Directional derivatives of the derived optical parameters. `direction` stores
37    /// changes of epsilon, mu and kappa, rather than another physical material.
38    pub(super) fn pushforward(self, direction: Self) -> MaterialTangent {
39        MaterialTangent {
40            index: (self.mu * direction.epsilon + self.epsilon * direction.mu)
41                / (2.0 * self.index()),
42            impedance: self.impedance()
43                * (direction.mu / self.mu - direction.epsilon / self.epsilon)
44                / 2.0,
45            kappa: direction.kappa,
46        }
47    }
48
49    /// The refractive index `sqrt(epsilon mu)` on the principal branch.
50    #[must_use]
51    pub fn index(self) -> Complex {
52        crate::numerics::complex_sqrt(self.epsilon * self.mu)
53    }
54    /// The negative and positive helicity indices `sqrt(epsilon mu) -+ kappa`, each
55    /// negated if its imaginary part is negative.
56    ///
57    /// Upstream: `treams.misc.refractive_index`.
58    #[must_use]
59    pub fn indices(self) -> [Complex; 2] {
60        let n = self.index();
61        [n - self.kappa, n + self.kappa].map(|v| if v.im < 0.0 { -v } else { v })
62    }
63    /// The relative impedance `sqrt(mu / epsilon)` on the principal branch.
64    ///
65    /// Upstream: `treams.Material.impedance`.
66    #[must_use]
67    pub fn impedance(self) -> Complex {
68        crate::numerics::complex_sqrt(self.mu / self.epsilon)
69    }
70    /// Derivatives of the index, impedance and chirality with respect to epsilon, mu
71    /// and kappa, in that order.
72    pub(super) fn tangents(self) -> [MaterialTangent; 3] {
73        let (n, impedance) = (self.index(), self.impedance());
74        let zero = Complex::default();
75        [
76            MaterialTangent {
77                index: self.mu / (2.0 * n),
78                impedance: -impedance / (2.0 * self.epsilon),
79                kappa: zero,
80            },
81            MaterialTangent {
82                index: self.epsilon / (2.0 * n),
83                impedance: impedance / (2.0 * self.mu),
84                kappa: zero,
85            },
86            MaterialTangent {
87                index: zero,
88                impedance: zero,
89                kappa: Complex::new(1.0, 0.0),
90            },
91        ]
92    }
93}
94
95/// Negative- and positive-helicity refractive indices `sqrt(epsilon mu) -+ kappa` of a
96/// real material.
97///
98/// A negative product `epsilon mu` gives NaN, as `numpy.sqrt` of a negative float does.
99/// Complex materials use [`Material::indices`].
100///
101/// Upstream: `treams.misc.refractive_index` with real arguments.
102#[inline]
103#[must_use]
104pub fn real_refractive_indices(epsilon: f64, mu: f64, kappa: f64) -> [f64; 2] {
105    let n = (epsilon * mu).sqrt();
106    [n - kappa, n + kappa]
107}
108
109/// Gradients of concentric layers: the radii of their boundaries and their media.
110#[derive(Clone, Debug)]
111pub struct LayerGradient {
112    /// Radius gradients, from inner to outer boundary.
113    pub radii: Vec<f64>,
114    /// Permittivity gradients, including the embedding medium.
115    pub epsilon: Vec<Complex>,
116    /// Permeability gradients, including the embedding medium.
117    pub mu: Vec<Complex>,
118    /// Chirality gradients, including the embedding medium.
119    pub kappa: Vec<Complex>,
120}
121
122impl LayerGradient {
123    /// Zero gradients of `boundaries` radii and `boundaries + 1` media.
124    pub(crate) fn zeros(boundaries: usize) -> Self {
125        let media = vec![Complex::default(); boundaries + 1];
126        Self {
127            radii: vec![0.0; boundaries],
128            epsilon: media.clone(),
129            mu: media.clone(),
130            kappa: media,
131        }
132    }
133
134    /// Add the gradient of the same layers along another cotangent.
135    pub(crate) fn accumulate(&mut self, other: &Self) {
136        add_entries(&mut self.radii, &other.radii);
137        add_entries(&mut self.epsilon, &other.epsilon);
138        add_entries(&mut self.mu, &other.mu);
139        add_entries(&mut self.kappa, &other.kappa);
140    }
141}
142
143/// `sums[i] += values[i]` for every entry.
144pub(super) fn add_entries<T: Copy + AddAssign>(sums: &mut [T], values: &[T]) {
145    for (sum, &value) in sums.iter_mut().zip(values) {
146        *sum += value;
147    }
148}
149
150/// Derivatives of a material's index, impedance and chirality along one parameter.
151#[derive(Clone, Copy, Debug)]
152pub(super) struct MaterialTangent {
153    pub(super) index: Complex,
154    pub(super) impedance: Complex,
155    pub(super) kappa: Complex,
156}
157
158/// Shared material and concentric-boundary validation.
159pub(crate) fn validate_layers(sizes: &[f64], materials: &[Material]) -> Result<()> {
160    if sizes.is_empty() || materials.len() != sizes.len() + 1 {
161        return Err(Error::InvalidInput(
162            "require one more material than radii".into(),
163        ));
164    }
165    if sizes.iter().any(|x| !x.is_finite() || *x <= 0.0) || sizes.windows(2).any(|x| x[0] >= x[1]) {
166        return Err(Error::InvalidInput(
167            "radii and size parameters must be finite, positive and strictly increasing".into(),
168        ));
169    }
170    for material in materials {
171        if !finite(material.epsilon)
172            || !finite(material.mu)
173            || !finite(material.kappa)
174            || material.epsilon.norm_sqr() == 0.0
175            || material.mu.norm_sqr() == 0.0
176        {
177            return Err(Error::InvalidInput(
178                "require finite materials with nonzero epsilon and mu".into(),
179            ));
180        }
181    }
182    Ok(())
183}
184
185/// Validate a direction of the layer inputs, without imposing physical material or
186/// positive-radius constraints on derivatives.
187pub(crate) fn validate_layer_tangents(
188    boundaries: usize,
189    sizes: &[f64],
190    materials: &[Material],
191) -> Result<()> {
192    if sizes.len() != boundaries
193        || materials.len() != boundaries + 1
194        || sizes.iter().any(|x| !x.is_finite())
195        || materials
196            .iter()
197            .any(|m| !finite(m.epsilon) || !finite(m.mu) || !finite(m.kappa))
198    {
199        return Err(Error::InvalidInput(
200            "layer tangents must be finite and match the recorded inputs".into(),
201        ));
202    }
203    Ok(())
204}
treams_core/coeffs/material.rs

treams_core/coeffs/
material.rs
