6.7 Bessel Functions of Fractional Order 243 Internal arithmetic in the routine is carried out in double precision.The complex arithmetic is carried out explicitly with real variables. #include <math.h> #include "nrutil.h" #define EPS 1.0e-10 #define FPMIN 1.0e-30 #define MAXIT 10000 #define XMIN 2.0 #def1nePI3.141592653589793 http://www. void bessjy(float x,float xnu,float *rj,float *ry,float *rjp,float *ryp) Returns the Bessel functions rj =Jv.ry =Yo and their derivatives rjp=J ryp =Y,for positive x and for xnu=>0.The relative accuracy is within one or two significant digits of EPS,except near a zero of one of the functions,where EPS controls its absolute accuracy. 8 FPMIN is a number close to the machine's smallest floating-point number.All internal arithmetic 虽 granted for 18881892 is in double precision.To convert the entire routine to double precision,change the float declarations above to double and decrease EPS to 10-16.Also convert the function beschb 1800 void beschb(double x,double *gam1,double *gam2,double *gampl, double *gammi); from NUMERICAL RECIPESI int i,isign,1,nl; 可 double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2, fact3,ff,gam,gami,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,rjl, (Nort rjl1,rjmu,rjpi,rjpl,rjtemp,ry1,rymu,rymup,rytemp,sum,sumi, temp,日，x2,x1,x12,Xmu,Xmu2; America server computer, one paper University Press. THE if (x <=0.0 II xnu 0.0)nrerror("bad arguments in bessjy"); ART nl=(x XMIN (int)(xnu+0.5):IMAX(0,(int)(xnu-x+1.5))) nl is the number of downward recurrences of the J's and upward recurrences of Y's.xmu Programs lies between -1/2 and 1/2 for x XMIN,while it is chosen so that x is greater than the turning point for x XMIN. xmuexnu-nl: xmu2=xmu*xmu; x1=1.0/x; x12=2.0*x1: to dir w=xi2/PI; The Wronskian. 1s1gn=1; Evaluate CF1 by modified Lentz's method($5.2). h=xnu*xi; isign keeps track of sign changes in the de- if (h FPMIN)h-FPMIN; nominator. b=xi2*xnu; d=0.0; c=h; for (i=1;i<=MAXIT;i++){ b+=x12; 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN 10-621 d=b-d; -43108 if (fabs(d)<FPMIN)d=FPMIN: c=b-1.0/c: if (fabs(c)<FPMIN)c=FPMIN; d=1.0/d; Software. del=c*d: h=del*h; (outside North Amer 1f(d<0.0)isig即=-1s1g即； if (fabs(del-1.0)<EPS)break; if (i MAXIT)nrerror("x too large in bessjy;try asymptotic expansion"); rjl-isign*FPMIN; InitializeJ and J for downward recurrence. riplsh*ril: rjli=rjl; Store values for later rescaling. ripi=ripl; fact=xnu*xi; for(1=n1;1>=1;1--){ rjtemp=fact*rjl+rjpl; fact -xi;6.7 Bessel Functions of Fractional Order 243 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). Internal arithmetic in the routine is carried out in double precision. The complex arithmetic is carried out explicitly with real variables. #include <math.h> #include "nrutil.h" #define EPS 1.0e-10 #define FPMIN 1.0e-30 #define MAXIT 10000 #define XMIN 2.0 #define PI 3.141592653589793 void bessjy(float x, float xnu, float *rj, float *ry, float *rjp, float *ryp) Returns the Bessel functions rj = Jν, ry = Yν and their derivatives rjp = J
ν, ryp = Y
ν , for positive x and for xnu = ν ≥ 0. The relative accuracy is within one or two significant digits of EPS, except near a zero of one of the functions, where EPS controls its absolute accuracy. FPMIN is a number close to the machine’s smallest floating-point number. All internal arithmetic is in double precision. To convert the entire routine to double precision, change the float declarations above to double and decrease EPS to 10−16. Also convert the function beschb. { void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi); int i,isign,l,nl; double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2, fact3,ff,gam,gam1,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,rjl, rjl1,rjmu,rjp1,rjpl,rjtemp,ry1,rymu,rymup,rytemp,sum,sum1, temp,w,x2,xi,xi2,xmu,xmu2; if (x <= 0.0 || xnu < 0.0) nrerror("bad arguments in bessjy"); nl=(x < XMIN ? (int)(xnu+0.5) : IMAX(0,(int)(xnu-x+1.5))); nl is the number of downward recurrences of the J’s and upward recurrences of Y ’s. xmu lies between −1/2 and 1/2 for x < XMIN, while it is chosen so that x is greater than the turning point for x ≥ XMIN. xmu=xnu-nl; xmu2=xmu*xmu; xi=1.0/x; xi2=2.0*xi; w=xi2/PI; The Wronskian. isign=1; Evaluate CF1 by modified Lentz’s method (§5.2). isign keeps track of sign changes in the de￾nominator. h=xnu*xi; if (h < FPMIN) h=FPMIN; b=xi2*xnu; d=0.0; c=h; for (i=1;i<=MAXIT;i++) { b += xi2; d=b-d; if (fabs(d) < FPMIN) d=FPMIN; c=b-1.0/c; if (fabs(c) < FPMIN) c=FPMIN; d=1.0/d; del=c*d; h=del*h; if (d < 0.0) isign = -isign; if (fabs(del-1.0) < EPS) break; } if (i > MAXIT) nrerror("x too large in bessjy; try asymptotic expansion"); rjl=isign*FPMIN; Initialize Jν and J
ν for downward recurrence. rjpl=h*rjl; rjl1=rjl; Store values for later rescaling. rjp1=rjpl; fact=xnu*xi; for (l=nl;l>=1;l--) { rjtemp=fact*rjl+rjpl; fact -= xi;