#include <cstdlib>
#include <ctime>
#include <cmath>
#include <iostream>
#include <fstream>

using namespace std;


class MonteCarloIntegrator
{
  private:
    // Private variables
    unsigned long iNumberOfPoints;
    double ixLeftLimit, ixRightLimit, iyLeftLimit, iyRightLimit, iYLimit;
    double iResult;
    double (*iIntegrandFunction)(double);
    double (*iDoubleIntegrandFunction)(double, double);
	
  public:
    // Class functions
    MonteCarloIntegrator(double aIntegrandFunction(double), double aLeftLimit,
			 double aRightLimit, double aYLimit, 
                         unsigned long aNumberOfPoints)
    {
      iIntegrandFunction = aIntegrandFunction;
      ixLeftLimit = aLeftLimit;
      ixRightLimit = aRightLimit;
      iYLimit = aYLimit;
      iNumberOfPoints = aNumberOfPoints;
    }

    MonteCarloIntegrator(double aIntegrandFunction(double, double), 
                         double axLeftLimit, double axRightLimit, 
                         double ayLeftLimit, double ayRightLimit, 
                         double aYLimit, unsigned long aNumberOfPoints)
    {
      iDoubleIntegrandFunction = aIntegrandFunction;
      ixLeftLimit = axLeftLimit;
      ixRightLimit = axRightLimit;
      iyLeftLimit = ayLeftLimit;
      iyRightLimit = ayRightLimit;
      iYLimit = aYLimit;
      iNumberOfPoints = aNumberOfPoints;
    }
	
    double integrateByMeanValue() 
    {
      double funcVal, sum, randX;

      for (int i = 0; i < iNumberOfPoints; i++)
      {
	randX = ixLeftLimit + (ixRightLimit - ixLeftLimit) * drand48();
	funcVal = iIntegrandFunction(randX);
	sum += funcVal;
      }
      iResult = ( (ixRightLimit - ixLeftLimit) / 
                          static_cast<double>(iNumberOfPoints) ) * sum;

      return iResult;
    }
		
    double integrateBySamples()
    {
      double boxArea, boxWidth = ixRightLimit - ixLeftLimit;
      double randX, randY;
      unsigned long count = 0;
		
      for (int i = 0; i < iNumberOfPoints; i++)
      {
	randX = ixLeftLimit + boxWidth * drand48();
	randY = iYLimit * drand48();

	if (randY <= iIntegrandFunction(randX))
	  count++;
      }
      boxArea = boxWidth * iYLimit;
      iResult = boxArea * ( static_cast<double>(count) /
                                 static_cast<double>(iNumberOfPoints));

      return iResult;
    }
	
    // New class to handle double integrals
    double doubleintegral()
    {
      double randX, randY;
      double sum = 0;

      for (int i = 0; i < iNumberOfPoints; i++)
      {
	randX = ixLeftLimit + (ixRightLimit - ixLeftLimit) * drand48();
	randY = iyLeftLimit + (iyRightLimit - iyLeftLimit) * drand48();
	sum += iDoubleIntegrandFunction(randX, randY);
      }
      iResult = (ixRightLimit - ixLeftLimit) * (iyRightLimit - iyLeftLimit)  
                           * sum / static_cast<double>(iNumberOfPoints);

      return iResult;
    }

    // Get Functions
    unsigned long numberOfPoints() { return iNumberOfPoints; }
    double xLeftLimit() { return ixLeftLimit; } 
    double xRightLimit() { return ixRightLimit; }
    double yLeftLimit() { return iyLeftLimit; }
    double yRightLimit() { return iyRightLimit; }
    double yLimit() { return iYLimit; }
    double result() { return iResult; }
	 
    // Set Functions
    void setNumberOfPoints(unsigned long aNumberOfPoints)
	 	{ iNumberOfPoints = aNumberOfPoints; }
    void setxLeftLimit(double aLeftLimit)
	 	{ ixLeftLimit = aLeftLimit; }
    void setxRightLimit(double aRightLimit)
		{ ixRightLimit = aRightLimit; }
    void setyLeftLimit(double aLeftLimit)
		{ iyLeftLimit = aLeftLimit; }
    void setyRightLimit(double aRightLimit)
		{ iyRightLimit = aRightLimit; }
    void setYLimit(double aYLimit)
		{ iYLimit = aYLimit; } 
    void setResult(double aResult)
		{ iResult = aResult; }
};


double f(double aX, double aY)
{
  return sin(2 * aX * aY);
}


int main()
{
  ofstream fout ("doubleintegral.dat"); // Set up for file output
  time_t t;
  srand48(t);  // Set up for random vars
	
  // Set up for the integrator
  double xLeftLimit = 0.0, xRightLimit= M_PI;
  double yLeftLimit = 0.0, yRightLimit = M_PI, YLimit = 1.0;

  // Output vars
  double error, result, actual = 1.76112766828;

  // Loop for calculations
  for (unsigned long numberOfPoints = 10; numberOfPoints <= 1e8;
                                                      numberOfPoints *= 10)
  {
    MonteCarloIntegrator MC1(f, xLeftLimit, xRightLimit, yLeftLimit,  
                                      yRightLimit, YLimit, numberOfPoints);
    cout << "The number of MC points is: " << numberOfPoints << endl;
    cout << "The result of the integral is: " 
	 << (result = MC1.doubleintegral() ) << endl;
    cout << "Error: " << (error = fabs(actual - result)) << endl;
    fout << numberOfPoints << " " << error << endl;
  }
}