/* test_nav1_deriv.java */

#include <stdlib.h>
#include <stdio.h>
#include "nav1.h"
#include "nav1_deriv.h"

int main(int argc, char* argv[])
{
  /* PDE definition */
  nav1();

  /* solution domain */
  double xmin = 0.0;
  double xmax = 1.0;
  double ymin = 0.0;
  double ymax = 1.0;
  double zmin = 0.0;
  double zmax = 1.0;
  double tmin = 0.0;
  double tmax = 1.0;
  int nx = 5; /* minimum size */
  int ny = 5;
  int nz = 5;
  int nt = 5;

  /* working storage */
  double xc[81];  /* first derivative coeficients */
  double yc[81];  /* including a[]/(b*h) */
  double zc[81];
  double tc[81];
  double xxc[81]; /* second derivative coeficients */
  double yyc[81]; /* including a[]/(b*h^2) */
  double zzc[81];
  double x, y, z, t;
  int i, j, k, m, ii;
  double hx, hy, hz, ht;
  double xd, yd, zd, td;
  double myd, mydd, navd, navdd;

  printf("test_nav1_deriv running \n");
  test(); /* from nav1_deriv */
  hx = (xmax-xmin)/(nx-1);
  hy = (ymax-ymin)/(ny-1);
  hz = (zmax-zmin)/(nz-1);
  ht = (tmax-tmin)/(nt-1);

  printf("nx=%d, xmin=%6.3f, xmax=%6.3f, hx=%8.6f \n",nx,xmin,xmax,hx);
  printf("ny=%d, ymin=%6.3f, ymax=%6.3f, hy=%8.6f \n",ny,ymin,ymax,hy);
  printf("nz=%d, zmin=%6.3f, zmax=%6.3f, hz=%8.6f \n",nz,zmin,zmax,hz);
  printf("nt=%d, tmin=%6.3f, tmax=%6.3f, ht=%8.6f \n",nt,tmin,tmax,ht);

  printf("\n");

  printf("dx direction \n");
  nav1_deriv(1,nx,xmin,xmax,xc);  /* set up descrete derivative */
  nav1_deriv(2,nx,xmin,xmax,xxc); /* first order, second order, x */
  printf("dy direction \n");
  nav1_deriv(1,ny,ymin,ymax,yc);
  nav1_deriv(2,ny,ymin,ymax,yyc);
  printf("dz direction \n");
  nav1_deriv(1,nz,zmin,zmax,zc);
  nav1_deriv(2,nz,zmin,zmax,zzc);
  printf("dt direction \n");
  nav1_deriv(1,nt,tmin,tmax,tc);
  /* end setup_deriv */

  printf("test du/dx  ddu/ddx \n");
  y=ymin+hy;
  z=zmin+hz;
  t=tmin+ht;
  for( i=0; i<nx; i++)
  {
    x = xmin+i*hx;
    myd = 0.0;
    mydd = 0.0;
    for( ii=0; ii<nx; ii++)
    {
      xd = xmin+ii*hx;
      myd = myd + xc[i*nx+ii]*u(xd,y,z,t); 
      mydd = mydd + xxc[i*nx+ii]*u(xd,y,z,t); 
    }
    navd = ux(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
    navdd = uxx(x,y,z,t);
    printf("mydd=%15.9f, navdd=%15.9f, error=%g \n",mydd,navdd,(navdd-mydd));
  } /* end du/dx */
  printf("\n");

  printf("test dv/dx  ddv/ddx \n");
  y=ymin+hy;
  z=zmin+hz;
  t=tmin+ht;
  for( i=0; i<nx; i++)
  {
    x = xmin+i*hx;
    myd = 0.0;
    mydd = 0.0;
    for( ii=0; ii<nx; ii++)
    {
      xd = xmin+ii*hx;
      myd = myd + xc[i*nx+ii]*v(xd,y,z,t); 
      mydd = mydd + xxc[i*nx+ii]*v(xd,y,z,t); 
    }
    navd = vx(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
    navdd = vxx(x,y,z,t);
    printf("mydd=%15.9f, navdd=%15.9f, error=%g \n",mydd,navdd,(navdd-mydd));
  } /* end dv/dx */
  printf("\n");

  printf("test dw/dx  ddw/ddx \n");
  y=ymin+hy;
  z=zmin+hz;
  t=tmin+ht;
  for( i=0; i<nx; i++)
  {
    x = xmin+i*hx;
    myd = 0.0;
    mydd = 0.0;
    for( ii=0; ii<nx; ii++)
    {
      xd = xmin+ii*hx;
      myd = myd + xc[i*nx+ii]*w(xd,y,z,t); 
      mydd = mydd + xxc[i*nx+ii]*w(xd,y,z,t); 
    }
    navd = wx(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
    navdd = wxx(x,y,z,t);
    printf("mydd=%15.9f, navdd=%15.9f, error=%g \n",mydd,navdd,(navdd-mydd));
  } /* end dw/dx */
  printf("\n");

  printf("test dp/dx \n");
  y=ymin+hy;
  z=zmin+hz;
  t=tmin+ht;
  for( i=0; i<nx; i++)
  {
    x = xmin+i*hx;
    myd = 0.0;
    for( ii=0; ii<nx; ii++)
    {
      xd = xmin+ii*hx;
      myd = myd + xc[i*nx+ii]*p(xd,y,z,t); 
    }
    navd = px(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
  } /* end dp/dx */
  printf("\n");

  printf("test du/dy  ddu/ddy \n");
  x=xmin+hx;
  z=zmin+hz;
  t=tmin+ht;
  for( i=0; i<ny; i++)
  {
    y = ymin+i*hy;
    myd = 0.0;
    mydd = 0.0;
    for( ii=0; ii<ny; ii++)
    {
      yd = ymin+ii*hy;
      myd = myd + yc[i*ny+ii]*u(x,yd,z,t); 
      mydd = mydd + yyc[i*ny+ii]*u(x,yd,z,t); 
    }
    navd = uy(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
    navdd = uyy(x,y,z,t);
    printf("mydd=%15.9f, navdd=%15.9f, error=%g \n",mydd,navdd,(navdd-mydd));
  } /* end du/dy */
  printf("\n");

  printf("test dv/dy  ddv/ddy \n");
  x=xmin+hx;
  z=zmin+hz;
  t=tmin+ht;
  for( i=0; i<ny; i++)
  {
    y = ymin+i*hy;
    myd = 0.0;
    mydd = 0.0;
    for( ii=0; ii<ny; ii++)
    {
      yd = ymin+ii*hy;
      myd = myd + yc[i*ny+ii]*v(x,yd,z,t); 
      mydd = mydd + yyc[i*ny+ii]*v(x,yd,z,t); 
    }
    navd = vy(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
    navdd = vyy(x,y,z,t);
    printf("mydd=%15.9f, navdd=%15.9f, error=%g \n",mydd,navdd,(navdd-mydd));
  } /* end dv/dy */
  printf("\n");

  printf("test dw/dy  ddw/ddy \n");
  x=xmin+hx;
  z=zmin+hz;
  t=tmin+ht;
  for( i=0; i<ny; i++)
  {
    y = ymin+i*hy;
    myd = 0.0;
    mydd = 0.0;
    for( ii=0; ii<ny; ii++)
    {
      yd = ymin+ii*hy;
      myd = myd + yc[i*ny+ii]*w(x,yd,z,t); 
      mydd = mydd + yyc[i*ny+ii]*w(x,yd,z,t); 
    }
    navd = wy(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
    navdd = wyy(x,y,z,t);
    printf("mydd=%15.9f, navdd=%15.9f, error=%g \n",mydd,navdd,(navdd-mydd));
  } /* end dw/dy */
  printf("\n");

  printf("test dp/dy \n");
  x=xmin+hx;
  z=zmin+hz;
  t=tmin+ht;
  for( i=0; i<ny; i++)
  {
    y = ymin+i*hy;
    myd = 0.0;
    for( ii=0; ii<ny; ii++)
    {
      yd = ymin+ii*hy;
      myd = myd + yc[i*ny+ii]*p(x,yd,z,t); 
    }
    navd = py(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
  } /* end dp/dy */
  printf("\n");

  printf("test du/dz  ddu/ddz \n");
  x=xmin+hx;
  y=ymin+hy;
  t=tmin+ht;
  for( i=0; i<nz; i++)
  {
    z = zmin+i*hz;
    myd = 0.0;
    mydd = 0.0;
    for( ii=0; ii<nz; ii++)
    {
      zd = zmin+ii*hz;
      myd = myd + zc[i*nz+ii]*u(x,y,zd,t); 
      mydd = mydd + zzc[i*nz+ii]*u(x,y,zd,t); 
    }
    navd = uz(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
    navdd = uzz(x,y,z,t);
    printf("mydd=%15.9f, navdd=%15.9f, error=%g \n",mydd,navdd,(navdd-mydd));
  } /* end du/dz */
  printf("\n");

  printf("test dv/dz  ddv/ddz \n");
  x=xmin+hx;
  y=ymin+hy;
  t=tmin+ht;
  for( i=0; i<nz; i++)
  {
    z = zmin+i*hz;
    myd = 0.0;
    mydd = 0.0;
    for( ii=0; ii<nz; ii++)
    {
      zd = zmin+ii*hz;
      myd = myd + zc[i*nz+ii]*v(x,y,zd,t); 
      mydd = mydd + zzc[i*nz+ii]*v(x,y,zd,t); 
    }
    navd = vz(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
    navdd = vzz(x,y,z,t);
    printf("mydd=%15.9f, navdd=%15.9f, error=%g \n",mydd,navdd,(navdd-mydd));
  } /* end dv/dz */
  printf("\n");

  printf("test dw/dz  ddw/ddz \n");
  x=xmin+hx;
  y=ymin+hy;
  t=tmin+ht;
  for( i=0; i<nz; i++)
  {
    z = zmin+i*hz;
    myd = 0.0;
    mydd = 0.0;
    for( ii=0; ii<nz; ii++)
    {
      zd = zmin+ii*hz;
      myd = myd + zc[i*nz+ii]*w(x,y,zd,t); 
      mydd = mydd + zzc[i*nz+ii]*w(x,y,zd,t); 
    }
    navd = wz(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
    navdd = wzz(x,y,z,t);
    printf("mydd=%15.9f, navdd=%15.9f, error=%g \n",mydd,navdd,(navdd-mydd));
  } /* end dw/dz */
  printf("\n");

  printf("test dp/dz \n");
  x=xmin+hx;
  y=ymin+hy;
  t=tmin+ht;
  for( i=0; i<nz; i++)
  {
    z = zmin+i*hz;
    myd = 0.0;
    for( ii=0; ii<nz; ii++)
    {
      zd = zmin+ii*hz;
      myd = myd + zc[i*nz+ii]*p(x,y,zd,t); 
    }
    navd = pz(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
  } /* end dp/dz */
  printf("\n");

  printf("test du/dt \n");
  x=xmin+hx;
  y=ymin+hy;
  z=zmin+hz;
  for( i=0; i<nt; i++)
  {
    t = tmin+i*ht;
    myd = 0.0;
    for( ii=0; ii<nt; ii++)
    {
      td = tmin+ii*ht;
      myd = myd + tc[i*nt+ii]*u(x,y,z,td); 
    }
    navd = ut(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
  } /* end du/dt */
  printf("\n");

  printf("test dv/dt \n");
  x=xmin+hx;
  y=ymin+hy;
  z=zmin+hz;
  for( i=0; i<nt; i++)
  {
    t = tmin+i*ht;
    myd = 0.0;
    for( ii=0; ii<nt; ii++)
    {
      td = tmin+ii*ht;
      myd = myd + tc[i*nt+ii]*v(x,y,z,td); 
    }
    navd = vt(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
  } /* end dv/dt */
  printf("\n");

  printf("test dw/dt \n");
  x=xmin+hx;
  y=ymin+hy;
  z=zmin+hz;
  for( i=0; i<nt; i++)
  {
    t = tmin+i*ht;
    myd = 0.0;
    for( ii=0; ii<nt; ii++)
    {
      td = tmin+ii*ht;
      myd = myd + tc[i*nt+ii]*w(x,y,z,td); 
    }
    navd = wt(x,y,z,t);
    printf("myd=%15.9f, navd=%15.9f, error=%g \n",myd,navd,(navd-myd));
  } /* end dw/dt */
  printf("\n");

  printf("end test_nav1_deriv \n");
  return 0;
} /* end test_nav1_deriv */