Branch data Line data Source code
1 : : /*
2 : : * cbesselj.cpp - Bessel function of first kind
3 : : *
4 : : * Copyright (C) 2007 Bastien Roucaries <roucaries.bastien@gmail.com>
5 : : *
6 : : * This is free software; you can redistribute it and/or modify
7 : : * it under the terms of the GNU General Public License as published by
8 : : * the Free Software Foundation; either version 2, or (at your option)
9 : : * any later version.
10 : : *
11 : : * This software is distributed in the hope that it will be useful,
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : : * GNU General Public License for more details.
15 : : *
16 : : * You should have received a copy of the GNU General Public License
17 : : * along with this package; see the file COPYING. If not, write to
18 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 : : * Boston, MA 02110-1301, USA.
20 : : *
21 : : * $Id$
22 : : *
23 : : */
24 : :
25 : : /*!\file cbesselj.cpp
26 : : \brief compute complex bessel J function
27 : :
28 : : Bibligraphy:
29 : :
30 : : [1] Bessel function of the first kind with complex argument
31 : : Yousif, Hashim A.; Melka, Richard
32 : : Computer Physics Communications, vol. 106, Issue 3, pp.199-206
33 : : 11/1997, ELSEVIER, doi://10.1016/S0010-4655(97)00087-8
34 : :
35 : : [2] Handbook of Mathematical Functions with
36 : : Formulas, Graphs, and Mathematical Tables
37 : : Milton Abramowitz and Irene A. Stegun
38 : : public domain (work of US government)
39 : : online http://www.math.sfu.ca/~cbm/aands/
40 : :
41 : : [3] Mathematica Manual
42 : : Bessel, Airy, Struve Functions> BesselJ[nu,z]
43 : : > General characteristics> Symmetries and periodicities
44 : : http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/04/02/01/
45 : :
46 : : [4] Mathematica Manual
47 : : Bessel, Airy, Struve Functions> BesselJ[nu,z]
48 : : Representations through equivalent functions
49 : : http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/27/ShowAll.html
50 : :
51 : : [5] Algorithms for the evaluation of Bessel functions
52 : : of complex argument and integer orders
53 : : G. D. C. Kuiken
54 : : Applied Mathematics Letters, Volume 2, Issue 4,
55 : : 1989, Pages 353-356
56 : : doi://10.1016/0893-9659(89)90086-4
57 : :
58 : : */
59 : :
60 : :
61 : : //#if HAVE_CONFIG_H
62 : : //# include <config.h>
63 : : //#endif
64 : : //
65 : : //#include <cmath>
66 : : //#include <assert.h>
67 : : //#include <errno.h>
68 : : //#include <stdio.h>
69 : : //#include <stdlib.h>
70 : : //
71 : : //#include "real.h"
72 : : //#include "complex.h"
73 : : //#include "constants.h"
74 : : //#include "precision.h"
75 : : //#include <limits>
76 : :
77 : : #define SMALL_J0_BOUND 1e6
78 : :
79 : : /*!\brief use ascending serie below this magnitude */
80 : : #define SMALL_JN_BOUND 5.0
81 : :
82 : : /*!\brief use assymptotic expression above this magnitude */
83 : : #define BIG_JN_BOUND 25.0
84 : :
85 : : /*! \brief Arbitrary value for iteration*/
86 : : #define MAX_SMALL_ITERATION 2048
87 : :
88 : :
89 : : /*!\brief Implement bessel J function for small arguments
90 : : according to [5]
91 : :
92 : : For small argument we use the following formulae:
93 : : \{align}
94 : : J_n(z)&=\sum_0^\infty R_k & & \\
95 : : R_{-1}&=1 & R_k &= a_k R_{k-1} \\
96 : : a_0&=\frac{\left(\frac{1}{2}z\right)^n}{n!} &
97 : : a_{+k}&=\frac{-\frac{1}{4} z^2}{k(n+k)}
98 : : \}
99 : :
100 : : \todo Not really adapted to high order
101 : : therefore we do not check overflow for n >> 1
102 : :
103 : : \param[in] n order
104 : : \param[in] arg arguments
105 : : */
106 : : static nr_complex_t
107 : 0 : cbesselj_smallarg (unsigned int n, nr_complex_t z)
108 : : {
109 : 0 : nr_complex_t ak, Rk;
110 : 0 : nr_complex_t R;
111 : 0 : nr_complex_t R0;
112 : : unsigned int k;
113 : :
114 : : /* a_0 */
115 [ # # ]: 0 : errno = 0;
116 [ # # ]: 0 : ak = pow (0.5 * z, n);
117 : : /* underflow */
118 [ # # ][ # # ]: 0 : if (errno == ERANGE)
119 : : {
120 [ # # ]: 0 : return n > 0 ? 0.0 : 1;
121 : : }
122 : :
123 [ # # ]: 0 : ak = ak / (nr_double_t) qucs::factorial (n);
124 : :
125 : : /* R_0 */
126 : 0 : R0 = ak * 1.0;
127 : 0 : Rk = R0;
128 : 0 : R = R0;
129 : :
130 [ # # ]: 0 : for (k = 1; k < MAX_SMALL_ITERATION; k++)
131 : : {
132 [ # # ]: 0 : ak = -(z * z) / (4.0 * k * (n + k));
133 [ # # ]: 0 : Rk = ak * Rk;
134 [ # # # # ]: 0 : if (fabs (real (Rk)) < fabs (real (R) * std::numeric_limits<nr_double_t>::epsilon()) &&
[ # # ]
135 : 0 : fabs (imag (Rk)) < fabs (imag (R) * std::numeric_limits<nr_double_t>::epsilon()))
136 : 0 : return R;
137 : :
138 : 0 : R += Rk;
139 : : }
140 : : /* impossible */
141 [ # # ]: 0 : assert (k < MAX_SMALL_ITERATION);
142 : 0 : abort ();
143 : : }
144 : :
145 : :
146 : : static nr_complex_t
147 : 0 : cbesselj_mediumarg_odd (unsigned int n, nr_complex_t z)
148 : : {
149 : 0 : nr_complex_t first, second;
150 : 0 : nr_complex_t ak;
151 : :
152 : : unsigned int m;
153 : : unsigned int k;
154 : : nr_double_t t;
155 : : nr_double_t m1pna2;
156 : :
157 : 0 : m = (2 * std::abs (z) + 0.25 * (n + std::abs (imag (z))));
158 : :
159 : : /* -1^(n/2) */
160 [ # # ]: 0 : m1pna2 = (n / 2) % 2 == 0 ? 1.0 : -1.0;
161 [ # # ]: 0 : first = (1.0 + m1pna2 * cos (z)) / (2.0 * m);
162 : :
163 : 0 : second = 0.0;
164 [ # # ]: 0 : for (k = 1; k <= m - 1; k++)
165 : : {
166 : 0 : t = (k * M_PI) / (2 * m);
167 [ # # ]: 0 : ak = cos (z * std::sin (t)) * std::cos (n * t);
168 : 0 : second += ak;
169 : : }
170 [ # # ]: 0 : return first + second / (nr_double_t) m;
171 : : }
172 : :
173 : : static nr_complex_t
174 : 0 : cbesselj_mediumarg_even (unsigned int n, nr_complex_t z)
175 : : {
176 : 0 : nr_complex_t first, second;
177 : 0 : nr_complex_t ak;
178 : :
179 : : unsigned int m;
180 : : unsigned int k;
181 : : nr_double_t t;
182 : : nr_double_t m1pn1a2;
183 : :
184 : 0 : m = (2 * std::abs (z) + 0.25 * (n + std::abs (imag (z))));
185 : :
186 : : /* -1^(n/2) */
187 [ # # ]: 0 : m1pn1a2 = ((n - 1) / 2) % 2 == 0 ? 1.0 : -1.0;
188 [ # # ]: 0 : first = (m1pn1a2 * sin (z)) / (2.0 * m);
189 : :
190 : 0 : second = 0.0;
191 [ # # ]: 0 : for (k = 1; k <= m - 1; k++)
192 : : {
193 : 0 : t = (k * M_PI) / (2 * m);
194 : 0 : ak = std::sin (z * std::sin (t)) * std::sin (n * t);
195 : 0 : second += ak;
196 : : }
197 [ # # ]: 0 : return first + second / (nr_double_t) m;
198 : : }
199 : :
200 : :
201 : : static nr_complex_t
202 : 0 : cbesselj_mediumarg (unsigned int n, nr_complex_t z)
203 : : {
204 [ # # ]: 0 : if (n % 2 == 0)
205 : 0 : return cbesselj_mediumarg_odd (n, z);
206 : : else
207 : 0 : return cbesselj_mediumarg_even (n, z);
208 : : }
209 : :
210 : :
211 : :
212 : : /*! \brief num of P(k) (n = 8) will overlow above this value */
213 : : #define MAX_LARGE_ITERATION 430
214 : :
215 : : /*!\brief besselj for large argument
216 : :
217 : : Based on [5] eq (5)
218 : : */
219 : : static nr_complex_t
220 : 0 : cbesselj_largearg (unsigned int n, nr_complex_t z)
221 : : {
222 : 0 : nr_complex_t Pk, P0, P, Qk, Q0, Q_;
223 : 0 : nr_complex_t tmp;
224 : : unsigned long num, denum;
225 : 0 : nr_complex_t m1a8z2;
226 : : unsigned int k;
227 : : nr_double_t l, m;
228 : :
229 : : /* P0 & Q0 */
230 : 0 : P0 = 1;
231 : 0 : P = P0;
232 : 0 : Pk = P0;
233 : :
234 [ # # ]: 0 : Q0 = (4.0 * n * n - 1) / (8.0 * z);
235 : 0 : Q_ = Q0;
236 : 0 : Qk = Q0;
237 : :
238 [ # # ]: 0 : m1a8z2 = (-1.0) / (8.0 * sqr (z));
239 : :
240 : : /* P */
241 : 0 : for (k = 1;; k++)
242 : : {
243 : : /* Usually converge before overflow */
244 [ # # ]: 0 : assert (k <= MAX_LARGE_ITERATION);
245 : :
246 [ # # ][ # # ]: 0 : num = (4 * sqr (n) - sqr (4 * k - 3)) * (4 * sqr (n) - (4 * k - 1));
[ # # ]
247 : 0 : denum = 2 * k * (2 * k - 1);
248 [ # # ]: 0 : Pk = Pk * ((nr_double_t) num * m1a8z2) / ((nr_double_t) denum);
249 : :
250 [ # # # # ]: 0 : if (real (Pk) < real (P0) * std::numeric_limits<nr_double_t>::epsilon() &&
[ # # ]
251 : 0 : imag (Pk) < imag (P0) * std::numeric_limits<nr_double_t>::epsilon())
252 : 0 : break;
253 : :
254 : 0 : P += Pk;
255 : : }
256 : :
257 : : /* P */
258 : 0 : for (k = 1;; k++)
259 : : {
260 : : /* Usually converge before overflow */
261 [ # # ]: 0 : assert (k <= MAX_LARGE_ITERATION);
262 : :
263 [ # # ][ # # ]: 0 : num = (4 * sqr (n) - sqr (4 * k - 1)) * (4 * sqr (n) - (4 * k - 1));
[ # # ]
264 : 0 : denum = 2 * k * (2 * k - 1);
265 [ # # ]: 0 : Qk = Qk * ((nr_double_t) num * m1a8z2) / ((nr_double_t) denum);
266 : :
267 [ # # # # ]: 0 : if (real (Qk) < real (Q0) * std::numeric_limits<nr_double_t>::epsilon() ||
[ # # ]
268 : 0 : imag (Qk) < imag (Q0) * std::numeric_limits<nr_double_t>::epsilon())
269 : 0 : break;
270 : :
271 : 0 : Q_ += Qk;
272 : : }
273 : :
274 : : /* l, m cf [5] eq (5) */
275 [ # # ][ # # ]: 0 : l = (n % 4 == 0) || (n % 4 == 3) ? 1.0 : -1.0;
276 [ # # ][ # # ]: 0 : m = (n % 4 == 0) || (n % 4 == 1) ? 1.0 : -1.0;
277 : :
278 : :
279 [ # # ][ # # ]: 0 : tmp = (l * P + m * Q_) * cos (z);
[ # # ]
280 [ # # ][ # # ]: 0 : tmp += (m * P - l * Q_) * sin (z);
[ # # ]
281 [ # # ][ # # ]: 0 : return tmp / sqrt (M_PI * z);
282 : : }
283 : :
284 : : /*!\brief Main entry point for besselj function
285 : :
286 : : */
287 : : nr_complex_t
288 : 0 : cbesselj (unsigned int n, nr_complex_t z)
289 : : {
290 : 0 : nr_double_t mul = 1.0;
291 : 0 : nr_complex_t ret;
292 : :
293 : : /* J_n(-z)=(-1)^n J_n(z) */
294 : : /*
295 : : if(real(z) < 0.0)
296 : : {
297 : : z = -z;
298 : : mul = (n % 2) == 0 ? 1.0 : -1.0 ;
299 : : }
300 : : */
301 [ # # ]: 0 : if (abs (z) < SMALL_JN_BOUND)
302 [ # # ]: 0 : ret = cbesselj_smallarg (n, z);
303 [ # # ]: 0 : else if (abs (z) > BIG_JN_BOUND)
304 [ # # ]: 0 : ret = cbesselj_largearg (n, z);
305 : : else
306 [ # # ]: 0 : ret = cbesselj_mediumarg (n, z);
307 : :
308 : 0 : return mul * ret;
309 : : }
|
