/* test_nuderiv3d.c   test  nuderiv3d  with 35 points        */
/*                    order 1  minimum 4 independent points  */
/*                    order 2  minimum 10 independent points */
/*                    order 3  minimum 20 independent points */
/*                    order 4  minimum 35 independent points */

#include <stdio.h>
#include <stdlib.h>
#include "nuderiv3d.h"

#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#undef  max
#define max(x,y) ((x)>(y)?(x):(y))

static double f1(double x, double y, double z)
{
  return x + 2.0*y + 3.0*z + 4.0;
}

static double fx1(double x, double y, double z)
{
  return 1.0;
}

static double fy1(double x, double y, double z)
{
  return 2.0;
}

static double fz1(double x, double y, double z)
{
  return 3.0;
}


static double f2(double x, double y, double z)
{
  return x*x + 2.0*x*y + 3.0*x*z + 4.0*y*y +
     5.0*y*z + 6.0*z*z + 7.0*x   + 8.0*y   +
     9.0*z   + 10.0;
}

static double fx2(double x, double y, double z)
{
  return 2.0*x + 2.0*y + 3.0*z + 7.0;
}
 
static double fy2(double x, double y, double z)
{
  return 2.0*x + 8.0*y + 5.0*z + 8.0;
}

static double fz2(double x, double y, double z)
{
  return 3.0*x + 5.0*y + 12.0*z + 9.0;
}

static double fxx2(double x, double y, double z)
{
  return 2.0;
}

static double fxy2(double x, double y, double z)
{
  return 2.0;
}

static double fxz2(double x, double y, double z)
{
  return 3.0;
}

static double fyy2(double x, double y, double z)
{
  return 8.0;
}

static double fyz2(double x, double y, double z)
{
  return 5.0;
}

static double fzz2(double x, double y, double z)
{
  return 12.0;
}

 
static double f3(double x, double y, double z)
{
  return x*x*x +  2.0*x*x*y +  3.0*x*x*z + 4.0*x*y*y + 
     5.0*x*y*z +  6.0*x*z*z +  7.0*y*y*y + 8.0*y*y*z + 
     9.0*y*z*z + 10.0*z*z*z + 11.0*x*x   + 12.0*x*y + 
    13.0*x*z   + 14.0*y*y   + 15.0*y*z   + 16.0*z*z  + 
    17.0*x     + 18.0*y     + 19.0*z     + 20.0;
}

static double fx3(double x, double y, double z)
{
  return 3.0*x*x +  4.0*x*y +  6.0*x*z +  4.0*y*y +
         5.0*y*z +  6.0*z*z +
                              22.0*x   + 12.0*y   +
        13.0*z  +
        17.0;
}

static double fy3(double x, double y, double z)
{
  return             2.0*x*x  +             8.0*x*y +
         5.0*x*z +              21.0*y*y + 16.0*y*z +
         9.0*z*z +                         12.0*x   +
                    28.0*y    + 15.0*z   +
                    18.0;
}

static double fz3(double x, double y, double z)
{
  return                       3.0*x*x +
         5.0*x*y + 12.0*x*z +             8.0*y*y + 
        18.0*y*z + 30*z*z   +
        13.0*x   +            15.0*y   + 32.0*z   +
                              19.0;
}

static double fxx3(double x, double y, double z)
{
  return 6.0*x + 4.0*y +  6.0*z + 22.0;
}

static double fxy3(double x, double y, double z)
{
  return 4.0*x + 8.0*y + 5.0*z + 12.0  ;
}

static double fxz3(double x, double y, double z)
{
  return 6.0*x + 5.0*y + 12.0*z + 13.0;
}

static double fyy3(double x, double y, double z)
{
  return 8.0*x + 42.0*y + 16.0*z + 28.0;
}

static double fyz3(double x, double y, double z)
{
  return 5.0*x + 16.0*y + 18.0*z + 15.0;
}

static double fzz3(double x, double y, double z)
{
  return 12.0*x + 18.0*y + 60.0*z + 32.0;
}

static double fxxx3(double x, double y, double z)
{
  return 6.0;
}

static double fxxy3(double x, double y, double z)
{
  return 4.0;
}

static double fxxz3(double x, double y, double z)
{
  return 6.0;
}

static double fxyy3(double x, double y, double z)
{
  return 8.0;
}

static double fxyz3(double x, double y, double z)
{
  return 5.0;
}

static double fxzz3(double x, double y, double z)
{
  return 12.0;
}

static double fyyy3(double x, double y, double z)
{
  return 42.0;
}

static double fyyz3(double x, double y, double z)
{
  return 16.0;
}

static double fyzz3(double x, double y, double z)
{
  return 18.0;
}

static double fzzz3(double x, double y, double z)
{
  return 60.0;
}


int main(int argc, char *argb[])
{
  int i, j;
  int debug = 0;
  int n;
  int order;
  double x[35] = {0.05, 0.14, 0.27, 0.42, 0.57, 0.71, 0.86,
                  0.04, 0.18, 0.31, 0.43, 0.63, 0.82, 0.95,
                  0.07, 0.16, 0.33, 0.41, 0.59, 0.82, 0.99,
                  0.11, 0.21, 0.35, 0.45, 0.65, 0.75, 0.97,
                  0.06, 0.17, 0.32, 0.46, 0.61, 0.79, 0.98};
  double y[35] = {0.06, 0.05, 0.07, 0.10, 0.13, 0.15, 0.17,
                  0.30, 0.33, 0.36, 0.39, 0.31, 0.32, 0.37,
                  0.48, 0.49, 0.52, 0.51, 0.47, 0.50, 0.53,
                  0.65, 0.77, 0.71, 0.67, 0.67, 0.75, 0.76,
                  0.92, 0.97, 0.93, 0.98, 0.87, 0.91, 0.99};
  double z[35] = {0.07, 0.28, 0.47, 0.69, 0.91,
                  0.15, 0.45, 0.68, 0.93, 0.05,
                  0.26, 0.50, 0.97, 0.08, 0.14,
                  0.39, 0.65, 0.76, 0.93, 0.27,
                  0.56, 0.82, 0.99, 0.11, 0.33,
                  0.72, 0.94, 0.09, 0.22, 0.65,
                  0.95, 0.10, 0.24, 0.42, 0.82};
  double U[35];
  double est, actual, dx, dy, dz;
  double cx[35], cy[35], cz[35], cxx[35], cxy[35], cxz[35];
  double cyy[35], cyz[35], czz[35], cxxx[35], cxxy[35];
  double cxxz[35], cxyy[35], cxyz[35], cxzz[35], cyyy[35], cyyz[35];
  double  cyzz[35], czzz[35];
  double error, maxerr;
  int point;

  printf("test_nuderiv3d.c running \n");

  printf("test f1 for dx, dy, dz \n");
  n = 4;
  order = 1;
  for(i=0; i<n; i++)
  {
    U[i] = f1(x[i],y[i],z[i]);
    dx = fx1(x[i],y[i],z[i]);
    dy = fy1(x[i],y[i],z[i]);
    dz = fz1(x[i],y[i],z[i]);
    if(debug) printf("point %3d  x=%5.3f, y=%5.3f, z=%5.3f\n      U=%7.3f, dx=%7.3f, dy=%7.3f, dz=%7.3f \n",
                     i, x[i], y[i], z[i], U[i], dx, dy, dz);
  } 
  maxerr = 0.0;
  for(point=0; point<n; point++) /* fx1 */
  {
    nu3dx(order, n, point, x, y, z, cx);
    if(debug)
      printf("f1, fx1, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cx[i]*U[i];
    actual = fx1(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f1 cx[%3d], est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end fx1 */
  printf("f1, fx1, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fy1 */
  {
    nu3dy(order, n, point, x, y, z, cy);
    if(debug) 
      printf("f1, fy1, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cy[i]*U[i];
    actual = fy1(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f1 cy[%3d], est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end fy1 */
  printf("f1, fy1, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fz1 */
  {
    nu3dz(order, n, point, x, y, z, cz);
    if(debug) 
      printf("f1, fz1, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cz[i]*U[i];
    actual = fz1(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f1 cz[%3d], est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end fz1 */
  printf("f1, fz1, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  printf("test f2 for dx, dy, dz, dxx, dxy. dxz, dyy, dyz, dzz \n");
  n = 10;
  order = 2;

  for(i=0; i<n; i++)
  {
    U[i] = f2(x[i],y[i],z[i]);
  } 
  maxerr = 0.0;
  for(point=0; point<n; point++) /* fx2 */
  {
    nu3dx(order, n, point, x, y, z, cx);
    if(debug)
      printf("f2, fx2, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cx[i]*U[i];
    actual = fx2(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 cx[%3d], est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end fx2 */
  printf("f2, fx2, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fy2 */
  {
    nu3dy(order, n, point, x, y, z, cy);
    if(debug) 
      printf("f2, fy2, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cy[i]*U[i];
    actual = fy2(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 cy[%3d], est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end fy2 */
  printf("f2, fy2, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fz2 */
  {
    nu3dz(order, n, point, x, y, z, cz);
    if(debug) 
      printf("f2, fz2, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cz[i]*U[i];
    actual = fz2(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 cz[%3d], est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end fz2 */
  printf("f2, fz2, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fxx2 */
  {
    nu3dxx(order, n, point, x, y, z, cxx);
    if(debug)
      printf("f2, fxx2, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cxx[i]*U[i];
    actual = fxx2(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 cxx[%3d], est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end fxx2 */
  printf("f2, fxx2, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 


  maxerr = 0.0;
  for(point=0; point<n; point++) /* fxy2 */
  {
    nu3dxy(order, n, point, x, y, z, cxy);
    if(debug)
      printf("f2, fxy2, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cxy[i]*U[i];
    actual = fxy2(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 cxy[%3d], est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end fxy2 */
  printf("f2, fxy2, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fxz2 */
  {
    nu3dxz(order, n, point, x, y, z, cxz);
    if(debug)
      printf("f2, fxz2, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cxz[i]*U[i];
    actual = fxz2(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 cxz[%3d], est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end fxz2 */
  printf("f2, fxz2, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fyy2 */
  {
    nu3dyy(order, n, point, x, y, z, cyy);
    if(debug)
      printf("f2, fyy2, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cyy[i]*U[i];
    actual = fyy2(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 cyy[%3d], est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end fyy2 */
  printf("f2, fyy2, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fyz2 */
  {
    nu3dyz(order, n, point, x, y, z, cyz);
    if(debug)
      printf("f2, fyz2, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cyz[i]*U[i];
    actual = fyz2(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 cyz[%3d], est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end fyz2 */
  printf("f2, fyz2, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fzz2 */
  {
    nu3dzz(order, n, point, x, y, z, czz);
    if(debug)
      printf("f2, fzz2, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += czz[i]*U[i];
    actual = fzz2(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 czz[%3d], est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end fzz2 */
  printf("f2, fzz2, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 


  printf("test f3 for dx, dy, dz, dxx, dxy. dxz, dyy, dyz, dzz \n");
  printf(" dxxx, dxxy, dxxz, dxyy, dxyz. dxzz, dyyy, dyyz, dyzz, dzzz \n");
  n = 20;
  order = 3;

  for(i=0; i<n; i++)
  {
    U[i] = f3(x[i],y[i],z[i]);
  } 
  maxerr = 0.0;
  for(point=0; point<n; point++) /* fx3 */
  {
    nu3dx(order, n, point, x, y, z, cx);
    if(debug)
      printf("f3, fx3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cx[i]*U[i];
    actual = fx3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("cx3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fx3 */
  printf("f3, fx3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fy3 */
  {
    nu3dy(order, n, point, x, y, z, cy);
    if(debug)
      printf("f3, fy3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cy[i]*U[i];
    actual = fy3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 cy3 %3d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end fy3 */
  printf("f3, fy3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fz3 */
  {
    nu3dz(order, n, point, x, y, z, cz);
    if(debug)
      printf("f3, fz3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cz[i]*U[i];
    actual = fz3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 cz3 %3d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end fz3 */
  printf("f3, fz3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fxx3 */
  {
    nu3dxx(order, n, point, x, y, z, cxx);
    if(debug)
      printf("f3, fxx3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cxx[i]*U[i];
    actual = fxx3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("cxx3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fxx3 */
  printf("f3, fxx3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fyy3 */
  {
    nu3dyy(order, n, point, x, y, z, cyy);
    if(debug)
      printf("f3, fyy3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cyy[i]*U[i];
    actual = fyy3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("cyy3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fyy3 */
  printf("f3, fyy3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fxxx3 */
  {
    nu3dxxx(order, n, point, x, y, z, cxxx);
    if(debug)
      printf("f3, fxxx3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cxxx[i]*U[i];
    actual = fxxx3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("cxxx3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fxxx3 */
  printf("f3, fxxx3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fxxy3 */
  {
    nu3dxxy(order, n, point, x, y, z, cxxy);
    if(debug)
      printf("f3, fxxy3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cxxy[i]*U[i];
    actual = fxxy3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("cxxy3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fxxy3 */
  printf("f3, fxxy3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fxxz3 */
  {
    nu3dxxz(order, n, point, x, y, z, cxxz);
    if(debug)
      printf("f3, fxxz3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cxxz[i]*U[i];
    actual = fxxz3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("cxxz3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fxxz3 */
  printf("f3, fxxz3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fxyy3 */
  {
    nu3dxyy(order, n, point, x, y, z, cxyy);
    if(debug)
      printf("f3, fxyy3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cxyy[i]*U[i];
    actual = fxyy3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("cxyy3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fxyy3 */
  printf("f3, fxyy3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fxyz3 */
  {
    nu3dxyz(order, n, point, x, y, z, cxyz);
    if(debug)
      printf("f3, fxyz3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cxyz[i]*U[i];
    actual = fxyz3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("cxyz3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fxyz3 */
  printf("f3, fxyz3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fxzz3 */
  {
    nu3dxzz(order, n, point, x, y, z, cxzz);
    if(debug)
      printf("f3, fxzz3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cxzz[i]*U[i];
    actual = fxzz3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("cxzz3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fxzz3 */
  printf("f3, fxzz3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fyyy3 */
  {
    nu3dyyy(order, n, point, x, y, z, cyyy);
    if(debug)
      printf("f3, fyyy3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cyyy[i]*U[i];
    actual = fyyy3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("cyyy3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fyyy3 */
  printf("f3, fyyy3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fyyz3 */
  {
    nu3dyyz(order, n, point, x, y, z, cyyz);
    if(debug)
      printf("f3, fyyz3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cyyz[i]*U[i];
    actual = fyyz3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("cyyz3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fyyz3 */
  printf("f3, fyyz3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fyzz3 */
  {
    nu3dyzz(order, n, point, x, y, z, cyzz);
    if(debug)
      printf("f3, fyzz3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += cyzz[i]*U[i];
    actual = fyzz3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("cyzz3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fyzz3 */
  printf("f3, fyzz3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 

  maxerr = 0.0;
  for(point=0; point<n; point++) /* fzzz3 */
  {
    nu3dzzz(order, n, point, x, y, z, czzz);
    if(debug)
      printf("f3, fzzz3, order=%1d, n=%3d, point=%3d \n", order, n, point);
    est = 0.0;
    for(i=0; i<n; i++) est += czzz[i]*U[i];
    actual = fzzz3(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("czzz3 %3d, est=%g, actual=%g, error=%g \n",
                      point, est, actual, error);
  } /* end fzzz3 */
  printf("f3, fzzz3, order=%1d, n=%3d, maxerr=%g\n", order, n, maxerr); 


  printf("test_nuderiv3d.c finished \n");
  return 0;
} /* end  test_nuderiv3d.c */