//Author: Amol Holkundkar, Department of Physics, BITS Pilani, Pilani, India.
//This code demonstrate the GSL integration routine on some array. In applied 
//scenario we need to integrate some set of data points. Here,
//the function DoIntegration(from,to,dx,fx,&result,&error) will directly
//return the result for a function which is passed  as vector<double>.
//$ g++ -std=c++11 -lgslcblas -lgsl this_file.cpp 

#include <iostream>
#include <fstream>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <vector>
#include <gsl/gsl_integration.h>
using namespace std;
typedef vector<double> VD;

double dx = 1e-5;       
double a = 0.0;	        /* lower limit a */
double b = M_PI;	/* upper limit b */

double Function(double x)
{
 return x*x*sin(x); //Analytical Integration Can be done, result = pi*pi - 4
}
double exact_result = M_PI*M_PI - 4.0 ; 

void DoIntegration(double from,double to,double dx0,VD fx,double &result,double &error)
{ 

  double abs_err = dx0*100; // to avoid round off errors
  double rel_err = dx0*100; // to avoid roung off errors
  
  VD alpha = fx;
  int N = fx.size();
  
  alpha.push_back(dx0);
     
  gsl_integration_workspace *w = gsl_integration_workspace_alloc (fx.size());
  gsl_function F;
  void *param = &alpha;
  auto integrand = [] (double x, void *p){VD al = *(VD *) p;return al[int(x/al[al.size()-1])];};  
 
  F.function = integrand;
  F.params = param;
     
  gsl_integration_qag(&F,from,to,abs_err,rel_err,N,GSL_INTEG_GAUSS41,w, &result,&error);
  
  gsl_integration_workspace_free (w);
}

 
int main ()
{
  VD fx;
   
  double res,err;
   
  for(double x = a;x <= b;x += dx)
  {
   fx.push_back(  Function(x) ) ;
  }
  
  DoIntegration(a,b,dx,fx,res,err); 
  
  cout << setprecision(15) << scientific << uppercase;   
  cout << endl << "Result           = " <<  res;
  cout << endl << "Exact Result     = " << exact_result;
  cout << endl << "Estimated Error  = " <<  err;
  cout << endl << "Actual Error     = " << abs(res-exact_result);
  cout << endl << endl;   
  
  return 0;
}