PROGRAM BESSELJD IMPLICIT NONE DOUBLE PRECISION a,b,pnu,eps,xcert(1000),xcmal(3) INTEGER i,numt,j,k,com,mode,itmal(3),ncc(1000), c itct(1000) c choose mode=1 if the values are to be tested for c possible loss of accuracy. This may be only c needed for negative orders smaller than -1 c For positive order mode=0 saves time mode=1 c write(*,*)'a,b,nu,eps' c read(*,*)a,b,pnu,eps a=0.d0 b=20.d0 pnu=-1.1171230d0 eps=1.d-14 CALL bessjd(a,b,pnu,eps,xcert,xcmal,itct, c itmal,ncc,numt,com,mode) write(*,30)'Zeros','Iterations','Number of zeros' DO i=1,numt write(*,31)xcert(i),itct(i),numt ENDDO write(*,*)' ' IF(com.gt.0) write(*,32)'Zeros with acc. loss', c 'Iterations','Digits lost' DO i=1,com write(*,31)xcmal(i),itmal(i),ncc(com) ENDDO 30 format(a8,16x,a10,5x,a16) 31 format(d22.16,5x,i5,15x,i5) 32 format(a20,4x,a10,2x,a16) END SUBROUTINE bessjd(a,b,pnu,eps, c xcert,xcmal,itct,itmal,ncc,numt,com,mode) c c Warning: this is a preliminary version of the program c Additional error handling has to be incorporated c IMPLICIT NONE DOUBLE PRECISION a,b,pnu,eps,xcer(1000),itc(1000) DOUBLE PRECISION xcert(1000),li(1:7),ui(1:7), c nuga,xcmal(3),hfu,ttt,desv INTEGER i,j,co,cot,numt,signo(1:7),num,com,itmal(3), c itct(1000),ncc(1000),mode CALL intervals(a,b,pnu,co,li,ui,signo) numt=0 do i=1,co num=0 if(li(i).lt.ui(i))then CALL besjz(signo(i),li(i),ui(i),pnu, c eps,xcer,itc,num) do j=numt+1,num+numt xcert(j)=xcer(j-numt) itct(j)=itc(j-numt) enddo endif numt=numt+num enddo com=0 i=0 if(mode.eq.1)then do while(i.lt.numt) i=i+1 c write(*,*)i,xcert(i), c c abs(xcert(i)**2-pnu*2)*hfu(pnu,xcert(i),eps) ttt=abs(xcert(i)**2-pnu*2)* c abs(hfu(pnu,xcert(i),eps)) if(ttt.gt.eps)then com=com+1 xcmal(com)=xcert(i) itmal(com)=itct(i) ncc(com)=int((log(ttt)-log(eps))/log(10.d0)) numt=numt do j=i,numt xcert(j)=xcert(j+1) itct(j)=itct(j+1) enddo i=i-1 numt=numt-1 endif enddo endif RETURN END SUBROUTINE besjz(j,a,b,pnu,eps,xcer,itc,num) IMPLICIT NONE DOUBLE PRECISION a,b,xm,xc,xcer(1000),pnu,eps,delta,omega DOUBLE PRECISION h,omega2,dpi,dpi2,itc(1000),atast,hm,ss,sa DOUBLE PRECISION dev,dest,omegam,omegab,ffc,fc,coef DOUBLE PRECISION omep,hfu,singu,l1,ht1,ht2,si,hfut,sip INTEGER j,i,num,k,iter,isal dpi=dacos(-1.d0) if(a.ge.b)then write(*,*)'a should be smaller than b' STOP endif if(int(pnu).eq.pnu)pnu=abs(pnu) dev=1.d0 if(j.eq.1)then xm=a xc=b else xm=b xc=a endif i=0 omega2=coef(pnu,xc) omegam=coef(pnu,xm) if (j.gt.0) then omegab=omega2 else omegab=omegam endif ss=omega2*omegam if (ss.lt.0) then k=0 elseif (omegab.le.0) then k=-1 else k=1 end if DO WHILE((j*(xc-xm).gt.0).AND.(k.GE.0)) iter=0 dev=eps+1 isal=0 hfut=1 sip=1 DO WHILE ((dev.gt.eps).and.(j*(xc-xm).gt.0).AND.(k.GE.0). c and.(abs(sip*hfut).gt.1.d-19)) iter=iter+1 IF(iter.GT.1000)THEN WRITE(*,*)'convergence failure' STOP ENDIF sip=xc**2-pnu**2 si=sip/abs(sip) hfut=hfu(pnu,xc,eps) h=dsqrt(omep(pnu,xc))*si*hfut omega=dsqrt(omega2) atast=datan(h) dest=atast/omega dev=ABS(dest)/xc if (isal.eq.1) dev=0.5*eps if (log(dev).lt.0.25d0*log(eps)) isal=1 if(j*dest.lt.0) then if(dev.gt.1.d-10) dest=dest+j*dpi/omega endif xc=xc-dest omega2=coef(pnu,xc) IF(omega2.LE.0.D0) k=-1 ENDDO IF ((j*(xc-xm).gt.0).AND.(k.GE.0)) THEN i=i+1 xcer(i)=xc itc(i)=iter xc=xc-j*dpi/omega omega2=coef(pnu,xc) IF(omega2.LE.0.D0) k=-1 ENDIF ENDDO if(k.lt.0) then ht1=hfu(pnu,xm,eps) ht2=hfu(pnu,xc,eps) si=xc**2-pnu**2 si=si/abs(si) h=dsqrt(-omep(pnu,xc))*si*ht2 if((j*h.gt.0).and.(ht1*ht2.lt.0))then iter=0 dev=eps+1 omega=sqrt(-omega2) isal=0 sip=1 hfut=1 DO WHILE ((dev.gt.eps).and.(abs(h).lt.1).and. c (j*(xc-xm).gt.0).and.(abs(sip*hfut).gt.1.d-17)) iter=iter+1 IF(iter.GT.1000)THEN WRITE(*,*)'convergence failure' STOP ENDIF sip=xc**2-pnu**2 si=si/abs(si) hfut=hfu(pnu,xc,eps) h=dsqrt(-omep(pnu,xc))*si*hfut omega=dsqrt(-omega2) if(abs(h).lt.1)then atast=0.5d0*log((1.d0+h)/(1.d0-h)) dest=atast/omega dev=ABS(dest)/xc if (isal.eq.1) dev=0.5*eps if (log(dev).lt.0.5d0*log(eps)) isal=1 xc=xc-dest omega2=coef(pnu,xc) c write(*,*)'neg',xc,pnu,xc*fc(xc,pnu+1,eps) endif ENDDO IF (dev.lt.eps) THEN i=i+1 xcer(i)=xc itc(i)=iter ENDIF endif endif num=i RETURN END SUBROUTINE intervals(ai,bi,nu,co,li,ui,signo) implicit none double precision pi,d,lambda,g,c,q,r,delta,Ad,x,the,sqq,nu2 double precision x1,x2,x3,nu,si,si2,xt(1:7),alam,limi double precision li(1:7),ui(1:7),xpp,a,b,al,bl,ai,bi double precision PP,QQ,RR,sl,sr integer signo(1:7) integer i,mod,j,k,ind,co,nogo,ja,jb nogo=0 mod=0 nu2=nu*nu lambda=nu2-g**2 if(lambda.ne.0)then si=lambda/abs(lambda) pi=acos(-1.d0) d=(3./4-2*g+nu2) a=si*(2*g-3*nu2+23./4)/d b=3.*(nu2-0.25)/d c=si*(0.25-nu2)/d q=(a**2-3*b)/9. r=(2.*a**3-9.*a*b+27.*c)/54. delta=r**2-q**3 limi=max(mod*si,0.d0) i=0 if (delta.gt.0) then Ad=-r/abs(r)*(abs(r)+sqrt(delta))**(1./3.) x=(Ad+Q/Ad)-a/3. if (x.gt.limi) then i=i+1 xt(i)=x endif else the=acos(r/q**(3./2.)) sqq=sqrt(q) x1=-2.*sqq*cos(the/3.)-a/3. if(x1.gt.limi)then i=i+1 xt(i)=x1 endif x2=-2.*sqq*cos((the+2*pi)/3.)-a/3. if(x2.gt.limi)then i=i+1 xt(i)=x2 endif x3=-2*sqq*cos((the-2*pi)/3.)-a/3. if(x3.gt.limi)then i=i+1 xt(i)=x3 endif endif alam=abs(lambda) al=ai**2/alam bl=bi**2/alam if(lambda.gt.0)then if(mod.eq.0)then i=i+1 xt(i)=1.d0 i=i+1 xt(i)=1.d0 elseif(al.lt.1.d0) then if(bl.gt.1.d0)then al=1.d0 else nogo=1 endif endif endif i=i+1 xt(i)=al i=i+1 xt(i)=bl if(nogo.eq.0)then do ind=i,2,-1 do k=1,ind-1 if(xt(k).gt.xt(k+1)) then xpp=xt(k) xt(k)=xt(k+1) xt(k+1)=xpp endif enddo enddo ja=1 if(xt(1).lt.al)then do while(xt(ja).lt.al) ja=ja+1 enddo else ja=1 endif jb=ja if(xt(i).gt.bl)then do while(xt(jb).lt.bl) jb=jb+1 enddo else jb=i endif co=0 do k=ja,jb-1 if(abs(1.d0-xt(k)/xt(k+1)).gt.1.d-8)then co=co+1 li(co)=xt(k) ui(co)=xt(k+1) endif enddo do j=1,co x=(li(j)+ui(j))/2 si2=(x-si)*d*(x**3+a*x**2+b*x+c) signo(j)=int(si2/abs(si2)) li(j)=sqrt(alam*li(j)) ui(j)=sqrt(alam*ui(j)) enddo PP=-3.d0*nu2+2.d0*g*(g+1)-0.75d0 QQ=(3.d0*nu2-g*(g+2.d0)-2.5d0)*lambda RR=(0.25d0-nu2)*lambda**2 do j=1,co x=li(j)**2 sl=x**3+PP*x**2+QQ*x+RR x=ui(j)**2 sr=x**3+PP*x**2+QQ*x+RR enddo endif else co=1 li(1)=ai ui(1)=bi signo(1)=0 if(nu.eq.g)then signo(1)=int((nu-1.d0)**2-0.25d0) else signo(1)=int((nu+1.d0)**2-0.25d0) endif if(signo(1).ne.0) then signo(1)=signo(1)/abs(signo(1)) else signo(1)=1 endif endif if(li(1).lt.1.d-6) li(1)=1.d-6 c write(*,*)'The intervals are' c do j=1,co c write(*,*)li(j),ui(j),signo(j) c enddo c write(*,*)'-------------' c write(*,*)'The zeros are' return end DOUBLE PRECISION FUNCTION hfu(nu,x,eps) implicit none double precision nu,x,x2,r,m,n,n2,fc, c eps x2=x*x n2=nu*nu m=-x2**2+nu*(2*nu-0.5d0)*x2-nu*n2*(nu+0.5d0) n=x*0.5d0*(x2+n2) r=fc(x,nu+1.d0,eps) hfu=(nu-x*r)/(m+n*r) RETURN END DOUBLE PRECISION FUNCTION coef(nu,x) implicit none double precision nu,x,n2,x2,omep x2=x*x n2=nu*nu coef=omep(nu,x)/x2/(x2-n2)**2 RETURN END DOUBLE PRECISION FUNCTION omep(nu,x) implicit none double precision nu,x,P,Q,R,x2,n2 n2=nu*nu x2=x*x P=(-3*n2-0.75d0) Q=(3*n2-2.5d0)*n2 R=(0.25d0-n2)*n2**2 omep=x2**3+P*x2**2+Q*x2+R RETURN END DOUBLE PRECISION FUNCTION fc(xc,pnu,eps) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cc Evaluation of the cf for the ratio Jpnu(x)/Jpnu-1(x) c cc We use Lentz-Thompson algorithm. c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc IMPLICIT NONE DOUBLE PRECISION xc,pnu,eps,tiny,tinysq,b,a,c0,d0, * delta INTEGER m PARAMETER(tiny=1.D-290) tinysq=DSQRT(tiny) m=0 b=2.D0*pnu/xc a=1.D0 fc=tinysq c0=fc d0=0.D0 delta=0.D0 DO WHILE(ABS(delta-1.D0).GT.eps) d0=b+a*d0 IF(ABS(d0).LT.tinysq) d0=tinysq c0=b+a/c0 IF(ABS(c0).LT.tinysq) c0=tinysq d0=1.D0/d0 delta=c0*d0 fc=fc*delta m=m+1 IF(m.GT.100000000)THEN WRITE(*,*)'convergence of the cf fails' STOP ENDIF a=-1.D0 b=2.D0*(pnu+m)/xc ENDDO RETURN END