treams_core/special/
mod.rs

1//! Special functions: Bessel, Legendre, Ferrers and Wigner functions, solid harmonics,
2//! incomplete gamma and Kambe integrals, and coordinate transforms.
3//!
4//! Upstream: `treams.special`.
5//!
6//! | File | Contents | treams source |
7//! |---|---|---|
8//! | `bessel` | Cylindrical and spherical Bessel functions and radial jets | `jv`, `yv`, `hankel1`, `hankel2` and `spherical_*` (scipy), `_bessel.pyx` |
9//! | `legendre` | Associated Legendre functions and the angular functions `pi` and `tau` | `lpmv`, `pi_fun`, `tau_fun` of `_waves.pyx` |
10//! | `ferrers` | Ferrers functions of real, non-integer degree | `treams.special.lpmv` |
11//! | `wigner` | Wigner 3j symbols and Wigner small-d and D functions | `_wigner3j.pyx`, `_wignerd.pyx` |
12//! | `harmonics` | Cartesian solid harmonics and the spherical-harmonic normalization | none: a treams-rs addition for translations and lattice sums |
13//! | [`coordinates`] | Point and vector-component transforms between charts | `_coord.pyx` |
14//! | `integrals` | Upper incomplete gamma function and Kambe integrals | `_integrals.pyx` |
15//! | `polarization` | Helicity and parity combinations of vector waves | the `vsw_A` convention of `_waves.pyx` |
16//!
17//! ## Label bounds
18//!
19//! treams accepts any degree. Here the multipole degree, the integer-order radial jets
20//! and the Kambe and Wigner labels have fixed bounds, so the caches and tables that
21//! serve them have a fixed size:
22//!
23//! - [`MAX_DEGREE`] = 128 bounds the multipole degree `l`.
24//! - `MAX_ORDER` = 256 = `2 MAX_DEGREE` bounds the integer-order radial jets and the
25//!   log-factorial cache, since translations couple degrees `l1 + l2`.
26//! - [`MAX_LABEL`] = 260 bounds Kambe orders, and Wigner labels at the Python
27//!   interface. Lattice sums and their derivatives take Kambe orders `n` up to
28//!   `|n| = 2 MAX_DEGREE + 3` = 259, which the derivatives of 1D spherical sums reach;
29//!   the bound keeps one order of margin.
30//! - `SERIES_RADIUS` = 0.5 is the argument modulus below which regular radial functions
31//!   use their power series instead of AMOS evaluations.
32//!
33//! ## Two kinds of jet
34//!
35//! A jet is a value together with some of its derivatives. [`RadialJet`] holds a radial
36//! function and its first two derivatives in its one complex argument, as the vector
37//! waves need them. `numerics::Jet<N>` holds a value and its first derivatives along `N`
38//! inputs; kernels generic in `N` use it to evaluate (`N = 0`) and to differentiate with
39//! one code path. [`spherical_radial`] and [`cylindrical_radial`] return a `RadialJet`,
40//! and `radial_jet` turns one into a `Jet<N>` by the chain rule.
41
42mod bessel;
43pub mod coordinates;
44mod ferrers;
45mod harmonics;
46mod integrals;
47mod legendre;
48mod polarization;
49mod wigner;
50
51pub use bessel::{
52    Bessel, BesselResidual, Radial, RadialJet, bessel, bessel_array, cylindrical_radial,
53    spherical_radial,
54};
55pub(crate) use bessel::{hankel_below, radial_jet, spherical_hankels, spherical_radial_sequence};
56#[cfg(test)]
57pub(crate) use ferrers::ferrers_real_degree;
58#[cfg(test)]
59pub(crate) use ferrers::ferrers_route_switches;
60pub(crate) use harmonics::{
61    Solid, SolidTable, direction, harmonic_normalization, on_sphere, solid, tangent,
62};
63pub(crate) use integrals::{
64    EvenKambe, ScaledGammaLadder, SmallSplitKambe, kambe, kambe_argument, kambe_with_bound,
65    upper_gamma,
66};
67pub use integrals::{
68    IncgammaResidual, IntkambeResidual, incgamma, incgamma_array, intkambe, intkambe_array,
69};
70pub use legendre::{Angular, AngularResidual, angular_array, angular_value, lpmv_real};
71pub(crate) use legendre::{legendre_factor, legendre_factors, polar_trig, polarized_angular};
72pub use polarization::pol_index;
73pub(crate) use polarization::{check_pol, helicity_sign, polarized_wave};
74pub(crate) use wigner::{Wigner3jRow, index, ladder, wigner_small_d_matrix};
75pub use wigner::{WignerDResidual, wigner_d, wigner_d_array, wigner_small_d, wigner3j};
76
77use crate::MAX_DEGREE;
78
79/// The largest order of the integer-order radial jets ([`spherical_radial`],
80/// `spherical_radial_sequence`, [`cylindrical_radial`]) and of the order difference of
81/// cylindrical translations: translations couple degrees `l1 + l2` up to `2 MAX_DEGREE`,
82/// and cylindrical translations reach orders `m1 - m2` of the same size. [`bessel()`]
83/// takes any finite order.
84pub(crate) const MAX_ORDER: i32 = 2 * MAX_DEGREE;
85
86/// The largest `|n|` of an integer special-function label (Kambe orders, Wigner labels).
87///
88/// Lattice sums of degree `l` up to [`MAX_DEGREE`] use Kambe orders `n` with `|n|` up
89/// to `2 l + 1`, and their derivatives up to `2 l + 3`; the reciprocal terms of 1D
90/// spherical sums reach both. Translations use Wigner labels up to
91/// `l1 + l2 = 2 MAX_DEGREE`. The bound adds one order of margin.
92pub const MAX_LABEL: i32 = 2 * MAX_DEGREE + 4;
93
94/// The radius `|z|` below which regular radial functions use their power series instead
95/// of AMOS Bessel evaluations.
96///
97/// Kernels that expand small arguments themselves (vector waves, polar translations,
98/// fields) switch at the same radius, so values and derivatives stay continuous there.
99pub(crate) const SERIES_RADIUS: f64 = 0.5;
100
101/// `ln n!`: libm's `lgamma(n + 1)`, cached for `0 <= n <= MAX_ORDER`. The factorials of
102/// wave labels recur in every normalization, lattice image and angular closed form, so
103/// the cache computes each one once.
104pub(crate) fn log_factorial(n: i32) -> f64 {
105    const CACHED: usize = MAX_ORDER.unsigned_abs() as usize + 1;
106    static LOG_FACTORIAL: std::sync::OnceLock<[f64; CACHED]> = std::sync::OnceLock::new();
107    let values = LOG_FACTORIAL.get_or_init(|| {
108        std::array::from_fn(|i| libm::lgamma(f64::from(u32::try_from(i).unwrap_or_default()) + 1.0))
109    });
110    usize::try_from(n)
111        .ok()
112        .and_then(|i| values.get(i))
113        .copied()
114        .unwrap_or_else(|| libm::lgamma(f64::from(n) + 1.0))
115}
116
117/// `sqrt((l - m)! / (l + m)!)`, the order factor of spherical-harmonic normalizations.
118pub(crate) fn factorial_ratio_sqrt(l: i32, m: i32) -> f64 {
119    (0.5 * (log_factorial(l - m) - log_factorial(l + m))).exp()
120}
treams_core/special/mod.rs

treams_core/special/
mod.rs
