#include"TH1D.h"
#include "TF1.h"
#include "TMath.h"
double X[100];
double Y[100];
//int n=19; //i nodi sono n+1 tra a e b [provare con 19 o 25 tra a e b]
int n=10; //i nodi sono n+1 tra 0 e 9 
double fteorica(double x)
{ 
  double y = x*x*x;
  y = y * cos(x);
  return y;
}
double Lagrange(double x, int i)
{
	double y=1.;
	for(int j=0; j<n+1; j++)
	{
		if(j==i) continue;
		y=y*(x-X[j])/(X[i]-X[j]);
	}
 	return y;
}
double ginterpol(double *x, double *par)
{
	double y=0;
	for(int j=0;j<n+1; j++)
	{
		y=y+Y[j]*Lagrange(x[0],j);
	}
 	return y;
}
double resto(double *x, double *par)
{
  return fteorica(x[0]) - ginterpol(x,par);
}

void interpolPolLagrange()
{

  // estremi del doinio della funzione teorica 
	double a=-10.;
	double b=10.;
	TF1 *ptfteorica= new TF1("tfteorica","fteorica(x)",a,b);
	TCanvas *c1=new TCanvas("c1","c1", 800,600);
	ptfteorica->SetLineColor(kRed);
	ptfteorica->SetLineWidth(5);
	ptfteorica->Draw();
	//return;

	//inizializzo i nodi (n+1 componenti degli array X e Y, con indici da 0 a n)
	X[0]=a;
	X[n]=b;
	double passo=(X[n]-X[0])/n;
	cout<< "al passo i=0, X="<<X[0]<<endl;
	Y[0]=fteorica(X[0])	;
	cout<< "al passo i=0, Y="<<Y[0]<<endl;  
	for(int i=1; i<n+1;i++)
	{
		X[i]=X[i-1]+passo;
		//cout<< "al passo i="<<i<<" ,Y="<<Y[i]<<endl; 
		Y[i]=fteorica(X[i])	;
		cout<< "al passo i="<<i<<" ,X="<<X[i]<<" , Y="<<Y[i]<<std::endl; 
	}
	//return;
	
	// uso il costruttore di TF1 a https://root.cern.ch/doc/master/classTF1.html#F4
	// si noti che questo costruttore richiede la segnatura della funzione secondo questa sintassi:
	// double ginterpol(double *x, double *par)
	TF1 *pfinterpol= new TF1("finterpol",ginterpol,X[0],X[n]);
	pfinterpol->SetLineColor(kBlue);
	pfinterpol->SetLineStyle(2);
	pfinterpol->SetLineWidth(5);
	pfinterpol->Draw("same");
	
	//costruisco un grafico a punti dei nodi 
	TGraph * p=new TGraph(n+1,X,Y);
	p->SetMarkerStyle(21);
	p->SetMarkerSize(1.2);
	p->SetMarkerColor(1);
	p->Draw("Psame");

	//costruisco un grafico della funzione resto  
	TF1 *pfresto= new TF1("resto",resto,X[0], X[n]);
	TCanvas *c2=new TCanvas("c2","c2", 800,600);
	pfresto->Draw();
	
	
	
		
}