!*********************************************************
!*       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.329402327209023               *
!*  Estimated Error   = -1.327032626597764E-010          *
!*  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_RATINT
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 RATINT(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 RATINT(XA,YA,N,X,Y,DY)
!*****************************************************
!*    Interpolation or Extrapolation of a Discrete   *
!*       Function By a Quotient of Polynomials       *
!* ------------------------------------------------- *
!* 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,TINY=1.E-25)
REAL*8 XA(N),YA(N), X,Y,DY
REAL*8 C(NMAX),D(NMAX)
REAL*8 DD,H,HH,T,W
NS=1
HH=DABS(X-XA(1))
DO I=1,N
  H=DABS(X-XA(I))
  IF(H.EQ.0.D0) THEN
    Y=YA(I)
    DY=0.D0                 
    RETURN
  ELSE IF (H.LT.HH) THEN
    NS=I
    HH=H
  ENDIF
  C(I)=YA(I)
  D(I)=YA(I)+TINY        !The TINY part is needed to prevent a rare 
END DO                   !zero-over-zero condition.
Y=YA(NS)
NS=NS-1
DO M=1,N-1
  DO I=1,N-M
    W=C(I+1)-D(I)
    H=XA(I+M)-X          !H<>0 (tested before)
    T=(XA(I)-X)*D(I)/H
    DD=T-C(I+1)
    IF(DD.EQ.0.) Pause ' Pole at requested value of X.'
    DD=W/DD
    D(I)=C(I+1)*DD
    C(I)=T*DD
  END DO
  IF(2*NS.LT.N-M) THEN
    DY=C(NS+1)
  ELSE 
    DY=D(NS) 
    NS=NS-1
  ENDIF
  Y=Y+DY
END DO 

RETURN
END

!end of file tratint.f90