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