laddu_core/utils/
functions.rs

1use factorial::Factorial;
2use nalgebra::ComplexField;
3use num::{complex::Complex, Integer};
4
5#[cfg(feature = "f32")]
6use std::f32::consts::PI;
7#[cfg(not(feature = "f32"))]
8use std::f64::consts::PI;
9
10use crate::Float;
11
12fn alp_pos_m(l: usize, m: usize, x: Float) -> Float {
13    let mut p = 1.0;
14    if l == 0 && m == 0 {
15        return p;
16    }
17    let y = Float::sqrt(1.0 - Float::powi(x, 2));
18    for m_p in 0..m {
19        p *= -((2 * m_p + 1) as Float) * y;
20    }
21    if l == m {
22        return p;
23    }
24    let mut p_min_2 = p;
25    let mut p_min_1 = (2 * m + 1) as Float * x * p_min_2;
26    if l == m + 1 {
27        return p_min_1;
28    }
29    for l_p in (m + 1)..l {
30        p = ((2 * l_p + 1) as Float * x * p_min_1 - (l_p + m) as Float * p_min_2)
31            / (l_p - m + 1) as Float;
32        p_min_2 = p_min_1;
33        p_min_1 = p;
34    }
35    p
36}
37
38/// Computes the spherical harmonic $`Y_{\ell}^m(\theta, \phi)`$ (given $`\cos\theta`$). Note that
39/// this formulation includes the Condon-Shortley phase.
40pub fn spherical_harmonic(l: usize, m: isize, costheta: Float, phi: Float) -> Complex<Float> {
41    let abs_m = isize::abs(m) as usize;
42    let mut res = alp_pos_m(l, abs_m, costheta); // Includes Condon-Shortley phase already
43    res *= Float::sqrt(
44        (2 * l + 1) as Float / (4.0 * PI) * ((l - abs_m).factorial()) as Float
45            / ((l + abs_m).factorial()) as Float,
46    );
47    if m < 0 {
48        res *= if abs_m.is_even() { 1.0 } else { -1.0 }; // divide out Condon-Shortley phase
49                                                         // (it's just +/-1 so division is the
50                                                         // same as multiplication here)
51    }
52    Complex::new(
53        res * Float::cos(m as Float * phi),
54        res * Float::sin(m as Float * phi),
55    )
56}
57
58/// Computes $`\chi_+(s, m_1, m_2) = 1 - \frac{(m_1 + m_2)^2}{s}`$.
59pub fn chi_plus(s: Float, m1: Float, m2: Float) -> Float {
60    1.0 - (m1 + m2) * (m1 + m2) / s
61}
62
63/// Computes $`\chi_-(s, m_1, m_2) = 1 - \frac{(m_1 - m_2)^2}{s}`$.
64pub fn chi_minus(s: Float, m1: Float, m2: Float) -> Float {
65    1.0 - (m1 - m2) * (m1 - m2) / s
66}
67
68/// Computes the phase-space factor $`\rho(s, m_1, m_2) = \sqrt(\chi_+(s, m_1, m_2)\chi_-(s, m_1, m_2))`$
69pub fn rho(s: Float, m1: Float, m2: Float) -> Complex<Float> {
70    let x: Complex<Float> = (chi_plus(s, m1, m2) * chi_minus(s, m1, m2)).into();
71    x.sqrt()
72}
73
74/// Computes the breakup momentum (often denoted $`q`$) for a particle with mass $`m_0`$ decaying
75/// into two particles with masses $`m_1`$ and $`m_2`$: $`\frac{m_0 \left|\rho(m_0^2, m_1, m_2)\right|}{2}`$.
76///
77/// Note that this version will always yield a real-valued number.
78pub fn breakup_momentum(m0: Float, m1: Float, m2: Float) -> Float {
79    let s = m0.powi(2);
80    let x: Complex<Float> = (chi_plus(s, m1, m2) * chi_minus(s, m1, m2)).into();
81    x.abs().sqrt() * m0 / 2.0
82}
83
84/// Computes the breakup momentum (often denoted $`q`$) for a particle with mass $`m_0`$ decaying
85/// into two particles with masses $`m_1`$ and $`m_2`$: $`\frac{m_0 \rho(m_0^2, m_1, m_2)}{2}`$.
86pub fn complex_breakup_momentum(m0: Float, m1: Float, m2: Float) -> Complex<Float> {
87    rho(Float::powi(m0, 2), m1, m2) * m0 / 2.0
88}
89
90/// Computes the Blatt-Weisskopf centrifugal barrier factor for a particle with mass $`m_0`$ and
91/// angular momentum $`\ell`$ decaying to two particles with masses $`m_1`$ and $`m_2`$.
92///
93/// Note that this version uses absolute form of the breakup momentum, so the results are
94/// real-valued. Additionally, this method uses an impact parameter of $`0.1973\text{GeV}^{-1}`$.
95/// Currently, only values of $`\ell <= 4`$ are implemented.
96pub fn blatt_weisskopf(m0: Float, m1: Float, m2: Float, l: usize) -> Float {
97    let q = breakup_momentum(m0, m1, m2);
98    let z = q.powi(2) / 0.1973.powi(2);
99    match l {
100        0 => 1.0,
101        1 => ((2.0 * z) / (z + 1.0)).sqrt(),
102        2 => ((13.0 * z.powi(2)) / ((z - 3.0).powi(2) + 9.0 * z)).sqrt(),
103        3 => {
104            ((277.0 * z.powi(3)) / (z * (z - 15.0).powi(2) + 9.0 * (2.0 * z - 5.0).powi(2))).sqrt()
105        }
106        4 => ((12746.0 * z.powi(4))
107            / ((z.powi(2) - 45.0 * z + 105.0).powi(2) + 25.0 * z * (2.0 * z - 21.0).powi(2)))
108        .sqrt(),
109        l => panic!("L = {l} is not yet implemented"),
110    }
111}
112
113/// Computes the Blatt-Weisskopf centrifugal barrier factor for a particle with mass $`m_0`$ and
114/// angular momentum $`\ell`$ decaying to two particles with masses $`m_1`$ and $`m_2`$.
115///
116/// Note that this method uses an impact parameter of $`0.1973\text{GeV}^{-1}`$. Currently, only
117/// values of $`\ell <= 4`$ are implemented.
118pub fn complex_blatt_weisskopf(m0: Float, m1: Float, m2: Float, l: usize) -> Complex<Float> {
119    let q = complex_breakup_momentum(m0, m1, m2);
120    let z = q.powi(2) / 0.1973.powi(2);
121    match l {
122        0 => Complex::ONE,
123        1 => ((z * 2.0) / (z + 1.0)).sqrt(),
124        2 => ((z.powi(2) * 13.0) / ((z - 3.0).powi(2) + z * 9.0)).sqrt(),
125        3 => {
126            ((z.powi(3) * 277.0) / (z * (z - 15.0).powi(2) + (z * 2.0 - 5.0).powi(2) * 9.0)).sqrt()
127        }
128        4 => ((z.powi(4) * 12746.0)
129            / ((z.powi(2) - z * 45.0 + 105.0).powi(2) + z * 25.0 * (z * 2.0 - 21.0).powi(2)))
130        .sqrt(),
131        l => panic!("L = {l} is not yet implemented"),
132    }
133}
134
135#[cfg(test)]
136mod test {
137    use approx::assert_relative_eq;
138    use num::complex::Complex;
139
140    use crate::{
141        utils::functions::{
142            blatt_weisskopf, breakup_momentum, chi_minus, complex_breakup_momentum, rho,
143            spherical_harmonic,
144        },
145        Float,
146    };
147
148    use super::{chi_plus, complex_blatt_weisskopf};
149
150    #[test]
151    fn test_spherical_harmonics() {
152        #[cfg(feature = "f32")]
153        use std::f32::consts::PI;
154        #[cfg(not(feature = "f32"))]
155        use std::f64::consts::PI;
156        let costhetas = [-1.0, -0.8, -0.3, 0.0, 0.3, 0.8, 1.0];
157        let phis = [0.0, 0.3, 0.5, 0.8, 1.0].map(|v| v * PI * 2.0);
158        for costheta in costhetas {
159            for phi in phis {
160                // L = 0
161                let y00 = spherical_harmonic(0, 0, costheta, phi);
162                let y00_true = Complex::from(Float::sqrt(1.0 / (4.0 * PI)));
163                assert_relative_eq!(y00.re, y00_true.re);
164                assert_relative_eq!(y00.im, y00_true.im);
165                // L = 1
166                let y1n1 = spherical_harmonic(1, -1, costheta, phi);
167                let y1n1_true = Complex::from_polar(
168                    Float::sqrt(3.0 / (8.0 * PI)) * Float::sin(Float::acos(costheta)),
169                    -phi,
170                );
171                assert_relative_eq!(y1n1.re, y1n1_true.re);
172                assert_relative_eq!(y1n1.im, y1n1_true.im);
173                let y10 = spherical_harmonic(1, 0, costheta, phi);
174                let y10_true = Complex::from(Float::sqrt(3.0 / (4.0 * PI)) * costheta);
175                assert_relative_eq!(y10.re, y10_true.re);
176                assert_relative_eq!(y10.im, y10_true.im);
177                let y11 = spherical_harmonic(1, 1, costheta, phi);
178                let y11_true = Complex::from_polar(
179                    -Float::sqrt(3.0 / (8.0 * PI)) * Float::sin(Float::acos(costheta)),
180                    phi,
181                );
182                assert_relative_eq!(y11.re, y11_true.re);
183                assert_relative_eq!(y11.im, y11_true.im);
184                // L = 2
185                let y2n2 = spherical_harmonic(2, -2, costheta, phi);
186                let y2n2_true = Complex::from_polar(
187                    Float::sqrt(15.0 / (32.0 * PI)) * Float::sin(Float::acos(costheta)).powi(2),
188                    -2.0 * phi,
189                );
190                assert_relative_eq!(y2n2.re, y2n2_true.re);
191                assert_relative_eq!(y2n2.im, y2n2_true.im);
192                let y2n1 = spherical_harmonic(2, -1, costheta, phi);
193                let y2n1_true = Complex::from_polar(
194                    Float::sqrt(15.0 / (8.0 * PI)) * Float::sin(Float::acos(costheta)) * costheta,
195                    -phi,
196                );
197                assert_relative_eq!(y2n1.re, y2n1_true.re);
198                assert_relative_eq!(y2n1.im, y2n1_true.im);
199                let y20 = spherical_harmonic(2, 0, costheta, phi);
200                let y20_true =
201                    Complex::from(Float::sqrt(5.0 / (16.0 * PI)) * (3.0 * costheta.powi(2) - 1.0));
202                assert_relative_eq!(y20.re, y20_true.re);
203                assert_relative_eq!(y20.im, y20_true.im);
204                let y21 = spherical_harmonic(2, 1, costheta, phi);
205                let y21_true = Complex::from_polar(
206                    -Float::sqrt(15.0 / (8.0 * PI)) * Float::sin(Float::acos(costheta)) * costheta,
207                    phi,
208                );
209                assert_relative_eq!(y21.re, y21_true.re);
210                assert_relative_eq!(y21.im, y21_true.im);
211                let y22 = spherical_harmonic(2, 2, costheta, phi);
212                let y22_true = Complex::from_polar(
213                    Float::sqrt(15.0 / (32.0 * PI)) * Float::sin(Float::acos(costheta)).powi(2),
214                    2.0 * phi,
215                );
216                assert_relative_eq!(y22.re, y22_true.re);
217                assert_relative_eq!(y22.im, y22_true.im);
218            }
219        }
220    }
221
222    #[test]
223    fn test_momentum_functions() {
224        assert_relative_eq!(
225            chi_plus(1.3, 0.51, 0.62),
226            0.01776923076923098,
227            epsilon = Float::EPSILON.sqrt()
228        );
229        assert_relative_eq!(
230            chi_minus(1.3, 0.51, 0.62),
231            0.9906923076923076,
232            epsilon = Float::EPSILON.sqrt()
233        );
234        let x0 = rho(1.3, 0.51, 0.62);
235        assert_relative_eq!(x0.re, 0.1326794642613792, epsilon = Float::EPSILON.sqrt());
236        assert_relative_eq!(x0.im, 0.0);
237        let x1 = rho(1.3, 1.23, 0.62);
238        assert_relative_eq!(x1.re, 0.0);
239        assert_relative_eq!(x1.im, 1.0795209736472833, epsilon = Float::EPSILON.sqrt());
240        let y0 = breakup_momentum(1.2, 0.4, 0.5);
241        assert_relative_eq!(y0, 0.3954823004889093, epsilon = Float::EPSILON.sqrt());
242        let y1 = breakup_momentum(1.2, 1.4, 1.5);
243        assert_relative_eq!(y1, 1.3154464282347478, epsilon = Float::EPSILON.sqrt());
244        let y2 = complex_breakup_momentum(1.2, 0.4, 0.5);
245        assert_relative_eq!(y2.re, 0.3954823004889093, epsilon = Float::EPSILON.sqrt());
246        assert_relative_eq!(y2.im, 0.0);
247        let y3 = complex_breakup_momentum(1.2, 1.4, 1.5);
248        assert_relative_eq!(y3.re, 0.0);
249        assert_relative_eq!(y3.im, 1.3154464282347478, epsilon = Float::EPSILON.sqrt());
250        assert_relative_eq!(y0, y2.re);
251        assert_relative_eq!(y1, y3.im);
252
253        let z0 = blatt_weisskopf(1.2, 0.4, 0.5, 0);
254        assert_relative_eq!(z0, 1.0);
255        let z1 = blatt_weisskopf(1.2, 0.4, 0.5, 1);
256        assert_relative_eq!(z1, 1.2654752018685698, epsilon = Float::EPSILON.sqrt());
257        let z2 = blatt_weisskopf(1.2, 0.4, 0.5, 2);
258        assert_relative_eq!(z2, 2.375285855793918, epsilon = Float::EPSILON.sqrt());
259        let z3 = blatt_weisskopf(1.2, 0.4, 0.5, 3);
260        assert_relative_eq!(z3, 5.6265876867850695, epsilon = Float::EPSILON.sqrt());
261        let z4 = blatt_weisskopf(1.2, 0.4, 0.5, 4);
262        assert_relative_eq!(z4, 12.747554064467208, epsilon = Float::EPSILON.sqrt());
263        let z0im = blatt_weisskopf(1.2, 1.4, 0.5, 0);
264        assert_relative_eq!(z0im, 1.0);
265        let z1im = blatt_weisskopf(1.2, 1.4, 1.5, 1);
266        assert_relative_eq!(z1im, 1.398569848337654, epsilon = Float::EPSILON.sqrt());
267        let z2im = blatt_weisskopf(1.2, 1.4, 1.5, 2);
268        assert_relative_eq!(z2im, 3.482294988252171, epsilon = Float::EPSILON.sqrt());
269        let z3im = blatt_weisskopf(1.2, 1.4, 1.5, 3);
270        assert_relative_eq!(z3im, 15.450855647831101, epsilon = Float::EPSILON.sqrt());
271        let z4im = blatt_weisskopf(1.2, 1.4, 1.5, 4);
272        assert_relative_eq!(z4im, 98.48799450927207, epsilon = Float::EPSILON.sqrt());
273
274        let w0 = complex_blatt_weisskopf(1.2, 0.4, 0.5, 0);
275        assert_relative_eq!(w0.re, 1.0);
276        assert_relative_eq!(w0.im, 0.0);
277        let w1 = complex_blatt_weisskopf(1.2, 0.4, 0.5, 1);
278        assert_relative_eq!(w1.re, 1.2654752018685698, epsilon = Float::EPSILON.sqrt());
279        assert_relative_eq!(w1.im, 0.0);
280        let w2 = complex_blatt_weisskopf(1.2, 0.4, 0.5, 2);
281        assert_relative_eq!(w2.re, 2.375285855793918, epsilon = Float::EPSILON.sqrt());
282        assert_relative_eq!(w2.im, 0.0);
283        let w3 = complex_blatt_weisskopf(1.2, 0.4, 0.5, 3);
284        assert_relative_eq!(w3.re, 5.62658768678507, epsilon = Float::EPSILON.sqrt());
285        assert_relative_eq!(w3.im, 0.0, epsilon = Float::EPSILON.sqrt());
286        let w4 = complex_blatt_weisskopf(1.2, 0.4, 0.5, 4);
287        assert_relative_eq!(w4.re, 12.747554064467208, epsilon = Float::EPSILON.sqrt());
288        assert_relative_eq!(w4.im, 0.0, epsilon = Float::EPSILON.sqrt());
289        let w0im = complex_blatt_weisskopf(1.2, 1.4, 1.5, 0);
290        assert_relative_eq!(w0im.re, 1.0);
291        assert_relative_eq!(w0im.im, 0.0);
292        let w1im = complex_blatt_weisskopf(1.2, 1.4, 1.5, 1);
293        assert_relative_eq!(w1im.re, 1.430394249144933, epsilon = Float::EPSILON.sqrt());
294        assert_relative_eq!(w1im.im, 0.0);
295        let w2im = complex_blatt_weisskopf(1.2, 1.4, 1.5, 2);
296        assert_relative_eq!(w2im.re, 3.724659004227952, epsilon = Float::EPSILON.sqrt());
297        assert_relative_eq!(w2im.im, 0.0, epsilon = Float::EPSILON.sqrt());
298        let w3im = complex_blatt_weisskopf(1.2, 1.4, 1.5, 3);
299        assert_relative_eq!(w3im.re, 17.689297320491015, epsilon = Float::EPSILON.sqrt());
300        assert_relative_eq!(w3im.im, 0.0, epsilon = Float::EPSILON.sqrt());
301        let w4im = complex_blatt_weisskopf(1.2, 1.4, 1.5, 4);
302        assert_relative_eq!(w4im.re, 124.0525841825899, epsilon = Float::EPSILON.sqrt());
303        assert_relative_eq!(w4im.im, 0.0, epsilon = Float::EPSILON.sqrt());
304
305        assert_relative_eq!(z0, w0.re);
306        assert_relative_eq!(z1, w1.re);
307        assert_relative_eq!(z2, w2.re);
308        assert_relative_eq!(z3, w3.re);
309        assert_relative_eq!(z4, w4.re);
310    }
311    #[test]
312    #[should_panic]
313    fn panicking_blatt_weisskopf() {
314        blatt_weisskopf(1.2, 0.4, 0.5, 5);
315    }
316    #[test]
317    #[should_panic]
318    fn panicking_complex_blatt_weisskopf() {
319        complex_blatt_weisskopf(1.2, 0.4, 0.5, 5);
320    }
321}
laddu_core/utils/functions.rs

laddu_core/utils/
functions.rs
