/************************************************************************
Returns the derivative of a tabulated function.
C.A. Bertulani        May/01/2000
*************************************************************************/

typedef double Number;

#include <math.h>
#include "butil.h"

#define n 20


int main(){
	void deriv(Number f[], Number [], Number );
	Number *f, *df;
	Number *x, h;
	int i;

	cout << "Here are the numerical derivative of sin(x) at your points: \n" ;
	f = vector(n);
	x = vector(n);
	df = vector(n);
	h = 0.1;
	x[1]=0.;
	f[1]=0.;
	for(i=2;i<=n;++i){
		x[i] = x[i-1] + h;
		f[i] = sin(x[i]);
	}
	deriv(f,df,h);
	for(i=1;i<=n;++i){
		cout << "x= " << x[i] << "\tdf/dx= " << df[i] << 
		"\t\tcos(x)= " << cos(x[i]) << endl;
	}
	return 0;
}


/************************************************************************
Returns the derivative of a tabulated function func (array func), of
equally spaced points with stepsize h using a five point differentation 
formula (Abramowitz, pag 914).
C.A. Bertulani        May/01/2000
*************************************************************************/

void deriv(Number func[], Number dfunc[], Number h)
{
	Number **a, sum, emfact;
	int i, j, jj, nmx,  k;
    a = matrix(1,5,1,5);
	Number b[5][5] = {{-50.,96.,-72.,32.,-6.}, {-6.,-20.,36.,-12.,2.},
	{2.,-16.,0.,16.,-2.},{-2.,12.,-36.,20.,6.}, {6.,-32.,72.,-96.,50.}};
	for(i=1;i<=5;i++){ 
		for(j=1;j<=5;j++){ 
			a[i][j]=b[i-1][j-1];
		}
	}
	emfact = 24.;
	nmx = n - 2;
	for(j=1;j<=n;j++){
		k = 3; 
		if(j<3) k = j;
		if(j>nmx) k = j-n+5;
		sum = 0.;
		for(i=1;i<=5;i++){
			jj = j + i - k;
			sum = sum + a[k][i]*func[jj];
		}
		dfunc[j] = sum/(h*emfact);
	}
}
#undef n

/*  End of routine deriv  */
/********************************************************************************************/