// test_nuderiv3dg.java   test  nuderiv3dg  with 28 points 
//          order 1  minimum independent points 4   Ux ... on first degree
//          order 2  minimum independent points 10   Uxx ... on second degree
//          order 3  minimum independent points 20  Uxxx ... on third degree
//          order 4  minimum independent points 35  Uxxxx ... on fourth degree
//          going beyond minimum points does not help much 
class test_nuderiv3dg
{
  test_nuderiv3dg()
  {
    int debug = 0;
    int m;
    int n = 35; // max 4th order this test 
    int ordpt[] = {0, 4, 10, 20, 35};
    double x[] = {0.05, 0.14, 0.27, 0.42, 0.57, 0.71, 0.86,
	          0.04, 0.18, 0.31, 0.43, 0.63, 0.82, 0.95,
	          0.07, 0.16, 0.33, 0.41, 0.59, 0.82, 0.99,
	          0.11, 0.21, 0.35, 0.45, 0.65, 0.75, 0.97,
	          0.06, 0.17, 0.32, 0.46, 0.61, 0.79, 0.98};
    double y[] = {0.06, 0.05, 0.07, 0.10, 0.13, 0.15, 0.17,
	          0.30, 0.33, 0.36, 0.39, 0.31, 0.32, 0.37,
	          0.48, 0.49, 0.52, 0.51, 0.47, 0.50, 0.53,
	          0.65, 0.77, 0.71, 0.67, 0.67, 0.75, 0.76,
	          0.92, 0.97, 0.93, 0.98, 0.87, 0.91, 0.99};
    double z[] = {0.07, 0.28, 0.47, 0.69, 0.91, 0.15, 0.45,
                  0.68, 0.93, 0.05, 0.26, 0.50, 0.97, 0.08,
                  0.14, 0.39, 0.65, 0.76, 0.93, 0.27, 0.56,
                  0.82, 0.99, 0.11, 0.33, 0.72, 0.94, 0.09,
                  0.22, 0.65, 0.95, 0.10, 0.24, 0.42, 0.82};
    double U[] = new double[35];
    double est, actual, dx, dy;
    double cx[] = new double[35];
    double cy[] = new double[35];
    double cz[] = new double[35];
    double cxx[] = new double[35];
    double cxy[] = new double[35];
    double cxz[] = new double[35];
    double cyy[] = new double[35];
    double cyz[] = new double[35];
    double czz[] = new double[35];
    double cxxx[] = new double[35];
    double cxxy[] = new double[35];
    double cxxz[] = new double[35];
    double cxyy[] = new double[35];
    double cxyz[] = new double[35];
    double cxzz[] = new double[35];
    double cyyy[] = new double[35];
    double cyyz[] = new double[35];
    double cyzz[] = new double[35];
    double czzz[] = new double[35];
    double cxxxx[] = new double[35];
    double cxxxy[] = new double[35];
    double cxxxz[] = new double[35];
    double cxxyy[] = new double[35];
    double cxxyz[] = new double[35];
    double cxyyy[] = new double[35];
    double cxyyz[] = new double[35];
    double cxyzz[] = new double[35];
    double cxzzz[] = new double[35];
    double cyyyy[] = new double[35];
    double cyyyz[] = new double[35];
    double cyyzz[] = new double[35];
    double cyzzz[] = new double[35];
    double czzzz[] = new double[35];
    int ip[] = new int[35];
    double error, maxerror;

    nuderiv3dg N = new nuderiv3dg();

    System.out.println("test_nuderiv3dg.java running");
    for(int i=0; i<n; i++)
    {
      U[i] = f3(x[i],y[i],z[i]);
      dx = f3x(x[i],y[i],z[i]);
      dy = f3y(x[i],y[i],z[i]);
      System.out.println("point "+i+"  x="+x[i]+", y="+y[i]+
                         ", U="+U[i]+", dx="+dx+", dy="+dy);
    } 

    // f3(x,y,z) is 3th order, thus 3 is minimum and maximum
    for(int order=3; order<4; order++)
    {
      n = ordpt[order]+1; // give extra point to prevent singular
      maxerror = 0.0;
      for(int point=0; point<n; point++) // Ux 
      {
        m = N.nu3dx(order, n, point, x, y, z, cx, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cx[i]*U[ip[i]];
        actual = f3x(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cx  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Ux 

      for(int point=0; point<n; point++) // Uy 
      {
        m = N.nu3dy(order, n, point, x, y, z, cy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cy[i]*U[ip[i]];
        actual = f3y(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uy 

      for(int point=0; point<n; point++) // Uz 
      {
        m = N.nu3dz(order, n, point, x, y, z, cz, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cz[i]*U[ip[i]];
        actual = f3z(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cz  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uz

      if(n<6) continue;
      for(int point=0; point<n; point++) // Uxx 
      {
        m = N.nu3dxx(order, n, point, x, y, z, cxx, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxx[i]*U[ip[i]];
        actual = f3xx(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxx  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxx 

      for(int point=0; point<n; point++) // Uxy 
      {
        m = N.nu3dxy(order, n, point, x, y, z, cxy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxy[i]*U[ip[i]];
        actual = f3xy(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxy 

      for(int point=0; point<n; point++) // Uxz 
      {
        m = N.nu3dxz(order, n, point, x, y, z, cxz, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxz[i]*U[ip[i]];
        actual = f3xz(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxz  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxz 

      for(int point=0; point<n; point++) // Uyy 
      {
        m = N.nu3dyy(order, n, point, x, y, z, cyy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cyy[i]*U[ip[i]];
        actual = f3yy(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cyy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uyy 

      for(int point=0; point<n; point++) // Uyz 
      {
        m = N.nu3dyz(order, n, point, x, y, z, cyz, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cyz[i]*U[ip[i]];
        actual = f3yz(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cyz  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uyz

      for(int point=0; point<n; point++) // Uzz 
      {
        m = N.nu3dzz(order, n, point, x, y, z, czz, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += czz[i]*U[ip[i]];
        actual = f3zz(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("czz  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uzz

      if(n<10) continue;
      for(int point=0; point<n; point++) // Uxxx 
      {
        m = N.nu3dxxx(order, n, point, x, y, z, cxxx, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxxx[i]*U[ip[i]];
        actual = f3xxx(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxxx  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxxx 

      for(int point=0; point<n; point++) // Uxxy 
      {
        m = N.nu3dxxy(order, n, point, x, y, z, cxxy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxxy[i]*U[ip[i]];
        actual = f3xxy(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxxy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxxy 

      for(int point=0; point<n; point++) // Uxxz 
      {
        m = N.nu3dxxz(order, n, point, x, y, z, cxxz, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxxz[i]*U[ip[i]];
        actual = f3xxz(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxxz  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxxz 

      for(int point=0; point<n; point++) // Uxyy 
      {
        m = N.nu3dxyy(order, n, point, x, y, z, cxyy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxyy[i]*U[ip[i]];
        actual = f3xyy(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxyy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxyy 

      for(int point=0; point<n; point++) // Uxyz 
      {
        m = N.nu3dxyz(order, n, point, x, y, z, cxyz, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxyz[i]*U[ip[i]];
        actual = f3xyz(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxyz  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxyz 

      for(int point=0; point<n; point++) // Uyyy 
      {
        m = N.nu3dyyy(order, n, point, x, y, z, cyyy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cyyy[i]*U[ip[i]];
        actual = f3yyy(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cyyy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uyyy 

      for(int point=0; point<n; point++) // Uyyz 
      {
        m = N.nu3dyyz(order, n, point, x, y, z, cyyz, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cyyz[i]*U[ip[i]];
        actual = f3yyz(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cyyz  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uyyz 

      for(int point=0; point<n; point++) // Uyzz 
      {
        m = N.nu3dyzz(order, n, point, x, y, z, cyzz, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cyzz[i]*U[ip[i]];
        actual = f3yzz(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cyzz  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uyzz 

      for(int point=0; point<n; point++) // Uzzz 
      {
        m = N.nu3dzzz(order, n, point, x, y, z, czzz, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += czzz[i]*U[ip[i]];
        actual = f3zzz(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cyzz  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uzzz 

      if(n<15) continue;
/*
      for(int point=0; point<n; point++) // Uxxxx 
      {
        m = N.nu3dxxxx(order, n, point, x, y, z, cxxxx, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxxxx[i]*U[ip[i]];
        actual = f3xxxx(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxxxx  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxxxx 

      for(int point=0; point<n; point++) // Uxxxy 
      {
        m = N.nu3dxxxy(order, n, point, x, y, z, cxxxy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxxxy[i]*U[ip[i]];
        actual = f3xxxy(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxxxy %2d,  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxxxy 

      for(int point=0; point<n; point++) // Uxxyy 
      {
        m = N.nu3dxxyy(order, n, point, x, y, z, cxxyy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxxyy[i]*U[ip[i]];
        actual = f3xxyy(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxxyy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxxyy 

      for(int point=0; point<n; point++) // Uxyyy 
      {
        m = N.nu3dxyyy(order, n, point, x, y, z, cxyyy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxyyy[i]*U[ip[i]];
        actual = f3xyyy(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxyyy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxyyy 

      for(int point=0; point<n; point++) // Uyyyy 
      {
        m = N.nu3dyyyy(order, n, point, x, y, z, cyyyy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cyyyy[i]*U[ip[i]];
        actual = f3yyyy(x[point],y[point],z[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cyyyy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uyyyy 

*/
      System.out.println("maxerror for order="+order+", n="+n+
                         ", is "+maxerror);
    } // end order loop 
    System.out.println("test_nuderiv3dg.java finished");
  }

  // test functions
  // 1st order test case 
  double f1(double x, double y, double z)
  {
    return x + 2.0*y + 3.0*z + 4.0;
  }
  double f1x(double x, double y, double z)
  {
    return 1.0;
  }
  double f1y(double x, double y, double z)
  {
    return 2.0;
  }
  double f1z(double x, double y, double z)
  {
    return 3.0;
  }

  // 2nd order test case 
  double f2(double x, double y, double z)
  {
    return x*x + 2.0*x*y + 3.0*x*z + 4.0*y*y +
       5.0*y*z + 6.0*z*z + 7.0*x   + 8.0*y   +
       9.0*z   + 10.0;
  }
  double f2x(double x, double y, double z)
  {
    return 2.0*x + 2.0*y + 3.0*z + 7.0;
  }
  double f2y(double x, double y, double z)
  {
    return 2.0*x + 8.0*y + 5.0*z + 8.0;
  }
  double f2z(double x, double y, double z)
  {
   return 3.0*x + 5.0*y + 12.0*z + 9.0;
  }
  double f2xx(double x, double y, double z)
  {
    return 2.0;
  }
  double f2xy(double x, double y, double z)
  {
    return 2.0;
  }
  double f2xz(double x, double y, double z)
  {
    return 3.0;
  }
  double f2yy(double x, double y, double z)
  {
    return 8.0;
  }
  double f2yz(double x, double y, double z)
  {
    return 5.0;
  }
  double f2zz(double x, double y, double z)
  {
    return 12.0;
  }

  // 3rd order test case 
  double f3(double x, double y, double z)
  {
    return x*x*x +  2.0*x*x*y +  3.0*x*x*z + 4.0*x*y*y +
       5.0*x*y*z +  6.0*x*z*z +  7.0*y*y*y + 8.0*y*y*z +
       9.0*y*z*z + 10.0*z*z*z + 11.0*x*x   + 12.0*x*y +
      13.0*x*z   + 14.0*y*y   + 15.0*y*z   + 16.0*z*z  +
      17.0*x     + 18.0*y     + 19.0*z     + 20.0;
  }
  double f3x(double x, double y, double z)
  {
    return 3.0*x*x +  4.0*x*y +  6.0*x*z +  4.0*y*y +
           5.0*y*z +  6.0*z*z +
                                22.0*x   + 12.0*y   +
          13.0*z  +
      17.0;
  }
  double f3y(double x, double y, double z)
  {
    return             2.0*x*x  +             8.0*x*y +
           5.0*x*z +              21.0*y*y + 16.0*y*z +
           9.0*z*z +                         12.0*x   +
                      28.0*y    + 15.0*z   +
      18.0;
  }
  double f3z(double x, double y, double z)
  {
    return                       3.0*x*x +
           5.0*x*y + 12.0*x*z +             8.0*y*y +
          18.0*y*z + 30*z*z   +
          13.0*x   +            15.0*y   + 32.0*z   +
      19.0;
  }
  double f3xx(double x, double y, double z)
  {
    return 6.0*x + 4.0*y +  6.0*z + 22.0;
  }
  double f3xy(double x, double y, double z)
  {
    return 4.0*x + 8.0*y + 5.0*z + 12.0  ;
  }
  double f3xz(double x, double y, double z)
  {
    return 6.0*x + 5.0*y + 12.0*z + 13.0;
  }
  double f3yy(double x, double y, double z)
  {
    return 8.0*x + 42.0*y + 16.0*z + 28.0;
  }
  double f3yz(double x, double y, double z)
  {
    return 5.0*x + 16.0*y + 18.0*z + 15.0;
  }
  double f3zz(double x, double y, double z)
  {
    return 12.0*x + 18.0*y + 60.0*z + 32.0;
  }
  double f3xxx(double x, double y, double z)
  {
    return 6.0;
  }
  double f3xxy(double x, double y, double z)
  {
    return 4.0;
  }
  double f3xxz(double x, double y, double z)
  {
    return 6.0;
  }
  double f3xyy(double x, double y, double z)
  {
    return 8.0;
  }
  double f3xyz(double x, double y, double z)
  {
    return 5.0;
  }
  double f3xzz(double x, double y, double z)
  {
    return 12.0;
  }
  double f3yyy(double x, double y, double z)
  {
    return 42.0;
  }
  double f3yyz(double x, double y, double z)
  {
    return 16.0;
  }
  double f3yzz(double x, double y, double z)
  {
    return 18.0;
  }
  double f3zzz(double x, double y, double z)
  {
    return 60.0;
  }

  // 4th order test case  TBD

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