!*********************************************************
!*       Polynomial Interpolation or Extrapolation       *
!*              of a Discrete Function F(x)              *
!* ----------------------------------------------------- *
!* SAMPLE RUN:                                           *
!* (Example: Function sin(x) - 2*cos(x) is given by 12   *
!*          points from x=0 to x=1.1.                    *
!*          Extrapolate for x=1.255).                    *
!*                                                       *
!*  For X             =   1.25500000000000               *
!*  Estimated Y value =  0.329402327224849               *
!*  Estimated Error   = -8.258969146682216E-011          *
!*  Exact Y value     =  0.329402327220005               *
!*                                                       *
!* ----------------------------------------------------- *
!* REFERENCE: "Numerical Recipes, The Art of Scientific  *
!*             Computing By W.H. Press, B.P. Flannery,   *
!*             S.A. Teukolsky and W.T. Vetterling,       *
!*             Cambridge University Press, 1986"         *
!*             [BIBLI 08].                               *
!*                                                       *
!*                    F90 Release By J-P Moreau, Paris.  *
!*                           (www.jpmoreau.fr)           *
!*********************************************************   
PROGRAM TEST_POLINT
PARAMETER(NMAX=25)
REAL*8 X(NMAX), Y(NMAX), XX, YY, XERR, FCT

  N = 12      !Number of points

  !define tables X and Y
  X(1) = 0.d0;  Y(1)=FCT(X(1))
  X(2) = 0.1d0; Y(2)=FCT(X(2))
  X(3) = 0.2d0; Y(3)=FCT(X(3))
  X(4) = 0.3d0; Y(4)=FCT(X(4))
  X(5) = 0.4d0; Y(5)=FCT(X(5))
  X(6) = 0.5d0; Y(6)=FCT(X(6))
  X(7) = 0.6d0; Y(7)=FCT(X(7))
  X(8) = 0.7d0; Y(8)=FCT(X(8))
  X(9) = 0.8d0; Y(9)=FCT(X(9))
  X(10) = 0.9d0; Y(10)=FCT(X(10))
  X(11) = 1.d0;  Y(11)=FCT(X(11))
  X(12) = 1.1d0; Y(12)=FCT(X(12))

  !define interpolation abscissa
  XX = 1.255d0

  !call interpolation subroutine
  CALL POLINT(X,Y,N,XX,YY,XERR)
  
  !print results
  print *,' '
  print *,' For X             =', XX
  print *,' Estimated Y value =', YY
  print *,' Estimated Error   =', XERR
  print *,' Exact Y value     =', FCT(XX)
  print *,' '  

END

FUNCTION FCT(X)
  REAL*8 X, FCT
  FCT=DSIN(X) - 2.D0*DCOS(X)
  RETURN
END

SUBROUTINE POLINT(XA,YA,N,X,Y,DY)
!*****************************************************
!*     Polynomial Interpolation or Extrapolation     *
!*            of a Discrete Function                 *
!* ------------------------------------------------- *
!* INPUTS:                                           *
!*   XA:    Table of abcissas  (N)                   *
!*   YA:    Table of ordinates (N)                   *
!*    N:    Number of points                         *
!*    X:    Interpolating abscissa value             *
!* OUTPUTS:                                          *
!*    Y:    Returned estimation of function for X    *
!*   DY:    Estimated error for Y                    *
!*****************************************************
PARAMETER(NMAX=25)
REAL*8 XA(N),YA(N), X,Y,DY
REAL*8 C(NMAX),D(NMAX)
REAL*8 DEN,DIF,DIFT,HO,HP,W
NS=1
DIF=DABS(X-XA(1))
DO I=1,N
  DIFT=DABS(X-XA(I))
  IF(DIFT.LT.DIF) THEN
    NS=I                 !index of closest table entry
	DIF=DIFT
  ENDIF
  C(I)=YA(I)             !Initialize the C's and D's
  D(I)=YA(I)
END DO
Y=YA(NS)                 !Initial approximation of Y
NS=NS-1
DO M=1,N-1
  DO I=1,N-M
    HO=XA(I)-X
	HP=XA(I+M)-X
    W=C(I+1)-D(I) 
    DEN=HO-HP
	IF(DEN.EQ.0.) Pause 'Error: two identical abcissas)'
	DEN=W/DEN
	D(I)=HP*DEN          !Update the C's and D's
	C(I)=HO*DEN
  END DO
  IF(2*NS.LT.N-M) THEN   !After each column in the tableau XA
    DY=C(NS+1)           !is completed, we decide which correction,
  ELSE                   !C or D, we want to add to our accumulating 
    DY=D(NS)             !value of Y, i.e. which path to take through 
    NS=NS-1              !the tableau, forking up or down. We do this
  ENDIF                  !in such a way as to take the most "straight 
  Y=Y+DY	         !line" route through the tableau to its apex,
END DO                   !updating NS accordingly to keep track of 
                         !where we are. This route keeps the partial
RETURN                   !approximations centered (insofar as possible)
END                      !on the target X.The last DY added is thus the
                         !error indication.

!end of file tpolint.f90