// root.cern.ch/function-interpolation
//
// This macro performs a simple poly-type interpolation throughout a set of datapoints
//
// See also: www.gnu.org/software/gsl/manual/html_node/Interpolation.html

#include <Math/Interpolator.h>

void macro_interpolation ()
{
  float xmin = -3, xmax = 2.5;
 
  int Npts = 10;
  double* xi = new double [Npts];
  double* yi = new double [Npts];
 
  TF1* funPoly = new TF1 ("funpoly", "[0]+[1]*x+[2]*x^2+[3]*x^3", xmin, xmax);
  funPoly->SetParameters (1, -1.5, 1, 1); 
 
  for ( int i = 0; i < Npts; ++i )
  {   
    xi[i] = (double) i * (xmax - xmin) / (Npts-1) + xmin;
    yi[i] = funPoly->Eval ( xi[i] );
  }

  ROOT::Math::Interpolator inter (Npts, ROOT::Math::Interpolation::kPOLYNOMIAL);
/*
  Interpolator (
    unsigned int ndata,                   : how many data points
    ROOT::Math::Interpolation::Type type  : kLINEAR, kPOLYNOMIAL, kCSPLINE, kAKIMA, kCSPLINE_PERIODIC, kAKIMA_PERIODIC
  )

  Setting up the interpolation tool. Choose among the following methods:
   
  o kLINEAR           : straight line between points
  o kPOLYNOMIAL       : to be used for small number of points since introduces large oscillations
  o kCSPLINE          : cubic spline with natural  boundary conditions
  o kCSPLINE_PERIODIC : cubic spline with periodic boundary conditions
  o kAKIMA            : Akima spline with natural  boundary conditions ( requires min. 5 points)
  o kAKIMA_PERIODIC   : Akima spline with periodic boundaries          ( requires min. 5 points)

  See: root.cern.ch/root/html/ROOT__Math__Interpolator.html#ROOT__Math__Interpolator:Interpolator
*/
 
  inter.SetData (Npts, xi, yi);

/*
  SetData is overloaded and allows to plug in:
  - 2x array of doubles or
  - 2x vector<double>
*/

  int Nprob = 100;
  double* Xprob  = new double [Nprob];
  double* Yinter = new double [Nprob];
  
  for ( int i = 0; i < Nprob; ++i ) {
    Xprob [i] = (double) i*(xmax-xmin)/(Nprob-1) + xmin;
    Yinter[i] = inter.Eval ( Xprob[i] );
  }
 
  TCanvas* can1 = new TCanvas ("c1", "Interpolation", 600, 400); 
 
  TGraph* gf = new TGraph (Npts, xi, yi);
  gf->SetLineColor(14);
  gf->SetLineWidth(3);
  gf->SetMarkerStyle(22);
  gf->SetTitle ("Triangles: input.  Curve: probed interpolator response");
  gf->Draw ("AP");
 
  TGraph* gi = new TGraph (Nprob, Xprob, Yinter);
  gi->SetLineColor (28);
  gi->SetLineWidth (1);
  gi->SetTitle ("Integral");
  gi->Draw ("SAME L"); 
}