treams_core/cw/
cartesian.rs

1//! Cartesian cylindrical-wave translation and analytic derivatives.
2//!
3//! treams-rs extension: treams evaluates translations at polar arguments only.
4#![allow(clippy::indexing_slicing)] // Fixed-size Cartesian axes and jet slots.
5
6use super::{Mode, transverse, transverse_wavenumber};
7use crate::{
8    Complex, Error, Result,
9    basis::ModeLabel,
10    numerics::{Jet, finite},
11    special::{Radial, SERIES_RADIUS, cylindrical_radial},
12};
13
14/// Translation coefficient and its continuous derivatives, with the axial wavenumbers of
15/// both modes moving together.
16#[derive(Clone, Copy, Debug, Default)]
17pub struct CartesianTranslation {
18    /// Complex coefficient.
19    pub value: Complex,
20    /// Derivatives of `value` with respect to the displacement components.
21    pub position: [Complex; 3],
22    /// Derivative of `value` with respect to the medium wavenumber `k`.
23    pub k: Complex,
24    /// Derivative of `value` with respect to the common axial wavenumber `kz`.
25    pub kz: Complex,
26}
27
28/// Translation coefficient of the source mode `from` to the destination mode `to`, with
29/// its derivatives.
30///
31/// `position` is the destination position minus the source position. Distinct axial
32/// wavenumbers and polarizations decouple, and singular waves give zero at the origin.
33///
34/// treams-rs extension: the Cartesian single-pair jet behind the expansion matrices and
35/// the tests of the polar coefficients; `treams.cw.translate` takes polar arguments.
36#[allow(clippy::float_cmp)] // kz labels and the self origin are exact mode identities.
37pub fn cartesian_translation(
38    to: Mode,
39    from: Mode,
40    k: Complex,
41    position: [f64; 3],
42    radial: Radial,
43) -> Result<CartesianTranslation> {
44    to.validate()?;
45    from.validate()?;
46    if !finite(k) || k.norm_sqr() == 0.0 || position.iter().any(|v| !v.is_finite()) {
47        return Err(Error::InvalidInput(
48            "require a finite nonzero wavenumber and finite displacement".into(),
49        ));
50    }
51    if to.kz != from.kz
52        || to.pol != from.pol
53        || (radial == Radial::Singular && position == [0.0; 3])
54    {
55        return Ok(CartesianTranslation::default());
56    }
57    let order = from.m - to.m;
58    let kz = from.kz;
59    let [x, y, z] = position;
60    let rho = x.hypot(y);
61    let krho = transverse_wavenumber(k, kz);
62    if krho.norm_sqr() == 0.0 {
63        return Err(Error::InvalidInput(
64            "cylindrical cutoff requires a limiting formulation".into(),
65        ));
66    }
67    if radial == Radial::Regular && (krho * rho).norm() < SERIES_RADIUS {
68        let kz = Jet::<5>::variable(kz, 4);
69        let wave = regular_harmonic(
70            order,
71            transverse(Jet::variable(k, 3), kz),
72            std::array::from_fn(|i| Jet::variable(position[i], i)),
73            kz,
74        );
75        return Ok(CartesianTranslation {
76            value: wave.value,
77            position: std::array::from_fn(|i| wave.derivative[i]),
78            k: wave.derivative[3],
79            kz: wave.derivative[4],
80        });
81    }
82    // Here rho > 0: regular waves on the axis took the series above, and outgoing
83    // waves are singular there.
84    let phi = y.atan2(x);
85    let phase = (Complex::i() * (f64::from(order) * phi + kz * z)).exp();
86    let radial_jet = cylindrical_radial(order, krho * rho, radial)?;
87    let value = radial_jet.value * phase;
88    let dr = radial_jet.first * krho * phase;
89    let azimuthal = Complex::i() * f64::from(order) * (value / rho);
90    Ok(CartesianTranslation {
91        value,
92        position: [
93            dr * (x / rho) - azimuthal * (y / rho),
94            dr * (y / rho) + azimuthal * (x / rho),
95            Complex::i() * kz * value,
96        ],
97        k: radial_jet.first * rho * k / krho * phase,
98        kz: -radial_jet.first * rho * kz / krho * phase + Complex::i() * z * value,
99    })
100}
101
102/// Regular `J_m` harmonic as a Cartesian power series, for `|k_rho rho|` below
103/// [`SERIES_RADIUS`].
104/// Factoring the azimuthal polynomial keeps derivatives regular on the cylinder axis.
105pub(crate) fn regular_harmonic<const N: usize>(
106    order: i32,
107    transverse: Jet<N>,
108    r: [Jet<N>; 3],
109    kz: Jet<N>,
110) -> Jet<N> {
111    let m = order.abs();
112    let azimuth = r[0] + Complex::i() * if order < 0 { -r[1] } else { r[1] };
113    let sign = if order < 0 && m % 2 == 1 { -1.0 } else { 1.0 };
114    let step = -0.25 * transverse * transverse * (r[0] * r[0] + r[1] * r[1]);
115    // The radial series depends on one scalar: differentiate it once, then
116    // apply the Cartesian chain rule after summation instead of at every term.
117    let step_norm = step.derivative.iter().map(|g| g.norm()).fold(0.0, f64::max);
118    let mut term = Complex::new((-crate::special::log_factorial(m)).exp(), 0.0);
119    let mut slope_term = Complex::default();
120    let mut sum = term;
121    let mut slope = Complex::default();
122    for q in 1..32 {
123        let denominator = f64::from(q * (m + q));
124        if N > 0 {
125            slope_term = (slope_term * step.value + term) / denominator;
126            slope += slope_term;
127        }
128        term = term * step.value / denominator;
129        sum += term;
130        let term_norm = term.norm().max(slope_term.norm() * step_norm);
131        let sum_norm = sum.norm().max(slope.norm() * step_norm);
132        if term_norm <= f64::EPSILON * sum_norm {
133            break;
134        }
135    }
136    sign * (0.5 * transverse * azimuth).powi(m)
137        * step.chain(sum, slope)
138        * (Complex::i() * kz * r[2]).exp()
139}
140
141#[cfg(test)]
142mod tests {
143    //! The axis series and the Cartesian derivatives against Bessel values and order
144    //! recurrences. Identities of translations are in `properties/waves.rs`.
145
146    use super::{Complex, Jet, Mode, Radial, cartesian_translation, regular_harmonic};
147    use crate::test_support::{DEFAULT_CASES, complex, prop_assert_close, radial};
148    use proptest::{prelude::*, test_runner::TestCaseError};
149
150    proptest! {
151        // Above the default, so that each of the 25 orders appears about four times.
152        #![proptest_config(ProptestConfig::with_cases(96))]
153        #[test]
154        fn regular_harmonic_value_and_cartesian_derivatives(
155            order in -12_i32..13, x in -0.25_f64..0.25, y in -0.25_f64..0.25,
156            z in -1.0_f64..1.0, kz in -0.5_f64..0.5, loss in -0.2_f64..0.2,
157        ) {
158            check_regular_harmonic(order, [x, y, z], kz, loss)?;
159        }
160    }
161
162    proptest! {
163        #![proptest_config(ProptestConfig::with_cases(DEFAULT_CASES))]
164
165        #[test]
166        fn translation_derivative_recurrences(
167            order in -8_i32..=8,
168            rho in prop_oneof![Just(0.0), 0.01_f64..0.35, 0.45_f64..2.0],
169            phi in -3.1_f64..3.1,
170            z in -1.0_f64..1.0,
171            kz in -0.5_f64..0.5,
172            k in complex(0.1).prop_map(|k| k + 1.2),
173            radial in radial(),
174        ) {
175            check_derivative_recurrences(order, [rho, phi, z], kz, k, radial)?;
176        }
177    }
178
179    /// The Cartesian power series of `J_m` and its position, `k_rho` and `kz`
180    /// derivatives match Bessel values and their order recurrences, also on the axis.
181    fn check_regular_harmonic(
182        order: i32,
183        [x, y, z]: [f64; 3],
184        kz: f64,
185        loss: f64,
186    ) -> Result<(), TestCaseError> {
187        let transverse = Complex::new(1.1, loss);
188        for position in [[x, y, z], [1e-100, -1e-100, z], [0.0, 0.0, z]] {
189            let wave = regular_harmonic(
190                order,
191                Jet::<5>::variable(transverse, 3),
192                std::array::from_fn(|i| Jet::variable(position[i], i)),
193                Jet::variable(kz, 4),
194            );
195            let rho = position[0].hypot(position[1]);
196            let phase = |m: i32| {
197                (Complex::i() * (f64::from(m) * position[1].atan2(position[0]) + kz * z)).exp()
198            };
199            let reference = |m: i32| {
200                let sign = if m < 0 && m.abs() % 2 == 1 { -1.0 } else { 1.0 };
201                sign * complex_bessel::besselj(f64::from(m.abs()), transverse * rho).unwrap()
202                    * phase(m)
203            };
204            let [lower, value, upper] = [order - 1, order, order + 1].map(reference);
205            let dx = 0.5 * transverse * (lower - upper);
206            let dy = 0.5 * Complex::i() * transverse * (lower + upper);
207            let expected = [
208                dx,
209                dy,
210                Complex::i() * kz * value,
211                (position[0] * dx + position[1] * dy) / transverse,
212                Complex::i() * z * value,
213            ];
214            prop_assert_close!(wave.value, value, 2e-13 * (1.0 + value.norm()));
215            for (actual, expected) in wave.derivative.into_iter().zip(expected) {
216                prop_assert_close!(actual, expected, 2e-13 * (1.0 + expected.norm()));
217            }
218        }
219        Ok(())
220    }
221
222    /// Exact derivative identities of `T_n = Z_n(k_rho rho) exp(i n phi + i kz z)` in
223    /// terms of the neighbouring orders, on both sides of the regular series
224    /// threshold and (for regular waves) on the axis:
225    /// `∂x T_n = k_rho (T_{n-1} - T_{n+1}) / 2`, `∂y T_n = i k_rho (T_{n-1} + T_{n+1}) / 2`,
226    /// `∂z T_n = i kz T_n`, `∂k T_n = k (x ∂x + y ∂y) T_n / k_rho²` and
227    /// `∂kz T_n = i z T_n - kz (x ∂x + y ∂y) T_n / k_rho²`.
228    fn check_derivative_recurrences(
229        order: i32,
230        [rho, phi, z]: [f64; 3],
231        kz: f64,
232        k: Complex,
233        radial: Radial,
234    ) -> Result<(), TestCaseError> {
235        if rho == 0.0 && radial == Radial::Singular {
236            return Ok(());
237        }
238        let position = [rho * phi.cos(), rho * phi.sin(), z];
239        let at = |n: i32| {
240            let mode = |m| Mode { kz, m, pol: 1 };
241            cartesian_translation(mode(0), mode(n), k, position, radial).unwrap()
242        };
243        let [lower, jet, upper] = [order - 1, order, order + 1].map(at);
244        let mut krho = (k * k - kz * kz).sqrt();
245        if krho.im < 0.0 {
246            krho = -krho;
247        }
248        let dx = 0.5 * krho * (lower.value - upper.value);
249        let dy = 0.5 * Complex::i() * krho * (lower.value + upper.value);
250        // (x ∂x + y ∂y) T_n / k_rho², the k_rho derivative divided by k_rho.
251        let spectral = (position[0] * dx + position[1] * dy) / krho.powi(2);
252        let expected = [
253            dx,
254            dy,
255            Complex::i() * kz * jet.value,
256            k * spectral,
257            Complex::i() * z * jet.value - kz * spectral,
258        ];
259        let actual = [
260            jet.position[0],
261            jet.position[1],
262            jet.position[2],
263            jet.k,
264            jet.kz,
265        ];
266        let scale = [lower, jet, upper]
267            .iter()
268            .map(|t| t.value.norm())
269            .sum::<f64>();
270        prop_assert_close!(
271            actual,
272            expected,
273            1e-12 * (1.0 + scale * (1.0 + rho + z.abs()))
274        );
275        Ok(())
276    }
277}
treams_core/cw/cartesian.rs

treams_core/cw/
cartesian.rs
