# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <cmath>

using namespace std;

# include "test_int_2d.H"

int main ( );
void test01 ( );

//****************************************************************************80

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    MAIN is the main program for TEST_INT_2D_PRB.
//
//  Discussion:
//
//    TEST_INT_2D_PRB demonstrates the TEST_INT_2D integration test functions.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    16 January 2009
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );

  cout << "\n";
  cout << "TEST_INT_2D_PRB\n";
  cout << "  C++ version\n";
  cout << "  Test the routines in the TEST_INT_2D library.\n";
 
  test01 ( );

  cout << "\n";
  cout << "TEST_INT_2D_PRB\n";
  cout << "  Normal end of execution.\n";

  cout << "\n";
  timestamp ( );

  return 0;
}
//****************************************************************************80

void test01 ( )

//****************************************************************************80
//
//  Purpose:
//
//    TEST01 applies a Monte Carlo rule.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    16 January 2009
//
//  Author:
//
//    John Burkardt
//
{
  double a[2];
  double b[2];
  int dim;
  double error;
  double exact;
  double *fx;
  int i;
  int j;
  int n;
  int problem;
  int problem_num;
  double quad;
  int seed;
  double volume;
  double *x;

  cout << "\n";
  cout << "TEST01\n";
  cout << "  Use a Monte Carlo rule.\n";
  cout << "\n";
  cout << "  Repeatedly multiply the number of points by 16.\n";

  problem_num = p00_problem_num ( );

  cout << "\n";
  cout << "   Problem      Points         Approx            Exact          Error\n";
  cout << "\n";

  for ( problem = 1; problem <= problem_num; problem++ )
  {
    n = 1;
    for ( i = 1; i <= 6; i++ )
    {
      seed = 123456789;

      x = new double[2*n];
      fx = new double[n];

      x = r8mat_uniform_01 ( 2, n, &seed );

      p00_lim ( problem, a, b );

      for ( dim = 0; dim < 2; dim++ )
      {
        for ( j = 0; j < n; j++ )
        {
          x[dim+j*2] = ( 1.0 - x[dim+j*2] ) * a[dim]
                     +         x[dim+j*2]   * b[dim];
        }
      }
      volume = ( b[1] - a[1] ) * ( b[0] - a[0] );

      p00_fun ( problem, n, x, fx );

      quad = volume * r8vec_sum ( n, fx ) / ( double ) ( n );
   

      exact = p00_exact ( problem );

      error = r8_abs ( quad - exact );

      cout << "  " << setw(8)  << problem
           << "  " << setw(10) << n
           << "  " << setw(14) << quad
           << "  " << setw(14) << exact
           << "  " << setw(14) << error << "\n";

      delete [] fx;
      delete [] x;

      n = n * 16;
    }
    cout << "\n";
  }

  return;
}