// Kelvin Cheung
// 1.12 PS7

#include<math.h>
#include<iostream.h>

#define ONE 0.1127016654
#define TWO 0.8872983345

float lower_f(float x, float y)
     // the lower bound for z
{
  return 0;
}


float upper_f(float x, float y)
     // the upper bound for z
{
  //  return 5*sin(pow(x*y*y/2500,3));
  return 1;
}

float f(float x, float y)
     // returns difference between upper_f and lower_f at a point
{
  return upper_f(x,y)-lower_f(x,y);
}

float three_term_gauss(float lowx, float highx, float y)
     // three-term, one-dimensional Gauss Quadrature integration to
     // determine the volume of the area defined by our bounds from
     // lower_f to upper_f
{
  return (highx-lowx)*(5/18.0*f(ONE*highx+TWO*lowx,y)+
		       4/9.0*f(0.5*highx+0.5*lowx,y)+
		       5/18.0*f(TWO*highx+ONE*lowx,y));
}

float double_three_term_gauss(float lowx, float highx, float lowy, float highy)
     // three-term, two-dimensional Gauss Quadrature integration to
     // determine the volume of the area defined by our bounds from
     // lower_f to upper_f
{
  return (highy-lowy)*(5/18.0*three_term_gauss(lowx,highx,ONE*highy+TWO*lowy)+
		       4/9.0*three_term_gauss(lowx,highx,0.5*highy+0.5*lowy)+
		       5/18.0*three_term_gauss(lowx,highx,TWO*highy+ONE*lowy));
}

main()
{
  cout << "The sand layer volume: " 
       << double_three_term_gauss(0.0,950.0,0.0,500.0) << endl;
}