/* rect_int.c  integrate a function over a rectangle            */
/*             this is numerical quadrature                     */
/*             basically compute a volume under surface function*/
/*             compare methods for accuracy                     */

/* first trapezoid rule, varying step size, int f(x,y) dx dx    */
/* then using Gauss-Legendre coordinates and weights            */

/* remember np points or nodes means np*np function evaluations */
/* or (np-1)*(np-1) rectangular elements                        */

#include <stdio.h>
#include <math.h>
#undef  abs
#define abs(x)  ((x)<0.0?(-(x)):(x))

static void gaulegf(double x1, double x2, double x[], double w[], int n);

static double f1(double x, double y)
{
  return x+2.0*y+3.0; /* linear, first order */
}

static double f2(double x, double y)
{
  return x*x+2.0*y*y+3.0*x*y+4.0*x+5.0*y+6.0; /* quadratic, second order */
}

static double f3(double x, double y)
{
  return x*x*x + 2.0*y*y*y + 3.0*x*x*y + 4.0*y*y*x +
         5.0*x*x + 6.0*y*y + 7.0*x*y + 8.0*x +
         9.0*y + 10.0; /* cubic, third order */
}

static double fe(double x, double y)
{
  return exp(x)*exp(y);
}

int main(int argc, char *argv[])
{
  /* simple trapezoidal integration from x=1 to x=3*/
  double volume;
  double x, y, hx, hy;
  double xmin = 1.0;
  double xmax = 3.0;
  double ymin = 0.5;
  double ymax = 2.0;
  double f1exact = 22.5;
  double f2exact = 106.75;
  double f3exact = 397.4375;
  double feexact = 99.69385929;
  double xx[200], wx[200]; /* for Gauss-Legendre */
  double yy[200], wy[200]; /* for Gauss-Legendre */
  int i, j, np;

  printf("rect_int.c integrate to find volume \n");
  printf("Numerical quadrature, x=%g to x=%g, y=%g to y=%g \n",
          xmin, xmax, ymin, ymax);
  printf("\n");

  printf("Using trapezoidal method \n");
  np = 4;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  printf("f1(x,y)=x+2y+3  np=%d,  hx=%g, hy=%g \n", np, hx, hy);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    volume = volume + (f1(x,ymin) + f1(x, ymax))/2.0;
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    volume = volume + (f1(xmin,y) + f1(xmax, y))/2.0;
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + f1(x,y);
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*volume/(double)((np-1)*(np-1));
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f1exact, volume-f1exact);

  np = 16;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  printf("f1(x,y)=x+2y+3  np=%d,  hx=%g, hy=%g \n", np, hx, hy);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    volume = volume + (f1(x,ymin) + f1(x, ymax))/2.0;
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    volume = volume + (f1(xmin,y) + f1(xmax, y))/2.0;
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + f1(x,y);
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*volume/(double)((np-1)*(np-1));
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f1exact, volume-f1exact);

  np = 64;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  printf("f1(x,y)=x+2y+3  np=%d,  hx=%g, hy=%g \n", np, hx, hy);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    volume = volume + (f1(x,ymin) + f1(x, ymax))/2.0;
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    volume = volume + (f1(xmin,y) + f1(xmax, y))/2.0;
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + f1(x,y);
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*volume/(double)((np-1)*(np-1));
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f1exact, volume-f1exact);
  printf("\n");

  np = 4;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  printf("f2(x,y)=x*x+2*y*y+3*x*y+4*x+5*y+6  np=%d,  hx=%g, hy=%g \n",
         np, hx, hy);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    volume = volume + (f2(x,ymin) + f2(x, ymax))/2.0;
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    volume = volume + (f2(xmin,y) + f2(xmax, y))/2.0;
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + f2(x,y);
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*volume/(double)((np-1)*(np-1));
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);

  np = 16;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  printf("f2(x,y)=x*x+2*y*y+3*x*y+4*x+5*y+6  np=%d,  hx=%g, hy=%g \n",
          np, hx, hy);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    volume = volume + (f2(x,ymin) + f2(x, ymax))/2.0;
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    volume = volume + (f2(xmin,y) + f2(xmax, y))/2.0;
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + f2(x,y);
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*volume/(double)((np-1)*(np-1));
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);

  np = 64;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  printf("f2(x,y)=x*x+2*y*y+3*x*y+4*x+5*y+6  np=%d,  hx=%g, hy=%g \n",
          np, hx, hy);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    volume = volume + (f2(x,ymin) + f2(x, ymax))/2.0;
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    volume = volume + (f2(xmin,y) + f2(xmax, y))/2.0;
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + f2(x,y);
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*volume/(double)((np-1)*(np-1));
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);
  printf("\n");

  np = 4;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  printf("f3(x,y)=x*x*x+2*y*y*y+3*x*x*y+4*y*y*x+5*x*x+6*y*y+7*x*y+8*x + \n");
  printf("        9*y+10   np=%d,  hx=%g, hy=%g \n", np, hx, hy);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    volume = volume + (f3(x,ymin) + f3(x, ymax))/2.0;
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    volume = volume + (f3(xmin,y) + f3(xmax, y))/2.0;
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + f3(x,y);
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*volume/(double)((np-1)*(np-1));
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);

  np = 16;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  printf("f3(x,y)=x*x*x+2*y*y*y+3*x*x*y+4*y*y*x+5*x*x+6*y*y+7*x*y+8*x + \n");
  printf("        9*y+10   np=%d,  hx=%g, hy=%g \n", np, hx, hy);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    volume = volume + (f3(x,ymin) + f3(x, ymax))/2.0;
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    volume = volume + (f3(xmin,y) + f3(xmax, y))/2.0;
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + f3(x,y);
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*volume/(double)((np-1)*(np-1));
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);

  np = 64;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  printf("f3(x.y)=x*x*x+2*y*y*y+3*x*x*y+4*y*y*x+5*x*x+6*y*y+7*x*y+8*x + \n");
  printf("        9*y+10   np=%d,  hx=%g, hy=%g \n", np, hx, hy);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    volume = volume + (f3(x,ymin) + f3(x, ymax))/2.0;
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    volume = volume + (f3(xmin,y) + f3(xmax, y))/2.0;
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + f3(x,y);
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*volume/(double)((np-1)*(np-1));
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);
  printf("\n");


  printf("Using Gauss Legendre method \n");

  np = 4;
  printf("f1(x,y)=x+2y+3 np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      volume = volume + wx[i]*wy[j]*f1(xx[i], yy[j]);
    }
  }
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f1exact, volume-f1exact);
  printf("\n");

  np = 4;
  printf("f2(x,y)=x*x+2y*y+3*x*y+4*x+5*y+6   np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      volume = volume + wx[i]*wy[j]*f2(xx[i], yy[j]);
    }
  }
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);
  printf("\n");

  np = 4;
  printf("f3(x,y)=x*x*x+2*y*y*y+3*x*x*y+4*y*y*x+5*x*x+6*y*y+7*x*y+8*x + \n");
  printf("        9*y+10   np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      volume = volume + wx[i]*wy[j]*f3(xx[i], yy[j]);
    }
  }
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);
  printf("\n");

  np = 4;
  printf("fe(x,y)=exp(x)*exp(y)  np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      volume = volume + wx[i]*wy[j]*fe(xx[i], yy[j]);
    }
  }
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);

  np = 8;
  printf("fe(x,y)=exp(x)*exp(y)  np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      volume = volume + wx[i]*wy[j]*fe(xx[i], yy[j]);
    }
  }
  printf("%d function evaluuations \n", np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);
  printf("\n");

  return 0;
} /* end main of rect_int.c */

static void gaulegf(double x1, double x2, double x[], double w[], int n)
{
  int i, j, m;
  double eps = 3.0E-14;
  double p1, p2, p3, pp, xl, xm, z, z1;
  
  m = (n+1)/2;
  xm = 0.5*(x2+x1);
  xl = 0.5*(x2-x1);
  for(i=1; i<=m; i++)
  {
    z = cos(3.141592654*((double)i-0.25)/((double)n+0.5));
    while(1)
    {
      p1 = 1.0;
      p2 = 0.0;
      for(j=1; j<=n; j++)
      {
        p3 = p2;
        p2 = p1;
        p1 = ((2.0*(double)j-1.0)*z*p2-((double)j-1.0)*p3)/
              (double)j;
      }
      pp = (double)n*(z*p1-p2)/(z*z-1.0);
      z1 = z;
      z = z1 - p1/pp;
      if(abs(z-z1) <= eps) break;
    }
    x[i] = xm - xl*z;
    x[n+1-i] = xm + xl*z;
    w[i] = 2.0*xl/((1.0-z*z)*pp*pp);
    w[n+1-i] = w[i];
  }
} /* end gaulegf */
/* end rect_int.c */