/* test_gaulegf3D.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 f13(double x, double y, double z)
{
  return x+2.0*y+3.0*z+0.4;
}

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

static double f33(double x, double y, double z)
{
  return x*x*x+2.0*y*y*y+3.0*z*z*z+4.0*x*x*y+5.0*x*y*y+
         7.0*x*y*z+7.0*x*x*z+8.0*x*z*z+9.0*y*y*z+10.0*y*z*z+
         11.0*x*x+12.0*y*y+13.0*z*z+14.0*x*y+15.0*x*z+
         16.0*y*z+17.0*x+18.0*y+19.0*z+20.0;
}

static double f43(double x, double y, double z)
{
  return x*x*x*x+2.0*y*y*y*y+3.0*z*z*z*z+4.0*x*x*x*y+
         5.0*x*x*y*y+6.0*x*x*x*y+7.0*x*x*x*z+8.0*x*x*z*z+
         9.0*x*z*z*z+10.0*y*y*y*z+
         11.0*x*x*x+12.0*y*y*y+13.0*x*x*y+
         14.0*x*y*y+15.0*x*x+16.0*y*y+17.0*x*y+18.0*x+
         19.0*y+20.0*z+21.0;
}

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


int main(int argc, char * argv[])
{
  int n;
  double xmin, xmax, ymin, ymax, zmin, zmax;
  double xx[100];
  double wx[100];
  double yy[100];
  double wy[100];
  double zz[100];
  double wz[100];
  double area, exact; /* area is 4D now */
  double dxg, xg[100];
  double dyg, yg[100];
  double dzg, zg[100];
  int i, j, k, m, mj;

  printf("test_gaulegf3D.c running \n");
  xmin = -1.0;
  xmax =  2.0;
  ymin = -2.0;
  ymax =  1.5;
  zmin = -1.5;
  zmax =  0.5;
  for(n=1; n<6; n++)
  {
    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("calling gaulegf(%f, %f, zz, wz, n) \n", zmin, zmax, n);
    printf("n=%d, equivalent h=%f \n", n, (zmax-zmin)/n);
    gaulegf(zmin, zmax, zz, wz, n);
    for(i=1; i<=n; i++)
      printf("zz[%d]=%f, wz[%d]=%f \n", i, zz[i], i, wz[i]);
    printf("\n");
  }

  printf("f13(x,y,z)=x+2.0*y+3.0 \n");
  for(n=1; n<6; n++)
  {
    gaulegf(xmin, xmax, xx, wx, n);
    gaulegf(ymin, ymax, yy, wy, n);
    gaulegf(zmin, zmax, zz, wz, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        for(k=1; k<=n; k++)
          area += wx[i]*wy[j]*wz[k]*f13(xx[i],yy[j],zz[k]);
    }
    printf("n=%d, n^3=%d, area=%f \n", n, n*n*n, area);
  }
  printf("\n");

  printf("f23(x,y,z)=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++)
  {
    gaulegf(xmin, xmax, xx, wx, n);
    gaulegf(ymin, ymax, yy, wy, n);
    gaulegf(zmin, zmax, zz, wz, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        for(k=1; k<=n; k++)
          area += wx[i]*wy[j]*wz[k]*f23(xx[i],yy[j],zz[k]);
    }
    printf("n=%d, n^3=%d, area=%f \n", n, n*n*n, area);
  }
  printf("\n");


  printf("f33(x,y,z)=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++)
  {
    gaulegf(xmin, xmax, xx, wx, n);
    gaulegf(ymin, ymax, yy, wy, n);
    gaulegf(zmin, zmax, zz, wz, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        for(k=1; k<=n; k++)
	  area += wx[i]*wy[j]*wz[k]*f33(xx[i],yy[j],zz[k]);
    }
    printf("n=%d, n^3=%d, area=%f \n", n, n*n*n, area);
  }

  printf("f43(x,y,z)=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++)
  {
    gaulegf(xmin, xmax, xx, wx, n);
    gaulegf(ymin, ymax, yy, wy, n);
    gaulegf(zmin, zmax, zz, wz, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        for(k=1; k<=n; k++)
          area += wx[i]*wy[j]*wz[k]*f43(xx[i],yy[j],zz[k]);
    }
    printf("n=%d, n^3=%d, area=%f \n", n, n*n*n, area);
  }
  printf("\n");

  printf("fe3(x,y,z)=exp(x+y+z) \n");
  for(n=1; n<6; n++)
  {
    gaulegf(xmin, xmax, xx, wx, n);
    gaulegf(ymin, ymax, yy, wy, n);
    gaulegf(zmin, zmax, zz, wz, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        for(k=1; k<=n; k++)
	  area += wx[i]*wy[j]*wz[k]*fe3(xx[i],yy[j],zz[k]);
    }
    printf("n=%d, n^3=%d, area=%f \n", 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)*phi(z) \n", mj, m);
    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]);
    }
    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]);
    }
    dzg = (zmax-zmin)/m;
    if(n==1) printf("phi order=%d, equivh=%d \n", m, dzg);
    for(j=0; j<=m; j++)
    {
      zg[j] = zmin+j*dzg;
      if(n==1) printf("j=%d, zg[j]=%f \n",j, zg[j]);
    }
    for(n=1; n<8; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      gaulegf(zmin, zmax, zz, wz, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          for(k=1; k<=n; k++)
            area += wx[i]*wy[j]*wz[k]*phi(xx[i],mj,m,xg)*
                    phi(yy[i],mj,m,yg)*phi(zz[i],mj,m,zg);
      }
      printf("n=%d, n^3=%d,  area=%f \n", n, n*n*n, area);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*phi(y)*phi(z)*f13(x,y,z) \n", mj, m);
    for(n=1; n<10; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      gaulegf(zmin, zmax, zz, wz, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          for(k=1; k<=n; k++)
            area += wx[i]*wy[j]*wz[k]*phi(xx[i],mj,m,xg)*phi(yy[i],mj,m,yg)*
	            phi(zz[k],mj,m,zg)*f13(xx[i],yy[i],zz[k]);
      }
      printf("n=%d, n^3=%d, area=%f \n", n, n*n*n, area);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*phi(y)*phi(z)*f23(x,y,z) \n", mj, m);
    for(n=1; n<10; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      gaulegf(zmin, zmax, zz, wz, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          for(k=1; k<=n; k++)
            area += wx[i]*wy[j]*wz[k]*phi(xx[i],mj,m,xg)*phi(yy[i],mj,m,yg)*
	            phi(zz[k],mj,m,zg)*f23(xx[i],yy[i],zz[k]);
      }
      printf("n=%d, n^3=%d, area=%f \n", n, n*n*n, area);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*phi(y)*phi(z)*f33(x,y,z) \n", mj, m);
    for(n=1; n<10; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      gaulegf(zmin, zmax, zz, wz, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          for(k=1; k<=n; k++)
            area += wx[i]*wy[j]*wz[k]*phi(xx[i],mj,m,xg)*phi(yy[i],mj,m,yg)*
	            phi(zz[k],mj,m,zg)*f33(xx[i],yy[i],zz[k]);
      }
      printf("n=%d, n^3=%d, area=%f \n", n, n*n*n, area);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*phi(y)*phi(z)*f43(x,y) \n", mj, m);
    for(n=1; n<10; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      gaulegf(zmin, zmax, zz, wz, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          for(k=1; k<=n; k++)
            area += wx[i]*wy[j]*wz[k]*phi(xx[i],mj,m,xg)*phi(yy[i],mj,m,yg)*
	            phi(zz[k],mj,m,zg)*f43(xx[i],yy[i],zz[i]);
      }
      printf("n=%d, n^3=%d, area=%f \n", n, n*n*n, area);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*phi(y)*phi(z)*fe3(x,y,z) \n", mj, m);
    for(n=1; n<14; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      gaulegf(zmin, zmax, zz, wz, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          for(k=1; k<=n; k++)
            area += wx[i]*wy[j]*wz[k]*phi(xx[i],mj,m,xg)*phi(yy[i],mj,m,yg)*
	            phi(zz[i],mj,m,zg)*fe3(xx[i],yy[i],zz[k]);
      }
      printf("n=%d, n^3=%d, area=%f \n", n, n*n*n, area);
    }
    printf("\n");
  }

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