c..this file contains routines to interpolate/extrapolate & table searchers 
c.. 
c..routine spline computes the natural cubic spline coefficients    
c..routine splinth uses the spline coefficients and hunting to interpolate 
c..routine hunt finds a table entry given a guess for a nearby entry 




      subroutine spline(x,y,n,yp1,ypn,y2)    
      include 'implno.dek' 
c..  
c..given arrays x(n) and y(n) containing a tabulated function y=f(x),  
c..this routine returns the array y2(n) of spline coefficients. 
c..x must be a monotonically increasing  
c.. 
c..if yp1 and/or ypn are set larger than 1e30, the routine sets the boundary 
c..condtion for a natural spline (zero second derivative) at that boundary, 
c..otherwise yp1 and ypn are the first derivatives to use at the endpoints 
c.. 
c..this routine is only called once for any given x and y.   
c..  
c..declare   
      integer            n,i,k,nmax  
      parameter          (nmax=1000)  
      double precision   x(n),y(n),y2(n),yp1,ypn,u(nmax),sig,p,qn,un 


c..the lower boundary condition is set to be either natural 
      if (yp1 .gt. 1.0d30) then  
       y2(1) = 0.0d0 
       u(1)  = 0.0d0 


c..or it has a specified first derivative    
      else   
       y2(1) = -0.5d0 
       u(1)  = (3.0d0/(x(2)-x(1))) * ((y(2)-y(1))/(x(2)-x(1)) - yp1)  
      end if 


c..this is the decomposition loop of the tridiagonal algorithm. y2 and u 
c..are used for temporary storage of the decomposition factors.  
      do i=2,n-1  
       sig   = (x(i) - x(i-1))/(x(i+1) - x(i-1))    
       p     = sig*y2(i-1) + 2.0d0 
       y2(i) = (sig - 1.0d0)/p 
       u(i)  = (6.0d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) / 
     1         (x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p    
      enddo

c..the upper boundary condition is set to be either natural  
      if (ypn .gt. 1.0d30) then  
       qn = 0.0d0 
       un = 0.0d0 

c..or it has a specified first derivative    
      else   
       qn = 0.5d0   
       un = (3.0d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))  
      end if 

c..this is the backsubstitution loop of the tridiagonal algorithm    
      y2(n) = (un-qn*u(n-1)) / (qn*y2(n-1) + 1.0d0)   
      do k=n-1,1,-1   
       y2(k) = y2(k)*y2(k+1) + u(k) 
      enddo
      return 
      end    




      subroutine splinth(xa,ya,y2a,n,x,y,klo) 
      include 'implno.dek' 
c..  
c..given arrays xa and ya of length n, which tablulate a monotonic function, 
c..and given the array y2a, which is the output of spline (above), and given 
c..a value of x this routine returns a cubic spline interpolated value of y. 
c..  
c..declare   
      integer             n,klo,khi 
      double precision    xa(n),ya(n),y2a(n),x,y,h,a,b   

c..assume the last value was close; hunt it down  
      call hunt(xa,n,x,klo) 
      khi = klo + 1 

c..klo and khi now bracket the input value of x; the xa's must be distinct 
      h = xa(khi) - xa(klo)  
      if (h .eq. 0.0) stop 'bad xa input in routine splint'  

c..evaluate the cubic spline 
      a = (xa(khi) - x)/h    
      b = (x - xa(klo))/h    
      y = a*ya(klo) + b*ya(khi) + 
     1    ((a*a*a-a)*y2a(klo) + (b*b*b-b)*y2a(khi)) * (h*h)/6.0d0 
      return     
      end    






      subroutine hunt(xx,n,x,jlo) 
      include 'implno.dek' 
c.. 
c..given an array xx of length n and a value x, this routine returns a value 
c..jlo such that x is between xx(jlo) and xx(jlo+1). the array xx must be  
c..monotonic. j=0 or j=n indicates that x is out of range. on input jlo is a  
c..guess for the table entry, and should be used as input for the next hunt.  
c..one assumes the next table entry is relatively close to the old one. 
c.. 
c..declare 
      logical          ascnd 
      integer          n,jlo,jhi,inc,jm 
      double precision xx(n),x 

c..logical function; true if ascending table, false if descending
      ascnd = xx(n) .ge. xx(1) 

c..horrid initial guess; go to bisection immediatly 
      if (jlo.le.0 .or. jlo.gt.n) then 
       jlo = 0 
       jhi = n+1 
       goto 3 
      end if 

c..set the initial hunting increment 
      inc = 1 

c..this is the hunt up section 
      if (x.ge.xx(jlo) .eqv. ascnd) then 
1      jhi = jlo + inc 

c..if we are the end of the table then we are done hunting 
       if (jhi .gt. n) then 
        jhi = n + 1 

c..otherwise replace the lower limit, double the increment and try again 
       else if (x.ge.xx(jhi) .eqv. ascnd) then 
        jlo = jhi 
        inc = inc + inc 
        goto 1 
       end if 

c..this is the hunt down section 
      else 
       jhi = jlo 
2      jlo = jhi - inc 

c..if we are the end of the table then we are done hunting 
       if (jlo.lt.1) then 
        jlo = 0 

c..otherwise replace the upper limit, double the increment and try again 
       else if (x.lt.xx(jlo) .eqv. ascnd) then 
        jhi = jlo 
        inc = inc + inc 
        goto 2 
       end if 
      end if 

c..this is the final bisection phase 
3     if (jhi-jlo .eq. 1) then 
       if (x .eq. xx(n)) jlo = n - 1 
       if (x .eq. xx(1)) jlo = 1 
       return 
      end if 
      jm = (jhi+jlo)/2 
      if (x .ge. xx(jm) .eqv. ascnd) then 
       jlo = jm 
      else 
       jhi = jm 
      end if 
      goto 3 
      end