// example of a double integral using GSL integrator
// must be run in compiled mode
// .x exampleDoubleIntegral.C+

#include "Math/GSLIntegrator.h"
#include "Math/WrappedFunction.h"
#include "Math/Functor.h"
#include <cmath>
#include <iostream>


using namespace ROOT::Math;

// starting integrand function

double f (double x, double y) { 
   return x*x + y*y;
} 


// first integrand function (to be integrated in y) 

struct IntegrandFunctionY { 

   IntegrandFunctionY(double x) : fX(x) {}

   // function of y and use x as data member
   double operator() (double y) const { 
      return f(fX,y); 
   }

   double fX; 
};

// second integrand function (to be integrated in x) 

class IntegrandFunctionX { 

public: 
   IntegrandFunctionX() :
      fIntegrator( new GSLIntegrator( 1.E-9,1.E-9,1000) )
   {}

   ~IntegrandFunctionX() { delete fIntegrator; }
 
   // function of y and use x as data member
   double operator() (double x) const { 

      IntegrandFunctionY fy(x);
      // wrapped the function object fFuncY in the required interface by GSLIntegrator
      // use reference
      WrappedFunction<IntegrandFunctionY > wf(fy); 
      double a = - std::sqrt(1. - x*x);
      double b = - a;

      return fIntegrator->Integral(wf,a,b); 
   }
private: 

   //avoid copying 
   IntegrandFunctionX(const IntegrandFunctionX &) {}
   IntegrandFunctionX & operator=(const IntegrandFunctionX &) {return *this;}
   
   mutable GSLIntegrator  * fIntegrator;
};


void exampleDoubleIntegral() { 

   // calculate the integral in x
   GSLIntegrator  ig(1.E-9,1.E-9,1000); 
   IntegrandFunctionX fx;
//   Functor1D func( new IntegrandFunctionX() );
   // use  reference as type to avoid copy
   WrappedFunction<const IntegrandFunctionX &> func(fx);
   double val = ig.Integral(func, -1.,1.);
   std::cout << "integral is " << val << std::endl;
}

int main() { 
   exampleDoubleIntegral();
}