laddu_core/utils/
functions.rs

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

laddu_core/utils/
functions.rs
