#include <ansi_c.h>
//---------------------------------------------------------------------------
//  This Unit contains routines for spline operations.
//  spline calculates spline coefficients y2[] for data arrays x[] and y[].
//  splint calculates spline interpolation points ys[], based upon y2[], at points xs[].
//  spli1 calculates first integral spline values using y2[].
//  spld1 calculates first derivative values.
//  spld2 calculates second derivative values.
//---------------------------------------------------------------------------

#define EPS 1e-3

//---------------------------------------------------------------------------
//  spline : Calculates spline coefficients from input data
//      Input:
//          double xa[], ya[] : input data arrays
//          int na : number of input data values
//      Output:
//          double y2[] : spline coeficient array
//---------------------------------------------------------------------------
void spline (double xa[], double ya[], int na, double y2[])
{
    int     i, k, n1;
    double  p, qn, sig, un, *u;

    n1 = na - 1;
    u = calloc (na, sizeof(double));
    y2[0] = 0.0;
    u[0] = 0.0;
    for ( i = 1; i < na-1; i++ ) {
        sig = (xa[i] - xa[i-1])/(xa[i+1] - xa[i-1]);
        p = sig*y2[i-1] + 2.0;
        y2[i] = (sig - 1.0)/p;
        u[i] = (ya[i+1] - ya[i])/(xa[i+1] - xa[i]) - (ya[i] - ya[i-1])/(xa[i] - xa[i-1]);
        u[i] = (6.0*u[i]/(xa[i+1] - xa[i-1]) - sig*u[i-1])/p;
    }
    qn = un = 0.0;
    y2[n1] = (un - qn*u[na-2])/(qn*y2[na-2] + 1.0);
    for (k = na-2; k >= 0; k--) {
        y2[k] = y2[k]*y2[k+1] + u[k];
    }

    free (u);

}
//---------------------------------------------------------------------------
//  splint :    Calculates spline interpolation from input data and spline coefficients.
//              Must call spline first.
//      Input :
//              double xa[], ya[] : input data
//              double y2[] : spline coefficients (calculated with call to spline)
//              int na : number of input data values
//              double xs[] : interpolation grid values
//              int ns : number of interpolation grid values
//      Output:
//              double ys[] : interpolation values
//---------------------------------------------------------------------------
int splint(double xa[], double ya[], double y2[], int na, double xs[], double ys[], int ns)
{
    int     klo, klo0, khi, i;
    double  a, b, c, d, x, dx;

    klo0 = 0;
    for (i = 0; i < ns; i++) {
        x = xs[i];
        // take care of degenerate cases
        if (x < xa[0] - EPS)
            ys[i] = 0.0;
        else if (x > xa[na-1] + EPS)
        	ys[i] = 0.0;
        else {
            for (klo = klo0; (xa[klo] < x) && (klo < na); klo++);       // find klo which falls just short of xs[i]
            if (klo > 0)
                klo--;
            if (klo == na-1)
                klo--;
            klo0 = klo;
            khi = klo + 1;              // xa[klo] and xa[khi] span xs[i]
            dx = xa[khi] - xa[klo];
            if (dx <= 0.0)
                return(1);                          // input data not monotonic
            a = (xa[khi] - x)/dx;
            b = (x - xa[klo])/dx;
            c = (a*a*a - a)*dx*dx/6.0;
            d = (b*b*b - b)*dx*dx/6.0;
            ys[i] = a*ya[klo] + b*ya[khi] + c*y2[klo] + d*y2[khi];
        }
    }

    return (0);

}

//---------------------------------------------------------------------------
//  spd1 :      Calculates the spline first derivative from input data and spline coefficients.
//              Must call spline first.
//      Input :
//              double xa[], ya[] : input data
//              doulbe y2[] : spline coefficients (calculated with call to spline)
//              int na : number of input data values
//              double xs[] : interpolation grid values
//              int ns : number of interpolation grid values
//      Output :
//              double  ys[] : derivative values
// ***** There is some error in this routine - values generated are off a bit
//---------------------------------------------------------------------------

int spld1(double xa[], double ya[], double y2[], int na, double xs[], double ys[], int ns)
{
    int     klo, klo0, khi, i;
    double  a, b, ap, bp, cp, dp, x, dx;

    klo0 = 0;
    for (i = 0; i < ns; i++) {
        x = xs[i];
        // take care of degenerate cases
        if (x < xa[0] - EPS)
            ys[i] = 0.0;
        else if (x > xa[na-1] + EPS)
        	ys[i] = 0.0;
        else {
            for (klo = klo0; (xa[klo] < x) && (klo < na); klo++);       // find klo which falls just short of xs[i]
            if (klo > 0)
                klo--;
            if (klo == na-1)
                klo--;
            klo0 = klo;
            khi = klo + 1;              // xa[klo] and xa[khi] span xs[i]
            dx = xa[khi] - xa[klo];
            if (dx <= 0.0)
                return(1);                          // input data not monotonic
            a = (xa[khi] - x)/dx;
            b = (x - xa[klo])/dx;
            ap = -1/dx;
            bp = 1/dx;
            cp = (3.0*a*a - 1.0)*ap*dx*dx/6.0;
            dp = (3.0*b*b - 1.0)*bp*dx*dx/6.0;
            ys[i] = ap*ya[klo] + bp*ya[khi] + cp*y2[klo] + dp*y2[khi];
        }
    }

    return (0);
}

//---------------------------------------------------------------------------
//  spli1 :    Calculates spline integration from input data and spline coefficients.
//              Must call spline first.
//              Calculate c0[] - continuity constants
//      Input :
//              double xa[], ya[] : input data
//              double y2[] : spline coefficients (calculated with call to spline)
//              int na : number of input data values
//              double xs[] : interpolation grid values
//              int ns : number of interpolation grid values
//      Output:
//              double ys[] : integration values
//---------------------------------------------------------------------------
int spli1(double xa[], double ya[], double y2[], int na, double xs[], double ys[], int ns)
{
    int     klo, klo0, khi, i;
    double  a, b, A, B, C, D, x, dx, dxo;
    double *alpha;

    // Calculate alpha
	alpha = calloc (na-1, sizeof(double));
    dx = xa[1] - xa[0];
    alpha[0] = ya[0]*dx/2.0;               // linear part
    alpha[0] -= y2[0]*dx*dx*dx/24.0;          // second derivative part
    dx = xa[1] - xa[0];
    for (i = 1; i < na-1; i++) {
        dxo = dx;
        dx = xa[i+1] - xa[i];
        alpha[i] = alpha[i-1] + ya[i]*(dxo + dx)/2.0;  // linear part
        alpha[i] -= y2[i]*(dxo*dxo*dxo + dx*dx*dx)/24.0;  // second derivative part
    }
    klo0 = 0;
    for (i = 0; i < ns; i++) {
        x = xs[i];
        // take care of degenerate cases
        if (x < xa[0] - EPS)
            ys[i] = 0.0;
        else if (x > xa[na-1] + EPS)
            ys[i] = ys[ns-1];
        else {
            for (klo = klo0; (xa[klo] < x) && (klo < na); klo++);       // find klo which falls just short of xs[i]
            if (klo > 0)
                klo--;
            if (klo == na-1)
                klo--;
            klo0 = klo;
            khi = klo + 1;
            dx = xa[khi] - xa[klo];
            if (dx <= 0.0)
                return(1);                          // input data not monotonic
            a = (xa[khi] - x)/dx;
            b = (x - xa[klo])/dx;
            A = -dx*a*a/2.0;
            B = dx*b*b/2.0;
            C = (-a*a*a*a/24.0 + a*a/12.0)*dx*dx*dx;
            D = (b*b*b*b/24.0 - b*b/12.0)*dx*dx*dx;
            ys[i] = alpha[klo] + A*ya[klo] +B*ya[khi];  // linear part
            ys[i] += C*y2[klo] + D*y2[khi];                 // second derivative part
        }
    }

    free (alpha);
    return (0);

}