/********************************************************
*       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.25500000e+000                  *
*  Estimated Y value = 3.29402327e-001                  *
*  Estimated Error   = -1.32471020e-010                 *
*  Exact Y value     = 3.29402327e-001                  *
*                                                       *
* ----------------------------------------------------- *
* 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].                               *
*                                                       *
*                    C++ Release By J-P Moreau, Paris.  *
*                          (www.jpmoreau.fr)            *
********************************************************/
#include <stdio.h>
#include <math.h>

#define  NMAX  25
#define  TINY  1e-15

double X[NMAX], Y[NMAX];
double XX, YY, XERR;
int    N;


double FCT(double X) {
  return (sin(X) - 2.0*cos(X));
}


void RATINT(double *XA, double *YA, int N, double X, double *Y, double *DY) {
/****************************************************
*    Interpolation or Extrapolation of a Discrete   *
*       Function By a Quotient of Polynomials       *
* ------------------------------------------------- *
* INPUTS:                                           *
*   XA:    Table of abscissas (N)                   *
*   YA:    Table of ordinates (N)                   *
*    N:    Number of points                         *
*    X:    Interpolation abscissa value             *
* OUTPUTS:                                          *
*    Y:    Returned estimation of function for X    *
*   DY:    Estimated error for Y                    *
****************************************************/
  double C[NMAX], D[NMAX];
  double DD,H,HH,T,W;
  int I,M,NS;

  NS=1;
  HH = fabs(X-XA[1]);
  for (I=1; I<=N; I++) {
    H = fabs(X-XA[I]);
    if (H == 0.0) {
      *Y=YA[I];
      *DY=0.0;
      return;
    }      
	else if (H < HH) {
      NS=I;               //index of closest table entry
      HH = H;
    }
    C[I]=YA[I];
    D[I]=YA[I] + TINY;    //TINY is to prevent a zero-over-zero  
	}                     //condition.
  *Y=YA[NS];
  NS--;
  for (M=1; M<N; M++) {
	for (I=1; I<=N-M; I++) {
      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 == 0.0) {
        printf("\n Error: Pole at requested value of X.\n\n");
        return;
      }
      DD=W/DD;
      D[I]=C[I+1]*DD;
      C[I]=T*DD;
    }
    if (2*NS < N-M)
      *DY=C[NS+1];   
    else {   
      *DY=D[NS];
      NS--;    
    }
    *Y = *Y + *DY;  
  } 
}


void main()  {

  N = 12;     //Number of points

  //define tables X and Y
  X[1] = 0.0; Y[1]=FCT(X[1]);
  X[2] = 0.1; Y[2]=FCT(X[2]);
  X[3] = 0.2; Y[3]=FCT(X[3]);
  X[4] = 0.3; Y[4]=FCT(X[4]);
  X[5] = 0.4; Y[5]=FCT(X[5]);
  X[6] = 0.5; Y[6]=FCT(X[6]);
  X[7] = 0.6; Y[7]=FCT(X[7]);
  X[8] = 0.7; Y[8]=FCT(X[8]);
  X[9] = 0.8; Y[9]=FCT(X[9]);
  X[10] = 0.9; Y[10]=FCT(X[10]);
  X[11] = 1.0; Y[11]=FCT(X[11]);
  X[12] = 1.1; Y[12]=FCT(X[12]);

  //define interpolation abscissa
  XX = 1.255;

  //call interpolation function
  RATINT(X,Y,N,XX,&YY,&XERR);
  
  //print results
  printf("\n");
  printf(" for X             = %15.8e\n", XX);
  printf(" Estimated Y value = %15.8e\n", YY);
  printf(" Estimated Error   = %15.8e\n", XERR);
  printf(" Exact Y value     = %15.8e\n", FCT(XX));
  printf("\n");
}

// end of file tratint.cpp