f6-7

234 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).6.7 Bessel Functions of Fractional Order, Airy Functions, Spherical Bessel Functions Many algorithmshavebeen proposed forcomputing Besselfunctions offractionalorder numerically. Most of them are, in fact, not very good in practice. The routines given here arerather complicated, but they can be recommended wholeheartedly. OrdinaryBessel Functions Thebasic idea is Steed’s method , which was originally developed [1]for Coulomb wave functions. The method calculates Jν,J/prime ν,Yν, andY/prime νsimultaneously, and so involves four relations among these functions. Three of the relations come from two continued fractions,one of which is complex. The fourth is provided by the Wronskian relation W≡J νY/prime ν−YνJ/prime ν=2 πx(6.7.1 ) The first continued fraction, CF1, is defined by fν≡J/prime ν Jν=ν x−Jν+1 Jν =ν x−1 2(ν+1 )/x−1 2(ν+2 )/x−···(6.7.2 ) Youcaneasilyderiveitfromthethree-termrecurrencerelationforBesselfunctions: Startwith equation (6.5.6) and use equation (5.5.18). Forward evaluation of the continued fraction byone of the methods of §5.2 is essentially equivalent to backward recurrence of the recurrence relation. The rate of convergence of CF1 is determined by the position of the turning point x tp=/radicalbig ν(ν+1 )≈ν, beyond which the Bessel functions become oscillatory. If x<∼xtp, convergence isveryrapid. If x>∼xtp,theneachiterationofthecontinued fractioneffectively increasesνby one untilx<∼xtp; thereafter rapid convergence sets in. Thus the number of iterations of CF1 is of order xfor largex. In the routine bessjywe set the maximum allowed number of iterations to 10,000. For larger x, you can use the usual asymptotic expressions for Bessel functions. One can show that the sign of Jνis the same as the sign of the denominator of CF1 once it has converged. The complex continued fraction CF2 is defined by p+iq≡J/prime ν+iY/prime ν Jν+iYν=−1 2x+i+i x(1/2)2−ν2 2(x+i)+(3/2)2−ν2 2(x+2i)+··· (6.7.3 ) (We sketch the derivation of CF2 in the analogous case of modified Bessel functions in the next subsection.) This continued fraction converges rapidly for x>∼xtp, while convergence failsasx→0. Wehave toadopt aspecial method forsmall x,which wedescribe below. For xnot too small, we can ensure that x>∼xtpby a stable recurrence of JνandJ/prime νdownwards to a valueν=µ<∼x, thus yielding the ratio fµat this lower value of ν. This is the stable direction for the recurrence relation. The initial values for the recurrence are Jν=arbitrary,J/prime ν=fνJν, (6.7.4 ) with the sign of the arbitrary initial value of Jνchosen to be the sign of the denominator of CF1. Choosing theinitialvalue of Jνvery smallminimizesthepossibility ofoverflowduring the recurrence. The recurrence relations are Jν−1=ν xJν+J/prime ν J/prime ν−1=ν−1 xJν−1−Jν(6.7.5 ) 6.7BesselFunctionsofFractionalOrder 235Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Once CF2 has been evaluated at ν=µ, then with the Wronskian (6.7.1) we have enough relationstosolveforallfourquantities. Theformulasaresimplifiedbyintroducingthequantity γ≡p−fµ q(6.7.6 ) Then Jµ=±/parenleftbiggW q+γ(p−fµ)/parenrightbigg1/2 (6.7.7 ) J/prime µ=fµJµ (6.7.8 ) Yµ=γJµ (6.7.9 ) Y/prime µ=Yµ/parenleftbigg p+q γ/parenrightbigg (6.7.10 ) The sign ofJµin (6.7.7) is chosen to be the same as the sign of the initial Jνin (6.7.4). Onceallfourfunctionshavebeendeterminedatthevalue ν=µ,wecanfindthematthe original value of ν.F o rJνandJ/prime ν, simply scale the values in (6.7.4) by the ratio of (6.7.7) to the value found after applying the recurrence (6.7.5). The quantities YνandY/prime νcan be found by starting with the values in (6.7.9) and (6.7.10) and using the stable upwards recurrence Yν+1=2ν xYν−Yν−1 (6.7.11 ) together with the relation Y/prime ν=ν xYν−Yν+1 (6.7.12 ) Now turn to the case of small x, when CF2 is not suitable. Temme [2]has given a good method of evaluating YνandYν+1, and henceY/prime νfrom (6.7.12), by series expansions that accurately handle the singularity as x→0. The expansions work only for |ν|≤1/2, and so now the recurrence (6.7.5) is used to evaluate fνat a valueν=µin this interval. Then one calculates Jµfrom Jµ=W Y/primeµ−Yµfµ(6.7.13 ) andJ/prime µfrom(6.7.8). Thevaluesattheoriginalvalueof νaredeterminedbyscalingasbefore, and theY’s are recurred up as before. Temme’s series are Yν=−∞/summationdisplay k=0ckgkYν+1=−2 x∞/summationdisplay k=0ckhk (6.7.14 ) Here ck=(−x2/4)k k!(6.7.15 ) while the coefficients gkandhkare defined in terms of quantities pk,qk, andfkthat can be found by recursion: gk=fk+2 νsin2/parenleftBigνπ 2/parenrightBig qk hk=−kgk+pk pk=pk−1 k−ν qk=qk−1 k+ν fk=kfk−1+pk−1+qk−1 k2−ν2(6.7.16 ) 236 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).The initial values for the recurrences are p0=1 π/parenleftBigx 2/parenrightBig−ν Γ(1 +ν) q0=1 π/parenleftBigx 2/parenrightBigν Γ(1−ν) f0=2 πνπ sinνπ/bracketleftbigg coshσΓ1(ν)+sinhσ σln/parenleftbigg2 x/parenrightbigg Γ2(ν)/bracketrightbigg(6.7.17 ) with σ=νln/parenleftbigg2 x/parenrightbigg Γ1(ν)=1 2ν/bracketleftbigg1 Γ(1−ν)−1 Γ(1 +ν)/bracketrightbigg Γ2(ν)=1 2/bracketleftbigg1 Γ(1−ν)+1 Γ(1 +ν)/bracketrightbigg(6.7.18 ) The whole point of writing the formulas in this way is that the potential problems as ν→0 canbecontrolledbyevaluating νπ/sinνπ,sinhσ/σ,and Γ1carefully. Inparticular,Temme gives Chebyshev expansions for Γ1(ν)andΓ2(ν). We have rearranged his expansion for Γ1 tobeexplicitlyanevenseriesin νsothatwecanuseourroutine chebevasexplained in §5.8. The routine assumes ν≥0. For negativeνyou can use the reflection formulas J−ν=c o sνπJ ν−sinνπY ν Y−ν=s i nνπJ ν+c o sνπY ν(6.7.19 ) The routine also assumes x> 0.F o rx< 0the functions are in general complex, but expressible in terms of functions with x> 0.F o rx=0,Yνis singular. Internal arithmetic in the routine is carried out in double precision. To maintain portability, complex arithmetic has been recoded with real variables. SUBROUTINE bessjy(x,xnu,rj,ry,rjp,ryp) INTEGER MAXIT REAL rj,rjp,ry,ryp,x,xnu,XMINDOUBLE PRECISION EPS,FPMIN,PI PARAMETER (EPS=1.e-10,FPMIN=1.e-30,MAXIT=10000,XMIN=2., * PI=3.141592653589793d0) C USES beschb Returns the Bessel functions rj=Jν,ry=Yνand their derivatives rjp=J/prime ν,ryp=Y/prime ν, for positive xand 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 REAL declaration above and decrease EPS to10−16. Also convert the subroutine beschb . INTEGER i,isign,l,nl DOUBLE PRECISION 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.le.0..or.xnu.lt.0.) pause ’bad arguments in bessjy’ if(x.lt.XMIN)then nl is the number of downward recurrences of the J’s and upward recurrences of Y’s.xmu lies between −1/2and 1/2 for x<XMIN , while it is chosen so that xis greater than the turning point for x≥XMIN .nl=int(xnu+.5d0) else nl=max(0,int(xnu-x+1.5d0)) endifxmu=xnu-nl xmu2=xmu*xmu xi=1.d0/xxi2=2.d0*xi w=xi2/PI The Wronskian. 6.7BesselFunctionsofFractionalOrder 237Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).isign=1 Evaluate CF1 by modified Lentz’s method ( §5.2). isign keeps track of sign changes in the denominator. h=xnu*xi if(h.lt.FPMIN)h=FPMIN b=xi2*xnud=0.d0 c=h do 11i=1,MAXIT b=b+xi2d=b-d if(abs(d).lt.FPMIN)d=FPMIN c=b-1.d0/cif(abs(c).lt.FPMIN)c=FPMINd=1.d0/d del=c*d h=del*hif(d.lt.0.d0)isign=-isign if(abs(del-1.d0).lt.EPS)goto 1 enddo 11 pause ’x too large in bessjy; try asymptotic expansion’ 1 continue rjl=isign*FPMIN Initialize JνandJ/prime νfor downward recurrence. rjpl=h*rjlrjl1=rjl Store values for later rescaling. rjp1=rjpl fact=xnu*xido 12l=nl,1,-1 rjtemp=fact*rjl+rjpl fact=fact-xi rjpl=fact*rjtemp-rjlrjl=rjtemp enddo 12 if(rjl.eq.0.d0)rjl=EPS f=rjpl/rjl Now have unnormalized JµandJ/prime µ. if(x.lt.XMIN) then Use series. x2=.5d0*x pimu=PI*xmuif(abs(pimu).lt.EPS)then fact=1.d0 else fact=pimu/sin(pimu) endif d=-log(x2) e=xmu*dif(abs(e).lt.EPS)then fact2=1.d0 else fact2=sinh(e)/e endif call beschb(xmu,gam1,gam2,gampl,gammi) Chebyshev evaluation of Γ 1andΓ2. ff=2.d0/PI*fact*(gam1*cosh(e)+gam2*fact2*d) f0. e=exp(e)p=e/(gampl*PI) p 0. q=1.d0/(e*PI*gammi) q0. pimu2=0.5d0*pimuif(abs(pimu2).lt.EPS)then fact3=1.d0 else fact3=sin(pimu2)/pimu2 endif r=PI*pimu2*fact3*fact3 c=1.d0d=-x2*x2sum=ff+r*q sum1=p 238 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).do13i=1,MAXIT ff=(i*ff+p+q)/(i*i-xmu2) c=c*d/i p=p/(i-xmu)q=q/(i+xmu) del=c*(ff+r*q) sum=sum+deldel1=c*p-i*delsum1=sum1+del1 if(abs(del).lt.(1.d0+abs(sum))*EPS)goto 2 enddo 13 pause ’bessy series failed to converge’ 2 continue rymu=-sum ry1=-sum1*xi2rymup=xmu*xi*rymu-ry1 rjmu=w/(rymup-f*rymu) Equation (6.7.13). else Evaluate CF2 by modified Lentz’s method (§5.2). a=.25d0-xmu2 p=-.5d0*xi q=1.d0 br=2.d0*xbi=2.d0 fact=a*xi/(p*p+q*q) cr=br+q*factci=bi+p*factden=br*br+bi*bi dr=br/den di=-bi/dendlr=cr*dr-ci*di dli=cr*di+ci*dr temp=p*dlr-q*dliq=p*dli+q*dlrp=temp do 14i=2,MAXIT a=a+2*(i-1)bi=bi+2.d0dr=a*dr+br di=a*di+bi if(abs(dr)+abs(di).lt.FPMIN)dr=FPMINfact=a/(cr*cr+ci*ci) cr=br+cr*fact ci=bi-ci*factif(abs(cr)+abs(ci).lt.FPMIN)cr=FPMINden=dr*dr+di*di dr=dr/den di=-di/dendlr=cr*dr-ci*di dli=cr*di+ci*dr temp=p*dlr-q*dliq=p*dli+q*dlrp=temp if(abs(dlr-1.d0)+abs(dli).lt.EPS)goto 3 enddo 14 pause ’cf2 failed in bessjy’ 3 continue gam=(p-f)/q Equations (6.7.6) – (6.7.10). rjmu=sqrt(w/((p-f)*gam+q))rjmu=sign(rjmu,rjl) rymu=rjmu*gam rymup=rymu*(p+q/gam)ry1=xmu*xi*rymu-rymup endif fact=rjmu/rjl 6.7BesselFunctionsofFractionalOrder 239Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).rj=rjl1*fact Scale original JνandJ/prime ν. rjp=rjp1*fact do15i=1,nl Upward recurrence of Yν. rytemp=(xmu+i)*xi2*ry1-rymurymu=ry1 ry1=rytemp enddo 15 ry=rymuryp=xnu*xi*rymu-ry1 return END SUBROUTINE beschb(x,gam1,gam2,gampl,gammi) INTEGER NUSE1,NUSE2DOUBLE PRECISION gam1,gam2,gammi,gampl,x PARAMETER (NUSE1=5,NUSE2=5) C USES chebev Evaluates Γ1andΓ2by Chebyshev expansion for |x|≤1/2. Also returns 1/Γ(1 + x)and 1/Γ(1−x). If converting to double precision, set NUSE1 =7,NUSE2 =8. REAL xx,c1(7),c2(8),chebevSAVE c1,c2DATA c1/-1.142022680371168d0,6.5165112670737d-3, * 3.087090173086d-4,-3.4706269649d-6,6.9437664d-9, * 3.67795d-11,-1.356d-13/ DATA c2/1.843740587300905d0,-7.68528408447867d-2, * 1.2719271366546d-3,-4.9717367042d-6,-3.31261198d-8, * 2.423096d-10,-1.702d-13,-1.49d-15/ xx=8.d0*x*x-1.d0 Multiply xby 2 to make range be −1to 1, and then apply transformation for evaluating even Cheby- shev series.gam1=chebev(-1.,1.,c1,NUSE1,xx) gam2=chebev(-1.,1.,c2,NUSE2,xx) gampl=gam2-x*gam1gammi=gam2+x*gam1return END ModifiedBessel Functions Steed’s method does not work for modified Bessel functions because in this case CF2 is purely imaginary and we have only three relations among the four functions. Temme [3]has given a normalization condition that provides the fourth relation. The Wronskian relation is W≡IνK/prime ν−KνI/prime ν=−1 x(6.7.20 ) The continued fraction CF1 becomes fν≡I/prime ν Iν=ν x+1 2(ν+1 )/x+1 2(ν+2 )/x+··· (6.7.21 ) TogetCF2and thenormalization condition inaconvenient form,consider thesequence of confluent hypergeometric functions zn(x)=U(ν+1/2+n,2ν+1,2x)( 6.7.22 ) for fixedν. Then Kν(x)=π1/2(2x)νe−xz0(x)( 6.7.23 ) Kν+1(x) Kν(x)=1 x/bracketleftbigg ν+1 2+x+/parenleftbigg ν2−1 4/parenrightbiggz1 z0/bracketrightbigg (6.7.24 ) 240 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Equation (6.7.23) is the standard expression for Kνin terms of a confluent hypergeometric function, while equation (6.7.24) follows from relations between contiguous confluent hy-pergeometric functions (equations 13.4.16 and 13.4.18 in Abramowitz and Stegun). Nowthe functionsz nsatisfy the three-term recurrence relation (equation 13.4.15 in Abramowitz and Stegun) zn−1(x)=bnzn(x)+an+1zn+1 (6.7.25 ) with bn=2 (n+x) an+1=−[(n+1/2)2−ν2](6.7.26 ) Following the steps leading to equation (5.5.18), we get the continued fraction CF2 z1 z0=1 b1+a2 b2+··· (6.7.27 ) from which (6.7.24) gives Kν+1/Kνand thusK/prime ν/Kν. Temme’s normalization condition is that ∞/summationdisplay n=0Cnzn=/parenleftbigg1 2x/parenrightbiggν+1/2 (6.7.28 ) where Cn=(−1)n n!Γ(ν+1/2+n) Γ(ν+1/2−n)(6.7.29 ) Note that theCn’s can be determined by recursion: C0=1,C n+1=−an+1 n+1Cn (6.7.30 ) We use the condition (6.7.28) by finding S=∞/summationdisplay n=1Cnzn z0(6.7.31 ) Then z0=/parenleftbigg1 2x/parenrightbiggν+1/21 1+S(6.7.32 ) and (6.7.23) gives Kν. Thompson and Barnett [4]have given a clever method of doing the sum (6.7.31) simultaneously with the forward evaluation of the continued fraction CF2. Suppose thecontinued fraction is being evaluated as z 1 z0=∞/summationdisplay n=0∆hn (6.7.33 ) wheretheincrements ∆hnarebeingfound by,e.g.,Steed’salgorithmorthemodified Lentz’s algorithm of §5.2. Then the approximation to Skeeping the first Nterms can be found as SN=N/summationdisplay n=1Qn∆hn (6.7.34 ) Here Qn=n/summationdisplay k=1Ckqk (6.7.35 ) andqkis found by recursion from qk+1=(qk−1−bkqk)/ak+1 (6.7.36 ) starting withq0=0,q1=1. For the case at hand, approximately three times as many terms are needed to get Sto converge as are needed simply for CF2 to converge. 6.7BesselFunctionsofFractionalOrder 241Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).To findKνandKν+1for smallxwe use series analogous to (6.7.14): Kν=∞/summationdisplay k=0ckfkKν+1=2 x∞/summationdisplay k=0ckhk (6.7.37 ) Here ck=(x2/4)k k! hk=−kfk+pk pk=pk−1 k−ν qk=qk−1 k+ν fk=kfk−1+pk−1+qk−1 k2−ν2(6.7.38 ) The initial values for the recurrences are p0=1 2/parenleftBigx 2/parenrightBig−ν Γ(1 +ν) q0=1 2/parenleftBigx 2/parenrightBigν Γ(1−ν) f0=νπ sinνπ/bracketleftbigg coshσΓ1(ν)+sinhσ σln/parenleftbigg2 x/parenrightbigg Γ2(ν)/bracketrightbigg(6.7.39 ) Both the series for small x, and CF2 and the normalization relation (6.7.28) require |ν|≤1/2. In both cases,therefore, werecurse Iνdown to avalue ν=µinthisinterval,find Kµthere, and recurse Kνback up to the original value of ν. The routine assumes ν≥0. For negativeνuse the reflection formulas I−ν=Iν+2 πsin(νπ)Kν K−ν=Kν(6.7.40 ) Note that for large x,Iν∼ex,Kν∼e−x, and so these functions will overflow or underflow. It is often desirable to be able to compute the scaled quantities e−xIνandexKν. Simply omitting the factor e−xin equation (6.7.23) will ensure that all four quantities will have the appropriate scaling. If you also want to scale the four quantities for small xwhen the series in equation (6.7.37) are used, you must multiply each series by ex. SUBROUTINE bessik(x,xnu,ri,rk,rip,rkp) INTEGER MAXIT REAL ri,rip,rk,rkp,x,xnu,XMINDOUBLE PRECISION EPS,FPMIN,PI PARAMETER (EPS=1.e-10,FPMIN=1.e-30,MAXIT=10000,XMIN=2., * PI=3.141592653589793d0) C USES beschb Returns the modified Bessel functions ri=Iν,rk=Kνand their derivatives rip=I/prime ν, rkp =K/prime ν,f o rp o s i t i v e xand for xnu =ν≥0. The relative accuracy is within one or two significant digits of EPS.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 REAL declaration above and decrease EPS to10−16.A l s o convert the subroutine beschb . INTEGER i,l,nlDOUBLE PRECISION a,a1,b,c,d,del,del1,delh,dels,e,f,fact, * fact2,ff,gam1,gam2,gammi,gampl,h,p,pimu,q,q1,q2, * qnew,ril,ril1,rimu,rip1,ripl,ritemp,rk1,rkmu,rkmup,* rktemp,s,sum,sum1,x2,xi,xi2,xmu,xmu2 if(x.le.0..or.xnu.lt.0.) pause ’bad arguments in bessik’ 242 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).nl=int(xnu+.5d0) nl is the number of downward recurrences of the I’s and upward recurrences ofK’s.xmu lies between −1/2and 1/2.xmu=xnu-nl xmu2=xmu*xmu xi=1.d0/xxi2=2.d0*xi h=xnu*xi Evaluate CF1 by modified Lentz’s method (§5.2). if(h.lt.FPMIN)h=FPMIN b=xi2*xnud=0.d0 c=h do 11i=1,MAXIT b=b+xi2d=1.d0/(b+d) Denominators cannot be zero here, so no need for special precautions. c=b+1.d0/c del=c*dh=del*h if(abs(del-1.d0).lt.EPS)goto 1 enddo 11 pause ’x too large in bessik; try asymptotic expansion’ 1 continue ril=FPMIN Initialize IνandI/prime νfor downward recur- rence. ripl=h*ril ril1=ril Store values for later rescaling. rip1=ripl fact=xnu*xido 12l=nl,1,-1 ritemp=fact*ril+ripl fact=fact-xi ripl=fact*ritemp+rilril=ritemp enddo 12 f=ripl/ril Now have unnormalized IµandI/prime µ. if(x.lt.XMIN) then Use series. x2=.5d0*x pimu=PI*xmu if(abs(pimu).lt.EPS)then fact=1.d0 else fact=pimu/sin(pimu) endifd=-log(x2) e=xmu*d if(abs(e).lt.EPS)then fact2=1.d0 else fact2=sinh(e)/e endifcall beschb(xmu,gam1,gam2,gampl,gammi) Chebyshev evaluation of Γ 1andΓ2. ff=fact*(gam1*cosh(e)+gam2*fact2*d) f0. sum=ffe=exp(e)p=0.5d0*e/gampl p 0. q=0.5d0/(e*gammi) q0. c=1.d0d=x2*x2 sum1=p do 13i=1,MAXIT ff=(i*ff+p+q)/(i*i-xmu2)c=c*d/i p=p/(i-xmu) q=q/(i+xmu)del=c*ffsum=sum+del del1=c*(p-i*ff) 6.7BesselFunctionsofFractionalOrder 243Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).sum1=sum1+del1 if(abs(del).lt.abs(sum)*EPS)goto 2 enddo 13 pause ’bessk series failed to converge’ 2 continue rkmu=sum rk1=sum1*xi2 else Evaluate CF2 by Steed’s algorithm ( §5.2), which is OK because there can be no zero denominators.b=2.d0*(1.d0+x) d=1.d0/b delh=dh=delhq1=0.d0 Initializations for recurrence (6.7.35). q2=1.d0 a1=.25d0-xmu2c=a1 q=c First term in equation (6.7.34). a=-a1s=1.d0+q*delhdo 14i=2,MAXIT a=a-2*(i-1) c=-a*c/iqnew=(q1-b*q2)/a q1=q2 q2=qnewq=q+c*qnewb=b+2.d0 d=1.d0/(b+a*d) delh=(b*d-1.d0)*delhh=h+delh dels=q*delh s=s+delsif(abs(dels/s).lt.EPS)goto 3 Need only test convergence of sum since CF2 itself converges more quickly. enddo 14 pause ’bessik: failure to converge in cf2’ 3 continue h=a1*hrkmu=sqrt(PI/(2.d0*x))*exp(-x)/s Omit the factor exp(−x)to scale all the returned functions by exp(x)forx≥ XMIN .rk1=rkmu*(xmu+x+.5d0-h)*xi endifrkmup=xmu*xi*rkmu-rk1 rimu=xi/(f*rkmu-rkmup) GetI µfrom Wronskian. ri=(rimu*ril1)/ril Scale original IνandI/prime ν. rip=(rimu*rip1)/rildo 15i=1,nl Upward recurrence of Kν. rktemp=(xmu+i)*xi2*rk1+rkmu rkmu=rk1rk1=rktemp enddo 15 rk=rkmu rkp=xnu*xi*rkmu-rk1return END AiryFunctions For positivex, the Airy functions are defined by Ai(x)=1 π/radicalbigg x 3K1/3(z)( 6.7.41 ) 244 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Bi(x)=/radicalbigg x 3[I1/3(z)+I−1/3(z)] ( 6.7.42 ) where z=2 3x3/2(6.7.43 ) By using the reflection formula (6.7.40), we can convert (6.7.42) into the computationally more useful form Bi(x)=√x/bracketleftbigg2√ 3I1/3(z)+1 πK1/3(z)/bracketrightbigg (6.7.44 ) so that AiandBican be evaluated with a single call to bessik. The derivatives should not be evaluated by simply differentiating the above expressions because of possible subtraction errors near x=0. Instead, use the equivalent expressions Ai/prime(x)=−x π√ 3K2/3(z) Bi/prime(x)=x/bracketleftbigg2√ 3I2/3(z)+1 πK2/3(z)/bracketrightbigg (6.7.45 ) The corresponding formulas for negative arguments are Ai(−x)=√x 2/bracketleftbigg J1/3(z)−1√ 3Y1/3(z)/bracketrightbigg Bi(−x)=−√x 2/bracketleftbigg1√ 3J1/3(z)+Y1/3(z)/bracketrightbigg Ai/prime(−x)=x 2/bracketleftbigg J2/3(z)+1√ 3Y2/3(z)/bracketrightbigg Bi/prime(−x)=x 2/bracketleftbigg1√ 3J2/3(z)−Y2/3(z)/bracketrightbigg(6.7.46 ) SUBROUTINE airy(x,ai,bi,aip,bip) REAL ai,aip,bi,bip,x C USES bessik,bessjy Returns Airy functions Ai(x),Bi(x), and their derivatives Ai/prime(x),Bi/prime(x). REAL absx,ri,rip,rj,rjp,rk,rkp,rootx,ry,ryp,z, * PI,THIRD,TWOTHR,ONOVRT PARAMETER (PI=3.1415927,THIRD=1./3.,TWOTHR=2.*THIRD, * ONOVRT=.57735027) absx=abs(x)rootx=sqrt(absx) z=TWOTHR*absx*rootx if(x.gt.0.)then call bessik(z,THIRD,ri,rk,rip,rkp) ai=rootx*ONOVRT*rk/PI bi=rootx*(rk/PI+2.*ONOVRT*ri)call bessik(z,TWOTHR,ri,rk,rip,rkp)aip=-x*ONOVRT*rk/PI bip=x*(rk/PI+2.*ONOVRT*ri) else if(x.lt.0.)then call bessjy(z,THIRD,rj,ry,rjp,ryp) ai=.5*rootx*(rj-ONOVRT*ry) bi=-.5*rootx*(ry+ONOVRT*rj)call bessjy(z,TWOTHR,rj,ry,rjp,ryp)aip=.5*absx*(ONOVRT*ry+rj) bip=.5*absx*(ONOVRT*rj-ry) else Case x=0. ai=.35502805 bi=ai/ONOVRT 6.7BesselFunctionsofFractionalOrder 245Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).aip=-.25881940 bip=-aip/ONOVRT endif returnEND Spherical Bessel Functions For integern, spherical Bessel functions are defined by jn(x)=/radicalbigg π 2xJn+(1/2)(x) yn(x)=/radicalbigg π 2xYn+(1/2)(x)(6.7.47 ) They can be evaluated by a call to bessjy, and the derivatives can safely be found from the derivatives of equation (6.7.47). NotethatinthecontinuedfractionCF2in(6.7.3)justthefirsttermsurvivesfor ν=1/2. Thus one can make a very simple algorithm for spherical Bessel functions along the lines ofbessjybyalwaysrecursing j ndownton=0,settingpandqfromthefirstterminCF2,and then recursingynup. No special series is required near x=0. However, bessjyis already so efficientthatwehavenot bothered toprovide an independent routine forspherical Bessels. SUBROUTINE sphbes(n,x,sj,sy,sjp,syp) INTEGER nREAL sj,sjp,sy,syp,x C USES bessjy Returns spherical Bessel functions jn(x),yn(x), and their derivatives j/prime n(x),y/prime n(x)for integer n. REAL factor,order,rj,rjp,ry,ryp,RTPIO2 PARAMETER (RTPIO2=1.2533141)if(n.lt.0.or.x.le.0.)pause ’bad arguments in sphbes’order=n+0.5 call bessjy(x,order,rj,ry,rjp,ryp) factor=RTPIO2/sqrt(x)sj=factor*rjsy=factor*ry sjp=factor*rjp-sj/(2.*x) syp=factor*ryp-sy/(2.*x)return END CITED REFERENCES AND FURTHER READING: Barnett, A.R., Feng, D.H., Steed, J.W., and Goldfarb, L.J.B. 1974, Computer Physics Commu- nications, vol. 8, pp. 377–395. [1] Temme, N.M. 1976, Journal of Computational Physics , vol. 21, pp. 343–350 [2]; 1975, op. cit., vol. 19, pp. 324–337. [3] Thompson, I.J., and Barnett, A.R. 1987, Computer Physics Communications , vol. 47, pp. 245– 257. [4] Barnett, A.R. 1981, Computer Physics Communications , vol. 21, pp. 297–314. Thompson,I.J.,andBarnett,A.R.1986, JournalofComputationalPhysics ,vol.64,pp.490–509. Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), Chapter 10.