treams_core/special/polarization.rs
1//! Polarization conventions of vector waves: how the polarization index `pol` combines
2//! the two multipole parts M and N of a wave.
3//!
4//! - Parity: `pol = 0` is M and `pol = 1` is N.
5//! - Helicity: `(N + (2 pol - 1) M) / sqrt(2)`, so `pol = 0` is negative and `pol = 1`
6//! positive helicity.
7//!
8//! Upstream: the convention of `treams.special.vsw_A` and `vcw_A`. [`polarized_wave`]
9//! applies it to vector waves, and `legendre::polarized_angular` applies the same sign
10//! `2 pol - 1` ([`helicity_sign`]) to the angular factors `tau + (2 pol - 1) pi` of
11//! plane-wave expansions.
12#![allow(clippy::indexing_slicing)] // Fixed Cartesian component arrays.
13
14use std::f64::consts::FRAC_1_SQRT_2;
15
16use crate::{Error, Result, numerics::Jet};
17
18/// The sign `2 pol - 1` of polarization index `pol` in the helicity convention: -1 for
19/// negative helicity (0), +1 for positive helicity (1).
20#[inline]
21pub(crate) fn helicity_sign(pol: u8) -> f64 {
22 2.0 * f64::from(pol) - 1.0
23}
24
25/// The polarization index of an integer label, which must be 0 or 1.
26#[inline]
27pub fn pol_index(label: i64) -> Result<u8> {
28 match label {
29 0 => Ok(0),
30 1 => Ok(1),
31 _ => Err(Error::InvalidInput("polarization must be 0 or 1".into())),
32 }
33}
34
35/// Checks that `pol` is a polarization index: 0 or 1.
36#[inline]
37pub(crate) fn check_pol(pol: u8) -> Result<()> {
38 if pol > 1 {
39 return Err(Error::InvalidInput("polarization must be 0 or 1".into()));
40 }
41 Ok(())
42}
43
44/// The vector wave of polarization `pol` from its M and N parts: M (0) or N (1) in
45/// parity polarization, `(N -+ M) / sqrt(2)` for negative (0) or positive (1) helicity.
46pub(crate) fn polarized_wave<const N: usize>(
47 m: [Jet<N>; 3],
48 n: [Jet<N>; 3],
49 pol: u8,
50 helicity: bool,
51) -> [Jet<N>; 3] {
52 if helicity {
53 std::array::from_fn(|i| (n[i] + helicity_sign(pol) * m[i]) * FRAC_1_SQRT_2)
54 } else if pol == 0 {
55 m
56 } else {
57 n
58 }
59}