6.5 Bessel Functions of Integer Order 235 an additive amount oforder [constant xn1/2,where the square root of the constant is,very roughly,the number of significant figures of accuracy The above considerations lead to the following function. #include <math.h> #define ACC 40.0 Make larger to increase accuracy. #define BIGNO 1.0e10 http://www.nr. readable files Permission is #define BIGNI 1.0e-10 float bessj(int n,float x) Returns the Bessel function n(x)for any real x and n 2. .com or call { (including this one) granted for i float bessio(float x); float bessj1(float x); void nrerror(char error_text []) 11-800-872 int j,jsum,m; float ax,bj,bim,bip,sum,tox,ans; from NUMERICAL RECIPES IN C: if (n <2)nrerror("Index n less than 2 in bessj"); ax-fabs(x); if(ax==0.0) return 0.0; 7423 (North America to any server computer,is tusers to make one paper 1988-1992 by Cambridge University Press. THE else if (ax >(float)n){ Upwards recurrence from Jo and J1. tox=2.0/ax; bim=bessi0(ax); only). bj=bessj1(ax); Programs for (j=1;j<n;i++){ bjp=j*tox*bj-bjm; send copy for their bjm=bj; st st bj=bjpi to dir Copyright (C) ans=bj; else Downwards recurrence from an even m here com- tox=2.0/ax: puted. m=2*((n+(int)sgrt(ACC*n))/2); rectcustser jsum=0; jsum will alternate between 0 and 1;when it is ART OF SCIENTIFIC COMPUTING(ISBN 0-521 bjp=ans=sum=0.0 1,we accumulate in sum the even terms in bj=1.0: (5.5.16) for (j=m;j>0;j--){ The downward recurrence. bjm=j*tox*bj-bjp; bjp-bj; bi=bim; idge.org 之 1988-1992 by Numerical Recipes -431085 if (fabs(bj)>BIGNO) Renormalize to prevent overflows. bj *BIGNI; bjp米=BIGNI: ans米=BIGNI; sum *BIGNI; (outside North Amer Software. if (jsum)sum +=bj; Accumulate the sum. jsum=!jsum; Change 0 to 1 or vice versa. if (j ==n)ans=bjp; Save the unnormalized answer sum=2.0+sum-bj; Compute (5.5.16) ans /sum; and use it to normalize the answer return x <0.0&(n&1)?-ans ans;6.5 Bessel Functions of Integer Order 235 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). an additive amount of order [constant × n]1/2, where the square root of the constant is, very roughly, the number of significant figures of accuracy. The above considerations lead to the following function. #include <math.h> #define ACC 40.0 Make larger to increase accuracy. #define BIGNO 1.0e10 #define BIGNI 1.0e-10 float bessj(int n, float x) Returns the Bessel function Jn(x) for any real x and n ≥ 2. { float bessj0(float x); float bessj1(float x); void nrerror(char error_text[]); int j,jsum,m; float ax,bj,bjm,bjp,sum,tox,ans; if (n < 2) nrerror("Index n less than 2 in bessj"); ax=fabs(x); if (ax == 0.0) return 0.0; else if (ax > (float) n) { Upwards recurrence from J0 and J1. tox=2.0/ax; bjm=bessj0(ax); bj=bessj1(ax); for (j=1;j<n;j++) { bjp=j*tox*bj-bjm; bjm=bj; bj=bjp; } ans=bj; } else { Downwards recurrence from an even m here com￾tox=2.0/ax; puted. m=2*((n+(int) sqrt(ACC*n))/2); jsum=0; jsum will alternate between 0 and 1; when it is 1, we accumulate in sum the even terms in (5.5.16). bjp=ans=sum=0.0; bj=1.0; for (j=m;j>0;j--) { The downward recurrence. bjm=j*tox*bj-bjp; bjp=bj; bj=bjm; if (fabs(bj) > BIGNO) { Renormalize to prevent overflows. bj *= BIGNI; bjp *= BIGNI; ans *= BIGNI; sum *= BIGNI; } if (jsum) sum += bj; Accumulate the sum. jsum=!jsum; Change 0 to 1 or vice versa. if (j == n) ans=bjp; Save the unnormalized answer. } sum=2.0*sum-bj; Compute (5.5.16) ans /= sum; and use it to normalize the answer. } return x < 0.0 && (n & 1) ? -ans : ans; }