// nuderiv2dg.java  non uniformly spaced derivative coefficients 
//    given  x0,y0, x1,y1, x2,y2, ... xn,yn  non uniformly spaced
//    find the coefficients  c???  of u0, u1, u2, ... un
//    in order to compute the derivatives:
//      Ux(x0,y0)  =  cx00*u0 +  cx01*u1 +  cx02*u2 + ...  cx0n*un
//      Ux(x1,y1)  =  cx10*u0 +  cx11*u1 +  cx12*u2 + ...  cx1n*un
//      Uy(xk,yk)  =  cyk0*u0 +  cyk1*u1 +  cyk2*u2 + ...  cykn*un
//      Uxx(x3,y3) = cxx30*u0 + cxx31*u1 + cxx32*u2 + ... cxx3n*un
//      ...
//
//           with u0, u1, u2, ...  un variables:
//           (Remember the  values are unknown at this time.
//            The u values are the unknowns in the PDE.
//            We are finding the c coefficients in order to
//            build a system of simultaneous equations, that
//            when solved, provide the solution to the PDE.)
//
//  Method:   fit data   (the coefficients are lined up in four columns)
//
//  U(x,y) = c0           +  c1*x         +  c2*y          +  c3*x^2    +
//           c4*xy        +  c5*y^2       +  c6*x^3        +  c7*x^2y   +
//           c8*xy^2      +  c9*y^3       + c10*x^4        + c11*x^3y   +
//          c12*x^2y^2    + c13*xy^3      + c14*y^4        + c15*x^5    +
//          c16*x^4y      + c17*x^3y^2    + c18*x^2y^3     + c19*xy^4   +
//          c20*y^5       + c21*x^6       + c22*x^5y       + c23*x^4y^2 +
//          c24*x^3y^3    + c25*x^2y^4    + c26*x^y^5      + c27*y^6    +
//          ...
//
//  then, taking the derivative of U(x,y) with respect to x :
//  Ux(x,y) =                c1                          +  2*c3*x       +
//           c4*y                       +  3*c6*x^2      +  2*c7*x*y     +
//           c8*y^2                     + 4*c10*x^3      + 3*c11*x^2*y   +
//        2*c12*x*y^2   +   c13*y^3                      + 5*c15*x^4     +
//        4*c16*x^3*y   + 3*c17*x^2*y^2 + 2*c18*x*y^3    +   c19*y^4     +
//                        6*c21*x^5     + 5*c22*x^4*y    + 4*c23*x^3*y^2 +
//        3*c24*x^2*y^3 + 2*c25*x*y^4   +   c26*y^5      +
//        ...
//  
//  then, taking the derivative of U(x,y) with respect to y:
//  Uy(x,y)  =                               c2                          +
//          c4*x       +  2*c5*y                         +    c7*x^2     +
//        2*c8*x*y     +  3*c9*y^2                       +   c11*x^3     +
//       2*c12*x^2*y   + 3*c13*x*y^2    + 4*c14*y^3                      +
//         c16*x^4     + 2*c17*x^3*y    + 3*c18*x^2*y^2  + 4*c19*x*y^3   +
//       5*c20*y^4                      +   c22*x^5      + 2*c23*x^4*y   +
//       3*c24*x^3*y^2  +4*c25*x^2*y^3  + 5*c26*x*y^4    + 6*c27*y^5     +
//       ...
//
//  ...
//
//  then, taking the forth derivative of U(x,y) with respect to x:
//  Uxxxx(x,y) =  ...
//           ...
//
//  Algorithm:
//                                               u0 u1 u2 u3 u4 u5 ... 
//  form:  | 1  x0  y0  x0^2  x0*y0 y0^2 x0^3 ... 1  0  0  0  0  0 ... | = c0
//         | 1  x1  y1  x1^2  x1*y1 y1^2 x1^3 ... 0  1  0  0  0  0 ... | = c1
//         | 1  x2  y2  x2^2  x2*y2 y2^2 x2^3 ... 0  0  1  0  0  0 ... | = c2
//         | 1  x3  y3  x3^2  x3*y3 y3^2 x3^3 ... 0  0  0  1  0  0 ... | = c3
//         | 1  x4  y4  x4^2  x4*y4 y4^2 x4^3 ... 0  0  0  0  1  0 ... | = c4
//         | 1  x5  y5  x5^2  x5*y5 y5^2 x5^3 ... 0  0  0  0  0  1 ... | = c5
//         | 1  x6  y6  x6^2  x6*y6 y6^2 x5^3 ... 0  0  0  0  0  0 ... | = c6
//         |    ...
//         | 1  x_(npoints-1)
//
//  reduce
//     c0 c1 c2 c3 c4 c5 c6 ...
//   | 1  0  0  0  0  0  0 ... c00  c01  c02  c03  c04  c05  c06 ... | = u0
//   | 0  1  0  0  0  0  0 ... c10  c11  c12  c13  c14  c15  c16 ... | = u1
//   | 0  0  1  0  0  0  0 ... c20  c21  c22  c23  c24  c25  c26 ... | = u2
//   | 0  0  0  1  0  0  0 ... c30  c31  c32  c33  c34  c35  c36 ... | = u3
//   | 0  0  0  0  1  0  0 ... c40  c41  c42  c43  c44  c45  c46 ... | = u4
//   | 0  0  0  0  0  1  0 ... c50  c51  c52  c53  c54  c55  c56 ... | = u5
//   | 0  0  0  0  0  0  1 ... c60  c61  c62  c63  c64  c65  c66 ... | = u6
//     ...
//                            |      this is just the inverse        |
//                            | add points that do not make singular |
// 
// thus, the derivative of U(x,y) with respect to x at point x0,y0 is:
// Ux(x0,y0) =
//       c10*u0       +   c11*u1       +   c12*u2       +   c13*u3       + ...
//     2*c30*u0*x0    + 2*c31*u1*x0    + 2*c32*u2*x0    + 2*c33*u3*x0    + ...
//       c40*u0*y0    +   c41*u1*y0    +   c42*u2*y0    +   c43*u3*y0    + ...
//     3*c60*u0*x0^2  + 3*c61*u1*x0^2  + 3*c62*u2*x0^2  + 3*c63*u3*x0^2  + ...
//     2*c70*u0*x0*y0 + 2*c71*u1*x0*y0 + 2*c72*u2*x0*y0 + 2*c73*u2*x0*y0 + ...
//           ...
//
// rewriting column coefficients as values to be computed:
//       cx00 = c10 + 2*c30*x0 + c40*y0 + 3*c60*x0^2 + 2*c70*u0*x0*y0  + ...
//       cx01 = c11 + 2*c31*x0 + c41*y0 + 3*c61*x0^2 + 2*c71*u0*x0*y0  + ...
//       cx02 = c12 + 2*c32*x0 + c42*y0 + 3*c62*x0^2 + 2*c72*u0*x0*y0  + ...
//       cx03 = c13 + 2*c33*x0 + c43*y0 + 3*c63*x0^2 + 2*c73*u0*x0*y0  + ...
//       ...
//
//
// To compute the first derivative with respect to x at point x0,y0,
// given cx00, cx01, cx02, cx03, ... and knowing u0 at x0,y0, u1 ...
//
// Ux(x0,y0) = cx00*u0 + cx01*u1 + cx02*u2 +cx03*u3 ...
//             hint: the cx00, cx01, cx03, .. cx0n are the values returned 
//                   for point zero
//
//
// For point one, x1,y1 compute coefficients:
//
//       cx10 = c10 + 2*c30*x1 + c40*y1 + 3*c60*x1^2 + 2*c70*u0*x1*y1  + ...
//       cx11 = c11 + 2*c31*x1 + c41*y1 + 3*c61*x1^2 + 2*c71*u0*x1*y1  + ...
//       cx12 = c12 + 2*c32*x1 + c42*y1 + 3*c62*x1^2 + 2*c72*u0*x1*y1  + ...
//       cx13 = c13 + 2*c33*x1 + c43*y1 + 3*c63*x1^2 + 2*c73*u0*x1*y1  + ...
//       ...
//
// Ux(x1,y1) = cx10*u0 + cx11*u1 + cx12*u2 +cx13*u3 + ...
//                hint: the cx10, cx11, cx13, .. cx1n are the values returned 
//                      for point one
//
// 
// the first derivative of U(x,y) with respect to y at point x1,y1 is:
// Uy(x1,y1) = c20*u0        + c21*u1      +   c22*u2      +   c23*u3      ...
//             c40*u0*x1     + c41*u1*x1   +   c42*u2*x1   +   c43*u3*x1   ...
//           2*c50*u0*y1   + 2*c51*u1*y1   + 2*c52*u2*y1   + 2*c53*u3*y1   ...
//             c70*u0*x1^2 +   c71*u1*x1^2 +   c72*u2*x1^2 +   c73*u3*x1^2 ...
//
// rewriting column coefficients as values to be computed:
//
//       cy10 = c20 + c40*x1 + 2*c50*y1 + c70*x1^2 + ...
//       cy11 = c21 + c41*x1 + 2*c51*y1 + c71*x1^2 + ...
//       cy12 = c22 + c42*x1 + 2*c52*y1 + c72*x1^2 + ...
//       cy13 = c23 + c43*x1 + 2*c53*y1 + c73*x1^2 + ...
//
// To compute the first derivative with respect to y at point x1,y1,
// given cy00, cy01, cy02, cy03, ... and knowing u0 at x0,y0, u1 ...
//
// Uy(x1,y1) = cy10*u0 + cy11*u1 + cy12*u2 + cy13*u3 ...
//
//
// ... Uxy, Uyy, Uxxx, Uxxy, Uxyy, Uyyy, Uxxxx, Uxxxy, Uxxyy, Uxyyy, Uyyyy
//     all partial derivatives up to fourth
//
// Now, we really get clever and use xpow and ypow to compute
// cx, cy, cxx, cxy, cyy, ... terms
//
// Then we are very careful to find the set of points that do
// not create a singular matrix, and return indices of the useful points
// Pick the point, as the first point. Add other x,y as long as they
// do not make the matrix signular. Check that enough points are able
// to be used to get the 'order' requested.

public class nuderiv2dg
{
  int debug = 0;
  double AI[][] = new double[37][37]; // P then inverse
  double P[][]  = new double[37][37]; 
  int ordpt[] = {0, 3, 6, 10, 15, 21, 28, 36};
  int n;      // matrix size is n by n 
  int norder; // actual maximum order 
  int maxnpwr = 36;
                     // max of 36 points used now, thus max of 7th order 
  int xpow[] = {0, 1, 0, 2, 1, 0, 3, 2, 1, 0, 4, 3, 2, 1, 0,
                5, 4, 3, 2, 1, 0, 6, 5, 4, 3, 2, 1, 0,
                7, 6, 5, 4, 3, 2, 1, 0};
  int ypow[] = {0, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4,
                0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 6,
                0, 1, 2, 3, 4, 5, 6, 7};
  double xpwr[] = new double[8];
  double ypwr[] = new double[8];
  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  3 npoints   only dx, dy
  // order 2  6 npoints   only dxx, dxy, dyy and lower orders
  // order 3 10 npoints   only dxxx, dxxy, dxyy, dyyy and lower orders
  // order 4 15 npoints   only dxxxx, dxxxy, dxxyy, dxyyy, dyyyy, lower orders
  // order 5 21 npoints   only dxxxx, dxxxy, dxxyy, dxyyy, dyyyy, lower orders
  // order 6 28 npoints   only dxxxx, dxxxy, dxxyy, dxyyy, dyyyy, lower orders
  // order 7 36 npoints   only dxxxx, dxxxy, dxxyy, dxyyy, dyyyy, lower orders
  // point is 0..npoint-1 the point about which the derivative is computed
  nuderiv2dg()
  {
    System.out.println("nuderiv2dg instantiated");
  }

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

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

    m = build(order, npoint, point, x, y, 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];
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
      } 
    }
    return m;
  } // end nu2dx 

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

    m = build(order, npoint, point, x, y, 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
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
      } 
    }
    return m;
  } // end nu2dy 

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

    m = build(order, npoint, point, x, y, 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];
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
      } 
    }
    return m;
  } // end nu2dxx 

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

    m = build(order, npoint, point, x, y, 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
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
      } 
    }
    return m;
  } // end nu2dxy 

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

    m = build(order, npoint, point, x, y, 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
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
      } 
    }
    return m;
  } // end nu2dyy 


  // get all cxx, cxy, cyy, cx, cy with same points
  int nu2d(int kind, int order, int npoint, int point,
           double x[], double y[],
           double cxx[], double cxy[], double cyy[], double cx[], double cy[],
           int ip[])
  {
    double ci;
    int xp, yp, m;

    okind = kind;
    m = build(order, npoint, point, x, y, ip);
    for(int j=0; j<m; j++) // for point[u] 
    {
      cxx[j] = 0.0;
      cxy[j] = 0.0;
      cyy[j] = 0.0;
      cx[j] = 0.0;
      cy[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];
          cxx[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
	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
          cxy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
	if(ypow[i]>=2) // cyy
	{
          ci = (double)ypow[i]*(double)(ypow[i]-1);
          xp = xpow[i];
          yp = ypow[i]-2; // second derivative
          cyy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
        }
	if(xpow[i]>=1) // cx
	{
          ci = (double)xpow[i];
          xp = xpow[i]-1; // first derivative
          yp = ypow[i];
          cx[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
        if(ypow[i]>=1) // cy
	{
          ci = (double)ypow[i];
          xp = xpow[i];
          yp = ypow[i]-1; // first derivative
          cy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
      } 
    }
    return m;
  } // end nu2d 

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

    m = build(order, npoint, point, x, y, 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];
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
      } 
    }
    return m;
  } // end nu2dxxx 

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

    m = build(order, npoint, point, x, y, 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; // third derivative
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
      } 
    }
    return m;
  } // end nu2dxxy 

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

    m = build(order, npoint, point, x, y, 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
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
      } 
    }
    return m;
  } // end nu2dxyy 

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

    m = build(order, npoint, point, x, y, 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
          c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
      } 
    }
    return m;
  } // end nu2dyyy 


  // get all cxxx, cxxy, cxyy, cyyy, cxx, cxy, cyy, cx, cy with same points
  int nu2d3(int kind, int order, int npoint, int point,
            double x[], double y[],
            double cxxx[], double cxxy[], double cxyy[], double cyyy[],
            double cxx[], double cxy[], double cyy[], double cx[], double cy[],
            int ip[])
  {
    double ci;
    int xp, yp, m;

    okind = kind;
    m = build(order, npoint, point, x, y, ip);
    for(int j=0; j<m; j++) // for point[u] 
    {
      cxxx[j] = 0.0;
      cxxy[j] = 0.0;
      cxyy[j] = 0.0;
      cyyy[j] = 0.0;
      cxx[j] = 0.0;
      cxy[j] = 0.0;
      cyy[j] = 0.0;
      cx[j] = 0.0;
      cy[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];
          cxxx[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
        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; // third derivative
          cxxy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
        } 
        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
          cxyy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
        } 
        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
          cyyy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
        } 
        if(xpow[i]>=2) // cxx
	{
          ci = (double)xpow[i]*(double)(xpow[i]-1);
          xp = xpow[i]-2; // second derivative
          yp = ypow[i];
          cxx[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
	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
          cxy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
	if(ypow[i]>=2) // cyy
	{
          ci = (double)ypow[i]*(double)(ypow[i]-1);
          xp = xpow[i];
          yp = ypow[i]-2; // second derivative
          cyy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
        }
	if(xpow[i]>=1) // cx
	{
          ci = (double)xpow[i];
          xp = xpow[i]-1; // first derivative
          yp = ypow[i];
          cx[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
        if(ypow[i]>=1) // cy
	{
          ci = (double)ypow[i];
          xp = xpow[i];
          yp = ypow[i]-1; // first derivative
          cy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
      } 
    }
    return m;
  } // end nu2d3

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

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

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

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

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

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

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

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

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

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

  // get all cxxx, cxxy, cxyy, cyyy, cxx, cxy, cyy, cx, cy with same points
  int nu2d4(int kind, int order, int npoint, int point,
            double x[], double y[],
            double cxxxx[], double cxxxy[], double cxxyy[], double cxyyy[],
            double cyyyy[],
            double cxxx[], double cxxy[], double cxyy[], double cyyy[],
            double cxx[], double cxy[], double cyy[], double cx[], double cy[],
            int ip[])
  {
    double ci;
    int xp, yp, m;

    okind = kind;
    m = build(order, npoint, point, x, y, ip);
    for(int j=0; j<m; j++) // for point[u] 
    {
      cxxxx[j] = 0.0;
      cxxxy[j] = 0.0;
      cxxyy[j] = 0.0;
      cxyyy[j] = 0.0;
      cyyyy[j] = 0.0;
      cxxx[j] = 0.0;
      cxxy[j] = 0.0;
      cxyy[j] = 0.0;
      cyyy[j] = 0.0;
      cxx[j] = 0.0;
      cxy[j] = 0.0;
      cyy[j] = 0.0;
      cx[j] = 0.0;
      cy[j] = 0.0;
      for(int i=0; i<m; i++)
      {
	if(xpow[i]>=4) // cxxxx
	{
          ci = (double)xpow[i]*(double)(xpow[i]-1)*
               (double)(xpow[i]-2)*(double)(xpow[i]-3);
          xp = xpow[i]-4; // fourth derivative
          yp = ypow[i];
          cxxxx[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
	if(xpow[i]>=3 && ypow[i]>=1) // cxxxy
	{
          ci = (double)xpow[i]*(double)(xpow[i]-1)*
               (double)(xpow[i]-2)*(double)(ypow[i]);
          xp = xpow[i]-3; // third derivative
          yp = ypow[i]-1; // first derivative
          cxxxy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
	if(xpow[i]>=2 && ypow[i]>=2) // cxxyy
	{
          ci = (double)xpow[i]*(double)(xpow[i]-1)*
               (double)(ypow[i])*(double)(ypow[i]-1);
          xp = xpow[i]-2; // second derivative
          yp = ypow[i]-2; // second derivative
          cxxyy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
	if(xpow[i]>=1 && ypow[i]>=3) // cxyyy
	{
          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; // third derivative
          cxyyy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
	if(ypow[i]>=4) // cyyyy
	{
          ci = (double)ypow[i]*(double)(ypow[i]-1)*
	       (double)(ypow[i]-2)*(double)(ypow[i]-3);
          xp = xpow[i];
          yp = ypow[i]-4; // fourth derivative
          cyyyy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
	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];
          cxxx[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
        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; // third derivative
          cxxy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
        } 
        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
          cxyy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
        } 
        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
          cyyy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
        } 
        if(xpow[i]>=2) // cxx
	{
          ci = (double)xpow[i]*(double)(xpow[i]-1);
          xp = xpow[i]-2; // second derivative
          yp = ypow[i];
          cxx[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
	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
          cxy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
	if(ypow[i]>=2) // cyy
	{
          ci = (double)ypow[i]*(double)(ypow[i]-1);
          xp = xpow[i];
          yp = ypow[i]-2; // second derivative
          cyy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
        }
	if(xpow[i]>=1) // cx
	{
          ci = (double)xpow[i];
          xp = xpow[i]-1; // first derivative
          yp = ypow[i];
          cx[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
        if(ypow[i]>=1) // cy
	{
          ci = (double)ypow[i];
          xp = xpow[i];
          yp = ypow[i]-1; // first derivative
          cy[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
	}
      } 
    }
    return m;
  } // end nu2d4


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

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

    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];
    }
    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;     
    }
    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;             
      }
    }
    ni = 2; // trying  ni  points for non singular 
    while(ni<=npoint)
    {
      sing = test_build(ni, xa, ya); // 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>=36) 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];
      }
    }
    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]);
    }
    // build xpwr, ypwr at point 
    xpwr[0] = 1.0;
    ypwr[0] = 1.0;
    for(int k=1; k<8; k++)
    {
      xpwr[k] = xpwr[k-1]*xa[0]; // 'point' moved to zero 
      ypwr[k] = ypwr[k-1]*ya[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[]) // 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;
      for(int k=1; k<8; k++)
      {
        xpwr[k] = xpwr[k-1]*xa[i];
        ypwr[k] = ypwr[k-1]*ya[i];
      }
      for(int j=0; j<ni; j++)
      {
        P[i][j] = xpwr[xpow[j]]*ypwr[ypow[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

} // end class nuderiv2dg