/* test_gaulegf2D.c  test  gaulegf.c integration accuracy */
#include "gaulegf.h"
#include "laphi.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* test functions */

static double f12(double x, double y)
{
  return x+2.0*y+3.0;
}

static double f22(double x, double y)
{
  return x*x+2.0*y*y+3.0*x*y+4.0*x+5.0*y+6.0;
}

static double f32(double x, double y)
{
  return x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+
         5.0*x*x+6.0*y*y+7.0*x*y+8.0*x+9.0*y+10.0;
}

static double f42(double x, double y)
{
  return x*x*x*x+2.0*y*y*y*y+3.0*x*x*x*y+4.0*x*x*y*y+
         5.0*x*x*x*y+6.0*x*x*x+7.0*y*y*y+8.0*x*x*y+
         9.0*x*y*y+10.0*x*x+11.0*y*y+12.0*x*y+13.0*x+
         14.0*y+15.0;
}

static double fe2(double x, double y)
{
  return exp(x+y);
}


int main(int argc, char * argv[])
{
  int n;
  double xmin, xmax, ymin, ymax;
  double xx[100];
  double wx[100];
  double yy[100];
  double wy[100];
  double area, exact; /* area is volume now */
  double dxg, xg[100];
  double dyg, yg[100];
  int i, j, m, mj;

  printf("test_gaulegf2D.c running \n");
  for(n=1; n<6; n++)
  {
    xmin = -1.0;
    xmax =  2.0;
    ymin = -2.0;
    ymax =  1.5;
    printf("calling gaulegf(%f, %f, xx, wx, n) \n", xmin, xmax, n);
    printf("n=%d, equivalent h=%f \n", n, (xmax-xmin)/n);
    gaulegf(xmin, xmax, xx, wx, n);
    for(i=1; i<=n; i++)
      printf("xx[%d]=%f, wx[%d]=%f \n", i, xx[i], i, wx[i]);
    printf("\n");
    printf("calling gaulegf(%f, %f, yy, wy, n) \n", ymin, ymax, n);
    printf("n=%d, equivalent h=%f \n", n, (ymax-ymin)/n);
    gaulegf(ymin, ymax, yy, wy, n);
    for(i=1; i<=n; i++)
      printf("yy[%d]=%f, wy[%d]=%f \n", i, yy[i], i, wy[i]);
    printf("\n");
  }

  printf("f12(x,y)=x+2.0*y+3.0 \n");
  for(n=1; n<6; n++)
  {
    xmin = -1.0;
    xmax =  2.0;
    ymin = -2.0;
    ymax =  1.5;
    gaulegf(xmin, xmax, xx, wx, n);
    gaulegf(ymin, ymax, yy, wy, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        area += wx[i]*wy[j]*f12(xx[i],yy[j]);
    }
    printf("n=%d, n^2=%d, area=%f \n", n, n*n, area);
  }
  printf("\n");

  printf("f22(x,y)=x*x+2.0*y*y+3.0*x*y+4.0*x+5.0*y+6.0; \n");
  for(n=1; n<6; n++)
  {
    xmin = -1.0;
    xmax =  2.0;
    ymin = -2.0;
    ymax =  1.5;
    gaulegf(xmin, xmax, xx, wx, n);
    gaulegf(ymin, ymax, yy, wy, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        area += wx[i]*wy[j]*f22(xx[i],yy[j]);
    }
    printf("n=%d, n^2=%d, area=%f \n", n, n*n, area);
  }
  printf("\n");


  printf("f32(x,y)=x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+ \n");
  printf("         5.0*x*x+6.0*y*y+7.0*x*y+8.0*x+9.0*y+10.0 \n");
  for(n=1; n<6; n++)
  {
    xmin = -1.0;
    xmax =  2.0;
    ymin = -2.0;
    ymax =  1.5;
    gaulegf(xmin, xmax, xx, wx, n);
    gaulegf(ymin, ymax, yy, wy, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        area += wx[i]*wy[j]*f32(xx[i],yy[j]);
    }
    printf("n=%d, n^2=%d, area=%f \n", n, n*n, area);
  }

  printf("f42(x,y)=x*x*x*x+2.0*y*y*y*y+3.0*x*x*x*y+4.0*x*x*y*y+ \n");
  printf("         5.0*x*x*x*y+6.0*x*x*x+7.0*y*y*y+8.0*x*x*y+ \n");
  printf("         9.0*x*y*y+10.0*x*x+11.0*y*y+12.0*x*y+13.0*x+ \n");
  printf("         14.0*y+15.0 \n");
  for(n=1; n<6; n++)
  {
    xmin = -1.0;
    xmax =  2.0;
    ymin = -2.0;
    ymax =  1.5;
    gaulegf(xmin, xmax, xx, wx, n);
    gaulegf(ymin, ymax, yy, wy, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        area += wx[i]*wy[j]*f42(xx[i],yy[j]);
    }
    printf("n=%d, n^2=%d, area=%f \n", n, n*n, area);
  }
  printf("\n");

  printf("fe2(x,y)=exp(x+y) \n");
  for(n=1; n<6; n++)
  {
    xmin = -1.0;
    xmax =  2.0;
    ymin = -2.0;
    ymax =  1.5;
    gaulegf(xmin, xmax, xx, wx, n);
    gaulegf(ymin, ymax, yy, wy, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        area += wx[i]*wy[j]*fe2(xx[i],yy[j]);
    }
    printf("n=%d, n^2=%d, area=%f \n", n, n*n, area);
  }
  printf("\n");

  for(m=3; m<6; m=m+2)
  {
    mj=(m+1)/2;
    printf("phi(x,%d,%d,xg)*phi(y,%d,%d,yg) \n", mj, m, mj, m);
    xmin = -1.0;
    xmax =  2.0;
    dxg = (xmax-xmin)/m;
    if(n==1) printf("phi order=%d, equivh=%d \n", m, dxg);
    for(j=0; j<=m; j++)
    {
      xg[j] = xmin+j*dxg;
      if(n==1) printf("j=%d, xg[j]=%f \n",j, xg[j]);
    }
    ymin = -2.0;
    ymax =  1.5;
    dyg = (ymax-ymin)/m;
    if(n==1) printf("phi order=%d, equivh=%d \n", m, dyg);
    for(j=0; j<=m; j++)
    {
      yg[j] = ymin+j*dyg;
      if(n==1) printf("j=%d, yg[j]=%f \n",j, yg[j]);
    }
    for(n=1; n<8; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          area += wx[i]*wy[j]*phi(xx[i],mj,m,xg)*phi(yy[i],mj,m,yg);
      }
      printf("n=%d, n^2=%d,  area=%f \n", n, n*n, area);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*phi(y,%d,%d,yg]*f12(x,y) \n", mj, m, mj, m);
    for(n=1; n<10; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          area += wx[i]*wy[j]*phi(xx[i],mj,m,xg)*phi(yy[i],mj,m,yg)*
                              f12(xx[i],yy[i]);
      }
      printf("n=%d, n^2=%d, area=%f \n", n, n*n, area);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*phi(y,%d,%d,yg)*f22(x,y) \n", mj, m, mj, m);
    for(n=1; n<10; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          area += wx[i]*wy[j]*phi(xx[i],mj,m,xg)*phi(yy[i],mj,m,yg)*
                              f22(xx[i],yy[i]);
      }
      printf("n=%d, n^2=%d, area=%f \n", n, n*n, area);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*phi(y,%d,%d,yg)*f32(x,y) \n", mj, m, mj, m);
    for(n=1; n<10; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          area += wx[i]*wy[j]*phi(xx[i],mj,m,xg)*phi(yy[i],mj,m,yg)*
                              f32(xx[i],yy[i]);
      }
      printf("n=%d, n^2=%d, area=%f \n", n, n*n, area);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*phi(y,%d,%d,yg)*f42(x,y) \n", mj, m, mj, m);
    for(n=1; n<10; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          area += wx[i]*wy[j]*phi(xx[i],mj,m,xg)*phi(yy[i],mj,m,yg)*
                              f42(xx[i],yy[i]);
      }
      printf("n=%d, n^2=%d, area=%f \n", n, n*n, area);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*phi(y,%d,%d,yg)*fe2(x,y) \n", mj, m, mj, m);
    for(n=1; n<12; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          area += wx[i]*wy[j]*phi(xx[i],mj,m,xg)*phi(yy[i],mj,m,yg)*
                              fe2(xx[i],yy[i]);
      }
      printf("n=%d, n^2=%d, area=%f \n", n, n*n, area);
    }
    printf("\n");
  }

  return 0;
} /* end test_gaulegf2D.c */