c $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/spline.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      subroutine spline(xd,yd,n,yp1,ypn,y2,mode)
c     author:  c seftor
c     date:    3 feb 1991
c     purpose: given arrays xd and yd of length n containing a tabulated
c      function, and given values yp1 and ypn for the first derivative of
c      the interpolating function at points 1 and n, respectively, this
c      routinereturns an array y2 of length n which contains the second
c      derivatives of the interpolating function at the tabulated value
c      points xd.  if yp1 and/or ypn are equal to 1 x 10**30 or larger,
c      routine is signalled to set the corresponding boundary condition for
c      a natural spline, with zero second derivative on the boundary.
c      (algorithm taken from numerical recipes (fortran), 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
c
      implicit integer*4(i-n), real*8(a-h,o-z)
      parameter (nmax=487)
      real*8 xd(n), yd(n), y2(n), u(nmax)
      real*8 x(nmax),y(nmax)
      real*8 yp1, ypn,dr
      integer*4 n,mode
c
      if(n.gt.nmax)then
	write(6,*)'spline:  n gt nmax. will stop now'
        stop
      endif
c
c
c check to see if can do the modes greater than 0
      if(mode.gt.0)then
        if(mode.ge.8.and.mode.le.10)dr=acos(-1.0d0)/180.0d0
        if(mode.eq.1.or.mode.eq.3)then
          do i=1,n
             if(xd(i).le.0.0)then
               mode=mode-1
               goto 100
             endif
          enddo
        endif
100     if((mode.eq.2).or.(mode.eq.3).or.(mode.eq.10))then
          do i=1,n
             if(yd(i).le.0.0)then
               if(mode.eq.2)then
                 mode=0
               else if(mode.eq.3)then
                 mode=1
               else if(mode.eq.10)then
                 mode=9
               endif
               goto 200
             endif
          enddo
        endif
200     continue
      endif
c
c load the x and y arrays for the computations
c
      if(mode.eq.0.or.mode.eq.2)then
        do i=1,n
           x(i)=xd(i)
        enddo
      else if(mode.eq.1.or.mode.eq.3)then
        do i=1,n
           x(i)=log(xd(i))
        enddo
      else if(mode.eq.8)then
        do i=1,n
           x(i)=sin(dr*xd(i))
        enddo
      else if(mode.eq.9.or.mode.eq.10)then
        do i=1,n
           x(i)=cos(dr*xd(i))
        enddo
      endif
      if(mode.lt.2.or.mode.eq.8.or.mode.eq.9)then
        do i=1,n
           y(i)=yd(i)
        enddo
      else
        do i=1,n
           y(i)=log(yd(i))
        enddo
      endif
      if (yp1 .gt. .99d30) then
         y2(1) = 0.
         u(1)  = 0.
      else
         y2(1) = -0.5
         u(1)  = (3./(x(2)-x(1))) * ((y(2)-y(1))/(x(2)-x(1)) - yp1)
      endif
      do 11 i = 2,n-1
         sig = (x(i)-x(i-1))/(x(i+1)-x(i-1))
         p = sig * y2(i-1) + 2.
         y2(i) = (sig-1.)/p
         u(i) = (6.*((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
11    continue
      if (ypn .gt. .99d30) then
         qn = 0.
         un = 0.
      else
         qn = 0.5
         un = (3./(x(n)-x(n-1))) * (ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n) = (un - qn*u(n-1))/(qn*y2(n-1)+1.)
      do 12 k = n-1,1,-1
         y2(k) = y2(k) * y2(k+1) + u(k)
12    continue
      return
      end