/* nuderiv3dg.c  non uniformly spaced derivative coefficients */
/*    given  x0,y0,z0 x1,y1,z1 x2,y2,z2 ... xn,yn,zn  non uniformly spaced
 *    find the coefficients  c???  of u0, u1, u2, ... un
 *    in order to compute the derivatives:
 *      Ux(x0,y0,z0)  =  cx00*u0 +  cx01*u1 +  cx02*u2 + ...  cx0n*un
 *      Ux(x1,y1,z1)  =  cx10*u0 +  cx11*u1 +  cx12*u2 + ...  cx1n*un
 *      Uy(xk,yk,zk)  =  cyk0*u0 +  cyk1*u1 +  cyk2*u2 + ...  cykn*un
 *      Uxx(x3,y3,z3) = 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,z) = c0         +  c1*x         +  c2*y          +  c3*z      +
 *           c4*x^2       +  c5*x*y       +  c6*x*z        +  c7*y^2    +
 *           c8*y*z       +  c9*z^2       + c10*x^3        + c11*x^2*y  +
 *          c12*x^2*z     + c13*x*y^2     + c14*x*y*z      + c15*x*z^2  +
 *          c16*y^3       + c17*y^2*z     + c18*y*z^2      + c19*z^z
 *          c20*x^3       + ...
 *
 *  then, taking the derivative of U(x,y,z) with respect to x :
 *  Ux(x,y,z) =              c1                                          +
 *         2*c4*x       +    c5*y       +    c6*z                        +
 *                                      + 3*c10*x^2      + 2*c11*x  *y   +
 *        2*c12*x*z     +   c13*y^2     +   c14*y*z      +   c15*z^2     +
 *        
 *        3*c20*x^2     + ...
 *  
 *  then, taking the derivative of U(x,y,z) with respect to y:
 *  Uy(x,y,z)  =                             c2                          +
 *                     +    c5*x                         +  2*c7*y       +
 *          c8*z       +                                 +   c11*x^2     +
 *                     + 2*c13*x*y      +   c14*x*z                      +
 *       3*c16*y^2     + 2*c17*y*z      +   c18*z^2                      +
 *       ...
 *
 *  ...
 *
 *  then, taking the forth derivative of U(x,y,z) with respect to z:
 *  Uzzzz(x,y,z) =  ...
 *           ...
 *
 *  Algorithm:
 *  from:                                     u0 u1 u2 u3 u4 u5 ... 
 *   | 1  x0  y0  z0 x0^2  x0*y0 y0^2 x0^3 ... 1  0  0  0  0  0 ... | = c0
 *   | 1  x1  y1  z1 x1^2  x1*y1 y1^2 x1^3 ... 0  1  0  0  0  0 ... | = c1
 *   | 1  x2  y2  z2 x2^2  x2*y2 y2^2 x2^3 ... 0  0  1  0  0  0 ... | = c2
 *   | 1  x3  y3  z3 x3^2  x3*y3 y3^2 x3^3 ... 0  0  0  1  0  0 ... | = c3
 *   | 1  x4  y4  z4 x4^2  x4*y4 y4^2 x4^3 ... 0  0  0  0  1  0 ... | = c4
 *   | 1  x5  y5  z5 x5^2  x5*y5 y5^2 x5^3 ... 0  0  0  0  0  1 ... | = c5
 *   | 1  x6  y6  z6 x6^2  x6*y6 y6^2 x5^3 ... 0  0  0  0  0  0 ... | = c6
 *   |    ...
 *   | 1  x_(npoints-1)
 *
 *  reduce (invert matrix)
 *     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,z) with respect to x at point x0,y0,z0 is:
 * Ux(x0,y0,z0) =
 *       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,z0) = 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,z1 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,z1) = 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,z) with respect to y at point x1,y1,z1 is:
 * Uy(x1,y1,z1) = 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,z1) = cy10*u0 + cy11*u1 + cy12*u2 + cy13*u3 ...
 *
 *
 * ... Uxx, Uxy, Uxz, Uyy, Uyz, Uzz, Uxxx, Uxxy, Uxxz, Uxyy, Uxyz, Uxzz, 
 *     Uyyy, Uyyz, Uyzz, Uzzz 
 *
 * Now, we really get clever and use xpow and ypow and zpow to compute
 * cx, cy, cz, cxx, cxy, cxz, cyy, ... , czzzz 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
 * Always make 'point' be the first, 0,  when solving PDE's
 * "order" is power of fit, not derivative order, as in nd3dxxyz
 * the "3" is three dimensiona, independent variables
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#undef  abs
#define abs(x) (((x)<0.0)?(-(x)):(x))

static int debug = 0;
/* static int stat = 0; not global, returned */
static double *P;  /* point matrix */
static double *AI; /* inverse */
static int ordpt[] = {1, 4, 10, 20, 35}; /* minimum points */
static int maxnpwr = 35;
                     /* max of 35 points used now, thus max of 4th order */
static 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};
static 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};
static 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};
static double xpwr[5], ypwr[5], zpwr[5];

static int build(int order, int npoint, int point,
                 double x[], double y[], double z[], int ip[]);

static int test_build(int ni, double xa[], double ya[], double za[]);

static int inverse(int n, double A[], double AA[]);

/* the number of x,y,z  npoints  must be at least:
 * here, order means maximum order of terms used in computing cx, cy, cz, ...
 * order 0  1 npoints   only constant
 * order 1  4 npoints   only dx, dy, dz
 * order 2 10 npoints   only dxx, dxy, dxz, dyy, dyz, dzz  and lower orders
 * order 3 20 npoints   only dxxx, dxxy, dxxz, dxyy, dxyz, dxzz,
 *                           dyyy, dyyz, dyzz, dzzz  and lower orders
 * order 4 35 npoints   only dxxxx, dxxxy, dxxxz, dxxyy, dxxyz,
 *                           dxyyy, dxyyz, dxyzz, dxzzz, dyyyy,
 *                           dyyyz, dyyzz, dyzzz, dzzzz  and lower orders
 * point is 0..npoint-1 the point about which the derivative is computed
 */

/* m is actual number of entries in c[] and ip[] */
/* xpow[], ypow[] 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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cx-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      if(xpow[i]<1) continue;
      ci = (double)xpow[i];
      xp = xpow[i]-1; /* taking first derivative */
      yp = ypow[i];
      zp = zpow[i];
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      if(ypow[i]<1) continue;
      ci = (double)ypow[i];
      xp = xpow[i];
      yp = ypow[i]-1; /* taking first derivative */
      zp = zpow[i];
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cz-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      if(zpow[i]<1) continue;
      ci = (double)zpow[i];
      xp = xpow[i];
      yp = ypow[i];
      zp = zpow[i]-1; /* taking first derivative */
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cxx-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      if(xpow[i]<2) continue;
      ci = (double)xpow[i]*(double)(xpow[i]-1);
      xp = xpow[i]-2; /* taking second derivative */
      yp = ypow[i];
      zp = zpow[i];
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cxy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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];
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cxz-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      if(xpow[i]<1 || zpow[i]<1) continue;
      ci = (double)xpow[i]*(double)zpow[i];
      xp = xpow[i]-1;
      yp = ypow[i];
      zp = zpow[i]-1;
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cyy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      if(ypow[i]<2) continue;
      ci = (double)ypow[i]*(double)(ypow[i]-1);
      xp = xpow[i];
      yp = ypow[i]-2;
      zp = zpow[i];
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cyz-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      if(ypow[i]<1 || zpow[i]<1) continue;
      ci = (double)ypow[i]*(double)zpow[i];
      xp = xpow[i];
      yp = ypow[i]-1;
      zp = zpow[i]-1;
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for czz-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      if(zpow[i]<2) continue;
      ci = (double)zpow[i]*(double)(zpow[i]-1);
      xp = xpow[i];
      yp = ypow[i];
      zp = zpow[i]-2;
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cxxx-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      if(xpow[i]<3) continue;
      ci = (double)xpow[i]*(double)(xpow[i]-1)*(double)(xpow[i]-2);
      xp = xpow[i]-3; /* taking third derivative */
      yp = ypow[i];
      zp = zpow[i];
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cxxy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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;
      zp = zpow[i];
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cxxz-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      if(xpow[i]<2 || zpow[i]<1) continue;
      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*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cxyy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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;
      zp = zpow[i];
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

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

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

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

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

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cyyy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      if(ypow[i]<3) continue;
      ci = (double)ypow[i]*(double)(ypow[i]-1)*(double)(ypow[i]-2);
      xp = xpow[i];
      yp = ypow[i]-3;
      zp = zpow[i];
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cyyz-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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;
      zp = zpow[i]-1;
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cyzz-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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;
      zp = zpow[i]-2;
      c[j] += ci*AI[i*m+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[])
{
  int i, j, m;
  double ci;
  int xp, yp, zp;

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


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

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cxxxx-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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; /* taking fourth derivative */
      yp = ypow[i];
      zp = zpow[i];
      c[j] += ci*AI[i*m+j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
    } 
  }
  return m;
} /* end nu3dxxxx */

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

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cxxxy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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; /* taking third derivative */
      yp = ypow[i]-1; /* taking first derivative */
      zp = zpow[i];
      c[j] += ci*AI[i*m+j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
    } 
  }
  return m;
} /* end nu3dxxxy */

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

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cxxxz-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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; /* taking third derivative */
      yp = ypow[i];
      zp = zpow[i]-1; /* taking first derivative */
      c[j] += ci*AI[i*m+j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
    } 
  }
  return m;
} /* end nu3dxxxz */

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

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

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

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

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

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

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

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cxyyy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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;
      zp = zpow[i];
      c[j] += ci*AI[i*m+j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
    } 
  }
  return m;
} /* end nu3dxyyy */

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

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

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

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

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

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cxzzz-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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;
      yp = ypow[i];
      zp = zpow[i]-3;
      c[j] += ci*AI[i*m+j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
    } 
  }
  return m;
} /* end nu3dxzzz */

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

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cyyyy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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; // taking 4th derivative
      zp = zpow[i];
      c[j] += ci*AI[i*m+j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
    } 
  }
  return m;
} /* end nu3dyyyy */

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

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


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

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

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

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for cyzzz-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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;
      zp = zpow[i]-3;
      c[j] += ci*AI[i*m+j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
    } 
  }
  return m;
} /* end nu3dyzzz */

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

  m = build(order, npoint, point, x, y, z, ip);
  for(j=0; j<m; j++) /* for czzzz-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<m; i++)
    {
      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; // taking 4th derivative
      c[j] += ci*AI[i*m+j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
    } 
  }
  return m;
} /* end nu3dzzzz */


static int build(int order, int npoint, int point,
		 double x[], double y[], double z[], int ip[])
{
  int i, j, k, n, sing, stat;
  int n2i, nd, ni, nj;
  double xa[1026], ya[1026], za[1026], tmp;

  /* checks  */
  if(point<0 || point>=npoint)
     printf("ERROR nuderiv3dg point outside range\n"); 
  n = npoint;
  // if(n>100) printf("npoint=%d, probably unstable \n", n);

  P = (double *)calloc(n*n, sizeof(double));
  if(P==NULL) printf("**** nuderiv3dg ERROR P not allocated \n");
  AI = (double *)calloc(n*n, sizeof(double));
  if(AI==NULL) printf("**** nuderiv3dg ERROR AI not allocated \n");

  n2i = 1; /* good points so far, this will become  n  */
  nd  = 0; /* number deleted */
  for(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;     
  }
  ni = 2; /* trying  ni  points for non singular */
  while(ni<=npoint)
  {
    if(debug) printf("test build with ni=%d points \n", ni);
    sing = test_build(ni, xa, ya, za); /* P and AI global */
    if(sing==0)
    {
      n2i = ni;
      if(debug) printf("test_build added point %d for index %d, nd=%d  \n",
	               ip[ni-1], ni-1, nd);
      if(ni>=99) break;
      ni++;
    }
    else /* sing != 0 */
    {
      nd++;
      if(debug)
        printf("test_build deleted point %d for index %d, nd=%d, ni+nd=%d \n",
	       ip[ni-1], ni-1, nd, ni+nd+1);
      if(ni+nd>=npoint) break;
      for(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;

  if(debug)
  {
    printf("\n\nactual points selected\n");
    for(i=0; i<n; i++) 
      printf("i=%2d, ip[i]=%2d, xa[i]=%f,  ya[i]=%f, za[i]=%f \n",
             i, ip[i], xa[i], ya[i], za[i]);
  }
  /* build xpwr, ypwr at point */
  xpwr[0] = 1.0;
  ypwr[0] = 1.0;
  zpwr[0] = 1.0;
  for(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])
    printf("probable loss of accuracy, order=%d, needs n=%d, has only %d \n",
           order, ordpt[order], n);
  return n;
} /* end build */


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

  for(i=0; i<ni; i++) /* build matrix that will be inverted */
  {
    xpwr[0] = 1.0;
    ypwr[0] = 1.0;
    zpwr[0] = 1.0;
    for(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(j=0; j<ni; j++)
    {
      P[i*ni+j] = xpwr[xpow[j]]*ypwr[ypow[j]]*zpwr[zpow[j]];
    }
  }

  sing = inverse(ni, P, AI); /* invert the matrix */

  return sing;
} /* end test_build */


static int inverse(int n, double A[], double AA[])
{
  int *ROW, *COL;
  double *TEMP;
  int HOLD , I_pivot , J_pivot;
  double pivot, abs_pivot;
  int i, j, k;

  ROW = (int *)calloc(n, sizeof(int));
  COL = (int *)calloc(n, sizeof(int));
  TEMP = (double *)calloc(n, sizeof(double));

  /*              set up row and column interchange vectors */
  for (k=0; k<n; k++)
  {
    ROW[k] = k ;
    COL[k] = k ;
  }

  /*               begin main reduction loop */
  for (k=0; k<n; k++)
  {
    /*             find largest element for pivot */
    pivot = A[ROW[k]*n+COL[k]] ;
    I_pivot = k;
    J_pivot = k;
    for(i=k; i<n; i++)
    {
      for(j=k; j<n; j++)
      {
        abs_pivot = abs(pivot) ;
        if( abs(A[ROW[i]*n+COL[j]]) > abs_pivot )
        {
          I_pivot = i ;
          J_pivot = j ;
          pivot = A[ROW[i]*n+COL[j]] ;
        }
      }
    }
    if(abs(pivot) < 1.0E-12)
    {
      free(ROW);
      free(COL);
      free(TEMP);
      return 1; /* bad */
    }

    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]*n+COL[k]] = 1.0 / pivot ;
    for (j=0; j<n; j++)
    {
      if( j != k )
      {
        A[ROW[k]*n+COL[j]] = A[ROW[k]*n+COL[j]] * A[ROW[k]*n+COL[k]];
      }
    }

    /*               inner reduction loop */
    for (i=0; i<n; i++)
    {
      if( k != i )
      {
        for (j=0; j<n; j++)
        {
          if( k != j )
          {
            A[ROW[i]*n+COL[j]] = A[ROW[i]*n+COL[j]] - A[ROW[i]*n+COL[k]] *
                                 A[ROW[k]*n+COL[j]] ;
          }
        }
        A[ROW[i]*n+COL [k]] = - A[ROW[i]*n+COL[k]] * A[ROW[k]*n+COL[k]] ;
      }
    }
  }
  /*      end of main reduction loop */

  /*      unscramble rows */
  for (j=0; j<n; j++) {
    for (i=0; i<n; i++) {
      TEMP[COL[i]] = A[ROW[i]*n+j];
    }
    for (i=0; i<n; i++) {
      A[i*n+j] = TEMP[i] ;
    }
  }

  /*      unscramble columns */
  for (i=0; i<n; i++)
  {
    for (j=0; j<n; j++)
    {
      TEMP[ROW[j]] = A[i*n+COL[j]] ;
    }
    for (j=0; j<n; j++)
    {
      AA[i*n+j] = TEMP[j]; /* final output */
    }
  }
  free(ROW);
  free(COL);
  free(TEMP);
  return 0; /* sing==0 is good */
} /* end inverse */

void self_test()
{
  int i, j;
  int error = 0;

  for(i=0; i<maxnpwr-1; i++)
  {
    for(j=i+1; j<maxnpwr; j++)
    {
      if(xpow[i]==xpow[j] && ypow[i]==ypow[j] && zpow[i]==zpow[j])
      {
	error++;
	printf("pow[%d]==pow[%d]\n", i, j);
      }
    }
    if(i==0)
    {
      if(xpow[i]+ypow[i]+zpow[i]!=0)
      {
        error++;
        printf("%d column not zero\n",i);
      }
    }
    else if(i<4)
    {
      if(xpow[i]+ypow[i]+zpow[i]!=1)
      {
        error++;
        printf("%d column not 1\n",i);
      }
    }
    else if(i<10)
    {
      if(xpow[i]+ypow[i]+zpow[i]!=2)
      {
        error++;
        printf("%d column not 2\n",i);
      }
    }
    else if(i<20)
    {
      if(xpow[i]+ypow[i]+zpow[i]!=3)
      {
	error++;
        printf("%d column not 3\n",i);
      }
    }
    else
    {
      if(xpow[i]+ypow[i]+zpow[i]!=4)
      {
        error++;
        printf("%d column not 4\n",i);
      }
    }
  }
  printf("nuderiv3dg.c self test completed with %d errors\n",error);
} // end self_test

/* end nuderiv3dg.c  with a few callable functions */