SUBROUTINE bessjy(x,xnu,rj,ry,rjp,ryp) INTEGER MAXIT REAL rj,rjp,ry,ryp,x,xnu,XMIN DOUBLE PRECISION EPS,FPMIN,PI PARAMETER (EPS=1.e-10,FPMIN=1.e-30,MAXIT=10000,XMIN=2., *PI=3.141592653589793d0) CU USES 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=int(xnu+.5d0) else nl=max(0,int(xnu-x+1.5d0)) endif xmu=xnu-nl xmu2=xmu*xmu xi=1.d0/x xi2=2.d0*xi w=xi2/PI isign=1 h=xnu*xi if(h.lt.FPMIN)h=FPMIN b=xi2*xnu d=0.d0 c=h do 11 i=1,MAXIT b=b+xi2 d=b-d if(abs(d).lt.FPMIN)d=FPMIN c=b-1.d0/c if(abs(c).lt.FPMIN)c=FPMIN d=1.d0/d del=c*d h=del*h if(d.lt.0.d0)isign=-isign if(abs(del-1.d0).lt.EPS)goto 1 11 continue pause 'x too large in bessjy; try asymptotic expansion' 1 continue rjl=isign*FPMIN rjpl=h*rjl rjl1=rjl rjp1=rjpl fact=xnu*xi do 12 l=nl,1,-1 rjtemp=fact*rjl+rjpl fact=fact-xi rjpl=fact*rjtemp-rjl rjl=rjtemp 12 continue if(rjl.eq.0.d0)rjl=EPS f=rjpl/rjl if(x.lt.XMIN) then x2=.5d0*x pimu=PI*xmu if(abs(pimu).lt.EPS)then fact=1.d0 else fact=pimu/sin(pimu) endif d=-log(x2) e=xmu*d if(abs(e).lt.EPS)then fact2=1.d0 else fact2=sinh(e)/e endif call beschb(xmu,gam1,gam2,gampl,gammi) ff=2.d0/PI*fact*(gam1*cosh(e)+gam2*fact2*d) e=exp(e) p=e/(gampl*PI) q=1.d0/(e*PI*gammi) pimu2=0.5d0*pimu if(abs(pimu2).lt.EPS)then fact3=1.d0 else fact3=sin(pimu2)/pimu2 endif r=PI*pimu2*fact3*fact3 c=1.d0 d=-x2*x2 sum=ff+r*q sum1=p do 13 i=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+del del1=c*p-i*del sum1=sum1+del1 if(abs(del).lt.(1.d0+abs(sum))*EPS)goto 2 13 continue pause 'bessy series failed to converge' 2 continue rymu=-sum ry1=-sum1*xi2 rymup=xmu*xi*rymu-ry1 rjmu=w/(rymup-f*rymu) else a=.25d0-xmu2 p=-.5d0*xi q=1.d0 br=2.d0*x bi=2.d0 fact=a*xi/(p*p+q*q) cr=br+q*fact ci=bi+p*fact den=br*br+bi*bi dr=br/den di=-bi/den dlr=cr*dr-ci*di dli=cr*di+ci*dr temp=p*dlr-q*dli q=p*dli+q*dlr p=temp do 14 i=2,MAXIT a=a+2*(i-1) bi=bi+2.d0 dr=a*dr+br di=a*di+bi if(abs(dr)+abs(di).lt.FPMIN)dr=FPMIN fact=a/(cr*cr+ci*ci) cr=br+cr*fact ci=bi-ci*fact if(abs(cr)+abs(ci).lt.FPMIN)cr=FPMIN den=dr*dr+di*di dr=dr/den di=-di/den dlr=cr*dr-ci*di dli=cr*di+ci*dr temp=p*dlr-q*dli q=p*dli+q*dlr p=temp if(abs(dlr-1.d0)+abs(dli).lt.EPS)goto 3 14 continue pause 'cf2 failed in bessjy' 3 continue gam=(p-f)/q 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 rj=rjl1*fact rjp=rjp1*fact do 15 i=1,nl rytemp=(xmu+i)*xi2*ry1-rymu rymu=ry1 ry1=rytemp 15 continue ry=rymu ryp=xnu*xi*rymu-ry1 return END C (C) Copr. 1986-92 Numerical Recipes Software .