// nuderiv4d.java  non uniformly spaced derivative coefficients 
//                 see comments in file  nuderiv4d.c
public class nuderiv4d
{
  int debug = 0;
  double AI[][]; // P then inverse 
  int ordpt[] = {0, 5, 15, 35, 70};
  int n;      // matrix size is n by n 
  int norder; // actual maximum order 
  int maxnpwr = 70;
                     // max of 70 points used now, thus max of 4th order 
    int xpow[] = {0, 1, 0, 0, 0, 2, 1, 1, 1, 0, 0, 0, 0, 0, 0,
                  3, 2, 2, 2, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                  4, 3, 3, 3, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};

    int ypow[] = {0, 0, 1, 0, 0, 0, 1, 0, 0, 2, 1, 1, 0, 0, 0,
                  0, 1, 0, 0, 2, 1, 1, 0, 0, 0, 3, 2, 2, 1, 1, 1, 0, 0, 0, 0,
                  0, 1, 0, 0, 2, 1, 1, 0, 0, 0, 3, 2, 2, 1, 1, 1, 0, 0, 0, 0,
                  4, 3, 3, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, 0, 0};

    int zpow[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 2, 1, 0,
                  0, 0, 1, 0, 0, 1, 0, 2, 1, 0, 0, 1, 0, 2, 1, 0, 3, 2, 1, 0,
                  0, 0, 1, 0, 0, 1, 0, 2, 1, 0, 0, 1, 0, 2, 1, 0, 3, 2, 1, 0,
                  0, 1, 0, 2, 1, 0, 3, 2, 1, 0, 4, 3, 2, 1, 0};

    int tpow[] = {0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 2,
                  0, 0, 0, 1, 0, 0, 1, 0, 1, 2, 0, 0, 1, 0, 1, 2, 0, 1, 2, 3,
                  0, 0, 0, 1, 0, 0, 1, 0, 1, 2, 0, 0, 1, 0, 1, 2, 0, 1, 2, 3,
                  0, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4};
  double xpwr[] = new double[6];
  double ypwr[] = new double[6];
  double zpwr[] = new double[6];
  double tpwr[] = new double[6];

// the number of x,y,z,t  npoints  must be at least sum i=1,order+1 of i
// here, order means maximum order of terms used in computing cx, cy, ...
// order 1  5 npoints   only dx, dy, dz, dt
// order 2 15 npoints   only dxx, dxy, dxz, dxt, dyy, dyz, dyt,
//                           dzz, dzt, dtt,     and lower orders
// order 3 35 npoints   only dxxx, dxxy, dxxz, dxxt, dxyy, dxyz,
//                           dxyt, dxzz, dxzt, dxtt, dyyy, dyyz,
//                           dyyt, dyzz, dyzt, dytt, dzzz, dzzt,
//                           dztt, dttt,        and lower orders
// order 4 70 npoints   only dxxxx, dxxxy, dxxxz, dxxxt, dxxyy,
//                           dxxyz, dxxyt, dxxzz, dxxzt, dxxtt,
//                           dxyyy, dxyyz, dxyyt, dxyzz, dxyzt,
//                           dxytt, dxzzz, dxzzt, dxztt, dxttt,
//                           dyyyy, dyyyz, dyyyt, dyyzz, dyyzt,
//                           dyytt, dyzzz, dyzzt, dyztt, dyttt,
//                           dzzzz, dzzzt, dzztt, dzttt, dtttt,
//                                              and lower orders
//
// point is 0..npoint-1 the point about which the derivative is computed
  nuderiv4d()
  {
    System.out.println("nuderiv4d instantiated");
  }

  nuderiv4d(int adebug)
  {
    System.out.println("nuderiv4d instantiated");
    debug = adebug;
    for(int i=0; i<maxnpwr; i++)
    {
      System.out.print(xpow[i]+"  "+ypow[i]+"  "+zpow[i]+"  "+tpow[i]+
                       "   sum="+(xpow[i]+ypow[i]+zpow[i]+tpow[i])+"   ");
      for(int j=0; j<xpow[i]; j++) System.out.print("x");
      for(int j=0; j<ypow[i]; j++) System.out.print("y");
      for(int j=0; j<zpow[i]; j++) System.out.print("z");
      for(int j=0; j<tpow[i]; j++) System.out.print("t");
      System.out.println(" ");
    }
  }

  void nu4dx(int order, int npoint, int point,
             double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cx-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1) continue;
        ci = (double)xpow[i];
        xp = xpow[i]-1; // first derivative
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dx

  void nu4dy(int order, int npoint, int point,
             double x[], double y[], double z[], double t[], double c[])

  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(ypow[i]<1) continue;
        ci = (double)ypow[i];
        xp = xpow[i];
        yp = ypow[i]-1; // first derivative
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dy 

  void nu4dz(int order, int npoint, int point,
             double x[], double y[], double z[], double t[], double c[])

  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(zpow[i]<1) continue;
        ci = (double)zpow[i];
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i]-1; // first derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dz 

  void nu4dt(int order, int npoint, int point,
             double x[], double y[], double z[], double t[], double c[])

  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for ct-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(tpow[i]<1) continue;
        ci = (double)tpow[i];
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dt 

  void nu4dxx(int order, int npoint, int point,
             double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxx-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<2) continue;
        ci = (double)xpow[i]*(double)(xpow[i]-1);
        xp = xpow[i]-2;
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxx 

  void nu4dxy(int order, int npoint, int point,
             double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || ypow[i]<1) continue;
        ci = (double)xpow[i]*(double)(ypow[i]);
        xp = xpow[i]-1;
        yp = ypow[i]-1;
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxy 

  void nu4dxz(int order, int npoint, int point,
              double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || zpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(zpow[i]);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i];
        zp = zpow[i]-1; // first derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxz

  void nu4dxt(int order, int npoint, int point,
              double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || tpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(tpow[i]);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxt

  void nu4dyy(int order, int npoint, int point,
             double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(ypow[i]<2) continue;
        ci = (double)ypow[i]*(double)(ypow[i]-1);
        xp = xpow[i];
        yp = ypow[i]-2;
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyy 

  void nu4dyz(int order, int npoint, int point,
             double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(ypow[i]<1 || zpow[i]<1) continue;
        ci = (double)ypow[i]*(double)(zpow[i]);
        xp = xpow[i];
        yp = ypow[i]-1; // first derivative
        zp = zpow[i]-1; // first derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyz

  void nu4dyt(int order, int npoint, int point,
              double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(ypow[i]<1 || tpow[i]<1) continue;
        ci = (double)ypow[i]*(double)(tpow[i]);
        xp = xpow[i];
        yp = ypow[i]-1; // first derivative
        zp = zpow[i];
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyt

  void nu4dzz(int order, int npoint, int point,
              double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for czz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(zpow[i]<2) continue;
        ci = (double)zpow[i]*(double)(zpow[i]-1);
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i]-2; // second derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dzz 

  void nu4dzt(int order, int npoint, int point,
             double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for czt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(zpow[i]<1 || tpow[i]<1) continue;
        ci = (double)zpow[i]*(double)(tpow[i]);
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i]-1; // first derivative
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dzt 

  void nu4dtt(int order, int npoint, int point,
              double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for ctt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(tpow[i]<2) continue;
        ci = (double)tpow[i]*(double)(tpow[i]-1);
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i]-2; // second derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dtt 

  void nu4dxxx(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxxx-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<3) continue;
        ci = (double)xpow[i]*(double)(xpow[i]-1)*(double)(xpow[i]-2);
        xp = xpow[i]-3;
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxxx 

  void nu4dxxy(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxxy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<2 || ypow[i]<1) continue;
        ci = (double)xpow[i]*(double)(xpow[i]-1)*(double)(ypow[i]);
        xp = xpow[i]-2; // second derivative
        yp = ypow[i]-1; // first derivative
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxxy 

  void nu4dxxz(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxxz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<2 || zpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(xpow[i]-1)*(double)(zpow[i]);
        xp = xpow[i]-2; // second derivative
        yp = ypow[i];
        zp = zpow[i]-1; // first derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxxz 

  void nu4dxxt(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxxt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<2 || tpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(xpow[i]-1)*(double)(tpow[i]);
        xp = xpow[i]-2; // second derivative
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxxt 


  void nu4dxyy(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxyy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || ypow[i]<2) continue;
        ci = (double)xpow[i]*(double)(ypow[i])*(double)(ypow[i]-1);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i]-2; // second derivative
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxyy 

  void nu4dxyz(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxyz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || ypow[i]<1 || zpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(ypow[i])*(double)(zpow[i]);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i]-1; // first derivative
        zp = zpow[i]-1; // first derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxyz

  void nu4dxyt(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxyt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || ypow[i]<1 || tpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(ypow[i])*(double)(tpow[i]);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i]-1; // first derivative
        zp = zpow[i];
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxyt

  void nu4dxzz(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxzz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || zpow[i]<2) continue;
        ci = (double)xpow[i]*(double)(zpow[i])*(double)(zpow[i]-1);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i];
        zp = zpow[i]-2; // second derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxzz

  void nu4dxzt(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxzt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || zpow[i]<1 || tpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(zpow[i])*(double)(tpow[i]);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i];
        zp = zpow[i]-1; // first derivative
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxzt

  void nu4dxtt(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxtt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || tpow[i]<2) continue;
        ci = (double)xpow[i]*(double)(tpow[i])*(double)(tpow[i]-1);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i]-2; // second derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxtt

  void nu4dyyy(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyyy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(ypow[i]<3) continue;
        ci = (double)ypow[i]*(double)(ypow[i]-1)*(double)(ypow[i]-2);
        xp = xpow[i];
        yp = ypow[i]-3; // third derivative
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyyy 

  void nu4dyyz(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyyz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(ypow[i]<2 || zpow[i]<1) continue;
        ci = (double)ypow[i]*(double)(ypow[i]-1)*(double)(zpow[i]);
        xp = xpow[i];
        yp = ypow[i]-2; // second derivative
        zp = zpow[i]-1; // third derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyyz

  void nu4dyyt(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyyt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(ypow[i]<2 || tpow[i]<1) continue;
        ci = (double)ypow[i]*(double)(ypow[i]-1)*(double)(tpow[i]);
        xp = xpow[i];
        yp = ypow[i]-2; // second derivative
        zp = zpow[i];
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyyt

  void nu4dyzz(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyzz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(ypow[i]<1 || zpow[i]<2) continue;
        ci = (double)ypow[i]*(double)(zpow[i])*(double)(zpow[i]-1);
        xp = xpow[i];
        yp = ypow[i]-1; // first derivative
        zp = zpow[i]-2; // second derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyzz

  void nu4dyzt(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyzt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(ypow[i]<1 || zpow[i]<1 || tpow[i]<1) continue;
        ci = (double)ypow[i]*(double)(zpow[i])*(double)(tpow[i]);
        xp = xpow[i];
        yp = ypow[i]-1; // first derivative
        zp = zpow[i]-1; // first derivative
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyzt

  void nu4dytt(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cytt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(ypow[i]<1 || tpow[i]<2) continue;
        ci = (double)ypow[i]*(double)(tpow[i])*(double)(tpow[i]-1);
        xp = xpow[i];
        yp = ypow[i]-1; // first derivative
        zp = zpow[i];
        tp = tpow[i]-2; // second derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dytt

  void nu4dzzz(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for czzz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(zpow[i]<3) continue;
        ci = (double)zpow[i]*(double)(zpow[i]-1)*(double)(zpow[i]-2);
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i]-3; // third derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dzzz 

  void nu4dzzt(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for czzt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(zpow[i]<2 || tpow[i]<1) continue;
        ci = (double)zpow[i]*(double)(zpow[i]-1)*(double)(tpow[i]);
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i]-2; // second derivative
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dzzt 

  void nu4dztt(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cztt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(zpow[i]<1 || tpow[i]<2) continue;
        ci = (double)zpow[i]*(double)(tpow[i])*(double)(tpow[i]-1);
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i]-1; // first derivative
        tp = tpow[i]-2; // second derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dztt 

  void nu4dttt(int order, int npoint, int point,
               double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cttt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(tpow[i]<3) continue;
        ci = (double)tpow[i]*(double)(tpow[i]-1)*(double)(tpow[i]-2);
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i]-3; // third derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dttt 

  void nu4dxxxx(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxxxx-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<4) continue;
        ci = (double)xpow[i]*(double)(xpow[i]-1)*
             (double)(xpow[i]-2)*(double)(xpow[i]-3);
        xp = xpow[i]-4;
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxxxx 

  void nu4dxxxy(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxxxy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<3 || ypow[i]<1) continue;
        ci = (double)xpow[i]*(double)(xpow[i]-1)*
             (double)(xpow[i]-2)*(double)(ypow[i]);
        xp = xpow[i]-3;
        yp = ypow[i]-1;
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxxxy 

  void nu4dxxxz(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxxxz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<3 || zpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(xpow[i]-1)*
             (double)(xpow[i]-2)*(double)(zpow[i]);
        xp = xpow[i]-3; // thrid derivative
        yp = ypow[i];
        zp = zpow[i]-1; // first derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxxxz

  void nu4dxxxt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxxxt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<3 || tpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(xpow[i]-1)*
             (double)(xpow[i]-2)*(double)(tpow[i]);
        xp = xpow[i]-3; // third derivative
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxxxt

  void nu4dxxyy(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxxyy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<2 || ypow[i]<2) continue;
        ci = (double)xpow[i]*(double)(xpow[i]-1)*
             (double)(ypow[i])*(double)(ypow[i]-1);
        xp = xpow[i]-2;
        yp = ypow[i]-2;
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxxyy 

  void nu4dxxyz(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxxyz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<2 || ypow[i]<1 || zpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(xpow[i]-1)*
             (double)(ypow[i])*(double)(zpow[i]);
        xp = xpow[i]-2; // second derivative
        yp = ypow[i]-1; // first derivative
        zp = zpow[i]-1; // first derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxxyz

  void nu4dxxyt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxxyt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<2 || ypow[i]<1 || tpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(xpow[i]-1)*
             (double)(ypow[i])*(double)(tpow[i]);
        xp = xpow[i]-2; // second derivative
        yp = ypow[i]-1; // first derivative
        zp = zpow[i];
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxxyt

  void nu4dxyyy(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxyyy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || ypow[i]<3) continue;
        ci = (double)xpow[i]*(double)(ypow[i])*
             (double)(ypow[i]-1)*(double)(ypow[i]-2);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i]-3; // second derivative
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxyyy 

  void nu4dxyyz(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxyyz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || ypow[i]<2 || zpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(ypow[i])*
             (double)(ypow[i]-1)*(double)(zpow[i]);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i]-2; // second derivative
        zp = zpow[i]-1; // frist derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxyyz

  void nu4dxyyt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxyyt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || ypow[i]<2 || tpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(ypow[i])*
             (double)(ypow[i]-1)*(double)(tpow[i]);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i]-2; // second derivative
        zp = zpow[i];
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxyyt

  void nu4dxyzz(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxyzz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || ypow[i]<1 || zpow[i]<2) continue;
        ci = (double)xpow[i]*(double)(ypow[i])*
             (double)(zpow[i])*(double)(zpow[i]-1);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i]-1; // first derivative
        zp = zpow[i]-2; // second derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxyzz

  void nu4dxyzt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxyzt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || ypow[i]<1 || zpow[i]<1 || tpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(ypow[i])*
             (double)(zpow[i])*(double)(tpow[i]);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i]-1; // first derivative
        zp = zpow[i]-1; // first derivative
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxyzt

  void nu4dxytt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxytt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || ypow[i]<1 || tpow[i]<2) continue;
        ci = (double)xpow[i]*(double)(ypow[i])*
             (double)(tpow[i])*(double)(tpow[i]-1);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i]-1; // first derivative
        zp = zpow[i];
        tp = tpow[i]-2; // second derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxytt

  void nu4dxzzz(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxzzz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || zpow[i]<3) continue;
        ci = (double)xpow[i]*(double)(zpow[i])*
             (double)(zpow[i]-1)*(double)(zpow[i]-2);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i];
        zp = zpow[i]-3; // third derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxzzz

  void nu4dxzzt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxzzt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || zpow[i]<2 || tpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(zpow[i])*
             (double)(zpow[i]-1)*(double)(tpow[i]);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i];
        zp = zpow[i]-2; // second derivative
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxzzt

  void nu4dxztt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxztt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || zpow[i]<1 || tpow[i]<1) continue;
        ci = (double)xpow[i]*(double)(zpow[i])*
             (double)(tpow[i])*(double)(tpow[i]-1);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i];
        zp = zpow[i]-1; // first derivative
        tp = tpow[i]-2; // second derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxztt

  void nu4dxttt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cxttt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      {
        if(i>=n) continue;
        if(xpow[i]<1 || tpow[i]<3) continue;
        ci = (double)xpow[i]*(double)(tpow[i])*
             (double)(tpow[i]-1)*(double)(tpow[i]-2);
        xp = xpow[i]-1; // first derivative
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i]-3; // third derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dxttt

  void nu4dyyyy(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyyyy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(ypow[i]<4) continue;
        ci = (double)ypow[i]*(double)(ypow[i]-1)*
	     (double)(ypow[i]-2)*(double)(ypow[i]-3);
        xp = xpow[i];
        yp = ypow[i]-4;
        zp = zpow[i];
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyyyy 

  void nu4dyyyz(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyyyz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(ypow[i]<3 || zpow[i]<1) continue;
        ci = (double)ypow[i]*(double)(ypow[i]-1)*
	     (double)(ypow[i]-2)*(double)(zpow[i]);
        xp = xpow[i];
        yp = ypow[i]-3; // third derivative
        zp = zpow[i]-1; // first derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyyyz

  void nu4dyyyt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyyyt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(ypow[i]<3 || tpow[i]<1) continue;
        ci = (double)ypow[i]*(double)(ypow[i]-1)*
	     (double)(ypow[i]-2)*(double)(tpow[i]);
        xp = xpow[i];
        yp = ypow[i]-3; // third derivative
        zp = zpow[i];
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyyyt

  void nu4dyyzz(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyyzz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(ypow[i]<2 || zpow[i]<2) continue;
        ci = (double)ypow[i]*(double)(ypow[i]-1)*
	     (double)(zpow[i])*(double)(zpow[i]-1);
        xp = xpow[i];
        yp = ypow[i]-2; // second derivative
        zp = zpow[i]-2; // second derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyyzz

  void nu4dyyzt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyyzt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(ypow[i]<2 || zpow[i]<1 || tpow[i]<1) continue;
        ci = (double)ypow[i]*(double)(ypow[i]-1)*
	     (double)(zpow[i])*(double)(tpow[i]);
        xp = xpow[i];
        yp = ypow[i]-2; // second derivative
        zp = zpow[i]-1; // first derivative
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyyzt

  void nu4dyytt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyytt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(ypow[i]<2 || tpow[i]<2) continue;
        ci = (double)ypow[i]*(double)(ypow[i]-1)*
	     (double)(tpow[i])*(double)(tpow[i]-1);
        xp = xpow[i];
        yp = ypow[i]-2; // second derivative
        zp = zpow[i];
        tp = tpow[i]-2; // second derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyytt

  void nu4dyzzz(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyzzz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(ypow[i]<1 || zpow[i]<3) continue;
        ci = (double)ypow[i]*(double)(zpow[i])*
	     (double)(zpow[i]-1)*(double)(zpow[i]-2);
        xp = xpow[i];
        yp = ypow[i]-1; // first derivative
        zp = zpow[i]-3; // second derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyzzz

  void nu4dyzzt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyzzt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(ypow[i]<1 || zpow[i]<2 || tpow[i]<1) continue;
        ci = (double)ypow[i]*(double)(zpow[i])*
	     (double)(zpow[i]-1)*(double)(tpow[i]);
        xp = xpow[i];
        yp = ypow[i]-1; // first derivative
        zp = zpow[i]-2; // second derivative
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyzzt

  void nu4dyztt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyztt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(ypow[i]<1 || zpow[i]<1 || tpow[i]<2) continue;
        ci = (double)ypow[i]*(double)(zpow[i])*
	     (double)(tpow[i])*(double)(tpow[i]-1);
        xp = xpow[i];
        yp = ypow[i]-1; // first derivative
        zp = zpow[i]-1; // first derivative
        tp = tpow[i]-2; // second derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyztt

  void nu4dyttt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for cyttt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(ypow[i]<1 || tpow[i]<3) continue;
        ci = (double)ypow[i]*(double)(tpow[i])*
	     (double)(tpow[i]-1)*(double)(tpow[i]-2);
        xp = xpow[i];
        yp = ypow[i]-1; // first derivative
        zp = zpow[i];
        tp = tpow[i]-3; // third derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dyttt

  void nu4dzzzz(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for czzzz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(zpow[i]<4) continue;
        ci = (double)zpow[i]*(double)(zpow[i]-1)*
	     (double)(zpow[i]-2)*(double)(zpow[i]-3);
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i]-4; // fourth derivative
        tp = tpow[i];
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dzzzz 

  void nu4dzzzt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for czzzt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(zpow[i]<3 || tpow[i]<1) continue;
        ci = (double)zpow[i]*(double)(zpow[i]-1)*
	     (double)(zpow[i]-2)*(double)(tpow[i]);
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i]-3; // third derivative
        tp = tpow[i]-1; // first derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dzzzt 

  void nu4dzztt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for czztt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(zpow[i]<2 || tpow[i]<2) continue;
        ci = (double)zpow[i]*(double)(zpow[i]-1)*
	     (double)(tpow[i])*(double)(tpow[i]-1);
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i]-2; // second derivative
        tp = tpow[i]-2; // second derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dzztt 

  void nu4dzttt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for czttt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(zpow[i]<1 || tpow[i]<3) continue;
        ci = (double)zpow[i]*(double)(tpow[i])*
	     (double)(tpow[i]-1)*(double)(tpow[i]-2);
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i]-1; // first derivative
        tp = tpow[i]-3; // third
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dzttt 

  void nu4dtttt(int order, int npoint, int point,
                double x[], double y[], double z[], double t[], double c[])
  {
    double ci;
    int xp, yp, zp, tp;

    build(order, npoint, point, x, y, z, t);
    for(int j=0; j<n; j++) // for ctttt-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<maxnpwr; i++)
      { 
        if(i>=n) continue;
        if(tpow[i]<4) continue;
        ci = (double)tpow[i]*(double)(tpow[i]-1)*
	     (double)(tpow[i]-2)*(double)(tpow[i]-3);
        xp = xpow[i];
        yp = ypow[i];
        zp = zpow[i];
        tp = tpow[i]-4; // fourth derivative
        c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp]*tpwr[tp];
      } 
    }
  } // end nu4dtttt 


  void build(int order, int npoint, int point, 
             double x[], double y[], double z[], double t[])
  {
    // check order for n, exact for now  
    n = npoint;
    norder = order;
    if(order>4)
    {
      System.out.println("nuderiv4d.java error, order="+order+" limit is 4");
      norder = 4;
    }
    if(n<ordpt[norder]) 
    {
      System.out.println("nuderiv4d.java error, n="+n+
                         " should be at least  "+ordpt[norder]);
      n = Math.min(n, maxnpwr); // limited here
    }
    AI = new double[n][n];

    for(int i=0; i<n; i++) // build matrix that will be inverted 
    {
      xpwr[0] = 1.0;
      ypwr[0] = 1.0;
      zpwr[0] = 1.0;
      tpwr[0] = 1.0;
      for(int k=1; k<5; k++)
      {
        xpwr[k] = xpwr[k-1]*x[i];
        ypwr[k] = ypwr[k-1]*y[i];
        zpwr[k] = zpwr[k-1]*z[i];
        tpwr[k] = tpwr[k-1]*t[i];
      }
      for(int j=0; j<n; j++)
      {
	AI[i][j] = xpwr[xpow[j]]*ypwr[ypow[j]]*zpwr[zpow[j]]*tpwr[tpow[j]];
      }
    }

    new invert(AI); // invert the matrix 

    // build xpwr, ypwr, zpwr, tpwr at point 
    xpwr[0] = 1.0;
    ypwr[0] = 1.0;
    zpwr[0] = 1.0;
    tpwr[0] = 1.0;
    for(int k=1; k<5; k++)
    {
      xpwr[k] = xpwr[k-1]*x[point];
      ypwr[k] = ypwr[k-1]*y[point];
      zpwr[k] = zpwr[k-1]*z[point];
      tpwr[k] = tpwr[k-1]*t[point];
    }
  } // end build 
} // end class nuderiv4d