/***********************************************/
/* T.Kouya's GSL sample program collection     */
/*            Interpolation of Discrete Points */
/*                   Written by Tomonori Kouya */
/*                                             */
/* Version 0.01: 2007-10-04                    */
/***********************************************/
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>

/* Number of Discrete Points */
#define DIM 10

/* exact func(x) */
/* sin(x) */
double func(double x)
{
	return sin(x);
}

int main(void)
{
	int i, times, status;
	gsl_interp *workspace;
	gsl_interp_accel *accel;
	double h;
	double x[DIM], y[DIM], xtmp, ytmp;

/* Polynomial Interpolation */

	/* Set discrete points */
	x[0] = 0.0;
	x[DIM - 1] = 2 * M_PI;
	h = (x[DIM - 1] - x[0]) / DIM;
	for(i = 0; i < DIM; i++)
	{
		x[i] = x[0] + (double)i * h;
		y[i] = func(x[i]);
	}

	printf("Discrete Points for Interpolation (sin(x))\n");
	printf(" i                 x                     y\n");
	for(i = 0; i < DIM; i++)
		printf("%3d: %25.17e %25.17e\n", i, x[i], y[i]);
	printf("\n");

	/* Define type of Interpolation */
	workspace = gsl_interp_alloc(gsl_interp_polynomial, DIM);
	accel = gsl_interp_accel_alloc();

	/* Initialize */
	gsl_interp_init(workspace, x, y, DIM);

	/* Evaluation */
	h /= 3.0;
	printf("               x                 interpolated y      abs_error     \n");
	for(i = 0; i < 3 * DIM; i++)
	{
		xtmp = x[0] + h * (double)i;
		ytmp = gsl_interp_eval(workspace, x, y, xtmp, NULL);
		printf("%25.17e %25.17e %10.3e\n", xtmp, ytmp, fabs(ytmp - func(xtmp)));
	}

	/* free */
	gsl_interp_free(workspace);

	return 0;
}