program main
!====================================================================
! Largange interpolation (lagint.f90)
! Comments: values of function f(x) are calculated in n base points
! then: Lagrange interpolation is computed in 2n-1 points, as well as
! a difference sum|f(u)-lagint(u)| 
!====================================================================
implicit none
integer, parameter :: n=11      ! base points for interpolation
integer, parameter :: nint=21   ! compute interpolation in nint points
double precision xmin, xmax     ! given interval of x()
double precision, dimension (n) :: xi(n), yi(n)
double precision x, y, step, ys, error, errav
integer i, order
double precision f, lagint

! open files
!open (unit=1, file='tablex1.dat')
!open (unit=2, file='tablex2.dat')

xmin = 0.0
xmax = 2.0

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

!  step 2: interpolation at nint points
!  where order - is the order if interpolation
order = 3
write(*,*) '    Lagrange Interpolation'
write(*,203) order
errav = 0.0
step = (xmax-xmin)/(nint-1)
write(*,201) 
do i=1, nint
  x = xmin + step*float(i-1) 
  y = f(x) 

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

  error = ys-y
  write (*,200) x, ys, error
! step 4: calculate quality of interpolation               
  errav = errav + abs(y-ys)/nint 
end do
write (*,202) errav
200 format (3f12.5)
201 format ('        x        Lagrange     error')    
202 format ('           Average error',f12.5)
203 format ('     order of interpolation = ',i2)

end program main

!
!  Function f(x)
!
  function f(x)
  double precision f, x
  f = sin(x) 
  end function f

  function lagint(xx, xi, yi, ni, n)
!====================================================================
! Lagrange interpolation (straightforward implementation)
! Alex G: January 2010
!--------------------------------------------------------------------
! input ...
! xx    - the abscissa at which the interpolation is to be evaluated
! xi()  - the arrays of data abscissas
! yi()  - the arrays of data ordinates
! ni    - size of the arrays xi() and yi()
! n     - number of points for interpolation (order of interp. = n-1)
! output ...
! lagint- interpolated value
! comments ...
! if (n > ni) n = ni
! program works for both equally and unequally spaced xi()
!=====================================================================
implicit none
double precision lagint, xx
integer ni, n
double precision xi(ni), yi(ni)
double precision lambda(ni)
integer i, j, k, js, jl
double precision y

! check order of interpolation
if (n > ni) n = ni

! if x is ouside the xi(1)-xi(ni) interval take a boundary value
if (xx <= xi(1)) then
  lagint = yi(1)
  return
end if
if (xx >= xi(ni)) then
  lagint = yi(ni)
  return
end if

! a binary (bisectional) search to find i so that xi(i) < x < 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 that will correspond to n-th order of interpolation
! the search point will be in the middle in x_i, x_i+1, x_i+2 ...
i = i + 1 - n/2

! check boundaries: if i is ouside of the range [1, ... n] -> shift i
if (i < 1) i=1
if (i + n > ni) i=ni-n+1

! Lagrange interpolation
y = 0.0
do js=i,i+n-1
  lambda(js)=1.0
  do jl=i,i+n-1
    if(jl /= js) lambda(js)=lambda(js)*(xx-xi(jl))/(xi(js)-xi(jl))
  end do
  y = y + yi(js)*lambda(js)
end do
lagint = y
end function lagint