treams_core/sw/
cartesian.rs

1//! Spherical translation coefficients of single mode pairs at Cartesian displacements,
2//! with their displacement and wavenumber derivatives.
3//!
4//! treams-rs extension: treams evaluates translations at polar arguments only.
5#![allow(clippy::indexing_slicing)] // Fixed indices of the three Cartesian axes.
6
7use super::{Mode, coupling::terms};
8use crate::special::{Radial, direction, on_sphere, spherical_radial};
9use crate::{Complex, Error, Result, basis::ModeLabel};
10
11/// Translation coefficient with its Cartesian and wavenumber derivatives.
12#[derive(Clone, Copy, Debug, Default)]
13pub struct CartesianTranslation {
14    /// Coefficient value.
15    pub value: Complex,
16    /// Derivatives of `value` with respect to the displacement components.
17    pub position: [Complex; 3],
18    /// Derivative of `value` with respect to the wavenumber `k`.
19    pub k: Complex,
20}
21
22/// Translation coefficient of the source mode `from` to the destination mode `to`, with
23/// its derivatives.
24///
25/// `displacement` is the destination position minus the source position. Singular waves
26/// give zero at the origin and at subnormal displacements.
27///
28/// treams-rs extension: the Cartesian single-pair jet behind the expansion matrices and
29/// the tests of the polar coefficients; `treams.sw.translate` takes polar arguments.
30pub fn cartesian_translation(
31    to: Mode,
32    from: Mode,
33    k: Complex,
34    displacement: [f64; 3],
35    helicity: bool,
36    radial: Radial,
37) -> Result<CartesianTranslation> {
38    to.validate()?;
39    from.validate()?;
40    if !crate::numerics::finite(k) || displacement.iter().any(|v| !v.is_finite()) {
41        return Err(Error::InvalidInput(
42            "wavenumber and displacement must be finite".into(),
43        ));
44    }
45    let mut result = CartesianTranslation::default();
46    let geometry = direction(displacement);
47    if geometry.is_none() && radial == Radial::Singular {
48        return Ok(result);
49    }
50    let r = geometry.map_or(0.0, |(r, _)| r);
51    for (p, m, weight) in terms(to, from, helicity) {
52        let radial_value = spherical_radial(
53            u32::try_from(p).map_err(|_| Error::InvalidInput("negative radial order".into()))?,
54            k * r,
55            radial,
56        )?;
57        let term = harmonic(p, m, k, geometry, radial_value);
58        result.value += weight * term.value;
59        result.k += weight * term.k;
60        for (g, derivative) in result.position.iter_mut().zip(term.position) {
61            *g += weight * derivative;
62        }
63    }
64    Ok(result)
65}
66
67/// `z_p(k r) P_p^m(cos theta) exp(i m phi)` and its displacement and wavenumber
68/// derivatives, from the radial jet at `k r` and the displacement's
69/// [`direction`] (`None` at the origin and at subnormal displacements).
70pub(crate) fn harmonic(
71    p: i32,
72    m: i32,
73    k: Complex,
74    geometry: Option<(f64, [f64; 3])>,
75    radial: crate::special::RadialJet,
76) -> CartesianTranslation {
77    let Some((r, unit)) = geometry else {
78        let mut result = CartesianTranslation::default();
79        if p == 0 {
80            result.value = radial.value;
81        }
82        // j_1(kr) P_1^m(z/r) exp(im phi) is linear at the origin.
83        if p == 1 {
84            let gradient = match m {
85                -1 => [
86                    Complex::new(0.5, 0.0),
87                    Complex::new(0.0, -0.5),
88                    Complex::default(),
89                ],
90                0 => [
91                    Complex::default(),
92                    Complex::default(),
93                    Complex::new(1.0, 0.0),
94                ],
95                1 => [
96                    Complex::new(-1.0, 0.0),
97                    Complex::new(0.0, -1.0),
98                    Complex::default(),
99                ],
100                _ => [Complex::default(); 3],
101            };
102            result.position = gradient.map(|g| k * radial.first * g);
103        }
104        return result;
105    };
106    let (value, tangent) = on_sphere(p, m, unit);
107    combine(k, (r, unit), radial, value, tangent)
108}
109
110/// `z_p(k r) Y` and its displacement and wavenumber derivatives away from the origin,
111/// from the unit-sphere harmonic `Y` at the direction `unit` and its tangential
112/// gradient (see [`crate::special::tangent`]).
113pub(crate) fn combine(
114    k: Complex,
115    (r, unit): (f64, [f64; 3]),
116    radial: crate::special::RadialJet,
117    value: Complex,
118    tangent: [Complex; 3],
119) -> CartesianTranslation {
120    let scaled = radial.value / r;
121    CartesianTranslation {
122        value: radial.value * value,
123        k: radial.first * r * value,
124        position: std::array::from_fn(|axis| {
125            radial.first * k * unit[axis] * value + scaled * tangent[axis]
126        }),
127    }
128}
treams_core/sw/cartesian.rs

treams_core/sw/
cartesian.rs
