PolyFEM
Loading...
Searching...
No Matches
Bessel.hpp
Go to the documentation of this file.
1#include <cmath>
2
3namespace polyfem
4{
5 namespace utils
6 {
7 namespace
8 {
9 template <typename T>
10 class AbsVal
11 {
12 public:
13 T operator()(const T x) const
14 {
15 if (x.getValue() >= 0)
16 return x;
17
18 return -x;
19 }
20 };
21
22 template <>
23 class AbsVal<double>
24 {
25 public:
26 double operator()(const double x) const
27 {
28 return fabs(x);
29 }
30 };
31 } // namespace
32
33 template <typename T>
34 static T bessj0(T x)
35 {
36 static const AbsVal<T> abs_val;
37 const T ax = abs_val(x);
38
39 if (ax < 8.0)
40 {
41 const T y = x * x;
42 const T ans1 = 57568490574.0 + y * (-13362590354.0 + y * (651619640.7 + y * (-11214424.18 + y * (77392.33017 + y * (-184.9052456)))));
43 const T ans2 = 57568490411.0 + y * (1029532985.0 + y * (9494680.718 + y * (59272.64853 + y * (267.8532712 + y * 1.0))));
44 return ans1 / ans2;
45 }
46 else
47 {
48 const T z = 8.0 / ax;
49 const T y = z * z;
50 const T xx = ax - 0.785398164;
51 const T ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4 + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
52 const T ans2 = -0.1562499995e-1 + y * (0.1430488765e-3 + y * (-0.6911147651e-5 + y * (0.7621095161e-6 - y * 0.934935152e-7)));
53 return sqrt(0.636619772 / ax) * (cos(xx) * ans1 - z * sin(xx) * ans2);
54 }
55 }
56
57 template <typename T>
58 static T bessj1(T x)
59 {
60 static const AbsVal<T> abs_val;
61 const T ax = abs_val(x);
62
63 if (ax < 8.0)
64 {
65 const T y = x * x;
66 const T ans1 = x * (72362614232.0 + y * (-7895059235.0 + y * (242396853.1 + y * (-2972611.439 + y * (15704.48260 + y * (-30.16036606))))));
67 const T ans2 = 144725228442.0 + y * (2300535178.0 + y * (18583304.74 + y * (99447.43394 + y * (376.9991397 + y * 1.0))));
68 return ans1 / ans2;
69 }
70 else
71 {
72 const T z = 8.0 / ax;
73 const T y = z * z;
74 const T xx = ax - 2.356194491;
75 const T ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4 + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
76 const T ans2 = 0.04687499995 + y * (-0.2002690873e-3 + y * (0.8449199096e-5 + y * (-0.88228987e-6 + y * 0.105787412e-6)));
77 T ans = sqrt(0.636619772 / ax) * (cos(xx) * ans1 - z * sin(xx) * ans2);
78 if (x < 0.0)
79 ans = -ans;
80
81 return ans;
82 }
83 }
84
85 template <typename T>
86 static T bessy0(T x)
87 {
88 if (x < 8.0)
89 {
90 const T y = x * x;
91 const T ans1 = -2957821389.0 + y * (7062834065.0 + y * (-512359803.6 + y * (10879881.29 + y * (-86327.92757 + y * 228.4622733))));
92 const T ans2 = 40076544269.0 + y * (745249964.8 + y * (7189466.438 + y * (47447.26470 + y * (226.1030244 + y * 1.0))));
93 return (ans1 / ans2) + 0.636619772 * bessj0(x) * log(x);
94 }
95 else
96 {
97 const T z = 8.0 / x;
98 const T y = z * z;
99 const T xx = x - 0.785398164;
100 const T ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4 + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
101 const T ans2 = -0.1562499995e-1 + y * (0.1430488765e-3 + y * (-0.6911147651e-5 + y * (0.7621095161e-6 + y * (-0.934945152e-7))));
102 return sqrt(0.636619772 / x) * (sin(xx) * ans1 + z * cos(xx) * ans2);
103 }
104 }
105
106 template <typename T>
107 static T bessy1(T x)
108 {
109 if (x < 8.0)
110 {
111 const T y = x * x;
112 const T ans1 = x * (-0.4900604943e13 + y * (0.1275274390e13 + y * (-0.5153438139e11 + y * (0.7349264551e9 + y * (-0.4237922726e7 + y * 0.8511937935e4)))));
113 const T ans2 = 0.2499580570e14 + y * (0.4244419664e12 + y * (0.3733650367e10 + y * (0.2245904002e8 + y * (0.1020426050e6 + y * (0.3549632885e3 + y)))));
114 return (ans1 / ans2) + 0.636619772 * (bessj1(x) * log(x) - 1.0 / x);
115 }
116 else
117 {
118 const T z = 8.0 / x;
119 const T y = z * z;
120 const T xx = x - 2.356194491;
121 const T ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4 + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
122 const T ans2 = 0.04687499995 + y * (-0.2002690873e-3 + y * (0.8449199096e-5 + y * (-0.88228987e-6 + y * 0.105787412e-6)));
123 return sqrt(0.636619772 / x) * (sin(xx) * ans1 + z * cos(xx) * ans2);
124 }
125 }
126
127 template <typename T>
128 static T bessi0(T x)
129 {
130 static const AbsVal<T> abs_val;
131 const T ax = abs_val(x);
132
133 if (ax < 3.75)
134 {
135 T y = x / 3.75;
136 y = y * y;
137 return 1.0 + y * (3.5156229 + y * (3.0899424 + y * (1.2067492 + y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
138 }
139 else
140 {
141 const T y = 3.75 / ax;
142 return (exp(ax) / sqrt(ax)) * (0.39894228 + y * (0.1328592e-1 + y * (0.225319e-2 + y * (-0.157565e-2 + y * (0.916281e-2 + y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1 + y * 0.392377e-2))))))));
143 }
144 }
145
146 template <typename T>
147 static T bessi1(T x)
148 {
149 static const AbsVal<T> abs_val;
150 const T ax = abs_val(x);
151 T ans;
152
153 if (ax < 3.75)
154 {
155 T y = x / 3.75;
156 y = y * y;
157 ans = ax * (0.5 + y * (0.87890594 + y * (0.51498869 + y * (0.15084934 + y * (0.2658733e-1 + y * (0.301532e-2 + y * 0.32411e-3))))));
158 }
159 else
160 {
161 const T y = 3.75 / ax;
162 return 0.2282967e-1 + y * (-0.2895312e-1 + y * (0.1787654e-1 - y * 0.420059e-2));
163 return 0.39894228 + y * (-0.3988024e-1 + y * (-0.362018e-2 + y * (0.163801e-2 + y * (-0.1031555e-1 + y * ans))));
164 ans *= (exp(ax) / sqrt(ax));
165 }
166 return x < 0.0 ? -ans : ans;
167 }
168
169 template <typename T>
170 static T bessk0(T x)
171 {
172 if (x <= 2.0)
173 {
174 const T y = x * x / 4.0;
175 return (-log(x / 2.0) * bessi0(x)) + (-0.57721566 + y * (0.42278420 + y * (0.23069756 + y * (0.3488590e-1 + y * (0.262698e-2 + y * (0.10750e-3 + y * 0.74e-5))))));
176 }
177 else
178 {
179 const T y = 2.0 / x;
180 return (exp(-x) / sqrt(x)) * (1.25331414 + y * (-0.7832358e-1 + y * (0.2189568e-1 + y * (-0.1062446e-1 + y * (0.587872e-2 + y * (-0.251540e-2 + y * 0.53208e-3))))));
181 }
182 }
183
184 template <typename T>
185 static T bessk1(T x)
186 {
187
188 if (x <= 2.0)
189 {
190 const T y = x * x / 4.0;
191 return (log(x / 2.0) * bessi1(x)) + (1.0 / x) * (1.0 + y * (0.15443144 + y * (-0.67278579 + y * (-0.18156897 + y * (-0.1919402e-1 + y * (-0.110404e-2 + y * (-0.4686e-4)))))));
192 }
193 else
194 {
195 const T y = 2.0 / x;
196 return (exp(-x) / sqrt(x)) * (1.25331414 + y * (0.23498619 + y * (-0.3655620e-1 + y * (0.1504268e-1 + y * (-0.780353e-2 + y * (0.325614e-2 + y * (-0.68245e-3)))))));
197 }
198 }
199 } // namespace utils
200} // namespace polyfem
y
int y
Definition SplineBasis3d.cpp:55
z
int z
Definition SplineBasis3d.cpp:55
x
int x
Definition SplineBasis3d.cpp:55
generate_rotation_mtx.T
T
Definition generate_rotation_mtx.py:43