#include "nr3.h"
#include <cmath>
#include <iostream>
#include "interp_1d.h"
#include "interp_linear.h"
//#include "TROOT.h"
#include "TFile.h"
#include "TTree.h"

// interpolate a function. test spline, linear, 3rd order and 10th order polynomial interpolation schemes.
// as an addition, write out the results as a root tree for plotting purposes
Doub f(Doub x) {
// the function x^2/(x-4.75)*exp(-1/(x-4.75)^2)
// catch underflows close to x=4.75
// if (x-4.75) is smaller than 0.1, the function is smaller than 10x^2e^-100 which is smaller than 10^-30. this is completely negligible. 
   if (fabs(x-4.75)<0.1) return 0;
   return x*x/(x-4.75)*exp(-1/(x-4.75)/(x-4.75));
}

int main() {
   VecDoub xx(31),yy(31); // 31 values
   for (int i=0; i<31; i++) {
      xx[i]=0.3*i;
      yy[i]=f(xx[i]);
   }
   Poly_interp pol1(xx,yy,2); // linear interpolation, use 2 points
   Poly_interp pol3(xx,yy,4); // third-order polynomial interpolation: for linear use 2, second-order 3, third-order 4
   Poly_interp pol10(xx,yy,11); // 10th-order polynomial interpolation
   Spline_interp spline(xx,yy,1e101,1e101); // create a natural cubic spline
   //
   // optional: for root plotting, create a root file containing a tree with the results
   TFile *rootfil = new TFile("interpol.root","RECREATE"); // creates a root file on disk
// create the structure that we want to fill in the root tree
   struct rootstruct {
      double x;
      double y;
      double lin;
      double spline;
      double pol3;
      double pol10;
   } polresults;
   TTree *tree = new TTree("tree",""); // creates a root tree with root name tree and title "" (no title) 
   tree->Branch("pol",&polresults.x,"x/D:y:ylin:yspline:ypol3:ypol10"); // create a branch from the structure. it will be filled with the entries in polresults
  
   // print output to stdout for the exercise
   cout << "results : " << endl;
   cout << "x\t\ty\t\tlinear-diff\tspline-diff\tpol3-diff\tpol10-diff"<<endl;
   cout.precision(9); // print 9 digits and a tab
   for (int i=0; i<181; i++) {
      double x=i*0.05;
      double y=f(x);
      double ylin=pol1.interp(x);
      double yspline=spline.interp(x);
      double ypol3=pol3.interp(x);
      double ypol10=pol10.interp(x);
      cout <<x<<"\t"<<y<<"\t"<<ylin-y<<"\t"<<yspline-y<<"\t"<<ypol3-y<<"\t"<<ypol10-y<<endl;
     
      // optional: Fill the root tree
      polresults.x=x;
      polresults.y=y;
      polresults.pol3=ypol3;
      polresults.pol10=ypol10;
      polresults.spline=yspline;
      polresults.lin=ylin;
      tree->Fill();
   }
   rootfil->Write();
   rootfil->Close();
   return 0;
}