!****************************************************
!*   Program to demonstrate Newton interpolation    *
!*     of Function SIN(X) in real*8 precision       *
!* ------------------------------------------------ *
!* Reference: BASIC Scientific Subroutines, Vol. II *
!* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
!*                                                  *
!*              F90 Version by J.-P. Moreau, Paris. *
!*                       (www.jpmoreau.fr)          *
!* ------------------------------------------------ *
!* SAMPLE RUN:                                      *
!*                                                  *
!*  Newton interpolation of SIN(X):                 *
!*                                                  *
!*    Input X: 0.5                                  *
!*                                                  *
!*    SIN(X) = 0.47961465 (exact value: 0.47942554) *
!*    Error estimate: -0.00000028                   *
!*                                                  *
!****************************************************
PROGRAM Newton

!Labels: 50,60,100

parameter(SIZE=14)


integer iv, n
real*8  X(0:SIZE),Y(0:SIZE)
real*8  e,xx,yy

  iv=14
  print *,' '
  print *,'Newton interpolation of SIN(X):'
  ! Input sine table
  !-----------------------------------------------------------------
  ! Sine table values from  Handbook of mathematical functions
  ! by M. Abramowitz and I.A. Stegun, NBS, june 1964
  !-----------------------------------------------------------------
  X(1)=0.000; Y(1)=0.00000000; X(2)=0.125; Y(2)=0.12467473;
  X(3)=0.217; Y(3)=0.21530095; X(4)=0.299; Y(4)=0.29456472;
  X(5)=0.376; Y(5)=0.36720285; X(6)=0.450; Y(6)=0.43496553;
  X(7)=0.520; Y(7)=0.49688014; X(8)=0.589; Y(8)=0.55552980;
  X(9)=0.656; Y(9)=0.60995199; X(10)=0.721; Y(10)=0.66013615;
  X(11)=0.7853981634; Y(11)=0.7071067812;
  X(12)=0.849; Y(12)=0.75062005; X(13)=0.911; Y(13)=0.79011709;
  X(14)=0.972; Y(14)=0.82601466;
  !-----------------------------------------------------------------

  !Input interpolation point
100 print *,' '

  write(*,"('   Input X: ')",advance ='no')
  read *, xx

  call Interpol_Newton(iv,n,e,xx,yy,X,Y)

  if (n.ne.0) then
    write(*,50) yy, dsin(xx)
    write(*,60) e
    goto 100
  end if
  print *,' '

50 format(/'   SIN(X) = ',F10.8,' (exact value: ',F10.8)
60 format('   Error estimate: ',F10.8)
  stop
end


!********************************************************
!* Newton divided differences interpolation subroutine  *
!* ---------------------------------------------------- *
!* Calculates cubic interpolation for a given table.    *
!* iv is the total number of table values. X(i), Y(i)   *
!* are the coordinate table values, Y(i) being the      *
!* dependant variable. The X(i) may be arbitrarily spa- *
!* ced. xx is the interpolation point which is assumed  *
!* to be in the interval  with at least one table value *
!* to the left, and 3 to the right. If this is violated,*
!* n will be set to zero. It is assumed that the table  *
!* values are in ascending X(i) order. e is the error   *
!* estimate. The returned value is yy.                  *
!********************************************************
Subroutine Interpol_Newton(iv,n,e,xx,yy,X,Y)  
  !Labels: 100,200,300
  parameter(SIZE=14)
  integer i
  real*8  a,b,c
  real*8  X(0:SIZE),Y(0:SIZE),Y1(0:SIZE),Y2(0:SIZE),Y3(0:SIZE)
  real*8  e,xx,yy
    ! Check to see if interpolation point is correct
    n=1
    if (xx>=X(1)) goto 100
    n=0 ; return
100 if (xx<=X(iv-3)) goto 200
    n=0 ; return
    ! Generate divided differences
200 do i = 1, iv-1 
      Y1(i)=(Y(i+1)-Y(i))/(X(i+1)-X(i))
    end do
    do i = 1, iv-2
      Y2(i)=(Y1(i+1)-Y1(i))/(X(i+1)-X(i));
    end do
    do i = 1, iv-3
      Y3(i)=(Y2(i+1)-Y2(i))/(X(i+1)-X(i));
    end do
    ! Find relevant table interval
    i=0
300 i = i + 1
    if (xx>X(i)) goto 300
    i = i - 1
    ! Begin interpolation
    a=xx-X(i)
    b=a*(xx-X(i+1))
    c=b*(xx-X(i+2))
    yy=Y(i)+a*Y1(i)+b*Y2(i)+c*Y3(i)
    ! Calculate next term in the expansion for an error estimate
    e=c*(xx-X(i+3))*yy/24.d0
  
    return
end

! End of file Newton.f90