/* check_test_4d.c  check every derivative function */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "deriv.h"
#include "test_4d.h"

#undef abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#undef max
#define max(x,y) ((x)>(y)?(x):(y))

static double xg[10];
static double yg[10];
static double zg[10];
static double tg[10];
static double hx, hy, hz, ht;
static int nx, ny, nz, nt;
static double xmin, xmax, ymin, ymax, zmin, zmax, tmin, tmax;
static double maxerror, gmaxerror;

void check(char name[], 
           double(* deriv)(double x, double y, double z, double t),
           int points, int xord, int yord, int zord, int tord)
{
  double cx[10], cy[10], cz[10], ct[10];
  double p1[10], p2[10], p3[10];
  double analytic, discrete;
  double error, sum_error;
  int i, ii, iii, iiii;
  int j, jj, jjj, jjjj;
  double x, y, z, t, xx;

  maxerror = 0.0;
  for(i=0; i<points; i++)
  {
    rderiv(xord, points, i, hx, cx);
    x = xg[i];
    for(ii=0; ii<points; ii++)
    {
      rderiv(yord, points, ii, hy, cy);
      y = yg[ii];
      for(iii=0; iii<points; iii++)
      {
        rderiv(zord, points, iii, hz, cz);
        z = zg[iii];
        for(iiii=0; iiii<points; iiii++)
        {
          rderiv(tord, points, iiii, ht, ct);
          t = tg[iiii];
          analytic = deriv(x, y, z, t);
          discrete = 0.0;

          /* four cases of one partial x, y, z, t to any order */
	  if(xord!=0 && yord==0 && zord==0 && tord==0)
            for(j=0; j<points; j++)
              discrete += cx[j]*u(xg[j],y,z,t);
          else if(xord==0 && yord!=0 && zord==0 && tord==0)
            for(jj=0; jj<points; jj++)
              discrete += cy[jj]*u(x,yg[jj],z,t);
	  else if(xord==0 && yord==0 && zord!=0 && tord==0)
            for(jjj=0; jjj<points; jjj++)
                discrete += cz[jjj]*u(x,y,zg[jjj],t);
	  else if(xord==0 && yord==0 && zord==0 && tord!=0)
            for(jjjj=0; jjjj<points; jjjj++)
                discrete += ct[jjjj]*u(x,y,z,tg[jjjj]);
          /* six cases of two partials xy, xz, xt, yz, yt, zt to any order */
	  else if(xord!=0 && yord!=0 && zord==0 && tord==0)
	  {
            for(j=0; j<points; j++)
	    {
              p1[j] = 0.0;
              for(jj=0; jj<points; jj++)
                p1[j] += cy[jj]*u(xg[j],yg[jj],z,t);
              discrete += cx[j]*p1[j];
	    }
	  }
	  else if(xord!=0 && yord==0 && zord!=0 && tord==0)
	  {
            for(j=0; j<points; j++)
	    {
              p1[j] = 0.0;
              for(jjj=0; jjj<points; jjj++)
                p1[j] += cz[jjj]*u(xg[j],y,zg[jjj],t);
              discrete += cx[j]*p1[j];
	    }
	  }
	  else if(xord!=0 && yord==0 && zord==0 && tord!=0)
	  {
            for(j=0; j<points; j++)
	    {
              p1[j] = 0.0;
              for(jjjj=0; jjjj<points; jjjj++)
                p1[j] += ct[jjjj]*u(xg[j],y,z,tg[jjjj]);
              discrete += cx[j]*p1[j];
	    }
	  }
	  else if(xord==0 && yord!=0 && zord!=0 && tord==0)
	  {
            for(jj=0; jj<points; jj++)
	    {
              p1[jj] = 0.0;
              for(jjj=0; jjj<points; jjj++)
                p1[jj] += cz[jjj]*u(x,yg[jj],zg[jjj],t);
              discrete += cy[jj]*p1[jj];
	    }
	  }
	  else if(xord==0 && yord!=0 && zord==0 && tord!=0)
	  {
            for(jj=0; jj<points; jj++)
	    {
              p1[jj] = 0.0;
              for(jjjj=0; jjjj<points; jjjj++)
                p1[jj] += ct[jjjj]*u(x,yg[jj],z,tg[jjjj]);
              discrete += cy[jj]*p1[jj];
	    }
	  }
	  else if(xord==0 && yord==0 && zord!=0 && tord!=0)
	  {
            for(jjj=0; jjj<points; jjj++)
	    {
              p1[jjj] = 0.0;
              for(jjjj=0; jjjj<points; jjjj++)
                p1[jjj] += ct[jjjj]*u(x,y,zg[jjj],tg[jjjj]);
              discrete += cz[jjj]*p1[jjj];
	    }
	  }
          /* four cases of three partials xyz, xyt, xzt, yzt to any order */
	  else if(xord!=0 && yord!=0 && zord!=0 && tord==0)
	  {
            for(j=0; j<points; j++)
	    {
              p1[j] = 0.0;
              for(jj=0; jj<points; jj++)
	      {
                p2[jj] = 0.0;
                for(jjj=0; jjj<points; jjj++)
                  p2[jj] += cz[jjj]*u(xg[j],yg[jj],zg[jjj],t);
                p1[j] += cy[jj]*p2[jj];
	      }
              discrete += cx[j]*p1[j];
	    }
	  }
	  else if(xord!=0 && yord!=0 && zord==0 && tord!=0)
	  {
            for(j=0; j<points; j++)
	    {
              p1[j] = 0.0;
              for(jj=0; jj<points; jj++)
	      {
                p2[jj] = 0.0;
                for(jjjj=0; jjjj<points; jjjj++)
                  p2[jj] += ct[jjjj]*u(xg[j],yg[jj],z,tg[jjjj]);
                p1[j] += cy[jj]*p2[jj];
	      }
              discrete += cx[j]*p1[j];
	    }
	  }
	  else if(xord!=0 && yord==0 && zord!=0 && tord!=0)
	  {
            for(j=0; j<points; j++)
	    {
              p1[j] = 0.0;
              for(jjj=0; jjj<points; jjj++)
	      {
                p2[jjj] = 0.0;
                for(jjjj=0; jjjj<points; jjjj++)
                  p2[jjj] += ct[jjjj]*u(xg[j],y,zg[jjj],tg[jjjj]);
                p1[j] += cz[jjj]*p2[jjj];
	      }
              discrete += cx[j]*p1[j];
	    }
	  }
	  else if(xord==0 && yord!=0 && zord!=0 && tord!=0)
	  {
            for(jj=0; jj<points; jj++)
	    {
              p1[jj] = 0.0;
              for(jjj=0; jjj<points; jjj++)
	      {
                p2[jjj] = 0.0;
                for(jjjj=0; jjjj<points; jjjj++)
                  p2[jjj] += ct[jjjj]*u(x,yg[jj],zg[jjj],tg[jjjj]);
                p1[jj] += cz[jjj]*p2[jjj];
	      }
              discrete += cy[jj]*p1[jj];
	    }
	  }
          /* one case of four partials */
	  else
	  {
            for(j=0; j<points; j++)
	    {
              p1[j] = 0.0;
              for(jj=0; jj<points; jj++)
	      {
                p2[jj] = 0.0;
                for(jjj=0; jjj<points; jjj++)
	        {
                  p3[jjj] = 0.0;
                  for(jjjj=0; jjjj<points; jjjj++)
                    p3[jjj] += ct[jjjj]*u(xg[j],yg[jj],zg[jjj],tg[jjjj]);
                  p2[jj] += cz[jjj]*p3[jjj];
		} /* end jjj */
                p1[j] += cy[jj]*p2[jj];
	      } /* end jj */
              discrete += cx[j]*p1[j];
	    } /* end j */
	  } /* end if */
          error = analytic-discrete;
          maxerror = max(abs(error),maxerror);
          if(maxerror>0.001 && i==0 && ii==0 && iii==0 && iiii==0)
            printf("analytic=%f, discrete=%f, error=%g \n",
                   analytic, discrete, error);

	} /* end iiii for t */
      } /* end iii for z */
    } /* end ii for y */
  } /* end i for x */
  printf("                                       %s   maxerror=%g\n",
         name, maxerror);
  gmaxerror = max(maxerror,gmaxerror);
} /* end check */


int main(int argc, char *argv[])
{
  int i;

  xmin = -1.0;
  xmax = 1.0;
  ymin = -1.0;
  ymax = 1.0;
  zmin = -1.0;
  zmax = 1.0;
  tmin = -1.0;
  tmax = 1.0;
  nx = 5;
  ny = 5;
  nz = 5;
  nt = 5;

  hx = (xmax-xmin)/(double)(nx-1);
  hy = (ymax-ymin)/(double)(ny-1);
  hz = (zmax-zmin)/(double)(nz-1);
  ht = (tmax-tmin)/(double)(nt-1);
  for(i=0; i<nx; i++) xg[i] = xmin + hx*(double)i;
  for(i=0; i<ny; i++) yg[i] = ymin + hy*(double)i;
  for(i=0; i<nz; i++) zg[i] = zmin + hz*(double)i;
  for(i=0; i<nt; i++) tg[i] = tmin + ht*(double)i;


  printf("check_test_4d.c running \n");
  gmaxerror = 0.0;
  check("Ux   ", Ux,    5, 1, 0, 0, 0);
  check("Uy   ", Uy,    5, 0, 1, 0, 0);
  check("Uz   ", Uz,    5, 0, 0, 1, 0);
  check("Ut   ", Ut,    5, 0, 0, 0, 1);
  check("Uxx  ", Uxx,   5, 2, 0, 0, 0);
  check("Uxy  ", Uxy,   5, 1, 1, 0, 0);
  check("Uxz  ", Uxz,   5, 1, 0, 1, 0);
  check("Uxt  ", Uxt,   5, 1, 0, 0, 1);
  check("Uyy  ", Uyy,   5, 0, 2, 0, 0);
  check("Uyz  ", Uyz,   5, 0, 1, 1, 0);
  check("Uyt  ", Uyt,   5, 0, 1, 0, 1);
  check("Uzz  ", Uzz,   5, 0, 0, 2, 0);
  check("Uzt  ", Uzt,   5, 0, 0, 1, 1);
  check("Utt  ", Utt,   5, 0, 0, 0, 2);
  check("Uxxx ", Uxxx,  5, 3, 0, 0, 0);
  check("Uxxy ", Uxxy,  5, 2, 1, 0, 0);
  check("Uxxz ", Uxxz,  5, 2, 0, 1, 0);
  check("Uxxt ", Uxxt,  5, 2, 0, 0, 1);
  check("Uxyy ", Uxyy,  5, 1, 2, 0, 0);
  check("Uxyz ", Uxyz,  5, 1, 1, 1, 0);
  check("Uxyt ", Uxyt,  5, 1, 1, 0, 1);
  check("Uxzz ", Uxzz,  5, 1, 0, 2, 0);
  check("Uxzt ", Uxzt,  5, 1, 0, 1, 1);
  check("Uxtt ", Uxtt,  5, 1, 0, 0, 2);
  check("Uyyy ", Uyyy,  5, 0, 3, 0, 0);
  check("Uyyz ", Uyyz,  5, 0, 2, 1, 0);
  check("Uyyt ", Uyyt,  5, 0, 2, 0, 1);
  check("Uyzz ", Uyzz,  5, 0, 1, 2, 0);
  check("Uyzt ", Uyzt,  5, 0, 1, 1, 1);
  check("Uytt ", Uytt,  5, 0, 1, 0, 2);
  check("Uzzz ", Uzzz,  5, 0, 0, 3, 0);
  check("Uzzt ", Uzzt,  5, 0, 0, 2, 1);
  check("Uztt ", Uztt,  5, 0, 0, 1, 2);
  check("Uttt ", Uttt,  5, 0, 0, 0, 3);
  check("Uxxxx", Uxxxx, 5, 4, 0, 0, 0);
  check("Uxxxy", Uxxxy, 5, 3, 1, 0, 0);
  check("Uxxxz", Uxxxz, 5, 3, 0, 1, 0);
  check("Uxxxt", Uxxxt, 5, 3, 0, 0, 1);
  check("Uxxyy", Uxxyy, 5, 2, 2, 0, 0);
  check("Uxxyz", Uxxyz, 5, 2, 1, 1, 0);
  check("Uxxyt", Uxxyt, 5, 2, 1, 0, 1);
  check("Uxxzz", Uxxzz, 5, 2, 0, 2, 0);
  check("Uxxzt", Uxxzt, 5, 2, 0, 1, 1);
  check("Uxxtt", Uxxtt, 5, 2, 0, 0, 2);
  check("Uxyyy", Uxyyy, 5, 1, 3, 0, 0);
  check("Uxyyz", Uxyyz, 5, 1, 2, 1, 0);
  check("Uxyyt", Uxyyt, 5, 1, 2, 0, 1);
  check("Uxyzz", Uxyzz, 5, 1, 1, 2, 0);
  check("Uxyzt", Uxyzt, 5, 1, 1, 1, 1);
  check("Uxytt", Uxytt, 5, 1, 1, 0, 2);
  check("Uxzzz", Uxzzz, 5, 1, 0, 3, 0);
  check("Uxzzt", Uxzzt, 5, 1, 0, 2, 1);
  check("Uxztt", Uxztt, 5, 1, 0, 1, 2);
  check("Uxttt", Uxttt, 5, 1, 0, 0, 3);
  check("Uyyyy", Uyyyy, 5, 0, 4, 0, 0);
  check("Uyyyz", Uyyyz, 5, 0, 3, 1, 0);
  check("Uyyyt", Uyyyt, 5, 0, 3, 0, 1);
  check("Uyyzz", Uyyzz, 5, 0, 2, 2, 0);
  check("Uyyzt", Uyyzt, 5, 0, 2, 1, 1);
  check("Uyytt", Uyytt, 5, 0, 2, 0, 2);
  check("Uyzzz", Uyzzz, 5, 0, 1, 3, 0);
  check("Uyzzt", Uyzzt, 5, 0, 1, 2, 1);
  check("Uyztt", Uyztt, 5, 0, 1, 1, 2);
  check("Uyttt", Uyttt, 5, 0, 1, 0, 3);
  check("Uzzzz", Uzzzz, 5, 0, 0, 4, 0);
  check("Uzzzt", Uzzzt, 5, 0, 0, 3, 1);
  check("Uzztt", Uzztt, 5, 0, 0, 2, 2);
  check("Uzttt", Uzttt, 5, 0, 0, 1, 3);
  check("Utttt", Utttt, 5, 0, 0, 0, 4);
  printf("check_test_4d.c finished, max maxerror=%g \n", gmaxerror);
  return 0;
}