/* test_gaulegf4D.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 f14(double x, double y, double z, double t)
{
  return x+2.0*y+3.0*z+4.0*t+5.0;
}

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

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

static double f44(double x, double y, double z, double t)
{
  return x*x*x*x+2.0*y*y*y*y+3.0*z*z*z*z+4.0*t*t*t*t+5.0*x*x*x*y+
         5.0*x*x*x*z+6.0*x*x*x*t+7.0*x*x*y*y+8.0*x*x*y*z+
         9.0*x*x*y*t+10.0*x*x*z*z+11.0*x*x*z*t+
         12.0*x*x*x+13.0*y*y*y+14.0*z*z*z+15.0*t*t*t+16.0*x*x*y+
         16.0*x*x*z+17.0*x*x*t+18.0*x*y*y+19.0*x*y*z+20.0*x*y*t+
         21.0*y*y*z+22.0*y*y*t+23.0*y*z*t+24.0*z*z*t+
         25.0*x*x+26.0*y*y+27.0*z*z+28.0*x*y+
         29.0*x*z+30.0*x*t+31.0*y*z+32.0*y*t+33.0*z*t+
         34.0*x+35.0*y+36.0*z+37.0*t+38.0;
}

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


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

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

  printf("f14(x,y,z,t)=x+2.0*y+3.0*z+4.0*t+5.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);
    gaulegf(tmin, tmax, tt, wt, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        for(k=1; k<=n; k++)
          for(p=1; p<=n; p++)
	    area += wx[i]*wy[j]*wz[k]*wt[p]*f14(xx[i],yy[j],zz[k],tt[p]);
    }
    printf("n=%d, n^4=%d, area=%f \n", n, n*n*n*n, area);
  }
  printf("\n");

  printf("f24(x,y,z,t)=x*x+2.0*y*y+3.0*z*z+4.0*t*t+5.0*x*y+6.0*x*z+7.0*x*t+\n");
  printf("        8.0*y*z+9.0*y*t+10.0*z*t+11.0*x+12.0*y+13.0*z+14.0*t+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);
    gaulegf(tmin, tmax, tt, wt, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        for(k=1; k<=n; k++)
          for(p=1; p<=n; p++)
            area += wx[i]*wy[j]*wz[k]*wt[p]*f24(xx[i],yy[j],zz[k],tt[p]);
    }
    printf("n=%d, n^4=%d, area=%f \n", n, n*n*n*n, area);
  }
  printf("\n");


  printf("f34(x,y,z,t)=x*x*x+2.0*y*y*y+3.0*z*z*z+4.0*t*t*t+5.0*x*x*y+6.0*x*x*z+ \n");
  printf("             7.0*x*x*t+8.0*x*y*y+9.0*x*y*z+10.0*x*y*t+11.0*x*z*z+12.0*x*z*t+\n");
  printf("             13.0*x*x+14.0*y*y+15.0*z*z+16.0*t*t+17.0*x*y+18.0*x*z+\n");
  printf("             19.0*x*t+20.0*y*z+21.0*x+22.0*y+23.0*z+24.0*t+25.0\n");
  for(n=1; n<7; n++)
  {
    gaulegf(xmin, xmax, xx, wx, n);
    gaulegf(ymin, ymax, yy, wy, n);
    gaulegf(zmin, zmax, zz, wz, n);
    gaulegf(tmin, tmax, tt, wt, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        for(k=1; k<=n; k++)
          for(p=1; p<=n; p++)
	    area += wx[i]*wy[j]*wz[k]*wt[p]*f34(xx[i],yy[j],zz[k],tt[p]);
    }
    printf("n=%d, n^4=%d, area=%f \n", n, n*n*n*n, area);
  }

  printf("f44(x,y,z,t)=x*x*x*x+2.0*y*y*y*y+3.0*z*z*z*z+4.0*t*t*t*t+5.0*x*x*x*y+\n");
  printf("             5.0*x*x*x*z+6.0*x*x*x*t+7.0*x*x*y*y+8.0*x*x*y*z+\n");
  printf("             9.0*x*x*y*t+10.0*x*x*z*z+11.0*x*x*z*t+\n");
  printf("             12.0*x*x*x+13.0*y*y*y+14.0*z*z*z+15.0*t*t*t+16.0*x*x*y+\n");
  printf("             16.0*x*x*z+17.0*x*x*t+18.0*x*y*y+19.0*x*y*z+20.0*x*y*t+\n");
  printf("             21.0*y*y*z+22.0*y*y*t+23.0*y*z*t+24.0*z*z*t+\n");
  printf("             25.0*x*x+26.0*y*y+27.0*z*z+28.0*x*y+\n");
  printf("             29.0*x*z+30.0*x*t+31.0*y*z+32.0*y*t+33.0*z*t+\n");
  printf("             34.0*x+35.0*y+36.0*z+37.0*t+38.0\n");
  for(n=1; n<7; n++)
  {
    gaulegf(xmin, xmax, xx, wx, n);
    gaulegf(ymin, ymax, yy, wy, n);
    gaulegf(zmin, zmax, zz, wz, n);
    gaulegf(tmin, tmax, tt, wt, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=n; j++)
        for(k=1; k<=n; k++)
          for(p=1; p<=n; p++)
            area += wx[i]*wy[j]*wz[k]*wt[p]*f44(xx[i],yy[j],zz[k],tt[p]);
    }
    printf("n=%d, n^4=%d, area=%f \n", n, n*n*n*n, area);
  }
  printf("\n");

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

    printf("phi(x,%d,%d,xg)*phi(y)*phi(z)*phi(t)*f14(x,y,z,t) \n", mj, m);
    for(n=1; n<12; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      gaulegf(zmin, zmax, zz, wz, n);
      gaulegf(tmin, tmax, tt, wt, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          for(k=1; k<=n; k++)
            for(p=1; p<=n; p++)
            area += wx[i]*wy[j]*wz[k]*wt[p]*phi(xx[i],mj,m,xg)*
                    phi(yy[i],mj,m,yg)*phi(zz[k],mj,m,zg)*
	      phi(tt[p],mj,m,tg)*f14(xx[i],yy[i],zz[k],tt[p]);
      }
      printf("n=%d, n^4=%d, area=%f \n", n, n*n*n*n, area);
    }
    printf("\n");

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

    printf("phi(x,%d,%d,xg)*phi(y)*phi(z)*phi(t)*f34(x,y,z,t) \n", mj, m);
    for(n=1; n<12; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      gaulegf(zmin, zmax, zz, wz, n);
      gaulegf(tmin, tmax, tt, wt, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          for(k=1; k<=n; k++)
            for(p=1; p<=n; p++)
              area += wx[i]*wy[j]*wz[k]*wt[p]*phi(xx[i],mj,m,xg)*
                      phi(yy[i],mj,m,yg)*phi(zz[k],mj,m,zg)*
		      phi(tt[p],mj,m,tg)*f34(xx[i],yy[i],zz[k],tt[p]);
      }
      printf("n=%d, n^4=%d, area=%f \n", n, n*n*n*n, area);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*phi(y)*phi(z)*phi(t)*f44(x,y,z,t) \n", mj, m);
    for(n=1; n<12; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      gaulegf(zmin, zmax, zz, wz, n);
      gaulegf(tmin, tmax, tt, wt, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          for(k=1; k<=n; k++)
            for(p=1; p<=n; p++)
              area += wx[i]*wy[j]*wz[k]*wt[p]*phi(xx[i],mj,m,xg)*
                      phi(yy[i],mj,m,yg)*phi(zz[k],mj,m,zg)*
                      phi(tt[p],mj,m,tg)*f44(xx[i],yy[i],zz[i],tt[p]);
      }
      printf("n=%d, n^4=%d, area=%f \n", n, n*n*n*n, area);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*phi(y)*phi(z)*phi(t)*fe4(x,y,z,t) \n", mj, m);
    for(n=1; n<16; n++)
    {
      gaulegf(xmin, xmax, xx, wx, n);
      gaulegf(ymin, ymax, yy, wy, n);
      gaulegf(zmin, zmax, zz, wz, n);
      gaulegf(tmin, tmax, tt, wt, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=n; j++)
          for(k=1; k<=n; k++)
            for(p=1; p<=n; p++)
              area += wx[i]*wy[j]*wz[k]*wt[p]*phi(xx[i],mj,m,xg)*
                      phi(yy[i],mj,m,yg)*phi(zz[i],mj,m,zg)*
                      phi(tt[p],mj,m,tg)*fe4(xx[i],yy[i],zz[k],tt[k]);
      }
      printf("n=%d, n^4=%d, area=%f \n", n, n*n*n*n, area);
    }
    printf("\n");
  }

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