// test_nuderiv.java  
//          order 1  minimum independent points 3   Ux ... on first degree
//          order 2  minimum independent points 6   Uxx ... on second degree
//          order 3  minimum independent points 10  Uxxx ... on third degree
//          order 4  minimum independent points 15  Uxxxx ... on fourth degree
class test_nuderiv
{
  test_nuderiv()
  {
    int debug = 0;
    int order, n;
    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;

    double est, actual, hx, hy, hz, ht;
    double xg[] = new double[15];
    double yg[] = new double[15];
    double zg[] = new double[15];
    double tg[] = new double[15];
    double cx[] = new double[15];
    double cy[] = new double[15];
    double cz[] = new double[15];
    double ct[] = new double[15];
    double cxx[] = new double[15];
    double cyy[] = new double[15];
    double czz[] = new double[15];
    double ctt[] = new double[15];
    double cxxx[] = new double[15];
    double cyyy[] = new double[15];
    double czzz[] = new double[15];
    double cttt[] = new double[15];
    double cxxxx[] = new double[15];
    double cyyyy[] = new double[15];
    double czzzz[] = new double[15];
    double ctttt[] = new double[15];
    double err, maxerr=0.0;

    System.out.println("test_nuderiv.java running, testing:");

    System.out.println("order 1  minimum independent points 3   Ux ... on first degree");
    System.out.println("order 2  minimum independent points 6   Uxx ... on second degree");
    System.out.println("order 3  minimum independent points 10  Uxxx ... on third degree");
    System.out.println("order 4  minimum independent points 15  Uxxxx ... on fourth degree");
    System.out.println("f(x,y,z,t):= Math.exp(t)*Math.sin(x + y + z)");
    System.out.println(" ");
    System.out.println("xg, yg, zg, tg are grid coordinates where derivative computed");
    System.out.println("cx, cxx, cxxx, cxxxx are first, second, third, fourth derivative coefficients");
    System.out.println("cy, cyy, cyyy, cyyyy are first, second, third, fourth derivative coefficients");
    System.out.println("cz, czz, czzz, czzzz are first, second, third, fourth derivative coefficients");
    System.out.println("ct, ctt, cttt, ctttt are first, second, third, fourth derivative coefficients");
    System.out.println("check derivative functions fx,fxx,fxxx,fxxxx, then y,z,t ");
    System.out.println(" ");
    System.out.println("derivative is computed at xg[i],yg[ii],zg[iii],tg[iiii]");
    System.out.println("  new nuderiv(order, n, i, xg, cx);");
    System.out.println("  est = 0.0;");
    System.out.println("    for(int j=0; j<n; j++"); 
    System.out.println("      est += cx[j]*f(xg[j],yg[ii],zg[iii],tg[iiii]");
    System.out.println(" ");
    System.out.println(" est approximately fx(xg[i],yg[ii],zg[iii],tg[iiii])");
    System.out.println(" ");
    System.out.println("see source code for order=1,2,3,4  x,y,z,t ");
    System.out.println(" ");

    System.out.println(" ");
    System.out.println(" ");
    System.out.println("Check 1st order x");
    order = 1;
    n = 5;
    maxerr = 0.0;
    System.out.println("order="+order+", n="+n);
    hx = (xmax-xmin)/(n-1);
    System.out.println("hx="+hx);
    hy = (ymax-ymin)/(n-1);
    hz = (zmax-zmin)/(n-1);
    ht = (tmax-tmin)/(n-1);
    for(int i=0; i<n; i++)
    {
      xg[i] = xmin + i*hx;
      System.out.println("xg["+i+"]="+xg[i]);
      yg[i] = ymin + i*hy;
      zg[i] = zmin + i*hz;
      tg[i] = tmin + i*ht;
    }
    for(int point=0; point<n; point++)
    {
      new nuderiv(order, n, point, xg, cx);
      System.out.println("nuderiv("+order+", "+n+", "+point+", xg, cx)");
      for(int j=0; j<n; j++)
	System.out.println("xg["+j+"]="+xg[j]+"  cx["+j+"]="+cx[j]);
      est = 0.0;
      for(int j=0; j<n; j++)
      {
        est += cx[j]*f(xg[j],0.0,0.0,0.0);
        System.out.println("f(xg["+j+"],0,0,0)="+f(xg[j],0.0,0.0,0.0));
      } // end j	
      actual = fx(xg[point],0.0,0.0,0.0);
      err = Math.abs(actual-est);
      maxerr = Math.max(maxerr, err);
      System.out.println("point="+point+", est="+est+", actual="+actual+", err="+err);
    } // point
    System.out.println("maxerr="+maxerr);


    System.out.println(" ");
    System.out.println(" ");
    System.out.println("Check 1st order y");
    order = 1;
    n = 5;
    maxerr = 0.0;
    System.out.println("order="+order+", n="+n);
    hx = (xmax-xmin)/(n-1);
    hy = (ymax-ymin)/(n-1);
    System.out.println("hy="+hy);
    hz = (zmax-zmin)/(n-1);
    ht = (tmax-tmin)/(n-1);
    for(int i=0; i<n; i++)
    {
      xg[i] = xmin + i*hx;
      yg[i] = ymin + i*hy;
      System.out.println("yg["+i+"]="+yg[i]);
      zg[i] = zmin + i*hz;
      tg[i] = tmin + i*ht;
    }
    for(int point=0; point<n; point++)
    {
      new nuderiv(order, n, point, yg, cy);
      System.out.println("nuderiv("+order+", "+n+", "+point+", yg, cy)");
      for(int j=0; j<n; j++)
	System.out.println("yg["+j+"]="+yg[j]+"  cy["+j+"]="+cy[j]);
      est = 0.0;
      for(int j=0; j<n; j++)
      {
	est += cy[j]*f(0.0,yg[j],0.0,0.0);
        System.out.println("f(0,yg["+j+"],0,0)="+f(0.0,yg[j],0.0,0.0));
      } // end j	
      actual = fy(0.0,yg[point],0.0,0.0);
      err = Math.abs(actual-est);
      maxerr = Math.max(maxerr, err);
      System.out.println("point="+point+", est="+est+", actual="+actual+", err="+err);
    } // point
    System.out.println("maxerr="+maxerr);

    System.out.println(" ");
    System.out.println(" ");
    System.out.println("Check 4th order x");
    order = 4;
    n = 15;
    maxerr = 0.0;
    System.out.println("order="+order+", n="+n);
    hx = (xmax-xmin)/(n-1);
    System.out.println("hx="+hx);
    for(int i=0; i<n; i++)
    {
      xg[i] = xmin + i*hx;
      System.out.println("xg["+i+"]="+xg[i]);
    }
    for(int point=0; point<n; point++)
    {
      System.out.println("nuderiv("+order+", "+n+", "+point+", xg, cxxxxx)");
      new nuderiv(order, n, point, xg, cxxxx);
      for(int j=0; j<n; j++)
	System.out.println("xg["+j+"]="+xg[j]+"  cxxxx["+j+"]="+cxxxx[j]);
      est = 0.0;
      for(int j=0; j<n; j++)
        est += cxxxx[j]*f(xg[j],0.0,0.0,0.0);
      actual = fxxxx(xg[point],0.0,0.0,0.0);
      err = Math.abs(actual-est);
      maxerr = Math.max(maxerr, err);
      System.out.println("point="+point+", est="+est+", actual="+actual+", err="+err);
    } // point
    System.out.println("maxerr="+maxerr);


    System.out.println(" ");
    System.out.println(" ");
    System.out.println("Check 1st order in 4 dimensions");
    order = 1;
    n = 5;
    maxerr = 0.0;
    System.out.println("order="+order+", n="+n);
    hx = (xmax-xmin)/(n-1);
    System.out.println("hx="+hx);
    hy = (ymax-ymin)/(n-1);
    hz = (zmax-zmin)/(n-1);
    ht = (tmax-tmin)/(n-1);
    for(int i=0; i<n; i++)
    {
      xg[i] = xmin + i*hx;
      System.out.println("xg["+i+"]="+xg[i]);
      yg[i] = ymin + i*hy;
      zg[i] = zmin + i*hz;
      tg[i] = tmin + i*ht;
    }
    for(int i=0; i<n; i++)
    {
      System.out.println("nuderiv("+order+", "+n+", "+i+", xg, cx)");
      new nuderiv(order, n, i, xg, cx);
      for(int j=0; j<n; j++) System.out.println("cx["+j+"]="+cx[j]);
      for(int ii=0; ii<n; ii++)
      {
        System.out.println("nuderiv("+order+", "+n+", "+ii+", yg, cy)");
        new nuderiv(order, n, ii, yg, cy);
        for(int iii=0; iii<n; iii++)
        {
          System.out.println("nuderiv("+order+", "+n+", "+iii+", zg, cz)");
          new nuderiv(order, n, iii, zg, cz);
          for(int iiii=0; iiii<n; iiii++)
          {
            System.out.println("nuderiv("+order+", "+n+", "+iiii+", tg, ct)");
            new nuderiv(order, n, iiii, tg, ct);
            est = 0.0;
            for(int j=0; j<n; j++) 
              est += cx[j]*f(xg[j],yg[ii],zg[iii],tg[iiii]);
            actual = fx(xg[i],yg[ii],zg[iii],tg[iiii]);
            err = Math.abs(actual-est);
            maxerr = Math.max(maxerr, err);
            System.out.println("i="+i+", ii="+ii+", iii="+iii+", iiii="+iiii+", err="+err);
	  } // iiii
          System.out.println(" ");
          System.out.println(" ");
	} // iii
      } // ii
    } // i
    System.out.println("maxerr="+maxerr);

    
    System.out.println(" ");
    System.out.println(" ");
    System.out.println("Check 4th order in 4 dimensions");
    order = 4;
    n = 15;
    maxerr = 0.0;
    System.out.println("order="+order+", n="+n);
    hx = (xmax-xmin)/(n-1);
    System.out.println("hx="+hx);
    hy = (ymax-ymin)/(n-1);
    hz = (zmax-zmin)/(n-1);
    ht = (tmax-tmin)/(n-1);
    for(int i=0; i<n; i++)
    {
      xg[i] = xmin + i*hx;
      System.out.println("xg["+i+"]="+xg[i]);
      yg[i] = ymin + i*hy;
      zg[i] = zmin + i*hz;
      tg[i] = tmin + i*ht;
    }
    for(int i=0; i<n; i++)
    {
      System.out.println("nuderiv("+order+", "+n+", "+i+", xg, cxxxx)");
      new nuderiv(order, n, i, xg, cxxxx);
      for(int j=0; j<n; j++) System.out.println("cxxxx["+j+"]="+cxxxx[j]);
      for(int ii=0; ii<n; ii++)
      {
        System.out.println("nuderiv("+order+", "+n+", "+ii+", yg, cyyyy)");
        new nuderiv(order, n, ii, yg, cyyyy);
        for(int iii=0; iii<n; iii++)
        {
          System.out.println("nuderiv("+order+", "+n+", "+iii+", zg, czzzz)");
          new nuderiv(order, n, iii, zg, czzzz);
          for(int iiii=0; iiii<n; iiii++)
          {
            System.out.println("nuderiv("+order+", "+n+", "+iiii+", tg, ctttt)");
            new nuderiv(order, n, iiii, tg, ctttt);
            est = 0.0;
            for(int j=0; j<n; j++) 
              est += cxxxx[j]*f(xg[j],yg[ii],zg[iii],tg[iiii]);
            actual = fxxxx(xg[i],yg[ii],zg[iii],tg[iiii]);
            err = Math.abs(actual-est);
            maxerr = Math.max(maxerr, err);
            System.out.println("i="+i+", ii="+ii+", iii="+iii+", iiii="+iiii+", err="+err);
	  } // iiii
	  System.out.println(" ");
	  System.out.println(" ");
	} // iii
      } // ii
    } // i
    System.out.println("maxerr="+maxerr);
    System.out.println(" ");
    System.out.println("test_nuderiv.java finished");
  }
    
  // test functions
  // f(x, y, z, t) := exp(t) sin(x + y + z)    from maple test_nuderiv.mws
  // fx(x, y, z, t) := exp(t) cos(x + y + z)
  // fy(x, y, z, t) := exp(t) cos(x + y + z)
  // fz(x, y, z, t) := exp(t) cos(x + y + z)
  // ft(x, y, z, t) := exp(t) sin(x + y + z)
  // fxx(x, y, z, t) := -exp(t) sin(x + y + z)
  // fyy(x, y, z, t) := -exp(t) sin(x + y + z)
  // fzz(x, y, z, t) := -exp(t) sin(x + y + z)
  // ftt(x, y, z, t) := exp(t) sin(x + y + z)
  // fxxx(x, y, z, t) := -exp(t) cos(x + y + z)
  // fyyy(x, y, z, t) := -exp(t) cos(x + y + z)
  // fzzz(x, y, z, t) := -exp(t) cos(x + y + z)
  // fttt(x, y, z, t) := exp(t) sin(x + y + z)
  // fxxxx(x, y, z, t) := exp(t) sin(x + y + z)
  // fyyyy(x, y, z, t) := exp(t) sin(x + y + z)
  // fzzzz(x, y, z, t) := exp(t) sin(x + y + z)
  // ftttt(x, y, z, t) := exp(t) sin(x + y + z)
  //

  double f(double x, double y, double z, double t)
  {
    return Math.exp(t)*Math.sin(x + y + z);
  }

  double fx(double x, double y, double z, double t) // df(x,y,z,t)/dx
  {
    return Math.exp(t)*Math.cos(x + y + z);
  }

  double fy(double x, double y, double z, double t)
  {
    return Math.exp(t)*Math.cos(x + y + z);
  }

  double fz(double x, double y, double z, double t)
  {
    return Math.exp(t)*Math.cos(x + y + z);
  }

  double ft(double x, double y, double z, double t)
  {
    return Math.exp(t)*Math.sin(x + y + z);
  }

  double fxx(double x, double y, double z, double t)
  {
    return -Math.exp(t)*Math.sin(x + y + z);
  }

  double fyy(double x, double y, double z, double t)
  {
    return -Math.exp(t)*Math.sin(x + y + z);
  }

  double fzz(double x, double y, double z, double t)
  {
    return -Math.exp(t)*Math.sin(x + y + z);
  }

  double ftt(double x, double y, double z, double t)
  {
    return Math.exp(t)*Math.sin(x + y + z);
  }

  double fxxx(double x, double y, double z, double t)
  {
    return -Math.exp(t)*Math.cos(x + y + z);
  }

  double fyyy(double x, double y, double z, double t)
  {
    return -Math.exp(t)*Math.cos(x + y + z);
  }

  double fzzz(double x, double y, double z, double t)
  {
    return -Math.exp(t)*Math.cos(x + y + z);
  }

  double fttt(double x, double y, double z, double t)
  {
    return Math.exp(t)*Math.sin(x + y + z);
  }

  double fxxxx(double x, double y, double z, double t)
  {
    return Math.exp(t)*Math.sin(x + y + z);
  }

  double fyyyy(double x, double y, double z, double t)
  {
    return Math.exp(t)*Math.sin(x + y + z);
  }

  double fzzzz(double x, double y, double z, double t)
  {
    return Math.exp(t)*Math.sin(x + y + z);
  }

  double ftttt(double x, double y, double z, double t)
  {
      return Math.exp(t)*Math.sin(x + y + z);
  }

  public static void main (String[] args)
  {
    new test_nuderiv();
  }
} // end  class test_nuderiv