/* nuderiv_test.c  non uniformly spaced derivative coefficients test */

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

double yf4(double x) /* test function 4 */
{
  return 3.0 -0.5* x + 0.2*x*x - 0.1*x*x*x; 
}

double ypf4(double x) /* derivative of test function 4*/
{
  return -0.5 + 0.4*x - 0.3*x*x; 
}

double yppf4(double x) /* second derivative of test function 4 */
{
  return 0.4 - 0.6*x; 
}

double ypppf4(double x) /* third derivative of test function 4 */
{
  return -0.6; 
}

double yf5(double x) /* test function 5 */
{
  return 3.0 -0.5* x + 0.1*x*x - 0.05*x*x*x + 0.002*x*x*x*x; 
}

double ypf5(double x) /* derivative of test function 5*/
{
  return -0.5 + 0.2*x - 0.15*x*x + 0.008*x*x*x; 
}

double yppf5(double x) /* second derivative of test function 5 */
{
  return 0.2 - 0.3*x + 0.024*x*x; 
}

double ypppf5(double x) /* third derivative of test function 5 */
{
  return -0.3 + 0.048*x; 
}

double yppppf5(double x) /* fourth derivative of test function 5 */
{
  return 0.048; 
}

int main(int argc, char *argv[])
{
  int order, npoint, point;
  double x[5] = {1.0, 2.0, 4.0, 7.0, 8.5};
  double y[5];
  double yp[5];
  double ypp[5];
  double yppp[5];
  double ypppp[5];
  double c[5];
  double nuyp; /* any order */
  int i, j;

  printf("nuderiv_test.c running \n");
  npoint = 4;
  printf("\n testing  npoint=%d \n", npoint);
  for(i=0; i<npoint; i++)
  {
    y[i] = yf4(x[i]);
    yp[i] = ypf4(x[i]);
    ypp[i] = yppf4(x[i]);
    yppp[i] = ypppf4(x[i]);
    printf("x[%d]=%f, y=%f, y'=%f, y''=%f, y'''=%f \n",
           i, x[i], y[i], yp[i], ypp[i], yppp[i]);
  }
  for(order=1; order<npoint; order++)
  {
    for(point=0; point<npoint; point++)
    {
      printf("testing order=%d, npoint=%d, point=%d \n", order, npoint, point);
      nuderiv(order, npoint, point, x, c);
      for(i=0; i<npoint; i++) printf("nu c[%d]=%f \n", i, c[i]);
      nuyp = 0.0;
      for(i=0; i<npoint; i++)
      {
        nuyp = nuyp + c[i]*y[i];
      }
      if(order==1)
        printf("order=%d, npoint=%d, point=%d, yp=%f, nuyp=%f, err=%g \n",
               order, npoint, point, yp[point], nuyp, yp[point]-nuyp);
      else if(order==2)
        printf("order=%d, npoint=%d, point=%d, ypp=%f, nuyp=%f, err=%g \n",
               order, npoint, point, ypp[point], nuyp, ypp[point]-nuyp);
      else if(order==3)
        printf("order=%d, npoint=%d, point=%d, yppp=%f, nuyp=%f, err=%g \n",
               order, npoint, point, yppp[point], nuyp, yppp[point]-nuyp);
    }
  }

  npoint = 5;
  printf("\n testing  npoint=%d \n", npoint);
  for(i=0; i<npoint; i++)
  {
    y[i] = yf5(x[i]);
    yp[i] = ypf5(x[i]);
    ypp[i] = yppf5(x[i]);
    yppp[i] = ypppf5(x[i]);
    ypppp[i] = yppppf5(x[i]);
    printf("x[%d]=%f, y=%f, y'=%f, y''=%f, y'''=%f, y''''=%f \n",
           i, x[i], y[i], yp[i], ypp[i], yppp[i], ypppp[i]);
  }
  for(order=1; order<npoint; order++)
  {
    for(point=0; point<npoint; point++)
    {
      printf("testing order=%d, npoint=%d, point=%d \n", order, npoint, point);
      nuderiv(order, npoint, point, x, c);
      for(i=0; i<npoint; i++) printf("nu c[%d]=%f \n", i, c[i]);
      nuyp = 0.0;
      for(i=0; i<npoint; i++)
      {
        nuyp = nuyp + c[i]*y[i];
      }
      if(order==1)
        printf("order=%d, npoint=%d, point=%d, yp=%f, nuyp=%f, err=%g \n",
               order, npoint, point, yp[point], nuyp, yp[point]-nuyp);
      else if(order==2)
        printf("order=%d, npoint=%d, point=%d, ypp=%f, nuyp=%f, err=%g \n",
               order, npoint, point, ypp[point], nuyp, ypp[point]-nuyp);
      else if(order==3)
        printf("order=%d, npoint=%d, point=%d, yppp=%f, nuyp=%f, err=%g \n",
               order, npoint, point, yppp[point], nuyp, yppp[point]-nuyp);
      else if(order==4)
        printf("order=%d, npoint=%d, point=%d, ypppp=%f, nuyp=%f, err=%g \n",
               order, npoint, point, ypppp[point], nuyp, ypppp[point]-nuyp);
    }
  }

  return 0;
} /* end main of nuderiv.c */