Home Manual Packages Global Index Keywords Quick Reference ``` /* * bessel.i -- \$Id\$ * A few Bessel functions. */ /* Copyright (c) 2002. The Regents of the University of California. * All rights reserved. */ /* Taken from Numerical Recipes, repaired bessj, bessi for x<<1. */ /* ------------------------------------------------------------------------ */ func bessj0 (x) /* DOCUMENT bessj0(x) returns Bessel function J0 at points X. SEE ALSO: bessj */ { return mergef(x, _bessj0_1, abs(x)<8.0, _bessj0_2); } func _bessj0_1 (x) { y = x*x; return poly(y, 57568490574.0, -13362590354.0, 651619640.7, -11214424.18, 77392.33017, -184.9052456) / poly(y, 57568490411.0, 1029532985.0, 9494680.718, 59272.64853, 267.8532712, 1.0); } func _bessj0_2 (x) { ax = abs(x); z = 8.0/ax; y = z*z; x = ax-0.785398164; /* pi/4, rounded incorrectly */ return sqrt(0.636619772/ax) * (cos(x)*poly(y, 1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6) - sin(x)*z*poly(y, -0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7)); } func bessj1 (x) /* DOCUMENT bessj1(x) returns Bessel function J1 at points X. SEE ALSO: bessj */ { return mergef(x, _bessj1_1, abs(x)<8.0, _bessj1_2); } func _bessj1_1 (x) { y = x*x; return x * poly(y, 72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606) / poly(y, 144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0); } func _bessj1_2 (x) { ax = abs(x); z = 8.0/ax; y = z*z; xx = ax-2.356194491; /* 3*pi/4 */ return sign(x) * sqrt(0.636619772/ax) * (cos(xx)*poly(y, 1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6) - sin(xx)*z*poly(y, 0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6)); } func bessj (n, x) /* DOCUMENT bessj(n, x) returns Bessel function Jn of order N at points X. N must be scalar. SEE ALSO: bessy, bessi, bessk, bessj0, bessj1 */ { if (n>1) { ax = abs(x); bj = mergef(ax, _bessj_0, ax<0.02*sqrt(n), _bessj_1, ax>n, _bessj_2); if (n%2) bj *= sign(x); return bj; } else if (n==1) { return bessj1(x); } else if (!n) { return bessj0(x); } } func _bessj_0 (x) { x *= 0.5; nn = double(n); rn1 = 1./(nn+1.); rn2 = rn1/(nn+2.); return (x^n/nn) * poly(x, 1., -rn1, rn2, -rn2/(nn+3.)); } func _bessj_1 (x) { /* upward recurrence abs(x)>n */ ax = abs(x); tox = 2.0/ax; bjm = bessj0(ax); bj = bessj1(ax); for (i=1 ; i0 ; i--) { bjm = i*tox*bj-bjp; bjp = bj; bj = bjm; list = where(abs(bj) > bess_big); if (numberof(list)) { bess_nrm = 1./bess_big; bj(list) *= bess_nrm; bjp(list) *= bess_nrm; ans(list) *= bess_nrm; add(list) *= bess_nrm; } if (jsum) add += bj; jsum = !jsum; if (i==n) ans = bjp; } bj = ans/(2.0*add-bj); return bj; } /* ------------------------------------------------------------------------ */ func bessy0 (x) /* DOCUMENT bessy0(x) returns Bessel function Y0 at points X. SEE ALSO: bessy */ { return mergef(x, _bessy0_1, abs(x)<8.0, _bessy0_2); } func _bessy0_1 (x) { y= x*x; return poly(y, -2957821389.0, 7062834065.0, -512359803.6, 10879881.29, -86327.92757, 228.4622733) / poly(y, 40076544269.0, 745249964.8, 7189466.438, 47447.26470, 226.1030244, 1.0) + 0.636619772*bessj0(x)*log(x); } func _bessy0_2 (x) { ax = abs(x); z = 8.0/ax; y = z*z; xx = ax-0.785398164; /* pi/4, rounded incorrectly */ return sqrt(0.636619772/ax) * (sin(xx)*poly(y, 1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6) + cos(xx)*z*poly(y, -0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7)); } func bessy1 (x) /* DOCUMENT bessy1(x) returns Bessel function Y1 at points X. SEE ALSO: bessy */ { return mergef(x, _bessy1_1, abs(x)<8.0, _bessy1_2); } func _bessy1_1 (x) { y = x*x; return x * poly(y, -0.4900604943e13, 0.1275274390e13, -0.5153438139e11, 0.7349264551e9, -0.4237922726e7, 0.8511937935e4) / poly(y, 0.2499580570e14, 0.4244419664e12, 0.3733650367e10, 0.2245904002e8, 0.1020426050e6, 0.3549632885e3, 1.0) + 0.636619772*(bessj1(x)*log(x)-1.0/x); } func _bessy1_2 (x) { ax = abs(x); z = 8.0/ax; y = z*z; xx = ax-2.356194491; /* 3*pi/4 */ return sqrt(0.636619772/x) * (sin(xx)*poly(y, 1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6) + cos(xx)*z*poly(y, 0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6)); } func bessy (n, x) /* DOCUMENT bessy(n, x) returns Bessel function Yn of order N at points X. N must be scalar. SEE ALSO: bessj, bessi, bessk, bessy0, bessy1 */ { if (n>1) { /* upward recurrence */ tox = 2.0/x; bym = bessy0(x); by = bessy1(x); for (i=1 ; i1) { ax = abs(x); bi = mergef(ax, _bessi_0, ax<0.02*sqrt(n), _bessi_1); if (n%2) bi *= sign(x); return bi; } else if (n==1) { return bessi1(x); } else if (!n) { return bessi0(x); } } func _bessi_0 (x) { x *= 0.5; nn = double(n); rn1 = 1./(nn+1.); rn2 = rn1/(nn+2.); return (x^n/nn) * poly(x, 1., rn1, rn2, rn2/(nn+3.)); } func _bessi_1 (x) { /* downward recurrence abs(x)<=n */ tox = 2.0/x; m = 2*(n+long(sqrt(bess_acc*n))); bip = ans = array(0.0, numberof(x)); bi = array(1.0, numberof(x)); for (i=m ; i>0 ; i--) { bim = i*tox*bi+bip; bip = bi; bi = bim; list = where(abs(bi) > bess_big); if (numberof(list)) { bess_nrm = 1./bess_big; ans(list) *= bess_nrm; bi(list) *= bess_nrm; bip(list) *= bess_nrm; } if (i==n) ans = bip; } bi = ans*bessi0(x)/bi; return bi; } /* ------------------------------------------------------------------------ */ func bessk0 (x) /* DOCUMENT bessk0(x) returns Bessel function K0 at points X. SEE ALSO: bessk */ { return mergef(x, _bessk0_1, x<=2.0, _bessk0_2); } func _bessk0_1 (x) { y= x*x/4.0; return (-log(x/2.0)*bessi0(x)) + poly(y, -0.57721566, 0.42278420, 0.23069756, 0.3488590e-1, 0.262698e-2, 0.10750e-3, 0.74e-5); } func _bessk0_2 (x) { y = 2.0/x; return (exp(-x)/sqrt(x)) * poly(y, 1.25331414, -0.7832358e-1, 0.2189568e-1, -0.1062446e-1, 0.587872e-2, -0.251540e-2, 0.53208e-3); } func bessk1 (x) /* DOCUMENT bessk1(x) returns Bessel function K1 at points X. SEE ALSO: bessk */ { return mergef(x, _bessk1_1, x<=2.0, _bessk1_2); } func _bessk1_1 (x) { y = x*x/4.0; return (log(x/2.0)*bessi1(x)) + (1.0/x) * poly(y, 1.0, 0.15443144, -0.67278579, -0.18156897, -0.1919402e-1, -0.110404e-2, -0.4686e-4); } func _bessk1_2 (x) { y = 2.0/x; return (exp(-x)/sqrt(x)) * poly(y, 1.25331414, 0.23498619, -0.3655620e-1, 0.1504268e-1, -0.780353e-2, 0.325614e-2, -0.68245e-3); } func bessk (n, x) /* DOCUMENT bessk(n, x) returns Bessel function Kn of order N at points X. N must be scalar. SEE ALSO: bessi, bessj, bessy, bessi0, bessi1 */ { if (n>1) { /* upward recurrence */ tox = 2.0/x; bkm = bessk0(x); bk = bessk1(x); for (i=1 ; i