1 /*
2 * $Id: bessel.i,v 1.2 2007-11-10 20:03:49 dhmunro Exp $
3 * A few Bessel functions.
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 /* Taken from Numerical Recipes, repaired bessj, bessi for x<<1. */
12
13 /* ------------------------------------------------------------------------ */
14
bessj0(x)15 func bessj0(x)
16 /* DOCUMENT bessj0(x)
17 returns Bessel function J0 at points X.
18 SEE ALSO: bessj
19 */
20 {
21 if (structof(x)==complex) error, "Bessel functions not valid for complex";
22 return mergef(x, _bessj0_1, abs(x)<8.0, _bessj0_2);
23 }
24
_bessj0_1(x)25 func _bessj0_1(x)
26 {
27 y = x*x;
28 return poly(y, 57568490574.0, -13362590354.0, 651619640.7, -11214424.18,
29 77392.33017, -184.9052456) /
30 poly(y, 57568490411.0, 1029532985.0, 9494680.718, 59272.64853,
31 267.8532712, 1.0);
32 }
33
_bessj0_2(x)34 func _bessj0_2(x)
35 {
36 ax = abs(x);
37 z = 8.0/ax;
38 y = z*z;
39 x = ax-0.785398164; /* pi/4, rounded incorrectly */
40 return sqrt(0.636619772/ax) *
41 (cos(x)*poly(y, 1.0, -0.1098628627e-2,
42 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6) -
43 sin(x)*z*poly(y, -0.1562499995e-1, 0.1430488765e-3,
44 -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7));
45 }
46
bessj1(x)47 func bessj1(x)
48 /* DOCUMENT bessj1(x)
49 returns Bessel function J1 at points X.
50 SEE ALSO: bessj
51 */
52 {
53 if (structof(x)==complex) error, "Bessel functions not valid for complex";
54 return mergef(x, _bessj1_1, abs(x)<8.0, _bessj1_2);
55 }
56
_bessj1_1(x)57 func _bessj1_1(x)
58 {
59 y = x*x;
60 return x * poly(y, 72362614232.0, -7895059235.0, 242396853.1, -2972611.439,
61 15704.48260, -30.16036606) /
62 poly(y, 144725228442.0, 2300535178.0, 18583304.74, 99447.43394,
63 376.9991397, 1.0);
64 }
65
_bessj1_2(x)66 func _bessj1_2(x)
67 {
68 ax = abs(x);
69 z = 8.0/ax;
70 y = z*z;
71 xx = ax-2.356194491; /* 3*pi/4 */
72 return sign(x) * sqrt(0.636619772/ax) *
73 (cos(xx)*poly(y, 1.0, 0.183105e-2, -0.3516396496e-4,
74 0.2457520174e-5, -0.240337019e-6) -
75 sin(xx)*z*poly(y, 0.04687499995, -0.2002690873e-3, 0.8449199096e-5,
76 -0.88228987e-6, 0.105787412e-6));
77 }
78
bessj(n,x)79 func bessj(n, x)
80 /* DOCUMENT bessj(n, x)
81 returns Bessel function Jn of order N at points X. N must be scalar.
82 SEE ALSO: bessy, bessi, bessk, bessj0, bessj1
83 */
84 {
85 if (n>1) {
86 ax = abs(x);
87 bj = mergef(ax, _bessj_0, ax<0.02*sqrt(n), _bessj_1, ax>n, _bessj_2);
88 if (n%2) bj *= sign(x);
89 return bj;
90 } else if (n==1) {
91 return bessj1(x);
92 } else if (!n) {
93 return bessj0(x);
94 }
95 }
96
_bessj_0(x)97 func _bessj_0(x)
98 {
99 x *= 0.5;
100 nn = double(n);
101 rnf = exp(-sum(log(indgen(n))));
102 rn1 = 1./(nn+1.); rn2 = 0.5*rn1/(nn+2.);
103 return (x^n*rnf) * poly(x*x, 1., -rn1, rn2, -rn2/(3.*(nn+3.)));
104 }
105
_bessj_1(x)106 func _bessj_1(x)
107 {
108 /* upward recurrence abs(x)>n */
109 ax = abs(x);
110 tox = 2.0/ax;
111 bjm = bessj0(ax);
112 bj = bessj1(ax);
113 for (i=1 ; i<n ; i++) {
114 bjp = i*tox*bj-bjm;
115 bjm = bj;
116 bj = bjp;
117 }
118 return bj;
119 }
120
_bessj_2(x)121 func _bessj_2(x)
122 {
123 /* downward recurrence abs(x)<=n */
124 ax = abs(x);
125 tox = 2.0/ax; /* < 100/sqrt(n) */
126 m = 2*((n+long(sqrt(bess_acc*n)))/2);
127 jsum = 0;
128 bjp = ans = add = array(0.0, numberof(ax));
129 bj = array(1.0, numberof(ax));
130 for (i=m ; i>0 ; i--) {
131 bjm = i*tox*bj-bjp;
132 bjp = bj;
133 bj = bjm;
134 list = where(abs(bj) > bess_big);
135 if (numberof(list)) {
136 bess_nrm = 1./bess_big;
137 bj(list) *= bess_nrm;
138 bjp(list) *= bess_nrm;
139 ans(list) *= bess_nrm;
140 add(list) *= bess_nrm;
141 }
142 if (jsum) add += bj;
143 jsum = !jsum;
144 if (i==n) ans = bjp;
145 }
146 bj = ans/(2.0*add-bj);
147 return bj;
148 }
149
150 /* ------------------------------------------------------------------------ */
151
bessy0(x)152 func bessy0(x)
153 /* DOCUMENT bessy0(x)
154 returns Bessel function Y0 at points X.
155 SEE ALSO: bessy
156 */
157 {
158 if (structof(x)==complex) error, "Bessel functions not valid for complex";
159 return mergef(x, _bessy0_1, abs(x)<8.0, _bessy0_2);
160 }
161
_bessy0_1(x)162 func _bessy0_1(x)
163 {
164 y= x*x;
165 return poly(y, -2957821389.0, 7062834065.0, -512359803.6, 10879881.29,
166 -86327.92757, 228.4622733) /
167 poly(y, 40076544269.0, 745249964.8, 7189466.438, 47447.26470,
168 226.1030244, 1.0) + 0.636619772*bessj0(x)*log(x);
169 }
170
_bessy0_2(x)171 func _bessy0_2(x)
172 {
173 ax = abs(x);
174 z = 8.0/ax;
175 y = z*z;
176 xx = ax-0.785398164; /* pi/4, rounded incorrectly */
177 return sqrt(0.636619772/ax) *
178 (sin(xx)*poly(y, 1.0, -0.1098628627e-2, 0.2734510407e-4,
179 -0.2073370639e-5, 0.2093887211e-6) +
180 cos(xx)*z*poly(y, -0.1562499995e-1, 0.1430488765e-3,
181 -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7));
182 }
183
bessy1(x)184 func bessy1(x)
185 /* DOCUMENT bessy1(x)
186 returns Bessel function Y1 at points X.
187 SEE ALSO: bessy
188 */
189 {
190 if (structof(x)==complex) error, "Bessel functions not valid for complex";
191 return mergef(x, _bessy1_1, abs(x)<8.0, _bessy1_2);
192 }
193
_bessy1_1(x)194 func _bessy1_1(x)
195 {
196 y = x*x;
197 return x * poly(y, -0.4900604943e13, 0.1275274390e13, -0.5153438139e11,
198 0.7349264551e9, -0.4237922726e7, 0.8511937935e4) /
199 poly(y, 0.2499580570e14, 0.4244419664e12, 0.3733650367e10,
200 0.2245904002e8, 0.1020426050e6, 0.3549632885e3, 1.0) +
201 0.636619772*(bessj1(x)*log(x)-1.0/x);
202 }
203
_bessy1_2(x)204 func _bessy1_2(x)
205 {
206 ax = abs(x);
207 z = 8.0/ax;
208 y = z*z;
209 xx = ax-2.356194491; /* 3*pi/4 */
210 return sqrt(0.636619772/x) *
211 (sin(xx)*poly(y, 1.0, 0.183105e-2, -0.3516396496e-4,
212 0.2457520174e-5, -0.240337019e-6) +
213 cos(xx)*z*poly(y, 0.04687499995, -0.2002690873e-3, 0.8449199096e-5,
214 -0.88228987e-6, 0.105787412e-6));
215 }
216
bessy(n,x)217 func bessy(n, x)
218 /* DOCUMENT bessy(n, x)
219 returns Bessel function Yn of order N at points X. N must be scalar.
220 SEE ALSO: bessj, bessi, bessk, bessy0, bessy1
221 */
222 {
223 if (n>1) {
224 /* upward recurrence */
225 tox = 2.0/x;
226 bym = bessy0(x);
227 by = bessy1(x);
228 for (i=1 ; i<n ; i++) {
229 byp = i*tox*by-bym;
230 bym = by;
231 by = byp;
232 }
233 return by;
234 } else if (n==1) {
235 return bessy1(x);
236 } else if (!n) {
237 return bessy0(x);
238 }
239 }
240
241 /* ------------------------------------------------------------------------ */
242
bessi0(x)243 func bessi0(x)
244 /* DOCUMENT bessi0(x)
245 returns Bessel function I0 at points X.
246 SEE ALSO: bessi
247 */
248 {
249 if (structof(x)==complex) error, "Bessel functions not valid for complex";
250 x = abs(x);
251 return mergef(x, _bessi0_1, x<3.75, _bessi0_2);
252 }
253
_bessi0_1(x)254 func _bessi0_1(x)
255 {
256 x = x/3.75;
257 return poly(x*x, 1.0, 3.5156229, 3.0899424, 1.2067492, 0.2659732,
258 0.360768e-1, 0.45813e-2);
259 }
260
_bessi0_2(x)261 func _bessi0_2(x)
262 {
263 y = 3.75/x;
264 return (exp(x)/sqrt(x)) * poly(y, 0.39894228, 0.1328592e-1, 0.225319e-2,
265 -0.157565e-2, 0.916281e-2, -0.2057706e-1,
266 0.2635537e-1, -0.1647633e-1, 0.392377e-2);
267 }
268
bessi1(x)269 func bessi1(x)
270 /* DOCUMENT bessi1(x)
271 returns Bessel function I1 at points X.
272 SEE ALSO: bessi
273 */
274 {
275 if (structof(x)==complex) error, "Bessel functions not valid for complex";
276 return mergef(x, _bessi1_1, abs(x)<3.75, _bessi1_2);
277 }
278
_bessi1_1(x)279 func _bessi1_1(x)
280 {
281 y = x/3.75;
282 y *= y;
283 return x * poly(y, 0.5, 0.87890594, 0.51498869, 0.15084934,
284 0.2658733e-1, 0.301532e-2, 0.32411e-3);
285 }
286
_bessi1_2(x)287 func _bessi1_2(x)
288 {
289 ax = abs(x);
290 y = 3.75/ax;
291 return sign(x) * (exp(ax)/sqrt(ax)) *
292 poly(y, 0.39894228, -0.3988024e-1, -0.362018e-2, 0.163801e-2,
293 -0.1031555e-1, 0.2282967e-1, -0.2895312e-1, 0.1787654e-1,
294 -0.420059e-2);
295 }
296
bessi(n,x)297 func bessi(n, x)
298 /* DOCUMENT bessi(n, x)
299 returns Bessel function In of order N at points X. N must be scalar.
300 SEE ALSO: bessk, bessj, bessy, bessi0, bessi1
301 */
302 {
303 if (n>1) {
304 ax = abs(x);
305 bi = mergef(ax, _bessi_0, ax<0.02*sqrt(n), _bessi_1);
306 if (n%2) bi *= sign(x);
307 return bi;
308 } else if (n==1) {
309 return bessi1(x);
310 } else if (!n) {
311 return bessi0(x);
312 }
313 }
314
_bessi_0(x)315 func _bessi_0(x)
316 {
317 x *= 0.5;
318 nn = double(n);
319 rnf = exp(-sum(log(indgen(n))));
320 rn1 = 1./(nn+1.); rn2 = 0.5*rn1/(nn+2.);
321 return (x^n*rnf) * poly(x*x, 1., rn1, rn2, rn2/(3.*(nn+3.)));
322 }
323
_bessi_1(x)324 func _bessi_1(x)
325 {
326 /* downward recurrence abs(x)<=n */
327 tox = 2.0/x;
328 m = 2*(n+long(sqrt(bess_acc*n)));
329 bip = ans = array(0.0, numberof(x));
330 bi = array(1.0, numberof(x));
331 for (i=m ; i>0 ; i--) {
332 bim = i*tox*bi+bip;
333 bip = bi;
334 bi = bim;
335 list = where(abs(bi) > bess_big);
336 if (numberof(list)) {
337 bess_nrm = 1./bess_big;
338 ans(list) *= bess_nrm;
339 bi(list) *= bess_nrm;
340 bip(list) *= bess_nrm;
341 }
342 if (i==n) ans = bip;
343 }
344 bi = ans*bessi0(x)/bi;
345 return bi;
346 }
347
348 /* ------------------------------------------------------------------------ */
349
bessk0(x)350 func bessk0(x)
351 /* DOCUMENT bessk0(x)
352 returns Bessel function K0 at points X.
353 SEE ALSO: bessk
354 */
355 {
356 if (structof(x)==complex) error, "Bessel functions not valid for complex";
357 return mergef(x, _bessk0_1, x<=2.0, _bessk0_2);
358 }
359
_bessk0_1(x)360 func _bessk0_1(x)
361 {
362 y= x*x/4.0;
363 return (-log(x/2.0)*bessi0(x)) +
364 poly(y, -0.57721566, 0.42278420, 0.23069756, 0.3488590e-1, 0.262698e-2,
365 0.10750e-3, 0.74e-5);
366 }
367
_bessk0_2(x)368 func _bessk0_2(x)
369 {
370 y = 2.0/x;
371 return (exp(-x)/sqrt(x)) *
372 poly(y, 1.25331414, -0.7832358e-1, 0.2189568e-1, -0.1062446e-1,
373 0.587872e-2, -0.251540e-2, 0.53208e-3);
374 }
375
bessk1(x)376 func bessk1(x)
377 /* DOCUMENT bessk1(x)
378 returns Bessel function K1 at points X.
379 SEE ALSO: bessk
380 */
381 {
382 if (structof(x)==complex) error, "Bessel functions not valid for complex";
383 return mergef(x, _bessk1_1, x<=2.0, _bessk1_2);
384 }
385
_bessk1_1(x)386 func _bessk1_1(x)
387 {
388 y = x*x/4.0;
389 return (log(x/2.0)*bessi1(x)) +
390 (1.0/x) * poly(y, 1.0, 0.15443144, -0.67278579, -0.18156897,
391 -0.1919402e-1, -0.110404e-2, -0.4686e-4);
392 }
393
_bessk1_2(x)394 func _bessk1_2(x)
395 {
396 y = 2.0/x;
397 return (exp(-x)/sqrt(x)) *
398 poly(y, 1.25331414, 0.23498619, -0.3655620e-1, 0.1504268e-1,
399 -0.780353e-2, 0.325614e-2, -0.68245e-3);
400 }
401
bessk(n,x)402 func bessk(n, x)
403 /* DOCUMENT bessk(n, x)
404 returns Bessel function Kn of order N at points X. N must be scalar.
405 SEE ALSO: bessi, bessj, bessy, bessi0, bessi1
406 */
407 {
408 if (n>1) {
409 /* upward recurrence */
410 tox = 2.0/x;
411 bkm = bessk0(x);
412 bk = bessk1(x);
413 for (i=1 ; i<n ; i++) {
414 bkp = i*tox*bk+bkm;
415 bkm = bk;
416 bk = bkp;
417 }
418 return bk;
419 } else if (n==1) {
420 return bessk1(x);
421 } else if (!n) {
422 return bessk0(x);
423 }
424 }
425
426 /* ------------------------------------------------------------------------ */
427
428 bess_acc= 40.0;
429 bess_big= 1.e10;
430
431 /* ------------------------------------------------------------------------ */
432
433 #if 0
434 func bess_check(void)
435 {
436 /* values copied from Abramowitz and Stegun tables */
437 eg = 0.5772156649;
438 eps = 0.5e-30;
439 x = [2.*eps, 0.6, 3.0, 17.0];
440 fac5 = 5.*4.*3.*2.; fac6 = 6.*fac5;
441 j0 = [1., 0.912004863497211, -0.260051954901933, -0.169854252151184];
442 j1 = [eps, 0.2867009881, 0.3390589585, -0.0976684928];
443 j5 = [eps^5/fac5, 1.9948e-5, 4.3028e-2, -0.18704];
444 j6 = [eps^6/fac6, 9.9956e-7, 1.1394e-2, 0.00071533];
445 y0 = [2./pi*(log(eps)+eg), -0.3085098701, 0.3768500100, -0.0926371984];
446 y1 = [-1./(pi*eps), -1.2603913472, 0.3246744248, 0.1672050361];
447 y5 = [-24./(pi*eps^5), -3.2156e3, -1.9059, 0.06455];
448 y6 = [-120./(pi*eps^6), -5.3351e4, -5.4365, 0.19996];
449 i0 = [1., 0.5993272031, 0.2430003542, 0.0974943005]*exp(x);
450 i1 = [eps, 0.1721644195, 0.1968267133, 0.0945819107]*exp(x);
451 i5 = [eps^5/fac5, 1.1281e-5, 4.5409e-3, 4.5951e-2]*exp(x);
452 i6 = [eps^6/fac6, 5.6286e-7, 1.0796e-3, 3.3128e-2]*exp(x);
453 k0 = [-(log(eps)+eg), 1.4167376214, 0.6977615980, 0.3018080193]*exp(-x);
454 k1 = [0.5/eps, 2.3739200376, 0.8065634800, 0.3105612340]*exp(-x);
455 k5 = [12./eps^5, 8.7987e3, 1.8836e1, 6.1420e-1]*exp(-x);
456 k6 = [60./eps^6, 1.4730e5, 6.8929e1, 8.3734e-1]*exp(-x);
457 y = x;
458 for (i=1 ; i<=numberof(x) ; i++) y(i) = bessj(5,x(i));
459 if (anyof(y != bessj(5,x))) write, "ERROR - problem with scalar args";
460 write, "j0:", max(abs(bessj(0,x)/j0-1.)), max(abs(bessj(0,-x)/j0-1.));
461 write, "j1:", max(abs(bessj(1,x)/j1-1.)), max(abs(bessj(1,-x)/j1+1.));
462 write, "j5:", max(abs(bessj(5,x)/j5-1.)), max(abs(bessj(5,-x)/j5+1.));
463 write, "j6:", max(abs(bessj(6,x)/j6-1.)), max(abs(bessj(6,-x)/j6-1.));
464 write, "y0:", max(abs(bessy(0,x)/y0-1.));
465 write, "y1:", max(abs(bessy(1,x)/y1-1.));
466 write, "y5:", max(abs(bessy(5,x)/y5-1.));
467 write, "y6:", max(abs(bessy(6,x)/y6-1.));
468 write, "i0:", max(abs(bessi(0,x)/i0-1.)), max(abs(bessi(0,-x)/i0-1.));
469 write, "i1:", max(abs(bessi(1,x)/i1-1.)), max(abs(bessi(1,-x)/i1+1.));
470 write, "i5:", max(abs(bessi(5,x)/i5-1.)), max(abs(bessi(5,-x)/i5+1.));
471 write, "i6:", max(abs(bessi(6,x)/i6-1.)), max(abs(bessi(6,-x)/i6-1.));
472 write, "k0:", max(abs(bessk(0,x)/k0-1.));
473 write, "k1:", max(abs(bessk(1,x)/k1-1.));
474 write, "k5:", max(abs(bessk(5,x)/k5-1.));
475 write, "k6:", max(abs(bessk(6,x)/k6-1.));
476 }
477 #endif
478