subroutine spline(x,y,n,yp1,ypn,y2)
c     author:  c seftor
c     date:    3 feb 1991
c     purpose: given arrays x and y 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 x.  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
      implicit integer*4(i-n), real*8(a-h,o-z)
      parameter (nmax=487)
      real*8 x(n), y(n), y2(n), u(nmax)
      real*8 yp1, ypn
      integer*4 n
      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