// nuderiv2d.java  non uniformly spaced derivative coefficients 
//                 see comments in file  nuderiv2d.c
public class nuderiv2d
{
  int debug = 0;
  double AI[][]; // P then inverse 
  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];

// 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 and lower orders
// order 5 21 npoints   only dxxxx, dxxxy, dxxyy, dxyyy, dyyyy and lower orders
// order 6 28 npoints   only dxxxx, dxxxy, dxxyy, dxyyy, dyyyy and lower orders
// order 7 36 npoints   only dxxxx, dxxxy, dxxyy, dxyyy, dyyyy and lower orders
// point is 0..npoint-1 the point about which the derivative is computed
  nuderiv2d()
  {
    System.out.println("nuderiv2d instantiated");
  }

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

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

void nu2dy(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cy-point[u] 
  {
    c[j] = 0.0;
    for(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;
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dy 

void nu2dxx(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cxx-point[u] 
  {
    c[j] = 0.0;
    for(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];
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dxx 

void nu2dxy(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cxy-point[u] 
  {
    c[j] = 0.0;
    for(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;
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dxy 

void nu2dyy(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cyy-point[u] 
  {
    c[j] = 0.0;
    for(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;
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dyy 

void nu2dxxx(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cxxx-point[u] 
  {
    c[j] = 0.0;
    for(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];
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dxxx 

void nu2dxxy(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cxxx-point[u] 
  {
    c[j] = 0.0;
    for(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;
      yp = ypow[i]-1;
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dxxy 

void nu2dxyy(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cxxx-point[u] 
  {
    c[j] = 0.0;
    for(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;
      yp = ypow[i]-2;
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dxyy 

void nu2dyyy(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cyyy-point[u] 
  {
    c[j] = 0.0;
    for(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;
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dyyy 

void nu2dxxxx(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cxxxx-point[u] 
  {
    c[j] = 0.0;
    for(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];
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dxxxx 

void nu2dxxxy(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cxxxy-point[u] 
  {
    c[j] = 0.0;
    for(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;
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dxxxy 

void nu2dxxyy(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cxxyy-point[u] 
  {
    c[j] = 0.0;
    for(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;
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dxxyy 

void nu2dxyyy(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cxyyy-point[u] 
  {
    c[j] = 0.0;
    for(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;
      yp = ypow[i]-3;
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dxyyy 

void nu2dyyyy(int order, int npoint, int point, double x[], double y[], double c[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y);
  for(j=0; j<n; j++) // for cyyyy-point[u] 
  {
    c[j] = 0.0;
    for(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;
      c[j] += ci*AI[i][j]*xpwr[xp]*ypwr[yp];
    } 
  }
} // end nu2dyyyy 


  void build(int order, int npoint, int point, double x[], double y[])
  {
    // check order for n, exact for now  
    n = npoint;
    norder = order;
    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;
      for(int k=1; k<8; k++)
      {
        xpwr[k] = xpwr[k-1]*x[i];
        ypwr[k] = ypwr[k-1]*y[i];
      }
      for(int j=0; j<n; j++)
      {
        AI[i][j] = xpwr[xpow[j]]*ypwr[ypow[j]];
      }
    }
    if(debug>0 && n>3)
    {
      System.out.println("P matrix in build n="+n);
      System.out.println("AI00="+AI[0][0]+", AI01="+AI[0][1]+
                         ", AI02="+AI[0][2]+", AI03="+AI[0][3]);            
      System.out.println("AI10="+AI[1][0]+", AI11="+AI[1][1]+
                         ", AI12="+AI[1][2]+", AI13="+AI[1][3]);
      System.out.println("AI20="+AI[2][0]+", AI21="+AI[2][1]+
                         ", AI22="+AI[2][2]+", AI23="+AI[2][3]);
      System.out.println("AI30="+AI[3][0]+", AI31="+AI[3][1]+
                         ", AI32="+AI[3][2]+", AI33="+AI[3][3]);
    }

    new invert(AI); // invert the matrix 

    if(debug>0 && n>3)
    {
      System.out.println("inverted matrix in build n="+n);
      System.out.println("AI00="+AI[0][0]+", AI01="+AI[0][1]+
                         ", AI02="+AI[0][2]+", AI03="+AI[0][3]);            
      System.out.println("AI10="+AI[1][0]+", AI11="+AI[1][1]+
                         ", AI12="+AI[1][2]+", AI13="+AI[1][3]);
      System.out.println("AI20="+AI[2][0]+", AI21="+AI[2][1]+
                         ", AI22="+AI[2][2]+", AI23="+AI[2][3]);
      System.out.println("AI30="+AI[3][0]+", AI31="+AI[3][1]+
                         ", AI32="+AI[3][2]+", AI33="+AI[3][3]);
      System.out.println("leaving build of nuderiv2d n="+n);
    }
    // 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]*x[point];
      ypwr[k] = ypwr[k-1]*y[point];
    }
  } // end build 
} // end class nuderiv2d