c     dtable verified 6-25-96 for corrected table answers using scifix
	integer*4 n,j,k,iscale
	real*8 scale,xr,xi
	complex*16 jn,y0,j0,x,yn,xscale
1	print *,'Verified 6-25-96'
    	PRINT *,'INPUT Final ORDER,ARGUMENT OF SPHERICAL BESSEL FUNCTION '
	READ *,N,XR,XI
	X=DCMPLX(XR,XI) 
	if(xi.gt.231) then 
c     dlog is the natural or naptherian log to base e, not 10
         scale=100.d0*dlog(10.0d0)
	   iscale=dint(xi/scale)
	   print *,'iscale,scale = ',iscale,scale
	   xi=xi-scale*iscale
	   iscale=iscale*100
c	   print *,'iscale= ',iscale  
	else
	   ISCALE=0
	   scale=0.0
	endif     
	xscale=DCMPLX(XR,XI)  
c********************************************************************
	if(x.eq.0) then
	    j0=1.0
	    y0=0.0
	  else    
	    j0=cdsin(xscale)/x
	    y0=-cdcos(xscale)/x
	  endif  
c	print *,xr,xi,x,xscale,cdsin(x)/x,-cdcos(x)/x    
	write(6,50)
50    format(1h ,'            First Kind     
     x Second Kind         Scale')
	j=0          
	iscale=0
	write(6,100) j,j0,y0,iscale
	do 77 j=1,n
	call cisub(j,xr,xi,jn,scale) 
	if(iscale.gt.dint(scale)) then 
c	    print *,dint(scale),iscale,' scale,iscale'  
	    j0=j0/10**(dint(scale)-iscale) 
	    y0=y0*10**(dint(scale)-iscale)  
c	    print *,j0,'  j0'
	endif    
c     the first time, scale j0 to prevent double scaling
      iscale=dint(scale)
	yn=jn/j0*y0-1.0d0/(x*x*j0)    
	if(j.le.20.or.j.eq.99) then
          write(6,100) j,jn,yn,iscale      
      else
	  k=j 
	  if(int((k/10)*10).eq.k) then 
	    write(6,100) j,jn,yn,iscale
	  end if      
	  k=j-1
        if(int(k/10)*10.eq.k) then
          write(6,100) j,jn,yn,iscale
        end if 
      end if
100	format(1h ,I3,2x,2(D16.10,1x),2(D16.10,1x),I4)
200	format(1h ,'YN=',2D18.10)
c	print *,'jn= ',jn,iscale
c	print *,'yn= ',yn,j 
	
	j0=jn
	y0=yn
77	continue
c	print *,'scale= ',scale
	goto 1
	end

c     scifix from cifix verified 6-25-96 to correct tables
	subroutine cisub(n,xr,xi,jn,scale)
c	CIFIX fixes the overflow for large imaginary arguments
c	caused by an exponent overflow
c       Program cipub.for=generalized continued fraction
c       for spherical Bessel functions of complex argument using 
c       algorithm improvement only for num = 0 or den=0 exactly
c       
c       TEST CASE #1 100,1 j=7.444727742E-190
c       TEST CASE #2 100,100 j=1.088047701
c	TEST CASE #3 2,SQRT(15) j=0.2868111794
c       ICOUNT  Counts number of iterations in CF product
c       N       Integer order of final Bessel function
c       NFLAG   Flag for error of numerator being zero
c       DFLAG   Flag for error of denominator being zero
c       V       Half integer order of spherical Bessel function
c       XINC    2/x  to avoid unnecessary recalculation
c       AN      Negative simple continued fraction An(x)
c       J0      Zero order spherical Bessel function
c       JN      Nth order spherical Bessel function
c       X       Argument of the Spherical Bessel function
c       NUM     Current numerator term multiplied into PDT
c       DEN     Current denominator term divided into PDT
c       SCALE   Power of 10 scale to prevent overflow of the 
c               exponent in the computer.  This will differ for 
c               various compilers
c***************************************************************************
	INTEGER*4 ICOUNT,N,ISCALE
	LOGICAL NFLAG,DFLAG
	REAL*8 V,SCALE,XR,XI
	COMPLEX*16 XINC,AN,J0,JN,X,xscale,NUM,DEN,PDT,ival
c***************************************************************************

	if(xi.lt.0) then 
	    xi=dabs(xi)
	    print *,'Complex conjugate taken to prevent overflow'  
      endif	    

c***************************************************************************
C	Power of 10 scaling to prevent overflow

	X=DCMPLX(XR,XI)
	if(xi.gt.231) then 
c     dlog is the natural or naptherian log to base e, not 10
         scale=100.d0*dlog(10.0d0)
	   iscale=dint(xi/scale)
	   print *,'iscale,scale = ',iscale,scale
	   xi=xi-scale*iscale
	   iscale=iscale*100
	   print *,'iscale= ',iscale
	else
	   ISCALE=0
	   scale=0.0
	endif
c
	xscale=DCMPLX(XR,XI)
c***************************************************************************
c
c       Setup the initial values of the running product,numerator
c       denominator, and counter for number of iterations
c
	PDT=1.0D0
c
c       if zero order is calculated, use sin(x)/x and skip rest
c
	IF(N.EQ.0) GOTO 20
c
	V=DFLOAT(N) + 0.5
	NUM=0.0
	DEN=0.0
	AN=1.0D0/X
	XINC=2.0/X
	ICOUNT=0
	NFLAG=.TRUE.
	DFLAG=.TRUE.
c
c       Each successive An = AN is generated by adding xinc
c
   17   AN=AN+XINC
	NUM=AN-NUM
c
c       If NFLAGis true, the normal loop is used since no subtraction 
c       errors have previously occurred
c
	IF (NFLAG) THEN
c
c          If  the numerator  is not zero,  then the  new continued 
c          fraction product is made by multiplying the numerator
c          PDT by the next term AN
c
	  IF(NUM.NE.0) THEN
	    PDT=PDT*NUM
c
c          Prepare to generate the next term An + 1/[previous]
c          by taking the reciprocal 1/[previous]
c
	    NUM=1.0D0/NUM
	  ELSE
c           If the numerator is zero, set a flag so that the next
c           numerator product will be Am+2, skipping Am+1 since Beta   
c           is zero in equation (11)
c
	    NFLAG=.FALSE.
	    NUM=0.0D0
c
c	   If you forget to set the product negative, the final
c	   answer will have the wrong sign when NUM=0.0
c
	    PDT=-PDT
	  ENDIF
	ELSE
	  NFLAG=.TRUE.
	  NUM=0.0D0
	ENDIF
c
c       When the continued fraction product gets too big, scale it
c       by a factor of 10^-200 and keep an integer SCALE to printout
c
	IF(CDABS(PDT).LT.1.0D100) then
	  continue
	ELSE
	  PDT=PDT*1.0D-100
	  scale=scale-100
	ENDIF  
	ICOUNT=ICOUNT+1
	IF(ICOUNT.LE.N) GOTO 19
	DEN=AN-DEN
	IF(DFLAG) THEN
	  IF(DEN.NE.0) THEN
	    PDT=PDT/DEN
	    DEN=1.0D0/DEN
  	  ELSE
	    DFLAG=.FALSE.
	    DEN=0.0
	    PDT=-PDT
	  ENDIF
	
	ELSE
	  DFLAG=.TRUE.
	  DEN=0.0
	ENDIF
c
c       Force at least five terms in the continued fraction to
c       avoid false convergences for small Am + 1/[previious]  
c
   19   IF(NFLAG.AND.DFLAG) THEN
   	  IF(NUM.NE.DEN.OR.ICOUNT.LT.5) GOTO 17
   	ELSE
   	  GOTO 17
   	ENDIF
c
c	initial normalizing value or n=0 value
c
   20	J0=CDSIN(xscale)/X
c	or one could use in some compilers
c	J0=(cdexp(ival*xscale)-cdexp(-ival*xscale))/(2.0d0*ival*x)
c	if xscale is replaced by x at this point, the system
c	will overflow for large imaginary arguments because of
c	an exponent error
	ival=dcmplx(0,1)
	JN=J0/PDT
      return
	END
	