/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
#include "cddefines.h"
#include "interpolate.h"

#define LIMSPLINE	1000
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] -> [n], */
	  u[LIMSPLINE]/*[NLIM] -> [n]*/, 
	  un;

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

	if( n > LIMSPLINE )
	{
		fprintf(ioQQQ," PROBLEM spline called with too many points.\n");
		TotalInsanity();
	}
	/* >>chng 02 sep 1529, put LIMSPLINE back in for speed up */
	/* >>chng 01 nov 29, removed NLIM dependence, PvH */
	/*if( ( u = (double*)MALLOC((size_t)(n*sizeof(double))) ) == NULL ) 
		BadMalloc();*/

	if( yp1 > .99e30 )
	{
		y2a[0] = 0.;
		u[0] = 0.;
	}
	else
	{
		y2a[0] = -0.5;
		u[0] = 3.e0/(x[1] - x[0])*((y[1] -  y[0])/(x[1] - x[0]) - 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 = (x[i] - x[i-1])/(x[i+1] - x[i-1]);
			p = sig*y2a[i-1] + 2.e0;
			y2a[i] = (sig - 1.e0)/p;
			u[i] = (6.e0*((y[i+1] - y[i])/(x[i+1] - x[i]) -
				      (y[i] - y[i-1])/(x[i] - x[i-1])) / 
				(x[i+1] - x[i-1]) - sig*u[i-1])/p;
		}
	}

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

	/* >>chng 02 sep 1529, put LIMSPLINE back in for speed up */
	/*free( u );*/

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

/*splint spline interpolation */
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 = 0;
	khi = n-1;
	while( khi - klo > 1 )
	{
		k = (khi + klo)/2;
		if( xa[k] > x )
		{
			khi = k;
		}
		else
		{
			klo = k;
		}
	}

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

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

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

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


/* wrapper routine for splint that checks whether x-value is within bounds
 * if the x-value is out of bounds, a flag will be raised and the function
 * will be evaluated at the nearest boundary */
/* >>chng 03 jan 15, added splint_safe, PvH */
void splint_safe(double xa[], 
		 double ya[], 
		 double y2a[], 
		 long int n, 
		 double x, 
		 double *y,
		 int *lgOutOfBounds)
{
	double xsafe;

#	ifdef DEBUG_FUN
	fputs( "<+>splint_safe()\n", debug_fp );
#	endif

	const double lo_bound = MIN2(xa[0],xa[n-1]);
	const double hi_bound = MAX2(xa[0],xa[n-1]);
	const double SAFETY = MAX2(hi_bound-lo_bound,1.)*10.*DBL_EPSILON;
	if( x < lo_bound-SAFETY )
	{
		xsafe = lo_bound;
		*lgOutOfBounds = TRUE;
	}
	else if( x > hi_bound+SAFETY )
	{
		xsafe = hi_bound;
		*lgOutOfBounds = TRUE;
	}
	else
	{
		xsafe = x;
		*lgOutOfBounds = FALSE;
	}

	splint(xa,ya,y2a,n,xsafe,y);

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


/*spldrv spline interpolation, returns derivative instead of function itself */
void spldrv(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( "<+>spldrv()\n", debug_fp );
#	endif
	/* uses XA,YA deduced by SPLINE */

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

	h = xa[khi] - xa[klo];
	if( h == 0. )
	{
		puts( "[Stop in spldrv]" );
		cdEXIT(EXIT_FAILURE);
	}

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

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

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


/* wrapper routine for spldrv that checks whether x-value is within bounds
 * if the x-value is out of bounds, a flag will be raised and the function
 * will be evaluated at the nearest boundary */
/* >>chng 03 jan 15, added spldrv_safe, PvH */
void spldrv_safe(double xa[], 
		 double ya[], 
		 double y2a[], 
		 long int n, 
		 double x, 
		 double *y,
		 int *lgOutOfBounds)
{
	double xsafe;

#	ifdef DEBUG_FUN
	fputs( "<+>spldrv_safe()\n", debug_fp );
#	endif

	const double lo_bound = MIN2(xa[0],xa[n-1]);
	const double hi_bound = MAX2(xa[0],xa[n-1]);
	const double SAFETY = MAX2(fabs(lo_bound),fabs(hi_bound))*10.*DBL_EPSILON;
	if( x < lo_bound-SAFETY )
	{
		xsafe = lo_bound;
		*lgOutOfBounds = TRUE;
	}
	else if( x > hi_bound+SAFETY )
	{
		xsafe = hi_bound;
		*lgOutOfBounds = TRUE;
	}
	else
	{
		xsafe = x;
		*lgOutOfBounds = FALSE;
	}

	spldrv(xa,ya,y2a,n,xsafe,y);

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


/*
 * series of interpolation routines from CUJ Feb 2000, by Kyle Loudon
 */

/***************************************************************
*                                                              *
*  -------- Interpolation Using Lagrange Polynomials --------  *
*                                                              *
*  On Call:                                                    *
*  x           Array of interpolation points                   *
*  fx          Array of known values of f(x) at each x         *
*  n           Count of values in x and fx                     *
*  X           Value of x at which to compute f(x)             *
*                                                              *
*  On Return:                                                  *
*  lagrange    Value of f(x) at X                              *
*                                                              *
***************************************************************/


double lagrange(const double *x, const double *fx, int n,
   double X) {

double             l, p = 0.0;
int                i, k;

for (k = 0; k < n; k++) 
{

   l = 1.0;

   /*
   ** Obtain lk(x) for the next term to add to the interpolating
   ** polynomial.
   */

   for (i = 0; i < n; i++) 
   {

      if (i != k)
         l = l * ((X - x[i]) / (x[k] - x[i]));

   }

   /*
   ** Add the next term computed from lk(x) to the interpolating
   ** polynomial.
   */

   p = p + (fx[k] * l);

}

return p;

}
/***************************************************************
*                                                              *
*  --------- Interpolation Using Newton Polynomials ---------  *
*                                                              *
*  On Call:                                                    *
*  x           Array of interpolation points                   *
*  fx          Array of known values of f(x) at each x         *
*  n           Count of values in x and fx                     *
*  X           Value of x at which to compute f(x)             *
*  tolerance   Delta for successive pk(x) at which to stop     *
*                                                              *
*  On Return:                                                  *
*  degree      Degree of the interpolating polynomial          *
*  newton      Value of f(x) at X, or DBL_MAX on errors        *
*                                                              *
***************************************************************/


double newton(const double *x, const double *fx, int n,
   double X, double tolerance, int *degree) {

double             *table, delta, m = 1.0, p;
int                i, k;

/*
** Allocate enough storage for the longest diagonal required in
** the divided-diffrence table.
*/

if ((table = (double *)MALLOC(sizeof(double) * (unsigned)n)) == NULL)
   return DBL_MAX;

/*
** Set the initial interpolating polynomial and first element of 
** the divided-difference table to f(x0).
*/

p = table[0] = fx[0];
*degree = 0;

for (k = 0; k < n - 1; k++) {

   (*degree)++;
   table[k + 1] = fx[k + 1];

   /*
   ** Compute the next diagonal placed in the divided-difference
   ** table to obtain the next coefficient.
   */

   for (i = k; i >= 0; i--)
      table[i] = (table[i + 1] - table[i]) / (x[k + 1] - x[i]);

   /*
   ** Compute the next pk(x), one degree larger than the last,
   ** and determine whether the specified tolerance will be met.
   */

   m = m * (X - x[k]);
   p = p + (delta = table[0] * m);

   if (fabs(delta) < tolerance) {

      free(table);
      return p;

   }

}

/*
** The specified tolerance was never met for differences between
** successive interpolating polynomials.
*/

free(table);
return DBL_MAX;

}
/***************************************************************
*                                                              *
*  --------- Piecewise-Cubic Hermite Interpolation ----------  *
*                                                              *
*  On Call:                                                    *
*  x           Array of strictly increasing points             *
*  fx          Array of known values of f(x) at each x         *
*  fpx         Array of known values of f'(x) at each x        *
*  n           Count of values in x, fx, and fpx               *
*  X           Value of x at which to compute f(x)             *
*                                                              *
*  On Return:                                                  *
*  pwhermite   Value of f(x) at X                              *
*                                                              *
***************************************************************/


double pwhermite(const double *x, const double *fx, const
   double *fpx, int n, double X) {

double             a[4], m[3], dd, dx;
int                k;

/*
** Locate the interval on which X lies. If this interval is not
** found by k = n - 3, [x[n - 2], x[n - 1]] is used.
*/

for (k = 0; k < n - 2; k++) {

   if (X < x[k + 1])
      break;

}

/*
** Compute values to be used in the computation of coefficients
** for the cubic interpolating polynomial.
*/

dx = x[k + 1] - x[k];
dd = (fx[k + 1] - fx[k]) / dx;

/*
** Compute the coefficients of the interpolating polynomial and
** determine multipliers based on center calculations involving
** X and the interpolation points.
*/

a[0] = fx[k];
a[1] = fpx[k];
a[2] = (dd - fpx[k]) / dx;
a[3] = (fpx[k + 1] + fpx[k] - (2.0 * dd)) / (dx * dx);
m[0] = X - x[k];
m[1] = m[0] * (X - x[k]);
m[2] = m[1] * (X - x[k + 1]);

/*
** Compute p3(x), the cubic interpolating polynomial, to obtain
** a value for f(x) at X.
*/

return a[0] + (a[1] * m[0]) + (a[2] * m[1]) + (a[3] * m[2]);

}
/***************************************************************
*                                                              *
*  ------------- Piecewise-Linear Interpolation -------------  *
*                                                              *
*  On Call:                                                    *
*  x           Array of strictly increasing points             *
*  fx          Array of known values of f(x) at each x         *
*  n           Count of values in x and fx                     *
*  X           Value of x at which to compute f(x)             *
*                                                              *
*  On Return:                                                  *
*  pwlinear    Value of f(x) at X                              *
*                                                              *
***************************************************************/


double pwlinear(const double *x, const double *fx, int n,
   double X) {

int                k;

/*
** Locate the interval on which X lies. If this interval is not
** found by k = n - 3, [x[n - 2], x[n - 1]] is used.
*/

for (k = 0; k < n - 2; k++) {

   if (X < x[k + 1])
      break;

}

/*
** Compute p1(x), the linear interpolating polynomial, to obtain
** a value for f(x) at X.     
*/

return fx[k] + (((fx[k + 1] - fx[k]) / (x[k + 1] - x[k])) * (X -
   x[k]));

}
