PROGRAM ZEROBJ
      IMPLICIT NONE
      DOUBLE PRECISION a,b,pnu,eps,xcer(1000)
      INTEGER i,num,j,k,itc(1000)
      write(*,*)'a,b,nu,eps'
      read(*,*)a,b,pnu,eps
      CALL besjz(a,b,pnu,eps,xcer,itc,num)
      write(*,30)'Zeros','Iterations','Number of zeros'
      DO i=1,num
         write(*,31)xcer(i),itc(i),num
      ENDDO
 30   format(a8,16x,a10,5x,a16)
 31   format(d22.16,5x,i5,15x,i5)
c     write(*,*)'         Zero             Iterations   Number of zeros'
c     do i=1,num
c        write(*,*)xcer(i),itc(i),num
c     enddo
      END                                                                              
         
      SUBROUTINE besjz(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,atast,hm,ss,sa
      DOUBLE PRECISION dev,dest,pnu2,omegam,omegab,ffc,fc,ncc(1000)
      INTEGER j,i,num,k,apnu2,iter,isal,itc(1000)
      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)
      if(a.le.1.d-8) a=1.d-8
      if(pnu.ge.0.d0)a=max(max(pnu,a),2.4d0)
      pnu2=pnu*pnu-0.25d0
      if (pnu2.ne.0)then 
         j=int(pnu2/abs(pnu2))
      else
         j=1
      endif 
         dev=1.D0
         if(j.eq.1)then
            xm=a
            xc=b
         else
            xm=b
            xc=a
         endif
         i=0
         omega2=1.d0-(pnu*pnu-0.25d0)/xc/xc
         omegam=1.d0-(pnu*pnu-0.25d0)/xm/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
            h=1
            DO WHILE ((dev.gt.eps).and.(j*(xc-xm).gt.0).AND.(k.GE.0).
     c        and.(abs(h)/xc.gt.1.d-19))
               IF(iter.GT.1000)THEN
                  WRITE(*,*)'convergence failure'
                  STOP
               ENDIF
               ffc=fc(xc,pnu+1.d0,eps)
               h=1/((pnu+0.5d0)/xc-ffc)
               omega=dsqrt(omega2)
               atast=datan(omega*h)
               dest=atast/omega
               dev=ABS(dest)/xc
               if(j*dest.lt.0) then
                  if(dev.gt.1.d-10) dest=dest+j*dpi/omega
               endif
               if (isal.eq.1) dev=0.5*eps
               if (log(dev).lt.0.25d0*log(eps)) isal=1
               xc=xc-dest
               iter=iter+1
               omega2=1.d0-(pnu*pnu-0.25d0)/xc/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
               ffc=fc(xc,pnu+1.d0,eps)
               h=1/((pnu+0.5d0)/xc-ffc)
               xc=xc-j*dpi/omega
               omega2=1.d0-(pnu*pnu-0.25d0)/xc/xc 
               IF(omega2.LE.0.D0) k=-1
            ENDIF
         ENDDO

         if(k.lt.0) then
            ffc=fc(xm,pnu+1.d0,eps)
            hm=1/((pnu+0.5d0)/xm-ffc)
            ffc=fc(xc,pnu+1.d0,eps)
            h=1/((pnu+0.5d0)/xc-ffc)
            if((j*h.gt.0).and.(h*hm.lt.0))then
               iter=0
               dev=eps+1
               omega=sqrt(-omega2)
               isal=0
               DO WHILE ((dev.gt.eps).and.(abs(h*omega).lt.1).and.
     c                                          (j*(xc-xm).gt.0))
                  iter=iter+1 
                  IF(iter.GT.1000)THEN
                     WRITE(*,*)'convergence failure'
                     STOP
                  ENDIF
                  if(abs(h*omega).lt.1)then
                     atast=0.5*log((1+omega*h)/(1-omega*h))
                     dest=atast/omega
                     dev=ABS(dest)/xc
                     if (isal.eq.1) dev=0.5d0*eps
                     if (log(dev).lt.0.25d0*log(eps)) isal=1
                     xc=xc-dest
                     ffc=fc(xc,pnu+1.d0,eps)
                     omega2=1.d0-(pnu*pnu-0.25d0)/xc/xc
                     h=1/((pnu+0.5d0)/xc-ffc)
                     omega=dsqrt(-omega2)
                  endif
               ENDDO
               IF (dev.lt.eps) THEN
                  i=i+1
                  xcer(i)=xc
                  itc(i)=iter
               ENDIF
            endif
         endif
         num=i

         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