#include <stdio.h>
#include <string.h>
#include <math.h>

#include "array_utility.h"

/**********************************************************************

Given n data points (not equally spaced) of a time function data(time), 
nderiv finds the dth order numerical derivative of the function at the 
kth data point.

***********************************************************************/
float nderiv_old(y,x,k,d,n)
     float *y, *x;
     int k,d,n;
{
  float u,v,w,dy;
  int i,j,d2=d/2;

  if (k<d/2 || k>=n-d/2) 
    { printf("index %d outside range (%d -- %d)...",k,d/2,n-d/2); pause(); }

  u=1.0;
  for (i=1; i<=d2; i++) 
    if (i) u*=(x[k]-x[k-i])*(x[k]-x[k+i]);
  dy=0;
  for (i=-d2; i<=d2; i++) {
    v=1.0;
    for (j=-d2; j<=d2; j++) 
      if (i!=j) v=v*(x[k+i]-x[k+j]);
    if (i) dy+=y[k+i]*u/(x[k]-x[k+i])/v;
    else 
      for (j=-d2; j<=d2; j++) 
	if (j) dy+=y[k]/(x[k]-x[k+j]);
  }
  return dy;
}

float nderiv(y,x,k,d,n) 
     float *x, *y;
     int k,d,n;
{ float u,v,dy;
  int i,j,d2=d/2;

  u=1.0;
  for (i=-d2; i<=d2; i++)        // common numerator
    if (k+i>=0 && k+i<n && i) 
      u*=(x[k]-x[k+i]);
  dy=0;
  for (i=-d2; i<=d2; i++)        // for all d terms
    if (k+i>=0 && k+i<n) 
      {
	v=1.0;
	for (j=-d2; j<=d2; j++)  // common denominator
	  if (k+j>=0 && k+j<n && i!=j) 
	    v=v*(x[k+i]-x[k+j]);
	if (i) dy+=y[k+i]*u/(x[k]-x[k+i])/v;  // term in the center
	else 
	  for (j=-d2; j<=d2; j++)             // all other terms
	    if (k+j>=0 && k+j<n && j)
	      dy+=y[k]/(x[k]-x[k+j]);
      }
  return dy;
}


main()
{
  float u,v,w,s,pi=3.1416;
  float *x,*y, dy0, dy1, dy2, dy3, dy4, dy5, dy6, dy7, dy8;
  float e1,e2,e3,e4,e5,e6,e7,e8;
  int i,j,k,l,m,n;
  int K,N;

  N=30;
  x=alloc1df(N); y=alloc1df(N);

  for (n=0; n<N; n++)
    { x[n]=n;  
      if (!n) x[0]=0;
      x[n]=x[n-1]+1+(rand()%10-5)*0.1;
      y[n]=pow(x[n],2.3)+pow(x[n],1.2)+10*sin(x[n]+2)-15*cos(x[n]-3);
      //      printf("(%.2f %.2f) ",x[n],y[n]);
    }

  e1=e2=e3=e4=e5=e6=e7=e8=0.0;
  //  for (n=5; n<N-5; n++)
  for (n=0; n<N; n++) {
    dy0=2.3*pow(x[n],1.3)+1.2*pow(x[n],0.2)+10*cos(x[n]+2)+15*sin(x[n]-3); // true derivative 

    dy1=(y[n+1]-y[n])/(x[n+1]-x[n]);  // forward  difference
    dy2=(y[n]-y[n-1])/(x[n]-x[n-1]);  // backward difference
    dy3=(y[n+1]-y[n-1])/(x[n+1]-x[n-1]);  // average of the two

    u=x[n]-x[n-1]; v=x[n+1]-x[n]; w=x[n+1]-x[n-1]; 
    dy4=-y[n-1]*v/u/w-y[n]*(u-v)/u/v+y[n+1]*u/v/w;

    dy4=nderiv(y,x,n,3,N);   // 3rd order estimate (same as dy3)
    dy5=nderiv(y,x,n,5,N);   // 5rd order estimate
    dy6=nderiv(y,x,n,7,N);   // 7rd order estimate
    dy7=nderiv(y,x,n,9,N);   // 9rd order estimate
    dy8=nderiv(y,x,n,11,N);   // 11rd order estimate
	
    e1+=fabs((dy0-dy1)/dy0);
    e2+=fabs((dy0-dy2)/dy0);
    e3+=fabs((dy0-dy3)/dy0);
    e4+=fabs((dy0-dy4)/dy0);
    e5+=fabs((dy0-dy5)/dy0);
    e6+=fabs((dy0-dy6)/dy0);
    e7+=fabs((dy0-dy7)/dy0);
    e8+=fabs((dy0-dy8)/dy0);

    printf("(%6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f)\n",
	   (dy0-dy1)/dy0,(dy0-dy2)/dy0,(dy0-dy3)/dy0,(dy0-dy4)/dy0,
	   (dy0-dy5)/dy0,(dy0-dy6)/dy0,(dy0-dy7)/dy0,(dy0-dy8)/dy0);
  }
  printf("Overall error:\n(%f %f %f %f %f %f %f %f)\n",
	 e1,e2,e3,e4,e5,e6,e7,e8);
}