num_dual/
bessel.rs

1#![allow(clippy::excessive_precision)]
2use crate::DualNum;
3use num_traits::{Float, Zero};
4use std::f64::consts::{FRAC_2_PI, FRAC_PI_4};
5
6/// Implementation of bessel functions for double precision (hyper) dual numbers.
7pub trait BesselDual: DualNum<f64> + Copy {
8    /// 0th order bessel function of the first kind
9    fn bessel_j0(mut self) -> Self {
10        if self.is_negative() {
11            self = -self;
12        }
13
14        if self.re() <= 5.0 {
15            let z = self * self;
16            if self.re() < 1.0e-5 {
17                return Self::one() - z / 4.0;
18            }
19
20            (z - DR1) * (z - DR2) * polevl(z, &RP0) / p1evl(z, &RQ0)
21        } else {
22            let w = self.recip() * 5.0;
23            let q = w * w;
24            let p = polevl(q, &PP0) / polevl(q, &PQ0);
25            let q = polevl(q, &QP0) / p1evl(q, &QQ0);
26            let (s, c) = (self - FRAC_PI_4).sin_cos();
27            let p = p * c - w * q * s;
28            p * (Self::from(FRAC_2_PI) / self).sqrt()
29        }
30    }
31
32    /// 1st order bessel function of the first kind
33    fn bessel_j1(self) -> Self {
34        let x = self.abs();
35        if x.re() <= 5.0 {
36            {
37                let z = self * self;
38                polevl(z, &RP1) / p1evl(z, &RQ1) * self * (z - Z1) * (z - Z2)
39            }
40        } else {
41            let x = self.abs();
42            let w = x.recip() * 5.0;
43            let z = w * w;
44            let p = polevl(z, &PP1) / polevl(z, &PQ1);
45            let q = polevl(z, &QP1) / p1evl(z, &QQ1);
46            let (s, c) = (x - 3.0 * FRAC_PI_4).sin_cos();
47            let p = p * c - w * q * s;
48            self.signum() * p * (Self::from(FRAC_2_PI) / x).sqrt()
49        }
50    }
51
52    /// 2nd order bessel function of the first kind
53    fn bessel_j2(self) -> Self {
54        if self.re().is_zero() {
55            self * self / 8.0 * (self * self / 24.0 + 1.0)
56        } else {
57            self.bessel_j1() * 2.0 / self - self.bessel_j0()
58        }
59    }
60}
61
62impl<T: DualNum<f64> + Copy> BesselDual for T {}
63
64const DR1: f64 = 5.78318596294678452118E0;
65const DR2: f64 = 3.04712623436620863991E1;
66
67const RP0: [f64; 4] = [
68    -4.79443220978201773821E9,
69    1.95617491946556577543E12,
70    -2.49248344360967716204E14,
71    9.70862251047306323952E15,
72];
73const RQ0: [f64; 8] = [
74    4.99563147152651017219E2,
75    1.73785401676374683123E5,
76    4.84409658339962045305E7,
77    1.11855537045356834862E10,
78    2.11277520115489217587E12,
79    3.10518229857422583814E14,
80    3.18121955943204943306E16,
81    1.71086294081043136091E18,
82];
83
84const PP0: [f64; 7] = [
85    7.96936729297347051624E-4,
86    8.28352392107440799803E-2,
87    1.23953371646414299388E0,
88    5.44725003058768775090E0,
89    8.74716500199817011941E0,
90    5.30324038235394892183E0,
91    9.99999999999999997821E-1,
92];
93const PQ0: [f64; 7] = [
94    9.24408810558863637013E-4,
95    8.56288474354474431428E-2,
96    1.25352743901058953537E0,
97    5.47097740330417105182E0,
98    8.76190883237069594232E0,
99    5.30605288235394617618E0,
100    1.00000000000000000218E0,
101];
102
103const QP0: [f64; 8] = [
104    -1.13663838898469149931E-2,
105    -1.28252718670509318512E0,
106    -1.95539544257735972385E1,
107    -9.32060152123768231369E1,
108    -1.77681167980488050595E2,
109    -1.47077505154951170175E2,
110    -5.14105326766599330220E1,
111    -6.05014350600728481186E0,
112];
113const QQ0: [f64; 7] = [
114    6.43178256118178023184E1,
115    8.56430025976980587198E2,
116    3.88240183605401609683E3,
117    7.24046774195652478189E3,
118    5.93072701187316984827E3,
119    2.06209331660327847417E3,
120    2.42005740240291393179E2,
121];
122
123const Z1: f64 = 1.46819706421238932572E1;
124const Z2: f64 = 4.92184563216946036703E1;
125
126const RP1: [f64; 4] = [
127    -8.99971225705559398224E8,
128    4.52228297998194034323E11,
129    -7.27494245221818276015E13,
130    3.68295732863852883286E15,
131];
132const RQ1: [f64; 8] = [
133    6.20836478118054335476E2,
134    2.56987256757748830383E5,
135    8.35146791431949253037E7,
136    2.21511595479792499675E10,
137    4.74914122079991414898E12,
138    7.84369607876235854894E14,
139    8.95222336184627338078E16,
140    5.32278620332680085395E18,
141];
142
143const PP1: [f64; 7] = [
144    7.62125616208173112003E-4,
145    7.31397056940917570436E-2,
146    1.12719608129684925192E0,
147    5.11207951146807644818E0,
148    8.42404590141772420927E0,
149    5.21451598682361504063E0,
150    1.00000000000000000254E0,
151];
152const PQ1: [f64; 7] = [
153    5.71323128072548699714E-4,
154    6.88455908754495404082E-2,
155    1.10514232634061696926E0,
156    5.07386386128601488557E0,
157    8.39985554327604159757E0,
158    5.20982848682361821619E0,
159    9.99999999999999997461E-1,
160];
161
162const QP1: [f64; 8] = [
163    5.10862594750176621635E-2,
164    4.98213872951233449420E0,
165    7.58238284132545283818E1,
166    3.66779609360150777800E2,
167    7.10856304998926107277E2,
168    5.97489612400613639965E2,
169    2.11688757100572135698E2,
170    2.52070205858023719784E1,
171];
172const QQ1: [f64; 7] = [
173    7.42373277035675149943E1,
174    1.05644886038262816351E3,
175    4.98641058337653607651E3,
176    9.56231892404756170795E3,
177    7.99704160447350683650E3,
178    2.82619278517639096600E3,
179    3.36093607810698293419E2,
180];
181
182fn polevl<T: DualNum<F> + Copy, F: Float>(x: T, coef: &[F]) -> T {
183    coef.iter()
184        .skip(1)
185        .fold(T::from(coef[0]), |acc, &c| acc * x + c)
186}
187
188fn p1evl<T: DualNum<F> + Copy, F: Float>(x: T, coef: &[F]) -> T {
189    coef.iter().fold(T::one(), |acc, &c| acc * x + c)
190}
num_dual/bessel.rs

num_dual/
bessel.rs
