/* int2d.c  two dimensional integration to find volume */

#include <stdio.h>
#include <math.h>
#include "gaulegf.h"

static double f(double x, double y)  /* surface */
{
  return sin(x+y)*cos(x*y);
}

int main(int argc, char *argv[])
{
  double x[9], w[9], volume, exact;
                     /* need xy[] , wy[] if different limits             */
                     /* gaulegf(ymin, ymax, xy, wy, n)                   */
                     /* then x[j] becomes xy[j] and w[j] becomes wy[j]   */
  int i, j, n;

  exact = 0.7218583831887; /* analytic solution */
  printf("int2d.c find volume under sin(x+y)*cos(x*y)  \n");
  printf("        for 0 < x < 1   0 < y < 1            \n");
  printf("        using two dimensional Gauss Legendre \n");
  for(n=2; n<=7; n++)
  {
    gaulegf(0.0, 1.0, x, w, n);
    volume = 0.0;
    for(i=1; i<=n; i++) /* index in x direction */
    {
      for(j=1; j<=n; j++) /* index in y direction */
      {
        volume = volume + w[i]*w[j]*f(x[i],x[j]);
      }
    }
    printf("n=%d, volume=%18.15f, exact=%18.15f, err=%g \n",
           n, volume, exact, volume-exact);
  }
  return 0;
} /* end main of int2d.c */