/* nuderiv3d.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(x0,y0,z0)  =  cy00*u0 +  cy01*u1 +  cy02*u2 + ...  cy0n*un
 *      Uy(xk,yk,zk)  =  cyk0*u0 +  cyk1*u1 +  cyk2*u2 + ...  cykn*un
 *      Uz(x0,y0,z0)  =  cz00*u0 +  cz01*u1 +  cz02*u2 + ...  cz0n*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      +
 *             c4x^2       +  c5*xy        +  c6*xz        +  c7*y^2    +
 *             c8*yz       +  c9*z^2       + c10*x^3       + c11*x^2y   +
 *            c12*x^2z     + c13*xz^2      + c14*xyz       + c15*xz^2   +
 *            c16*y^3      + c17*y^2z      + c18*yz^2      + c19*z^3    +
 *            c20*x^4      + c21*x^3y      + c22*x^3z      + c23*x^2y^2 +
 *            c24*x^2yz    + c25*x^2z^2    + c26*xy^3      + c27*xy^2z  +
 *            c28*xyz^2    + c29*xz^3      + c30*y^4       + c31*y^3z   +
 *            c32*y^2z^2   + c33*y*z^3     + c34*z^4       +
 *          ...
 *
 *  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*z^2     +   c14*y*z     +   c15*z^2     +
 *          3*c20*x^2    + 3*c21*x^2*y   + 3*c22*x^2*z   + 2*c23*x*y^2   +
 *          2*c24*x*yz   + 2*c25*x*z^2   +   c26*y^3     +   c27*y^2*z   +
 *            c28*y*z^2  +   c29*z^3     +
 *          ...
 *  
 *
 *  Algorithm:
 *                                               u0 u1 u2 u3 u4 u5 ... 
 *  form:  | 1  x0  y0  z0  x0^2  x0*y0 x0*z0 ... 1  0  0  0  0  0 ... | = c0
 *         | 1  x1  y1  z0  x1^2  x1*y1 x0*z0 ... 0  1  0  0  0  0 ... | = c1
 *         | 1  x2  y2  z0  x2^2  x2*y2 x0*z0 ... 0  0  1  0  0  0 ... | = c2
 *         | 1  x3  y3  z0  x3^2  x3*y3 x0*z0 ... 0  0  0  1  0  0 ... | = c3
 *         | 1  x4  y4  z0  x4^2  x4*y4 x0*z0 ... 0  0  0  0  1  0 ... | = c4
 *         | 1  x5  y5  z0  x5^2  x5*y5 x0*z0 ... 0  0  0  0  0  1 ... | = c5
 *         | 1  x6  y6  z0  x6^2  x6*y6 x0*z0 ... 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      |
 * 
 * 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,z0
 * given cx00, cx01, cx02, cx03, ... and knowing u0 at x0,y0, u1 ...
 *
 * Ux(x0,y0,y0) = cx00*u0 + cx01*u1 + cx02*u2 +cx03*u3 ...
 *             hint: the cx00, cx01, cx03, .. cx0n are the values returned 
 *                   for point zero
 *
 * ... 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, ypow and zpow to compute
 * cx, cy, cxx, cxy, cyy, ... terms
 */

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

static int debug = 0;
static double *P = NULL;  /* point matrix */
static double *AI = NULL; /* inverse */
static int ordpt[] = {1, 4, 10, 20, 35};
static int n;      /* matrix size is n by n */
static int norder; /* actual maximum order */
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, 3, 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, 9, 1, 2, 3, 0, 1, 2, 3, 4};
static double xpwr[5], ypwr[5], zpwr[5];

static void inverse(int n, double A[], double AA[]);
static void build(int order, int npoint, int point,
                  double x[], double y[], double z[]);

/* the number of x,y,z  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  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
 */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  build(order, npoint, point, x, y, z);
  for(j=0; j<n; j++) /* for cxxy-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;
      zp = zpow[i];
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
    } 
  }
} /* end nu3dxxy */

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

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

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

  build(order, npoint, point, x, y, z);
  for(j=0; j<n; j++) /* for cxyy-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;
      zp = zpow[i];
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp]*zpwr[zp];
    } 
  }
} /* end nu3dxyy */

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

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

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

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

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

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

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

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

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

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

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

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


static void build(int order, int npoint, int point,
                  double x[], double y[], double z[])
{
  int i, j, k;

  /* check order for n, exact for now  */
  n = npoint;
  norder = order;
  if(P!=NULL) free(P);
  P = (double *)calloc(n*n, sizeof(double));
  if(P==NULL) printf("**** nuderiv3d ERROR P not allocated \n");
  if(AI!=NULL) free(AI);
  AI = (double *)calloc(n*n, sizeof(double));
  if(AI==NULL) printf("**** nuderiv3d ERROR AI not allocated \n");

  for(i=0; i<n; 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]*x[i];
      ypwr[k] = ypwr[k-1]*y[i];
      zpwr[k] = zpwr[k-1]*z[i];
    }
    for(j=0; j<n; j++)
    {
      P[i*n+j] = xpwr[xpow[j]]*ypwr[ypow[j]]*zpwr[zpow[j]];
    }
  }
  if(debug)
  {
    printf("in build n=%d \n", n);
    printf("P00=%g, P01=%g, P02=%g, P03=%g \n",
            P[0], P[1], P[2], P[3]);
    printf("P10=%g, P11=%g, P12=%g, P13=%g \n",
            P[n], P[n+1], P[n+2], P[n+3]);
    printf("P20=%g, P21=%g, P22=%g, P23=%g \n",
            P[2*n], P[2*n+1], P[2*n+2],P[2*n+3]);
    printf("P30=%g, P31=%g, P32=%g, P33=%g \n",
            P[3*n], P[3*n+1], P[3*n+2], P[3*n+3]);
  }
  inverse(n, P, AI); /* invert the matrix */

  if(debug)
  {
    printf("AI00=%g, AI01=%g, AI02=%g, AI03=%g \n",
            AI[0], AI[1], AI[2], AI[3]);
    printf("AI10=%g, AI11=%g, AI12=%g, AI13=%g \n",
            AI[n], AI[n+1], AI[n+2], AI[n+3]);
    printf("AI20=%g, AI21=%g, AI22=%g, AI23=%g \n",
            AI[2*n], AI[2*n+1], AI[2*n+2], AI[2*n+3]);
    printf("AI30=%g, AI31=%g, AI32=%g, AI33=%g \n",
            AI[3*n], AI[3*n+1], AI[3*n+2], AI[3*n+3]);
    printf("leaving build of nuderiv3d n=%d \n\n", n);
  }
  /* 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]*x[point];
    ypwr[k] = ypwr[k-1]*y[point];
    zpwr[k] = zpwr[k-1]*z[point];
  }
} /* end build */

static void 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));
    memcpy(AA,A,n*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 = AA[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(AA[ROW[i]*n+COL[j]]) > ABS_PIVOT ){
            I_PIVOT = i ;
            J_PIVOT = j ;
            PIVOT = AA[ROW[i]*n+COL[j]] ;
          }
        }
      }
      if(abs(PIVOT) < 1.0E-12){
        free(ROW);
        free(COL);
        free(TEMP);
        printf("nuderiv3d MATRIX is SINGULAR !!! \n");
        return;
      }

      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 */
      AA[ROW[k]*n+COL[k]] = 1.0 / PIVOT ;
      for (j=0; j<n; j++) {
        if( j != k ){
          AA[ROW[k]*n+COL[j]] = AA[ROW[k]*n+COL[j]] * AA[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 ){
              AA[ROW[i]*n+COL[j]] = AA[ROW[i]*n+COL[j]] - AA[ROW[i]*n+COL[k]] *
                                   AA[ROW[k]*n+COL[j]] ;
            }
          }
          AA[ROW[i]*n+COL [k]] = - AA[ROW[i]*n+COL[k]] * AA[ROW[k]*n+COL[k]] ;
        }
      }
    }
    /*      end main reduction loop */

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

    /*    unscramble columns */
    for (i=0; i<n; i++) {
      for (j=0; j<n; j++) {
        TEMP[ROW[j]] = AA[i*n+COL[j]] ;
      }
      for (j=0; j<n; j++) {
        AA[i*n+j] = TEMP[j] ;
      }
    }
  /*
  free(ROW);
  free(COL);
  free(TEMP);
  */
} /* end inverse */

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