#include <math.h>
#include <alloc.h>

#define   TNNearlyZero  1E-15
#define   ArraySize 30
#define   ArrayIntSize 300

typedef double *TNvector;
typedef struct  {
		     TNvector A, B, C, D;
		} TNcoefficients;
double _abs(x)
double x;
{
  return (x < 0) ? -x : x;
}

void CubicSplineFree();
/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-    Input:  NumPoints, XData, YData, NumInter, XInter                     -*/
/*-    Output: Coef0, Coef1, Coef2, Coef3, YInter, Error                     -*/
/*-                                                                          -*/
/*-            Purpose: To construct a cubic spline interpolant              -*/
/*-                     to a set of points. The second derivative            -*/
/*-                     of the interpolant is assumed to be zero             -*/
/*-                     at the endpoints (free cubic spline)                 -*/
/*-                     The spline is of the form                            -*/
/*-                                                                          -*/
/*-                        3               2                                 -*/
/*-            D[I](X-X[I])  + C[I](X-X[I])  + B[I](X-X[I]) + A[I]           -*/
/*-                                                                          -*/
/*-                     where 1 < I < NumPoints.                             -*/
/*-                                                                          -*/
/*-  User-defined Types: TNvector = ARRAY[0..TNArraySize] OF real;           -*/
/*-                                                                          -*/
/*- Global Variables: NumPoints : int;  Number of points                 -*/
/*-                   XData     : TNvector; X-coordinate data points         -*/
/*-                                         (must be in increasing order)    -*/
/*-                   YData     : TNvector; Y-coordinate data points         -*/
/*-                   Coef0     : TNvector; 0th coefficient of spline        -*/
/*-                   Coef1     : TNvector; 1st coefficient of spline        -*/
/*-                   Coef2     : TNvector; 2nd coefficient of spline        -*/
/*-                   Coef3     : TNvector; 3rd coefficient of spline        -*/
/*-                   NumInter  : int;  Number of interpolated points    -*/
/*-                   XInter    : TNvector; Points at which to interpolate   -*/
/*-                   YInter    : TNvector; Interpolated values at XInter    -*/
/*-                   Error     : byte;     0: flags if something goes wrong -*/
/*-                                                                          -*/
/*-              Errors: 0: No error                                         -*/
/*-                      1: X-coordinate points not                          -*/
/*-                         unique                                           -*/
/*-                      2: X-coordinate points not                          -*/
/*-                         in increasing order                              -*/
/*-                      3: NumPoints < 2                                    -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

void CalculateIntervals(NumPoints, XData, Interval, Error )
int NumPoints, *Error;
TNvector XData, *Interval;
/*----------------------------------------------------------*/
/*- Input: NumPoints, XData                                -*/
/*- Output: Interval, Error                                -*/
/*-                                                        -*/
/*- This function  calculates the length of the interval   -*/
/*- between two adjacent X values, contained in XData. If  -*/
/*- the X values are not sequential, Error = 2 is returned -*/
/*- and if ( the X values are not unique,  ) Error = 1 is  -*/
/*- returned.                                              -*/
/*----------------------------------------------------------*/
{
  int Index;

  *Interval = calloc(ArraySize,sizeof(double));
  *Error = 0;
  for ( Index = 1 ; Index <= NumPoints - 1; Index++ )
  {
    (*Interval)[Index] = XData[Index+1] - XData[Index];
    if ( _abs((*Interval)[Index]) < TNNearlyZero  )
      *Error = 1;     /* Data not unique */
    if ( (*Interval)[Index] < 0  )
      *Error = 2;     /* Data not in increasing order */
  }
} /* function CalculateIntervals */

void CalculateAs(NumPoints, YData, Spline )
int NumPoints;
TNvector YData;
TNcoefficients *Spline;
/*------------------------------------------*/
/*- Input: NumPoints, YData                -*/
/*- Ouput: Spline                          -*/
/*-                                        -*/
/*- This void calculates the constant -*/
/*- Term in each cubic.                    -*/
/*------------------------------------------*/
{
  int Index;

  for ( Index = 1 ; Index <= NumPoints ; Index++ )
    Spline->A[Index] = YData[Index];
} /* function CalculateAs */

void CalculateCs(NumPoints, XData, Interval, Spline )
int NumPoints;
TNvector XData, Interval;
TNcoefficients *Spline;
/*---------------------------------------------*/
/*- Input: NumPoints, XData, Interval         -*/
/*- Ouput: Spline                             -*/
/*-                                           -*/
/*- This void calculates the coefficient -*/
/*- of the squared Term in each cubic.        -*/
/*---------------------------------------------*/
{
  TNvector Alpha, L, Mu, Z;
  int Index;

  Alpha = calloc(ArraySize,sizeof(double));
  L = calloc(ArraySize,sizeof(double));
  Mu = calloc(ArraySize,sizeof(double));
  Z = calloc(ArraySize,sizeof(double));

  /* The next few lines solve a tridiagonal matrix */
  for ( Index = 2 ; Index <= NumPoints - 1;Index++ )
    Alpha[Index] = 3 * ((Spline->A[Index+1] * Interval[Index-1])
		       - (Spline->A[Index] * (XData[Index+1] - XData[Index-1]))
		       + (Spline->A[Index-1] * Interval[Index]))
		       / (Interval[Index-1] * Interval[Index]);
  L[1] = 0.0;
  Mu[1] = 0.0;
  Z[1] = 0.0;
  for ( Index = 2 ; Index <= NumPoints - 1; Index++ )
  {
    L[Index] = 2.0 * (XData[Index+1] - XData[Index-1])
		  - Interval[Index-1] * Mu[Index-1];
    Mu[Index] = Interval[Index]/L[Index];
    Z[Index] = (Alpha[Index] - Interval[Index-1] * Z[Index-1]) / L[Index];
  }
  /* Now calculate the C's */
  Spline->C[NumPoints] = 0.0;
  for ( Index = NumPoints - 1; Index >= 1; Index-- )
    Spline->C[Index] = Z[Index] - Mu[Index] * Spline->C[Index+1];
  free(Alpha);free(L);free(Mu);free(Z);
} /* function CalculateCs */

void CalculateBandDs(NumPoints, Interval, Spline )
int NumPoints;
TNvector Interval;
TNcoefficients *Spline;
/*------------------------------------------------*/
/*- Input: NumPoints, Interval                   -*/
/*- Ouput: Spline                                -*/
/*-                                              -*/
/*- This void calculates the coefficient of -*/
/*- the linear and cubic terms in each cubic.    -*/
/*------------------------------------------------*/
{
  int Index;

  for ( Index = NumPoints - 1;Index >= 1;Index-- )
  {
    Spline->B[Index] = (Spline->A[Index+1] - Spline->A[Index]) / Interval[Index]
		 - Interval[Index] * (Spline->C[Index+1] + 2.0 * Spline->C[Index]) / 3.0;
    Spline->D[Index] = (Spline->C[Index+1] - Spline->C[Index]) / (3.0 * Interval[Index]);
  }
} /* function CalculateDs */


void CalculateCoefficients(NumPoints, XData, YData, Interval, Spline )
int NumPoints;
TNvector XData, YData, Interval;
TNcoefficients *Spline;
/*---------------------------------------------------------------*/
/*- Input: NumPoints, XData, YData, Interval                    -*/
/*- Output: Spline                                              -*/
/*-                                                             -*/
/*- This void calculates the coefficients of each cubic    -*/
/*- in the interpolating spline. A separate cubic is calculated -*/
/*- for every interval between data points.  The coefficients   -*/
/*- are returned in the variable Spline.                        -*/
/*---------------------------------------------------------------*/
{ /* function CalculateCoefficients */
  Spline->A = calloc(ArraySize,sizeof(double));
  Spline->B = calloc(ArraySize,sizeof(double));
  Spline->C = calloc(ArraySize,sizeof(double));
  Spline->D = calloc(ArraySize,sizeof(double));

  CalculateAs(NumPoints, YData, Spline);
  CalculateCs(NumPoints, XData, Interval, Spline);
  CalculateBandDs(NumPoints, Interval, Spline);
} /* function CalculateCoefficients */

void Interpolate(NumPoints, XData, Spline, NumInter, XInter, YInter )
int NumPoints, NumInter;
TNvector XData, XInter, *YInter;
TNcoefficients Spline;
/*---------------------------------------------------------------*/
/*- Input: NumPoints, XData, Spline, NumInter, XInter           -*/
/*- Output: YInter                                              -*/
/*-                                                             -*/
/*- This function  uses the interpolating cubic spline (Spline) -*/
/*- to interpolate values for the X positions contained         -*/
/*- in XInter.  The interpolated values are returned in YInter. -*/
/*---------------------------------------------------------------*/
{
  int Index, Location, Term;
  double X;

  *YInter = calloc(ArrayIntSize,sizeof(double));

  for ( Index = 1 ; Index <= NumInter ; Index++ )
  {
    Location = 1;
    for ( Term = 1 ; Term <= NumPoints - 1 ; Term++)
      if ( XInter[Index] > XData[Term]  )
	Location = Term;
    X = XInter[Index] - XData[Location];
    (*YInter)[Index] = ((Spline.D[Location] * X + Spline.C[Location]) * X +
			 Spline.B[Location]) * X + Spline.A[Location];
  }
} /* function Interpolate */


void CubicSplineFree(NumPoints, XData, YData, NumInter, XInter,
		     Coef0, Coef1, Coef2, Coef3, YInter, Error )
int NumPoints, NumInter, *Error;
TNvector XData, YData, XInter, *Coef0, *Coef1, *Coef2, *Coef3, *YInter;
{ /* function CubicSplineFree */
  TNvector Interval;          /* Intervals between adjacent points */
  TNcoefficients Spline;      /* All the cubics */

  *Error = 0;
  if ( NumPoints < 2  )
    *Error = 3;
  else
    CalculateIntervals(NumPoints, XData, &Interval, Error);
  if ( *Error == 0  )
  {
    CalculateCoefficients(NumPoints, XData, YData, Interval, &Spline);
    Interpolate(NumPoints, XData, Spline, NumInter, XInter, YInter);
  }
  *Coef0 = Spline.A;
  *Coef1 = Spline.B;
  *Coef2 = Spline.C;
  *Coef3 = Spline.D;
} /* function CubicSplineFree */
