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

typedef double Number;

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

int main(){
	Number dfsafe(Number (*f)(Number), Number, Number, Number*);
	Number f(Number);
	Number x, h, error, res;
	int j;
	j=1;
    while(x!=0){
		if(j>10){ cout << "\n My patience is over. Try another function\n";
			break;
		}
		cout << "\n Enter x. To stop, enter x=0, h=anything.\n";
		cin >>x;
		cout << "OK. Now enter stepsize: an increment in x \n" <<
                 "over which function changes substantially\n";
		cin >> h;
		res = dfsafe(&f,x,h,&error);
		cout << "Here is the numerical derivative of sin(x): " << res
			<< "\n The analytical value is:  " << cos(x) << endl;
		j=j+1;
	}
	return 0;
}

Number f( Number x){
	return sin(x);
}


/************************************************************************
Returns the derivative of a function func at a point x by Ridders' method
of polynomial extrapolation. The value h is input as an estimated initial 
stepsize; it need not to be small, but rather should be an increment in x 
over which func changes "substantially". An estimate of the error in the
derivative is returned as err.
C.A. Bertulani        May/01/2000
*************************************************************************/
#define CON 1.4		/* Stepsize is increased by CON at each interaction */
#define CON2 (CON*CON)	
#define BIG 1.0e30
#define NTAB 10		/* Sets the maximum size of tableau */
#define SAFE 2.0	/* Return when eror is SAFE worse than the best so far */

Number dfsafe(Number (*func)(Number), Number x, Number h, Number *err)
{
	int i,j;
	Number errt,fac,hh,**a,ans;

	try {if (h == 0.0){ throw "h must be nonzero in dfsafe.\n";
	}
	}
	catch(char* message){
		cout << message;
	}
	a=matrix(1,NTAB,1,NTAB);
	hh=h;
	a[1][1]=((*func)(x+hh)-(*func)(x-hh))/(2.0*hh);
	*err=BIG;
	for (i=2;i<=NTAB;i++) {
		/* Successive columns in the Neville tableau go to smaller stepsizes and
		higher orders of extrapolation. */
		hh /= CON;
		a[1][i]=((*func)(x+hh)-(*func)(x-hh))/(2.0*hh);		/* Try new, smaller stepsize */
		fac=CON2;
		for (j=2;j<=i;j++) {		/* Compute extrapolations of various orders, requiring */
			a[j][i]=(a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0);  /* no new function evaluations */
			fac=CON2*fac;
			errt=FMAX(fabs(a[j][i]-a[j-1][i]),fabs(a[j][i]-a[j-1][i-1]));
			/* The error strategy is to compare each new extrapolation to one order lower, */
			/* both at the present stepsize and the previous one.  */
			if (errt <= *err) {	/* If error is decreased, save the improved answer  */
				*err=errt;
				ans=a[j][i];
			}
		}
		if (fabs(a[i][i]-a[i-1][i-1]) >= SAFE*(*err)) {
			/* If higher order is worse by a significant factor SAFE then quit early */
			free_matrix(a,1,NTAB,1,NTAB);
			return ans;
		}
	}
	free_matrix(a,1,NTAB,1,NTAB);
	return ans;
}
#undef CON
#undef CON2
#undef BIG
#undef NTAB
#undef SAFE
/*  End of routine dfsafe  */
/********************************************************************************************/