program main
!====================================================================
! Interpolation by divided difference polinomials (ddint.f90)
! Comments: values of function f(x) are calculated in n base points
! then: the DD interpolation is computed in 2n-1 points, as well as
! a difference sum|f(u)-ddint(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, ddint

! 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(*,*) ' Interpolation by Divided Differences'
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 = ddint(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        DDint       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 ddint(xx, xi, yi, ni, n)
!====================================================================
! Interpolation based on Divided Difference Polynomials:
! 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 ...
! ddint   - interpolated value
! comments ...
! if (n > ni) n = ni
! program works for both equally and unequally spaced xi()
!=====================================================================

implicit none
double precision ddint, xx
integer ni, n
double precision xi(ni), yi(ni)
double precision d(n,n), x(n)
integer i, j, k, ix
double precision c, pn

! 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
  ddint = yi(1)
  return
end if
if (xx >= xi(ni)) then
  ddint = 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

!  old output to test i
!  write(*,100) xx, i
!  100 format (f10.5, I5)

! just wanted to use index i later for d coefficients
ix = i

! initialization of d(n,1) and x(n)
do i=1,n
  d(i,1) = yi(ix+i-1)
  x(i)   = xi(ix+i-1)
end do

! calculations for the divided difference coefficients
do j=2,n
  do i=1,n-j+1
    d(i,j)=(d(i+1,j-1)-d(i,j-1))/(x(i+1+j-2)-x(i))
  end do
end do

! print results for the d coeff.
!  do i=1,n
!    write(*,200) (d(i,j),j=1,n-i+1)
!  end do
!200 format (4f10.6)

! divided difference interpolation
Pn = d(1,1)
do i=1,n-1
  c = 1.0
  do j=1,i
    c = c*(xx - x(j))
  end do
  Pn = Pn + c*d(1,i+1);
end do

ddint = Pn
end function ddint