/***********************************************/
/* T.Kouya's GSL sample program collection     */
/*                   Numerical Differentiation */
/*                   Written by Tomonori Kouya */
/*                                             */
/* Version 0.01: 2007-08-13                    */
/***********************************************/
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>

/* Return relative and absolute errors */
double rel_abs_error(double approx_val, double true_val, double *abserr)
{
	static double abs_error, rel_error;

	abs_error = fabs(approx_val - true_val);
	rel_error = abs_error;
	if(true_val >= 1.0e-15)
		rel_error /= fabs(true_val);

	if(abserr != NULL)
		*abserr = abs_error;

	return rel_error;
}

/* Original function: f(x) = sin(x) */
double func(double x, void *param)
{
	return sin(x);
}

/* Main function */
int main(void)
{
	int i;
	double x, h;
	double results[3], rel_errors[3], abs_errors[3], gsl_abs_errors[3];

	/* Definition of gsl_function:                         */
	/*   struct gsl_function_struct {                      */
	/*      double (* function) (double x, void * params); */
	/*      void * params;                                 */
	/*   };                                                */
	gsl_function f;

	/* Set the original function */
	f.function = &func;
	f.params = NULL;

	/* Get values of f'(x) by Central diff, absolute and relative errors in [0, 2*PI] */
	printf("    x           Central Diff               True f'(x)       gsl_abs Abs_err Rel_err\n");

	/* Initial guess for seeking the optimal stepsize */
	h = 1;

	for(i = 0; i <= 20; i++)
	{
		x = 0.0 + (double)i * (2.0 * M_PI) / 20.0;

		/* Central Difference with 5-points */
		gsl_deriv_central(&f, x, h, &results[0], &gsl_abs_errors[0]);

		rel_errors[0] = rel_abs_error(results[0], cos(x), &abs_errors[0]);
		printf("%9.2e %24.17e %24.17e %7.1e %7.1e %7.1e\n", x, results[0], cos(x), gsl_abs_errors[0], abs_errors[0], rel_errors[0]);
	}
	printf("\n");

	/* Get values of f'(x) by Forward diff, absolute and relative errors in [0, 2*PI] */
	printf("    x           Forward Diff               True f'(x)       gsl_abs Abs_err Rel_err\n");

	/* Initial guess for seeking the optimal stepsize */
	h = 1;

	for(i = 0; i <= 20; i++)
	{
		x = 0.0 + (double)i * (2.0 * M_PI) / 20.0;

		/* Forward Difference with 5-points */
		gsl_deriv_forward(&f, x, h, &results[1], &gsl_abs_errors[1]);

		rel_errors[1] = rel_abs_error(results[1], cos(x), &abs_errors[1]);
		printf("%9.2e %24.17e %24.17e %7.1e %7.1e %7.1e\n", x, results[1], cos(x), gsl_abs_errors[1], abs_errors[1], rel_errors[1]);
	}
	printf("\n");

	/* Get values of f'(x) by Backward diff, absolute and relative errors in [0, 2*PI] */
	printf("    x           Backward Diff              True f'(x)       gsl_abs Abs_err Rel_err\n");

	/* Initial guess for seeking the optimal stepsize */
	h = 1;

	for(i = 0; i <= 20; i++)
	{
		x = 0.0 + (double)i * (2.0 * M_PI) / 20.0;

		/* Backward Difference with 5-points */
		gsl_deriv_backward(&f, x, h, &results[2], &gsl_abs_errors[2]);

		rel_errors[2] = rel_abs_error(results[2], cos(x), &abs_errors[2]);
		printf("%9.2e %24.17e %24.17e %7.1e %7.1e %7.1e\n", x, results[2], cos(x), gsl_abs_errors[2], abs_errors[2], rel_errors[2]);
	}

	return 0;
}