c $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/splint2.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      subroutine splint2(xa,ya,y2a,n,x,y,y1,mode)
c     author:  c. seftor
c     date:    3 feb 1991
c     purpose: given the arrays xa and ya of length n, which tabulate a
c      function (with the xa's in order), and given the array y2a,
c      which is the output from the subroutine spline, and given a value
c      of x, this routine returns a cubic-spline interpolated value y.
c      (algorithm taken from numerical recipes (fortran) by press et. al.)
c
c added mode swith for doing interpolation in;
c mode=0, x-y space
c     =1, ln(x)-y space
c     =2, x-ln(y) space
c     =3, ln(x)-ln(y) space
c If mode=4, perform linear interpolation instead of cubic in x-y space.
c If mode=5, perform linear interpolation instead of cubic in ln(x)-y space.
c If mode=6, perform linear interpolation instead of cubic in x-ln(y) space.
c If mode=7, perform linear interpolation instead of cubic in ln(x)-ln(y) space.
c     =8, sin(x)-y space
c     =9, cos(x)-y space
c     =10, cos(x)-ln(y) space
c
      implicit integer*4(i-n), real*8(a-h,o-z)
      real*8 xa(n), ya(n), y2a(n), x, y,y1
      real*8 a, b,dr
      integer*4 n,mode
c 
      if(mode.ge.8.and.mode.le.10)dr=asin(-1.0d0)/180.0d0
c check to see if the data is increasing in the independent variable
c or decreasing
      klo = 1
      khi = n
      if(xa(khi).gt.xa(klo))then
10      continue
        if (khi-klo .gt. 1) then
           k = (khi+klo)/2
           if (xa(k) .eq. x) then
             y=ya(k)
             goto 20
           else if (xa(k) .gt. x) then
              khi = k
           else
              klo = k
           endif
           goto 10
        endif
      else
15      continue
        if (khi-klo.gt. 1) then
           k = (khi+klo)/2
           if (xa(k) .eq. x) then
             y=ya(k)
             goto 20
           else if (xa(k) .gt. x) then
              klo = k
           else
              khi = k
           endif
           goto 15
        endif
      endif
      if((mode.eq.0).or.(mode.eq.2).or.(mode.eq.4).or.
     &(mode.eq.6))then
        h = xa(khi) - xa(klo)
      else if((mode.eq.1).or.(mode.eq.3).or.(mode.eq.5).or.
     &(mode.eq.7))then
        h = log(xa(khi)) - log(xa(klo))
      else if(mode.eq.8)then
        h = sin(xa(khi)*dr) - sin(xa(klo)*dr)
      else if((mode.eq.9).or.(mode.eq.10))then
        h = cos(xa(khi)*dr) - cos(xa(klo)*dr)
      endif
      if (h .eq. 0.) then
         write (23,*)' bad xa input'
         stop
      endif
      if((mode.eq.0).or.(mode.eq.2).or.(mode.eq.4).or.
     &(mode.eq.6))then
        a = (xa(khi)-x)/h
        b = (x-xa(klo))/h
      else if((mode.eq.1).or.(mode.eq.3).or.(mode.eq.5).or.
     &(mode.eq.7))then
        a = log(xa(khi)/x)/h
        b = log(x/xa(klo))/h
      else if(mode.eq.8)then
        a = (sin(xa(khi)*dr)-sin(x*dr))/h
        b = (sin(x*dr)-sin(xa(klo)*dr))/h
      else if((mode.eq.9).or.(mode.eq.10))then
        a = (cos(xa(khi)*dr)-cos(x*dr))/h
        b = (cos(x*dr)-cos(xa(klo)*dr))/h
      endif
      if(mode.eq.4.or.mode.eq.5)then
        y = a*ya(klo) + b*ya(khi)
      else if(mode.eq.6.or.mode.eq.7)then
        y = a*log(ya(klo)) + b*log(ya(khi))
      else if((mode.eq.2).or.(mode.eq.3).or.(mode.eq.10))then
        y = a*log(ya(klo)) + b*log(ya(khi)) +
     1 ((a**3-a)*y2a(klo) + (b**3-b)*y2a(khi)) * (h**2)/6.
      else
        y = a*ya(klo) + b*ya(khi) +
     1 ((a**3-a)*y2a(klo) + (b**3-b)*y2a(khi)) * (h**2)/6.
        y1 = (ya(khi)-ya(klo))/h-(3.*a*a-1.)/6.*h*y2a(klo)+
     1 (3.*b*b-1.)/6.*h*y2a(khi)
      endif
      if((mode.eq.2).or.(mode.eq.3).or.(mode.eq.6).or.
     &(mode.eq.7).or.(mode.eq.10))then
        y=exp(y)
      endif
20    return
      end