treams_core/lattice/
spectral.rs

1//! The spectral series of 1D spherical sums off the axis, and when to take it.
2//!
3//! treams-rs extension: treams has only the Ewald sum.
4
5use std::f64::consts::PI;
6
7use super::{
8    accuracy::{MAX_LOSS, Rounding, derivative_scale},
9    reciprocal::THRESHOLD_DISTANCE,
10    shells::largest_modulus,
11    wave::azimuthal,
12};
13use crate::{
14    Complex, Error, Result,
15    numerics::Jet,
16    special::{self, harmonic_normalization},
17};
18
19/// Diffraction orders the spectral series of a 1D spherical sum may take on each side
20/// of zero before it is not used.
21const SERIES_SHELLS: f64 = 512.0;
22
23/// [`spectral_sw1d`] is used only where its terms fall by `e^-SERIES_DECAY` past their
24/// peak within [`SERIES_SHELLS`] orders.
25const SERIES_DECAY: f64 = 45.0;
26
27/// Part of the summed moduli below which a shell of [`spectral_sw1d`] counts as quiet;
28/// the series stops after two quiet shells past the peak.
29const QUIET_SHELL: f64 = 1e-18;
30
31/// The spectral series of a 1D spherical sum along z off the axis, with first-order
32/// bounds on its rounding errors in units of epsilon.
33///
34/// # Formula
35///
36/// By Poisson's formula over the chain,
37/// `S = (-i)^l i^|m| (-1)^m N_lm pi / (k a) ((x ± i y) / (k rho^2))^|m|
38/// sum_q F_l^m(q / k) u_q^|m| H_|m|(u_q) exp(-i q z)` over the diffraction orders
39/// `q = kpar + g reciprocal`, with `u_q = k_q rho`, `k_q = sqrt(k^2 - q^2)` and the
40/// Legendre factor `F_l^m = P_l^m / sin^|m|` continued with `sin = k_q / k`.
41///
42/// # Branch
43///
44/// `Im k_q >= 0`: the root of the Ewald reciprocal sums at every split (see
45/// [`split_sheet`]).
46///
47/// # Why
48///
49/// Evanescent orders fall like `(|q| / |k|)^l exp(-|q| rho)`, so the series converges
50/// fast off the axis, where the Ewald sum cancels. The value stops on its own moduli
51/// and the derivatives on the largest of all, both at [`QUIET_SHELL`] of them, so
52/// value-only and jet evaluations return the same value and bound.
53///
54/// `None` on the axis, for `Re k <= 0`, next to a diffraction threshold, beyond
55/// [`SERIES_SHELLS`] orders on either side or at a non-finite term.
56///
57/// [`split_sheet`]: super::sheets::split_sheet
58pub(super) fn spectral_sw1d<const N: usize>(
59    (l, m): (i32, i32),
60    k: Jet<N>,
61    r: [Jet<N>; 3],
62    kpar: Jet<N>,
63    reciprocal: Jet<N>,
64    period: Jet<N>,
65) -> Option<(Jet<N>, Rounding<N>)> {
66    let rho = (r[0] * r[0] + r[1] * r[1]).sqrt();
67    let spacing = reciprocal.value.norm();
68    if rho.value.re <= 0.0 || k.value.re <= 0.0 || spacing <= 0.0 {
69        return None;
70    }
71    let order = m.unsigned_abs();
72    // Terms fall monotonically once |q| exceeds `peak`; the orders to reach e^-45 past it.
73    let peak = k.value.norm() + f64::from(l + 1) / rho.value.re;
74    let reach = kpar.value.norm() + peak + SERIES_DECAY / rho.value.re;
75    if reach / spacing > SERIES_SHELLS {
76        return None;
77    }
78    let i = Complex::i();
79    let mut sum = Jet::default();
80    let mut bound = Rounding::default();
81    // Summed moduli and consecutive quiet shells of the value and of every component,
82    // and the value with its bound once its own terms are quiet.
83    let (mut magnitude, mut quiet) = ([0.0; 2], [0; 2]);
84    let mut value = None;
85    for shell in 0_i32.. {
86        let mut shell_magnitude = [0.0; 2];
87        for g in [shell, -shell]
88            .into_iter()
89            .take(if shell == 0 { 1 } else { 2 })
90        {
91            let q = kpar + f64::from(g) * reciprocal;
92            let square = k * k - q * q;
93            let mut root = square.value;
94            if root.im == 0.0 {
95                // An evanescent order of real k takes the +0 side.
96                root.im = 0.0;
97            }
98            // Next to a diffraction threshold the Ewald sum reports it instead.
99            if root.norm() <= THRESHOLD_DISTANCE * (k.value * k.value).norm() {
100                return None;
101            }
102            let mut root = crate::numerics::complex_sqrt(root);
103            if root.im < 0.0 {
104                root = -root;
105            }
106            let argument = square.chain(root, 0.5 / root) * rho;
107            // u^|m| H_|m|(u), whose derivative u^|m| H_(|m|-1)(u) the product rule
108            // would form by cancellation next to thresholds.
109            let [below, hankel] = special::hankel_below(order, argument.value).ok()?;
110            let power = argument.value.powu(order);
111            let radial = argument.chain(power * hankel, power * below);
112            let axial = q * r[2];
113            let term = special::legendre_factor(l, m, q / k) * radial * (-i * axial).exp();
114            if !term.finite() {
115                return None;
116            }
117            // The Hankel function and the phase lose |argument| and |q z| ulps, the
118            // Legendre recurrence about l.
119            let ulps = 8.0 + f64::from(l) + argument.value.norm() + axial.value.norm();
120            bound += Rounding::relative(&term, ulps);
121            sum += term;
122            shell_magnitude[0] += term.value.norm();
123            if N > 0 {
124                shell_magnitude[1] += largest_modulus(&term);
125            }
126        }
127        let past = f64::from(shell) * spacing - kpar.value.norm() > peak;
128        for ((total, count), added) in magnitude.iter_mut().zip(&mut quiet).zip(shell_magnitude) {
129            *total += added;
130            *count = if past && added <= QUIET_SHELL * *total {
131                *count + 1
132            } else {
133                0
134            };
135        }
136        if value.is_none() && quiet[0] >= 2 {
137            value = Some((sum.value, bound.value));
138        }
139        if let Some((complete, complete_bound)) = value
140            && (N == 0 || quiet[1] >= 2)
141        {
142            sum.value = complete;
143            bound.value = complete_bound;
144            break;
145        }
146        if f64::from(shell) > SERIES_SHELLS {
147            return None;
148        }
149    }
150    let sign = if order % 2 == 0 { 1.0 } else { -1.0 };
151    let order = m.abs();
152    let outer = (-i).powi(l) * i.powi(order) * (sign * harmonic_normalization(l, m) * PI)
153        / (k * period)
154        * (azimuthal(r[0], r[1], m) / (k * rho * rho)).powi(order);
155    Some((outer * sum, bound.times(&outer)))
156}
157
158/// The `w = |k| rho |eta|` from which a 1D spherical sum tries its spectral series
159/// before the Ewald sum, whose reciprocal part then cancels by `e^(w^2)` or more.
160pub(super) const SPECTRAL_FIRST_W: f64 = 2.5;
161
162/// Rounding error of the real-space part of a 1D spherical sum relative to its moduli,
163/// in ulps: its Kambe integrals of even order keep their relative accuracy.
164pub(super) const REAL_ULPS: f64 = 64.0;
165
166/// Relative rounding error, on the scale of `max(|S|, 1)`, from which a component of
167/// a 1D spherical Ewald sum is taken from its spectral series instead, if that bounds
168/// it more tightly. Derivatives count on the scale `|dS| + max(|S|, 1) s` with
169/// `s = |k|` for length derivatives and `1 / |k|` for inverse-length ones.
170pub(super) const SPECTRAL_SWITCH_LOSS: f64 = 1e-12;
171
172/// A 1D spherical Ewald sum `ewald` with the rounding bounds `rounding` of its components
173/// (the value's with the predicted loss of the [`cancellation`] of its parts), completed
174/// from the spectral series: each component whose bound exceeds [`SPECTRAL_SWITCH_LOSS`]
175/// comes from the series where that bounds it more tightly, and a failed sum is replaced
176/// by the series. The choice of the value reads only values and value bounds, as in
177/// value-only evaluations. Fails where the value's bound still exceeds [`MAX_LOSS`].
178///
179/// [`cancellation`]: super::accuracy::cancellation
180pub(super) fn prefer_spectral_sw1d<const N: usize>(
181    ewald: Result<Jet<N>>,
182    mut rounding: Rounding<N>,
183    series: impl FnOnce() -> Option<(Jet<N>, Rounding<N>)>,
184    k: Complex,
185) -> Result<Jet<N>> {
186    let accurate =
187        |value: Complex, bound: f64| f64::EPSILON * bound <= MAX_LOSS * value.norm().max(1.0);
188    let mut sum = match ewald {
189        Ok(sum) => sum,
190        // Non-finite summands, failed special functions or no convergence of either part.
191        Err(error) if error.is_numerical() => {
192            return match series() {
193                Some((series, bound)) if accurate(series.value, bound.value) => Ok(series),
194                _ => Err(error),
195            };
196        }
197        Err(error) => return Err(error),
198    };
199    // A bound that overflowed to infinity or NaN counts as lost.
200    let exceeds = |bound: f64, tolerance: f64| bound.is_nan() || f64::EPSILON * bound > tolerance;
201    let tighter =
202        |series: f64, ewald: f64| series.is_finite() && (ewald.is_nan() || series < ewald);
203    let scale = sum.value.norm().max(1.0);
204    let value_lost = exceeds(rounding.value, SPECTRAL_SWITCH_LOSS * scale);
205    let lost: [bool; N] = std::array::from_fn(|i| {
206        exceeds(
207            rounding.derivative[i],
208            SPECTRAL_SWITCH_LOSS
209                * scale.mul_add(derivative_scale(i, 1, k), sum.derivative[i].norm()),
210        )
211    });
212    if (value_lost || lost.contains(&true))
213        && let Some((series, bound)) = series()
214    {
215        if value_lost && tighter(bound.value, rounding.value) {
216            sum.value = series.value;
217            rounding.value = bound.value;
218        }
219        for (i, lost) in lost.into_iter().enumerate() {
220            if lost && tighter(bound.derivative[i], rounding.derivative[i]) {
221                sum.derivative[i] = series.derivative[i];
222            }
223        }
224    }
225    if !accurate(sum.value, rounding.value) {
226        return Err(Error::NotConverged(
227            "Ewald sum lost its accuracy to cancelling Kambe integrals; reduce the split parameter"
228                .into(),
229        ));
230    }
231    Ok(sum)
232}
treams_core/lattice/spectral.rs

treams_core/lattice/
spectral.rs
