#include <stdio.h>
#include <math.h>
#include "TGraph.h"


Double_t funz(Double_t x)
{
  /*
  if(x==0)
    return 0;
  else
    return (exp(x)-1)/sqrt(x);
  */

  return 100.-4*x*x +exp(x)*x;
  //return exp(x)*x;
  //return 3. - (x)*(x)+2.*x;
  //return x*x*x*x*x;
  
  // x*exp(x) -> primitiva exp(x)*(x-1) 
  
}
Double_t funzR(Double_t *x, Double_t *par)
{
  return funz(x[0]);
}
float pfunz(float x)
{
  return 100.*x- 4.*x*x*x/3. +exp(x)*(x-1);
  
  //return exp(x)*(x-1);
  //return 3.*x - x*x*x/3. +x*x;
  //return x*funz(x)/6.;
}
  
double cintRett(float a, float b, int N) {
	double h, integ, x;
	h = (b-a)/N;
	//metodo dei rettangoli
	integ=0;
	for(x=a; x<b-h/4.; x+=h)
	{
	integ+=funz(x);
	}
	integ*=h;
	return integ;
}
double cintRettPM(float a, float b, int N) {
	double h, integ, x;
	h = (b-a)/N;
	//metodo dei punti medi
	integ=0;
	for(x=a+h/2; x<b-h/4.; x+=h)
	{
	integ+=funz(x);
	}
	integ*=h;
	return integ;
}
double cintTrap(float a, float b, int N) {
	double h, integ, x;
	h = (b-a)/N;
	//metodo dei trapezi
	integ=(funz(a)+funz(b))/2 ;
	for(x=a+h; x<b-h/4.; x+=h)
	{
	integ+=funz(x);
	}
	integ*=h;
	return integ;
}

double cintSim(float a, float b, int N) {
	double h, integ, x;
	h = (b-a)/N;
	//metodo di Cavalieri-Simpson
	float y[10000];
	int i=0;
	//for(x=a; !(x>b); x+=h)
	for(x=a; x<b+h/4.; x+=h)
	{
		y[i]=funz(x);
		i++;
	}
	integ=0;
	for(int j=1; j<i-1; j+=2)
	{
	  if (j+1 > i-1)
	    {
	      std::cout<<" a problem here ..... "<<std::endl;
	    }
	  integ+=y[j-1]+4*y[j]+y[j+1];
	}
	integ*=h/3;
	return integ;
}


void confrontaInt_r5()
{

  double intRett = 0.;
  double intRettPM = 0.;
  double intTrap = 0.;
  double intSim = 0.;

  double a = 0.;
  double b = 3.;
  double h;

  double integraleEsatto = pfunz(b)-pfunz(a);
  std::cout<<" integraleEsatto "<<integraleEsatto<<std::endl;


  double fRett = 0.;
  double fRettPM = 0.;
  double fRettPM2 = 0.;
  double fTrap = 0.;
  double fSim = 0.;
  int NP = 21;
  double xvec[NP];
  double xvecSim[NP];
  double yvecRett[NP];
  double yvecRettPM[NP];
  double yvecTrap[NP];
  double yvecSim[NP];
  int countP = 0;
  int countPs = 0;
  for (int i=3; i<21; ++i)
    {
      h = (b-a)/i;
      xvec[countP] = h;

      intRett = cintRett(a, b, i);
      yvecRett[countP] =100.*fabs(integraleEsatto-intRett)/integraleEsatto;
      intRettPM = cintRettPM(a, b, i);
      yvecRettPM[countP] =100.*fabs(integraleEsatto-intRettPM)/integraleEsatto;
      intTrap = cintTrap(a, b, i);
      yvecTrap[countP] =100.*fabs(integraleEsatto-intTrap)/integraleEsatto;
      if (i%2 == 0) {
	intSim = cintSim(a, b, i);
	xvecSim[countPs] = h;
	yvecSim[countPs] =100.*fabs(integraleEsatto-intSim)/integraleEsatto;
	countPs++;
      }
      else intSim=0.;
      std::cout<<" n. passi =  "<<i<<" h "<<h<<" intRett,PM,Trap,Sim "<<intRett<<" "<<intRettPM<<" "<<intTrap<<" "<<intSim<<std::endl;
      if (i==4)
	{
	  fRett = 100.*fabs(integraleEsatto-intRett)/(h*integraleEsatto);
	  fRettPM  = 100.*fabs(integraleEsatto-intRettPM)/(h*integraleEsatto);
	  fRettPM2 = 100.*fabs(integraleEsatto-intRettPM)/(h*h*integraleEsatto);
	  fTrap = 100.*fabs(integraleEsatto-intTrap)/(h*h*integraleEsatto);
	  fSim = 100.*fabs(integraleEsatto-intSim)/(h*h*h*h*integraleEsatto);
	}
      countP++;
    }

  TGraph * gIntRett = new TGraph(countP,   xvec, yvecRett);
  TGraph * gIntRettPM = new TGraph(countP, xvec, yvecRettPM);
  TGraph * gIntTrap = new TGraph(countP,   xvec, yvecTrap);
  TGraph * gIntSim = new TGraph(countPs,    xvecSim, yvecSim);

  
  TCanvas* C = new TCanvas("c","c",800,600);
  C->Divide(2,2);
  C->cd(1);
  gIntRett->SetMarkerStyle(20);
  gIntRett->SetMarkerColor(kBlack);
  gIntRett->SetMarkerSize(0.5);
  gIntRett->Draw("AP");
  TF1 *fa = new TF1("fRett","[0]*x",0.,(b-a)/3.);
  fa->SetParameter(0,fRett);
  fa->Draw("SAME");
  C->cd(2);
  gIntRettPM->SetMarkerStyle(20);
  gIntRettPM->SetMarkerColor(kRed);
  gIntRettPM->SetMarkerSize(0.5);
  gIntRettPM->Draw("AP");
  TF1 *fb = new TF1("fRettPM","[0]*x",0.,(b-a)/3.);
  fb->SetParameter(0,fRettPM);
  //fb->Draw("SAME");
  TF1 *fb2 = new TF1("fRettPM2","[0]*x*x",0.,(b-a)/3.);
  fb2->SetParameter(0,fRettPM2);
  fb2->Draw("SAME");
  //  C->cd(3);
  gIntTrap->SetMarkerStyle(20);
  gIntTrap->SetMarkerColor(kBlue);
  gIntTrap->SetMarkerSize(0.5);
  //gIntTrap->Draw("AP");
  gIntTrap->Draw("SAMEP");
  TF1 *fc = new TF1("fTrap","[0]*x*x",0.,(b-a)/3.);
  fc->SetParameter(0,fTrap);
  fc->Draw("SAME");
  C->cd(3);  
  gIntSim->SetMarkerStyle(20);
  gIntSim->SetMarkerColor(kMagenta);
  gIntSim->SetMarkerSize(0.5);
  gIntSim->Draw("AP");
  TF1 *fd = new TF1("fSim","[0]*x*x*x*x",0.,(b-a)/3.);
  fd->SetParameter(0,fSim);
  fd->Draw("SAME");
  C->cd(4);
   
  
  TF1 *ftoInt = new TF1("function",funzR,a,b,0);
  ftoInt->Draw("");
 

  
}