program main
!====================================================================
! General polynomial interpolation
!
! Function values f(x) are calculated at n base points. A local
! polynomial of selected order is then used to interpolate the
! function at nint points across the interval. The program also
! calculates the average absolute interpolation error.
!
! The interpolating polynomial is evaluated in Lagrange form.
! For order = 1, 2, 3, ... the routine uses 2, 3, 4, ...
! neighboring data points, respectively.
!
! Original implementation 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, yi
double precision x, y, step, ys, error, errav
integer i, order
double precision f, polynomial_interp

xmin = 0.0
xmax = 2.0

! Step 1: generate xi and yi from f(x), xmin, xmax, and n
step = (xmax-xmin)/dble(n-1)

do i=1,n
  xi(i) = xmin + step*dble(i-1)
  yi(i) = f(xi(i))
end do

! Step 2: choose the interpolation order
! order = 1 -> linear
! order = 2 -> quadratic
! order = 3 -> cubic
order = 3

write(*,*) '    General Polynomial Interpolation'
write(*,203) order

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

write(*,201)

do i=1,nint
  x = xmin + step*dble(i-1)
  y = f(x)

  ys = polynomial_interp(x, xi, yi, n, order+1)

  error = ys-y
  write(*,200) x, ys, error

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

write(*,202) errav

200 format (3f12.5)
201 format ('        x     interpolated    error')
202 format ('           Average error',f12.5)
203 format ('     order of interpolation = ',i2)

end program main


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

f = sin(x)

end function f


function polynomial_interp(xx, xi, yi, ni, npts)
!====================================================================
! Local polynomial interpolation of selected order.
!
! The interpolating polynomial is evaluated in Lagrange form using
! npts neighboring data points. The interpolation order is npts-1.
!
! Input:
!   xx   = abscissa at which the interpolation is evaluated
!   xi   = array of data abscissas
!   yi   = array of data ordinates
!   ni   = number of data points in xi and yi
!   npts = number of points used for interpolation
!          (interpolation order = npts-1)
!
! Output:
!   polynomial_interp = interpolated value at xx
!
! Comments:
!   The routine works for both equally and unequally spaced xi.
!   A binary search locates the interval containing xx, and a local
!   group of neighboring data points is selected for interpolation.
!   If xx lies outside the data interval, the nearest boundary value
!   is returned.
!
! Original implementation by Alexander Godunov, January 2010.
! Revised for the companion website, 2026.
!====================================================================
implicit none

double precision polynomial_interp, xx
integer ni, npts
double precision xi(ni), yi(ni)
double precision lambda(ni)

integer i, j, k, js, jl, nuse
double precision y

! Check the number of interpolation points
nuse = npts
if (nuse > ni) nuse = ni

! If xx is outside the xi(1)-xi(ni) interval, return a boundary value
if (xx <= xi(1)) then
  polynomial_interp = yi(1)
  return
end if

if (xx >= xi(ni)) then
  polynomial_interp = yi(ni)
  return
end if

! Binary search to find i such that xi(i) < xx < xi(i+1)
i = 1
j = ni

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

! Shift i so that xx lies near the middle of the selected data points
i = i + 1 - nuse/2

! Keep the interpolation points inside the available data range
if (i < 1) i = 1
if (i + nuse - 1 > ni) i = ni-nuse+1

! Evaluate the interpolating polynomial in Lagrange form
y = 0.0

do js=i,i+nuse-1
  lambda(js) = 1.0

  do jl=i,i+nuse-1
    if (jl /= js) then
      lambda(js) = lambda(js)*(xx-xi(jl))/(xi(js)-xi(jl))
    end if
  end do

  y = y + yi(js)*lambda(js)
end do

polynomial_interp = y

end function polynomial_interp