/*splint spline interpolation */
#include "cddefines.h"
#include "splint.h"

void splint(double xa[], 
	  double ya[], 
	  double y2a[], 
	  long int n, 
	  double x, 
	  double *y)
{
	long int k, 
	  khi, 
	  klo;
	double a, 
	  b, 
	  h;

#	ifdef DEBUG_FUN
	fputs( "<+>splint()\n", debug_fp );
#	endif
	/* uses XA,YA deduced by SPLINE */

	klo = 1;
	khi = n;
	while( khi - klo > 1 )
	{
		k = (khi + klo)/2;
		if( xa[k-1] > x )
		{
			khi = k;
		}
		else
		{
			klo = k;
		}
	}

	h = xa[khi-1] - xa[klo-1];
	if( h == 0. )
	{
		puts( "[Stop in splint]" );
		exit(1);
	}

	a = (xa[khi-1] - x)/h;
	b = (x - xa[klo-1])/h;

	*y = a*ya[klo-1] + b*ya[khi-1] + ((POW3(a) - a)*
		y2a[klo-1] + (POW3(b) - b)*y2a[khi-1])*(POW2(h))/6.e0;

#	ifdef DEBUG_FUN
	fputs( " <->splint()\n", debug_fp );
#	endif
	return;
}