#include "cddefines.h"
#include "spline.h"
#	ifdef NLIM
#undef NLIM
#	endif
#define	NLIM	160

void spline(double x[], 
  double y[], 
  long int n, 
  double yp1, 
  double ypn, 
  double y2a[])
{
	long int i, 
	  k;
	double p, 
	  qn, 
	  sig, 
	  u[NLIM], 
	  un;

#	ifdef DEBUG_FUN
	fputs( "<+>spline()\n", debug_fp );
#	endif
	/*spline interpolation routine */
	/* returns array Y2 needed for SPLINT */

	/* check that not too many points */
	if( n > NLIM )
	{
		fprintf( ioQQQ, " Too many points for SPLINE.  Increase NLIM\n" );
		puts( "[Stop in spline]" );
		exit(1);
	}

	if( yp1 > .99e30 )
	{
		y2a[0] = 0.;
		u[0] = 0.;
	}
	else
	{
		y2a[0] = -0.5;
		u[0] = (3.e0/((double)(x[1]) - (double)(x[0])))*(((double)(y[1]) - 
		  (double)(y[0]))/((double)(x[1]) - (double)(x[0])) - (double)(yp1));
	}

	for( i=1; i < (n - 1); i++ )
	{
		if( (x[i+1] == x[i-1] || x[i+1] == x[i]) || x[i] == x[i-1] )
		{
			u[i] = u[i-1];
		}
		else
		{
			sig = ((double)(x[i]) - (double)(x[i-1]))/((double)(x[i+1]) - 
			  (double)(x[i-1]));
			p = sig*y2a[i-1] + 2.e0;
			y2a[i] = (sig - 1.e0)/p;
			u[i] = (6.e0*(((double)(y[i+1]) - (double)(y[i]))/((double)(x[i+1]) - 
			  (double)(x[i])) - ((double)(y[i]) - (double)(y[i-1]))/
			  ((double)(x[i]) - (double)(x[i-1])))/((double)(x[i+1]) - 
			  (double)(x[i-1])) - sig*u[i-1])/p;
		}
	}

	if( ypn > .99e30 )
	{
		qn = 0.;
		un = 0.;
	}
	else
	{
		qn = 0.5;
		un = (3.e0/((double)(x[n-1]) - (double)(x[n-2])))*((double)(ypn) - 
		  ((double)(y[n-1]) - (double)(y[n-2]))/((double)(x[n-1]) - 
		  (double)(x[n-2])));
	}

	y2a[n-1] = (un - qn*u[n-2])/(qn*y2a[n-2] + 1.e0);
	for( k=n - 2; k >= 0; k-- )
	{
		y2a[k] = y2a[k]*y2a[k+1] + u[k];
	}


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