// nuderiv3dg.java  non uniformly spaced derivative coefficients 
//                  see nuderiv2dg.java for comments

public class nuderiv3dg
{
  int debug = 0;
  double AI[][] = new double[36][36]; // P then inverse
  double P[][]  = new double[36][36]; 
  int ordpt[] = {1, 4, 10, 20, 35};
  int maxnpwr = 35;
                     // max of 35 points used now, thus max of 4th order 
  int xpow[] = {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 ypow[] = {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 zpow[] = {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[5];
  double ypwr[] = new double[5];
  double zpwr[] = new double[5];
  int okind = 0;

  // the number of x,y  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  6 npoints   only dx, dy
  // order 2 10 npoints   only dxx, dxy, dyy and lower orders
  // order 3 20 npoints   only dxxx, dxxy, dxyy, dyyy and lower orders
  // order 4 35 npoints   only dxxxx, dxxxy, dxxyy, dxyyy, dyyyy, lower orders

  nuderiv3dg()
  {
    System.out.println("nuderiv3dg instantiated");
  }

  // m is actual number of entries in c[] and ip[] 
  // xpow[] and ypow[] are correct for the entries in c[] 

  int nu3dx(int order, int npoint, int point,
	    double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

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

  int nu3dy(int order, int npoint, int point,
            double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

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

  int nu3dz(int order, int npoint, int point,
            double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

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

  int nu3dxx(int order, int npoint, int point,
             double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

    m = build(order, npoint, point, x, y, z, ip);
    for(int j=0; j<m; j++) // for cxx-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<m; i++)
      {
	if(xpow[i]>=2) // cxx
	{
          ci = (double)xpow[i]*(double)(xpow[i]-1);
          xp = xpow[i]-2; // second derivative
          yp = ypow[i];
          zp = zpow[i];
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
	}
      } 
    }
    return m;
  } // end nu3dxx 

  int nu3dxy(int order, int npoint, int point,
             double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

    m = build(order, npoint, point, x, y, z, ip);
    for(int j=0; j<m; j++) // for cxy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<m; i++)
      {
	if(xpow[i]>=1 && ypow[i]>=1) // cxy
	{
          ci = (double)xpow[i]*(double)(ypow[i]);
          xp = xpow[i]-1; // first derivative
          yp = ypow[i]-1; // first derivative
          zp = zpow[i];
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
	}
      } 
    }
    return m;
  } // end nu3dxy 

  int nu3dxz(int order, int npoint, int point,
             double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

    m = build(order, npoint, point, x, y, z, ip);
    for(int j=0; j<m; j++) // for cxz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<m; i++)
      {
	if(xpow[i]>=1 && zpow[i]>=1) // cxz
	{
          ci = (double)xpow[i]*(double)(zpow[i]);
          xp = xpow[i]-1; // first derivative
          yp = ypow[i];
          zp = zpow[i]-1; // first derivative
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
	}
      } 
    }
    return m;
  } // end nu3dxz

  int nu3dyy(int order, int npoint, int point,
             double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

    m = build(order, npoint, point, x, y, z, ip);
    for(int j=0; j<m; j++) // for cyy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<m; i++)
      {
	if(ypow[i]>=2)  // cyy
	{
          ci = (double)ypow[i]*(double)(ypow[i]-1);
          xp = xpow[i];
          yp = ypow[i]-2; // second derivative
          zp = zpow[i];
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
	}
      } 
    }
    return m;
  } // end nu3dyy 

  int nu3dyz(int order, int npoint, int point,
             double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

    m = build(order, npoint, point, x, y, z, ip);
    for(int j=0; j<m; j++) // for cyz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<m; i++)
      {
	if(ypow[i]>=1 && zpow[i]>=1)  // cyz
	{
          ci = (double)ypow[i]*(double)zpow[i];
          xp = xpow[i];
          yp = ypow[i]-1;
          zp = zpow[i]-1;
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
	}
      } 
    }
    return m;
  } // end nu3dyz 

  int nu3dzz(int order, int npoint, int point,
             double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

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

  int nu3dxxx(int order, int npoint, int point,
              double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

    m = build(order, npoint, point, x, y, z, ip);
    for(int j=0; j<m; j++) // for cxxx-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<m; i++)
      {
	if(xpow[i]>=3) // cxxx
	{
          ci = (double)xpow[i]*(double)(xpow[i]-1)*(double)(xpow[i]-2);
          xp = xpow[i]-3; // third derivative
          yp = ypow[i];
          zp = zpow[i];
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
	}
      } 
    }
    return m;
  } // end nu3dxxx 

  int nu3dxxy(int order, int npoint, int point,
              double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

    m = build(order, npoint, point, x, y, z, ip);
    for(int j=0; j<m; j++) // for cxxy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<m; i++)
      {
        if(xpow[i]>=2 && ypow[i]>=1) // cxxy
	{
          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];
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
	}
      } 
    }
    return m;
  } // end nu3dxxy 

  int nu3dxxz(int order, int npoint, int point,
              double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

    m = build(order, npoint, point, x, y, z, ip);
    for(int j=0; j<m; j++) // for cxxz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<m; i++)
      {
        if(xpow[i]>=2 && zpow[i]>=1) // cxxz
	{
          ci = (double)xpow[i]*(double)(xpow[i]-1)*(double)(zpow[i]);
          xp = xpow[i]-2;
          yp = ypow[i];
          zp = zpow[i]-1;
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
	}
      } 
    }
    return m;
  } // end nu3dxxz

  int nu3dxyy(int order, int npoint, int point,
              double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

    m = build(order, npoint, point, x, y, z, ip);
    for(int j=0; j<m; j++) // for cxyy-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<m; i++)
      {
	if(xpow[i]>=1 && ypow[i]>=2) // cxyy
	{
          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];
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
	}
      } 
    }
    return m;
  } // end nu3dxyy 

  int nu3dxyz(int order, int npoint, int point,
              double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

    m = build(order, npoint, point, x, y, z, ip);
    for(int j=0; j<m; j++) // for cxyz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<m; i++)
      {
	if(xpow[i]>=1 && ypow[i]>=1 && zpow[i]>=1) // cxyz
	{
          ci = (double)xpow[i]*(double)(ypow[i])*(double)zpow[i];
          xp = xpow[i]-1;
          yp = ypow[i]-1;
          zp = zpow[i]-1;
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
	}
      } 
    }
    return m;
  } // end nu3dxyz

  int nu3dyyy(int order, int npoint, int point,
              double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

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

  int nu3dyyz(int order, int npoint, int point,
              double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

    m = build(order, npoint, point, x, y, z, ip);
    for(int j=0; j<m; j++) // for cyyz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<m; i++)
      {
	if(ypow[i]>=2 && zpow[i]>=1)
	{
          ci = (double)ypow[i]*(double)(ypow[i]-1)*(double)zpow[i];
          xp = xpow[i];
          yp = ypow[i]-2;
          zp = zpow[i]-1;
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
	}
      } 
    }
    return m;
  } // end nu3dyyz 

  int nu3dyzz(int order, int npoint, int point,
              double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

    m = build(order, npoint, point, x, y, z, ip);
    for(int j=0; j<m; j++) // for cyzz-point[u] 
    {
      c[j] = 0.0;
      for(int i=0; i<m; i++)
      {
	if(ypow[i]>=1 && zpow[i]>=2)
	{
	  ci = (double)ypow[i]*(double)zpow[i]*(double)(zpow[i]-1);
          xp = xpow[i];
          yp = ypow[i]-1;
          zp = zpow[i]-2;
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
	}
      } 
    }
    return m;
  } // end nu3dyzz 

  int nu3dzzz(int order, int npoint, int point,
              double x[], double y[], double z[], double c[], int ip[])
  {
    double ci;
    int xp, yp, zp, m;

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





  int build(int order, int npoint, int point,
	    double x[], double y[], double z[], int ip[])
  {
    int sing;
    int n, n2i, nd, ni;
    double xa[] = new double[35];
    double ya[] = new double[35];
    double za[] = new double[35];
    double tmp;
    debug = 0;

    // checks  
    if(point<0 || point>=npoint)
       System.out.println("ERROR nuderiv3dg point outside range\n"); 

    n2i = 1; // good points so far, this will become  n  
    nd  = 0; // number deleted 
    for(int i=0; i<npoint; i++)
    {
      ip[i] = i;
      xa[i] = x[i];
      ya[i] = y[i];
      za[i] = z[i];
    }
    if(point != 0)
    {
      ip[0] = point; // make 'point' be first point, must be included 
      ip[point] = 0;
      tmp = xa[0];
      xa[0] = xa[point];
      xa[point] = tmp;
      tmp = ya[0];
      ya[0] = ya[point];
      ya[point] = tmp;     
      tmp = za[0];
      za[0] = za[point];
      za[point] = tmp;     
    }
    if(okind==1) // reverse order, generally for triangles with bound last
    {
      for(int i=1; i<npoint/2; i++)
      {
	ni = ip[i];
        ip[i] = ip[npoint-i];
        ip[npoint-i] = ni;
        tmp = xa[i];
        xa[i] = xa[npoint-i];
        xa[npoint-i] = tmp;
        tmp = ya[i];
        ya[i] = ya[npoint-i];
        ya[npoint-i] = tmp;             
        tmp = za[i];
        za[i] = za[npoint-i];
        za[npoint-i] = tmp;             
      }
    }
    ni = 2; // trying  ni  points for non singular 
    while(ni<=npoint)
    {
      sing = test_build(ni, xa, ya, za); // P and AI global 
      if(sing==0)
      {
        n2i = ni;
        if(debug>0)
          System.out.println("test_build added point "+point+" for index "+
                             (ni-1)+", nd="+nd);
        if(ni>=35) break;
        ni++;
        continue;
      }
      nd++;
      if(debug>0)
        System.out.println("test_build deleted point "+point+" for index "+
                           (ni-1)+", nd="+nd+", ni+nd="+(ni+nd));
      if(ni+nd>=npoint) break;
      for(int nj=ni-1; nj<npoint-1; nj++)
      {
        ip[nj] = ip[nj+1];
        xa[nj]  = xa[nj+1];
        ya[nj]  = ya[nj+1];
        za[nj]  = za[nj+1];
      }
    }
    n = n2i; // global
    // debug print 
    if(debug>0)
    {
      System.out.println("\n\nactual points selected\n"); 
      for(int i=0; i<n; i++)
        System.out.println("i="+i+", ip[i]="+ip[i]+", xa[i]="+xa[i]+
                           ", ya[i]="+ya[i]+", za[i]="+za[i]);
    }
    // build xpwr, ypwr at point 
    xpwr[0] = 1.0;
    ypwr[0] = 1.0;
    zpwr[0] = 1.0;
    for(int k=1; k<5; k++)
    {
      xpwr[k] = xpwr[k-1]*xa[0]; // 'point' moved to zero 
      ypwr[k] = ypwr[k-1]*ya[0];
      zpwr[k] = zpwr[k-1]*za[0];
    }
    if(n<ordpt[order])
      System.out.println("probable loss of accuracy for order="+order+
                         ", needs "+ordpt[order]+" points, have only "+n);
    return n;
  } // end build 


  int test_build(int ni, double xa[], double ya[], double za[])
                   // P and AI are global 
  {
    int sing;

    for(int i=0; i<ni; i++) // build matrix that will be inverted 
    {
      xpwr[0] = 1.0;
      ypwr[0] = 1.0;
      zpwr[0] = 1.0;
      for(int k=1; k<5; k++)
      {
        xpwr[k] = xpwr[k-1]*xa[i];
        ypwr[k] = ypwr[k-1]*ya[i];
        zpwr[k] = zpwr[k-1]*za[i];
      }
      for(int j=0; j<ni; j++)
      {
        P[i][j] = xpwr[xpow[j]]*ypwr[ypow[j]]*zpwr[zpow[j]];
      }
    }

    sing = inverse(ni, P); // invert the matrix P global 
                           // result in AI[][] global
    return sing;
  } // end test_build 

  int inverse(int n, double A[][])
  {
    int row[] = new int[n];
    int col[] = new int[n];
    double temp[] = new double[n];
    int hold , I_pivot , J_pivot;
    double pivot, abs_pivot;

    // set up row and column interchange vectors
    for(int k=0; k<n; k++)
    {
      row[k] = k ;
      col[k] = k ;
    }
    // begin main reduction loop
    for(int k=0; k<n; k++)
    {
      // find largest element for pivot
      pivot = A[row[k]][col[k]] ;
      I_pivot = k;
      J_pivot = k;
      for(int i=k; i<n; i++)
      {
        for(int j=k; j<n; j++)
        {
          abs_pivot = Math.abs(pivot) ;
          if(Math.abs(A[row[i]][col[j]]) > abs_pivot)
          {
            I_pivot = i ;
            J_pivot = j ;
            pivot = A[row[i]][col[j]] ;
          }
        }
      }
      if(Math.abs(pivot) < 1.0E-12)
      {
        return 1;
      }
      hold = row[k];
      row[k]= row[I_pivot];
      row[I_pivot] = hold ;
      hold = col[k];
      col[k]= col[J_pivot];
      col[J_pivot] = hold ;

      // reduce about pivot
      A[row[k]][col[k]] = 1.0 / pivot ;
      for(int j=0; j<n; j++)
      {
        if(j != k)
        {
          A[row[k]][col[j]] = A[row[k]][col[j]] * A[row[k]][col[k]];
        }
      }
      // inner reduction loop
      for(int i=0; i<n; i++)
      {
        if(k != i)
        {
          for(int j=0; j<n; j++)
          {
            if( k != j )
            {
              A[row[i]][col[j]] = A[row[i]][col[j]] - A[row[i]][col[k]] *
                                   A[row[k]][col[j]] ;
            }
          }
          A[row[i]][col [k]] = - A[row[i]][col[k]] * A[row[k]][col[k]] ;
        }
      }
    } // k
    // end main reduction loop

    // unscramble rows
    for(int j=0; j<n; j++)
    {
      for(int i=0; i<n; i++)
      {
        temp[col[i]] = A[row[i]][j];
      }
      for(int i=0; i<n; i++)
      {
	A[i][j] = temp[i];
      }
    }
    // unscramble columns
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        temp[row[j]] = A[i][col[j]] ;
      }
      for(int j=0; j<n; j++)
      {
	AI[i][j] = temp[j]; // global
      }
    }
    return 0;
  } // end inverse

  public void self_test()
  {
    int error = 0;

    for(int i=0; i<maxnpwr-1; i++)
    {
      for(int j=i+1; j<maxnpwr; j++)
	if(xpow[i]==xpow[j] && ypow[i]==ypow[j] && zpow[i]==zpow[j])
	{
	  error++;
          System.out.println("pow["+i+"]==pow["+j+"]");
	}
      if(i==0)
      {
	if(xpow[i]+ypow[i]+zpow[i]!=0)
	{
	  error++;
          System.out.println(i+" column not zero");
	}
      }
      else if(i<4)
      {
	if(xpow[i]+ypow[i]+zpow[i]!=1)
	{
	  error++;
          System.out.println(i+" column not 1");
	}
      }
      else if(i<10)
      {
	if(xpow[i]+ypow[i]+zpow[i]!=2)
	{
	  error++;
          System.out.println(i+" column not 2");
	}
      }
      else if(i<20)
      {
	if(xpow[i]+ypow[i]+zpow[i]!=3)
	{
	  error++;
          System.out.println(i+" column not 3");
	}
      }
      else
      {
	if(xpow[i]+ypow[i]+zpow[i]!=4)
	{
	  error++;
          System.out.println(i+" column not 4");
	}
      }
    }
    System.out.println("nuderiv3dg.java self test completed with "+
                        error+" errors");
  } // end self_test

} // end class nuderiv3dg