program main
!====================================================================
! Cubic spline interpolation
!
! Function values f(x) are calculated at n base points. The spline
! coefficients are then computed, and the spline is evaluated at
! nint points across the interval. The program also calculates the
! average absolute interpolation error.
!
! Origin:
!   Based on the SPLINE and SEVAL routines presented by
!   G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
!   Computer Methods for Mathematical Computations,
!   Prentice-Hall, 1977.
!
! This version was rewritten and adapted for modern Fortran
! by Alexander Godunov, January 2010.
! Revised for the companion website, 2026.
!====================================================================
implicit none
integer, parameter :: n=11      ! base points for interpolation
integer, parameter :: nint=21   ! compute interpolation in nint points
double precision xmin, xmax     ! interpolation interval
double precision, dimension (n) :: xi(n), yi(n), b(n), c(n), d(n)
double precision x, y, step, ys, error, errav
integer i
double precision f, spline_eval

xmin = 0.0
xmax = 2.0

! Step 1: generate xi and yi from f(x), xmin, xmax, and n
step = (xmax-xmin)/(n-1)
do i=1,n
  xi(i) = xmin + step*float(i-1)
  yi(i) = f(xi(i))
end do

! Step 2: calculate spline coefficients
call spline (xi, yi, b, c, d, n)

! Step 3: evaluate the spline at nint points
errav = 0.0
step = (xmax-xmin)/(nint-1)

write(*,201)
do i=1,nint
  x = xmin + step*float(i-1)
  y = f(x)
  ys = spline_eval(x, xi, yi, b, c, d, n)
  error = ys-y
  write (*,200) x, ys, error

! Step 4: calculate the average absolute interpolation error
  errav = errav + abs(y-ys)/nint
end do

write (*,202) errav

200 format (3f12.5)
201 format ('        x        spline      error')
202 format ('           Average error',f12.5)

end program main


!
! Function f(x)
!
function f(x)
implicit none
double precision f, x

f = sin(x)

end function f


subroutine spline (x, y, b, c, d, n)
!======================================================================
! Calculate the coefficients b(i), c(i), and d(i), i=1,2,...,n,
! for cubic spline interpolation
!
!   s(x) = y(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))**2
!                + d(i)*(x-x(i))**3
!
! for x(i) <= x <= x(i+1).
!
! Input:
!   x = array of data abscissas in strictly increasing order
!   y = array of data ordinates
!   n = number of data points (n >= 2)
!
! Output:
!   b, c, d = arrays of spline coefficients
!
! Origin:
!   Based on the SPLINE routine presented by
!   G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
!   Computer Methods for Mathematical Computations,
!   Prentice-Hall, 1977.
!
! This version was rewritten and adapted for modern Fortran
! by Alexander Godunov, January 2010.
!======================================================================
implicit none
integer n
double precision x(n), y(n), b(n), c(n), d(n)
integer i, j, gap
double precision h

gap = n-1

! Check input
if (n < 2) return

if (n < 3) then
  b(1) = (y(2)-y(1))/(x(2)-x(1))   ! linear interpolation
  c(1) = 0.0
  d(1) = 0.0
  b(2) = b(1)
  c(2) = 0.0
  d(2) = 0.0
  return
end if

!
! Step 1: preparation
!
d(1) = x(2) - x(1)
c(2) = (y(2) - y(1))/d(1)

do i = 2, gap
  d(i) = x(i+1) - x(i)
  b(i) = 2.0*(d(i-1) + d(i))
  c(i+1) = (y(i+1) - y(i))/d(i)
  c(i) = c(i+1) - c(i)
end do

!
! Step 2: end conditions
!
b(1) = -d(1)
b(n) = -d(n-1)
c(1) = 0.0
c(n) = 0.0

if (n /= 3) then
  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))
end if

!
! Step 3: forward elimination
!
do i = 2, n
  h = d(i-1)/b(i-1)
  b(i) = b(i) - h*d(i-1)
  c(i) = c(i) - h*c(i-1)
end do

!
! Step 4: back substitution
!
c(n) = c(n)/b(n)

do j = 1, gap
  i = n-j
  c(i) = (c(i) - d(i)*c(i+1))/b(i)
end do

!
! Step 5: compute spline coefficients
!
b(n) = (y(n) - y(gap))/d(gap) + d(gap)*(c(gap) + 2.0*c(n))

do i = 1, gap
  b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.0*c(i))
  d(i) = (c(i+1) - c(i))/d(i)
  c(i) = 3.0*c(i)
end do

c(n) = 3.0*c(n)
d(n) = d(n-1)

end subroutine spline


function spline_eval(u, x, y, b, c, d, n)
!======================================================================
! Evaluate the cubic spline at point u:
!
!   spline_eval = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))**2
!                        + d(i)*(u-x(i))**3
!
! where x(i) <= u <= x(i+1).
!
! Input:
!   u       = abscissa at which the spline is evaluated
!   x, y    = arrays of given data points
!   b, c, d = spline coefficients computed by spline
!   n       = number of data points
!
! Output:
!   spline_eval = interpolated value at u
!
! Origin:
!   Based on the SEVAL routine presented by
!   G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
!   Computer Methods for Mathematical Computations,
!   Prentice-Hall, 1977.
!
! This version was rewritten and adapted for modern Fortran
! by Alexander Godunov, January 2010.
!======================================================================
implicit none
double precision spline_eval
integer n
double precision u, x(n), y(n), b(n), c(n), d(n)
integer i, j, k
double precision dx

! If u is outside the x() interval, return the nearest boundary value.
if (u <= x(1)) then
  spline_eval = y(1)
  return
end if

if (u >= x(n)) then
  spline_eval = y(n)
  return
end if

!
! Binary search for i such that x(i) <= u <= x(i+1)
!
i = 1
j = n+1

do while (j > i+1)
  k = (i+j)/2
  if (u < x(k)) then
    j = k
  else
    i = k
  end if
end do

!
! Evaluate spline interpolation
!
dx = u - x(i)
spline_eval = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)))

end function spline_eval