C Routines from the book Computer Methods for Mathematical
C Computations, by Forsythe, Malcolm, and Moler.
C Developed by George Forsythe, Mike Malcolm, and Cleve Moler.

C Taken from files spline.f and seval.f at https://netlib.org/fmm/,
C August 22th, 2022.

C Natural spline adapted from here: https://www.pdas.com/fmm2.html
      
C =======================================================================
C =======================================================================
      
      subroutine nspline(n, x, y, b, c, d)
      
      integer n
      double precision x(n), y(n), b(n), c(n), d(n)
c
c     Construct the natural spline thru a set of points (i.e. zero
c     second derivative at both endpoints.

      integer :: k

c     Straight line for the special case n < 3
      
      if ( n .lt. 3 ) then           
         b(1) = 0.0d0
         if ( n .eq. 2 ) then
            b(1) = ( y(2) - y(1) ) / ( x(2) - x(1) )
         end if
         c(1) = 0.0d0
         d(1) = 0.0d0
         b(2) = b(1)
         c(2) = 0.0d0
         d(2) = 0.0d0
         return
      end if

c     Put the h-array of the text into array d
      do k = 1,n-1
         d(k) = x(k+1) - x(k)
      end do

c     Set up the upper triangular system in locations 2 thru n-1 of
c     arrays b and c. B holds the diagonal and c the right hand side.

      b(2) = 2.0d0 * ( d(1) + d(2) )
      c(2) = ( y(3) - y(2) ) / d(2) - ( y(2) - y(1) ) / d(1)
      do k = 3,n-1
         b(k) = 2.0d0 * ( d(k-1) + d(k) ) - d(k-1) * d(k-1) / b(k-1)
         c(k) = ( y(k+1) - y(k) ) / d(k) - ( y(k) - y(k-1) ) / d(k-1)
     *          - d(k-1) * c(k-1) / b(k-1)
      end do

c     Back substitute to get c-array
      
      c(n-1) = c(n-1) / b(n-1)      
      do k = n-2,2,-1
         c(k) = ( c(k) - d(k) * c(k+1) ) / b(k)
      end do
c     c now holds the sigma array of the text      
      c(1) = 0.0d0
      c(n) = 0.0d0  

c     Compute polynomial coefficients
      
      b(n) = ( y(n) - y(n-1) ) / d(n-1)
     *     + d(n-1) * ( c(n-1) + c(n) + c(n) )
      do k = 1,n-1
         b(k) = ( y(k+1) - y(k) ) / d(k)
     *        - d(k) * ( c(k+1) + c(k) + c(k) )
         d(k) = ( c(k+1) - c(k) ) / d(k)
         c(k) = 3.0d0 * c(k)
      end do
      c(n) = 3.0d0 * c(n)
      d(n) = d(n-1)

      return
      end 

C =======================================================================
C =======================================================================
      
      subroutine spline (n, x, y, b, c, d)
      integer n
      double precision x(n), y(n), b(n), c(n), d(n)
c
c  the coefficients b(i), c(i), and d(i), i=1,2,...,n are computed
c  for a cubic interpolating spline
c
c    s(x) = y(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))**2 + d(i)*(x-x(i))**3
c
c    for  x(i) .le. x .le. x(i+1)
c
c  input..
c
c    n = the number of data points or knots (n.ge.2)
c    x = the abscissas of the knots in strictly increasing order
c    y = the ordinates of the knots
c
c  output..
c
c    b, c, d  = arrays of spline coefficients as defined above.
c
c  using  p  to denote differentiation,
c
c    y(i) = s(x(i))
c    b(i) = sp(x(i))
c    c(i) = spp(x(i))/2
c    d(i) = sppp(x(i))/6  (derivative from the right)
c
c  the accompanying function subprogram  seval  can be used
c  to evaluate the spline.
c
c
      integer nm1, ib, i
      double precision t
c
      nm1 = n-1
      if ( n .lt. 2 ) return
      if ( n .lt. 3 ) go to 50
c
c  set up tridiagonal system
c
c  b = diagonal, d = offdiagonal, c = right hand side.
c
      d(1) = x(2) - x(1)
      c(2) = (y(2) - y(1))/d(1)
      do 10 i = 2, nm1
         d(i) = x(i+1) - x(i)
         b(i) = 2.*(d(i-1) + d(i))
         c(i+1) = (y(i+1) - y(i))/d(i)
         c(i) = c(i+1) - c(i)
   10 continue
c
c  end conditions.  third derivatives at  x(1)  and  x(n)
c  obtained from divided differences
c
      b(1) = -d(1)
      b(n) = -d(n-1)
      c(1) = 0.
      c(n) = 0.
      if ( n .eq. 3 ) go to 15
      c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1))
      c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3))
      c(1) = c(1)*d(1)**2/(x(4)-x(1))
      c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3))
c
c  forward elimination
c
   15 do 20 i = 2, n
         t = d(i-1)/b(i-1)
         b(i) = b(i) - t*d(i-1)
         c(i) = c(i) - t*c(i-1)
   20 continue
c
c  back substitution
c
      c(n) = c(n)/b(n)
      do 30 ib = 1, nm1
         i = n-ib
         c(i) = (c(i) - d(i)*c(i+1))/b(i)
   30 continue
c
c  c(i) is now the sigma(i) of the text
c
c  compute polynomial coefficients
c
      b(n) = (y(n) - y(nm1))/d(nm1) + d(nm1)*(c(nm1) + 2.*c(n))
      do 40 i = 1, nm1
         b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.*c(i))
         d(i) = (c(i+1) - c(i))/d(i)
         c(i) = 3.*c(i)
   40 continue
      c(n) = 3.*c(n)
      d(n) = d(n-1)
      return
c
   50 b(1) = (y(2)-y(1))/(x(2)-x(1))
      c(1) = 0.
      d(1) = 0.
      b(2) = b(1)
      c(2) = 0.
      d(2) = 0.
      return
      end

C =======================================================================
C =======================================================================

      double precision function seval(n, u, x, y, b, c, d)
      integer n
      double precision  u, x(n), y(n), b(n), c(n), d(n)
c
c  this subroutine evaluates the cubic spline function
c
c    seval = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))**2 + d(i)*(u-x(i))**3
c
c    where  x(i) .lt. u .lt. x(i+1), using horner's rule
c
c  if  u .lt. x(1) then  i = 1  is used.
c  if  u .ge. x(n) then  i = n  is used.
c
c  input..
c
c    n = the number of data points
c    u = the abscissa at which the spline is to be evaluated
c    x,y = the arrays of data abscissas and ordinates
c    b,c,d = arrays of spline coefficients computed by spline
c
c  if  u  is not in the same interval as the previous call, then a
c  binary search is performed to determine the proper interval.
c
      integer i, j, k
      double precision dx
      data i/1/
      if ( i .ge. n ) i = 1
      if ( u .lt. x(i) ) go to 10
      if ( u .le. x(i+1) ) go to 30
c
c  binary search
c
   10 i = 1
      j = n+1
   20 k = (i+j)/2
      if ( u .lt. x(k) ) j = k
      if ( u .ge. x(k) ) i = k
      if ( j .gt. i+1 ) go to 20
c
c  evaluate spline
c
   30 dx = u - x(i)
      seval = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)))
      return
      end

C =======================================================================
C =======================================================================

      double precision function dseval(n, u, x, y, b, c, d)
      integer n
      double precision  u, x(n), y(n), b(n), c(n), d(n)
c
c  this subroutine evaluates the DERIVATVE OF THE cubic spline function
c
c    seval = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))**2 + d(i)*(u-x(i))**3
c
c    where  x(i) .lt. u .lt. x(i+1), using horner's rule
c
c  if  u .lt. x(1) then  i = 1  is used.
c  if  u .ge. x(n) then  i = n  is used.
c
c  input..
c
c    n = the number of data points
c    u = the abscissa at which the spline is to be evaluated
c    x,y = the arrays of data abscissas and ordinates
c    b,c,d = arrays of spline coefficients computed by spline
c
c  if  u  is not in the same interval as the previous call, then a
c  binary search is performed to determine the proper interval.
c
      integer i, j, k
      double precision dx
      data i/1/
      if ( i .ge. n ) i = 1
      if ( u .lt. x(i) ) go to 10
      if ( u .le. x(i+1) ) go to 30
c
c  binary search
c
   10 i = 1
      j = n+1
   20 k = (i+j)/2
      if ( u .lt. x(k) ) j = k
      if ( u .ge. x(k) ) i = k
      if ( j .gt. i+1 ) go to 20
c
c  evaluate spline first-order derivative

   30 dx = u - x(i)
      dseval = b(i) + dx*(2.*c(i) + dx*3.*d(i))
      return
      end

C =======================================================================
C =======================================================================

      double precision function ddseval(n, u, x, y, b, c, d)
      integer n
      double precision  u, x(n), y(n), b(n), c(n), d(n)
c
c  this subroutine evaluates the SECOND DERIVATVE OF THE cubic spline function
c
c    seval = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))**2 + d(i)*(u-x(i))**3
c
c    where  x(i) .lt. u .lt. x(i+1), using horner's rule
c
c  if  u .lt. x(1) then  i = 1  is used.
c  if  u .ge. x(n) then  i = n  is used.
c
c  input..
c
c    n = the number of data points
c    u = the abscissa at which the spline is to be evaluated
c    x,y = the arrays of data abscissas and ordinates
c    b,c,d = arrays of spline coefficients computed by spline
c
c  if  u  is not in the same interval as the previous call, then a
c  binary search is performed to determine the proper interval.
c
      integer i, j, k
      double precision dx
      data i/1/
      if ( i .ge. n ) i = 1
      if ( u .lt. x(i) ) go to 10
      if ( u .le. x(i+1) ) go to 30
c
c  binary search
c
   10 i = 1
      j = n+1
   20 k = (i+j)/2
      if ( u .lt. x(k) ) j = k
      if ( u .ge. x(k) ) i = k
      if ( j .gt. i+1 ) go to 20
c
c  evaluate spline second-order derivative
c
   30 dx = u - x(i)
      ddseval = 2.*c(i) + dx*6.*d(i)
      return
      end
