/* test_deriv.c  high order, many points */
#include "deriv.h" /* basic deriv and rderiv */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.141592653589793238462643
#define nx 17

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

double f(double x)
{
  return sin(x);
}
double fp(double x) /* first derivative */
{
  return cos(x);
}
double fpp(double x) /* second derivative */
{
  return -sin(x);
}
double fppp(double x) /* third derivative */
{
  return -cos(x);
}
double fpppp(double x) /* fourth derivative */
{
  return sin(x);
}

int main(int argc, char *argv[])
{
  double cx[nx];
  double cxx[nx];
  double cxxx[nx];
  double cxxxx[nx];
  int i, j, ip, ij, iPI;
  int a[nx];
  int b[1];
  double xg[nx]; /* x values */
  double u[nx];  /* sine wave, one cycle, 9 points */
  double dx, x, xmin, xmax, hx;
  double up, upp, uppp, upppp; /* computed */
  double yp, ypp, yppp, ypppp; /* exact */
  int np, mid;
  double err, avgerr, maxerr;

  printf("test_deriv.c sine wave, to 4th order nx=%d \n", nx);
  printf("np is number of points used to compute derivative\n");
  printf("np must be one larger than derivative and 11 too big for third\n");

  xmin = 0.0;
  for(iPI=1; iPI<=2; iPI++)
  {
    xmax = iPI*PI; /* run pi or 2*pi */
    hx = (xmax-xmin)/(double)(nx-1);
    printf("nx=%d, xmin=%f, xmax=%f, hx=%f \n", nx, xmin, xmax, hx);

    for(i=0; i<nx; i++)
    {
      x = xmin+hx*i;
      xg[i] = x;
      u[i] = f(x);
      printf("i=%d, x=%f, u(x)=%f \n", i, x, u[i]);
    }
    printf("\n");

    for(np=3; np<12; np++)
    {
      avgerr = 0.0;
      maxerr = 0.0;
      mid = np/2;
      for(i=0; i<nx; i++)
      {
        ip = i; /* small */
        ij = 0;
        if(i>nx-1-mid) {ip = np-1-(nx-1-i); ij = nx-np;}
        else if(i>=mid) {ip = mid; ij = i-mid;}
        rderiv(1, np, ip, hx, cx);
        up = 0.0;
        for(j=0; j<np; j++)
        {
          /* printf("i=%d, ip=%d, ij=%d, j=%d, ij+j=%d \n",i,ip,ij,j,ij+j); */
          up += cx[j]*u[ij+j]; /* u[j] = f(xg[j]) */
        }
        yp = fp(xg[i]); /* exact derivative at xg[i] */
        err = abs(yp-up);
        avgerr += err;
        maxerr = max(maxerr,err);
        /* printf("np=%d, pt=%d, up=%f, yp=%f,    error=%e \n",
                   np,i,up,yp,err);*/
      } /* i point */
      printf("first,np=%d, maxerr=%e, avgerr=%e \n", np, maxerr, avgerr/nx);
    } /* np */
    printf("\n\n");

    for(np=4; np<12; np++)
    {
      avgerr = 0.0;
      maxerr = 0.0;
      mid = np/2;
      for(i=0; i<nx; i++)
      {
        ip = i; /* small */
        ij = 0;
        if(i>nx-1-mid) {ip = np-1-(nx-1-i); ij = nx-np;}
        else if(i>=mid) {ip = mid; ij = i-mid;}
        rderiv(2, np, ip, hx, cxx);
        upp = 0.0;
        for(j=0; j<np; j++)
        {
          /* printf("i=%d, ip=%d, ij=%d, j=%d, ij+j=%d \n",i,ip,ij,j,ij+j); */
          upp += cxx[j]*u[ij+j]; /* u[j] = f(xg[j]) */
        }
        ypp = fpp(xg[i]); /* exact derivative at xg[i] */
        err = abs(ypp-upp);
        avgerr += err;
        maxerr = max(maxerr,err);
        /* printf("np=%d, pt=%d, upp=%f, ypp=%f,    error=%e \n",
	           np,i,upp,ypp,err); */
      } /* i point */
      printf("second,np=%d, maxerr=%e, avgerr=%e \n", np, maxerr, avgerr/nx);
    } /* np */
    printf("\n\n");


    for(np=5; np<12; np++)
    {
      avgerr = 0.0;
      maxerr = 0.0;
      mid = np/2;
      for(i=0; i<nx; i++)
      {
        ip = i; /* small */
        ij = 0;
        if(i>nx-1-mid) {ip = np-1-(nx-1-i); ij = nx-np;}
        else if(i>=mid) {ip = mid; ij = i-mid;}
        rderiv(3, np, ip, hx, cxxx);
        uppp = 0.0;
        for(j=0; j<np; j++)
        {
          /* printf("i=%d, ip=%d, ij=%d, j=%d, ij+j=%d \n",i,ip,ij,j,ij+j); */
          uppp += cxxx[j]*u[ij+j]; /* u[j] = f(xg[j]) */
        }
        yppp = fppp(xg[i]); /* exact derivative at xg[i] */
        err = abs(yppp-uppp);
        avgerr += err;
        maxerr = max(maxerr,err);
        /* printf("np=%d, pt=%d, uppp=%f, yppp=%f,    error=%e \n",
                   np,i,uppp,yppp,err);*/
      } /* i point */
      printf("third,np=%d, maxerr=%e, avgerr=%e \n", np, maxerr, avgerr/nx);
    } /* np */
    printf("\n\n");


    for(np=5; np<11; np++) /* 11 too big */
    {
      avgerr = 0.0;
      maxerr = 0.0;
      mid = np/2;
      for(i=0; i<nx; i++)
      {
        ip = i; /* small */
        ij = 0;
        if(i>nx-1-mid) {ip = np-1-(nx-1-i); ij = nx-np;}
        else if(i>=mid) {ip = mid; ij = i-mid;}
        rderiv(4, np, ip, hx, cxxxx);
        upppp = 0.0;
        for(j=0; j<np; j++)
        {
          /* printf("i=%d, ip=%d, ij=%d, j=%d, ij+j=%d \n",i,ip,ij,j,ij+j); */
          upppp += cxxxx[j]*u[ij+j]; /* u[j] = f(xg[j]) */
        }
        ypppp = fpppp(xg[i]); /* exact derivative at xg[i] */
        err = abs(ypppp-upppp);
        avgerr += err;
        maxerr = max(maxerr,err);
        /* printf("np=%d, pt=%d, upppp=%f, ypppp=%f,    error=%e \n",
                   np,i,upppp,ypppp,err); */
      } /* i point */
      printf("fourth,np=%d, maxerr=%e, avgerr=%e \n", np, maxerr, avgerr/nx);
    } /* np */
    printf("\n\n");
  } /* iPI */
  printf("test_deriv.c finished \n");
  return 0;
} /* end main of test_deriv.c */