/*  spline.c  19 November 2001 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <sys/types.h>
#include <time.h>

#include "spline.h"

/**********************************************************************/

void parabola_val2 ( int ndata, float tdata[], float ydata[], int left, 
  float tval, float *yval ) {

/**********************************************************************/

/*
  Purpose:

    PARABOLA_VAL2 evaluates a parabolic function through 3 points in a table.

  Discussion:

    This routine is a utility routine used by OVERHAUSER_SPLINE_VAL.
    It constructs the parabolic interpolant through the data in
    3 consecutive entries of a table and evaluates this interpolant
    at a given abscissa value.

  Modified:

    28 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, int NDATA, the number of data points.
    NDATA must be at least 3.

    Input, float TDATA[NDATA], the abscissas of the data points.  The
    values in TDATA must be in strictly ascending order.

    Input, float YDATA[NDATA], the data points corresponding to
    the abscissas.

    Input, int LEFT, the location of the first of the three
    consecutive data points through which the parabolic interpolant
    must pass.  0 <= LEFT <= NDATA - 3.

    Input, float TVAL, the value of T at which the parabolic interpolant
    is to be evaluated.  Normally, TDATA[0] <= TVAL <= T[NDATA-1], and 
    the data will be interpolated.  For TVAL outside this range, 
    extrapolation will be used.

    Output, float *YVAL, the value of the parabolic interpolant 
    at TVAL.
*/
  float dif1;
  float dif2;
  float t1;
  float t2;
  float t3;
  float y1;
  float y2;
  float y3;

/* 
  Check. 
*/

  if ( left < 0 || left > ndata-3 ) {
    printf ( "\n" );
    printf ( "PARABOLA_VAL2 - Fatal error!\n" );
    printf ( "  LEFT < 0 or LEFT > NDATA-3.\n" );
    exit ( EXIT_FAILURE );
  }

/* 
  Copy out the three abscissas. 
*/

  t1 = tdata[left];
  t2 = tdata[left+1];
  t3 = tdata[left+2];

  if ( t1 >= t2 || t2 >= t3 ) {
    printf ( "\n" );
    printf ( "PARABOLA_VAL2 - Fatal error!\n" );
    printf ( "  T1 >= T2 or T2 >= T3.\n" );
    printf ( "  T1 = %f, T2 = %f, T3 = %f.\n", t1, t2, t3 );
    exit ( EXIT_FAILURE );
  }

/* 
  Construct and evaluate a parabolic interpolant for the data. 
*/

  y1 = ydata[left];
  y2 = ydata[left+1];
  y3 = ydata[left+2];

  dif1 = ( y2 - y1 ) / ( t2 - t1 );
  dif2 =
    ( ( y3 - y1 ) / ( t3 - t1 )
    - ( y2 - y1 ) / ( t2 - t1 ) ) / ( t3 - t2 );

  *yval = y1 + ( tval - t1 ) * ( dif1 + ( tval - t2 ) * dif2 );

  return;
}
/**********************************************************************/

void rvec_bracket ( int n, float x[], float xval, int *left, int *right ) {

/**********************************************************************/

/*
  Purpose:

    RVEC_BRACKET searches a sorted array for successive brackets of a value.

  Discussion:

    If the values in the vector are thought of as defining intervals
    on the real line, then this routine searches for the interval
    nearest to or containing the given value.

  Modified:

    08 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, int N, length of input array.

    Input, float X[N], an array that has been sorted into ascending order.

    Input, float XVAL, a value to be bracketed.

    Output, int *LEFT, *RIGHT, the results of the search.
    Either:
      XVAL < X[0], when *LEFT = 0, *RIGHT = 1;
      XVAL > X[N-1], when *LEFT = N-2, *RIGHT = N-1;
    or
      X[*LEFT] <= XVAL <= X[*RIGHT].
*/
  int i;

  for ( i = 1; i <= n - 2; i++ ) {

    if ( xval < x[i] ) {
      *left = i - 1;
      *right = i;
      return;
    }

   }

  *left = n - 2;
  *right = n - 1;

  return;
}
/**********************************************************************/

void rvec_bracket3 ( int n, float t[], float tval, int *left ) {

/**********************************************************************/

/*
  Purpose:

    RVEC_BRACKET3 finds the interval containing or nearest a given value.

  Discussion:

    The routine always returns the index LEFT of the sorted array
    T with the property that either
    *  T is contained in the interval [ T[LEFT], T[LEFT+1] ], or
    *  T < T[LEFT] = T[0], or
    *  T > T[LEFT+1] = T[N-1].

    The routine is useful for interpolation problems, where
    the abscissa must be located within an interval of data
    abscissas for interpolation, or the "nearest" interval
    to the (extreme) abscissa must be found so that extrapolation
    can be carried out.

  Modified:

    05 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, int N, length of the input array.

    Input, float T[N], an array that has been sorted into ascending order.

    Input, float TVAL, a value to be bracketed by entries of T.

    Input/output, int *LEFT.

    On input, if 0 <= LEFT <= N-2, LEFT is taken as a suggestion for the
    interval [ T[LEFT] T[LEFT+1] ] in which TVAL lies.  This interval
    is searched first, followed by the appropriate interval to the left
    or right.  After that, a binary search is used.

    On output, LEFT is set so that the interval [ T[LEFT], T[LEFT+1] ]
    is the closest to TVAL; it either contains TVAL, or else TVAL
    lies outside the interval [ T[0], T[N-1] ].
*/
  int high;
  int low;
  int mid;

/*  
  Check the input data.
*/
  if ( n < 2 ) {
    printf ( "\n" );
    printf ( "RVEC_BRACKET3 - Fatal error!\n" );
    printf ( "  N must be at least 2.\n" );
    exit ( EXIT_FAILURE );
  }
/*
  If *LEFT is not between 0 and N-2, set it to the middle value.
*/
  if ( *left < 0 || *left > n - 2 ) {
    *left = ( n - 1 ) / 2;
  }

/*
  CASE 1: TVAL < T[*LEFT]:
  Search for TVAL in (T[I],T[I+1]), for I = 0 to *LEFT-1.
*/
  if ( tval < t[*left] ) {

    if ( *left == 0 ) {
      return;
    }
    else if ( *left == 1 ) {
      *left = 0;
      return;
    }
    else if ( tval >= t[*left-1] ) {
      *left = *left - 1;
      return;
    }
    else if ( tval <= t[1] ) {
      *left = 0;
      return;
    }

/* 
  ...Binary search for TVAL in (T[I],T[I+1]), for I = 1 to *LEFT-2.
*/
    low = 1;
    high = *left - 2;

    for (;;) {

      if ( low == high ) {
        *left = low;
        return;
      }

      mid = ( low + high + 1 ) / 2;

      if ( tval >= t[mid] ) {
        low = mid;
      }
      else {
        high = mid - 1;
      }

    }
  }
/* 
  CASE 2: T[*LEFT+1] < TVAL:
  Search for TVAL in (T[I},T[I+1]) for intervals I = *LEFT+1 to N-2.
*/
  else if ( tval > t[*left+1] ) {

    if ( *left == n - 2 ) {
      return;
    }
    else if ( *left == n - 3 ) {
      *left = *left + 1;
      return;
    }
    else if ( tval <= t[*left+2] ) {
      *left = *left + 1;
      return;
    }
    else if ( tval >= t[n-2] ) {
      *left = n - 2;
      return;
    }
/* 
  ...Binary search for TVAL in (T[I],T[I+1]) for intervals I = *LEFT+2 to N-3.
*/
    low = *left + 2;
    high = n - 3;

    for (;;) {

      if ( low == high ) {
        *left = low;
        return;
      }

      mid = ( low + high + 1 ) / 2;

      if ( tval >= t[mid] ) {
        low = mid;
      }
      else {
        high = mid - 1;
      }
    }
  }
/*
  CASE 3: T[*LEFT] <= TVAL <= T[*LEFT+1]:
  T is just where the user said it might be.
*/
  else {
  }

  return;
}
/**********************************************************************/

void rvec_order_type ( int n, float x[], int *order ) {

/**********************************************************************/

/*
  Purpose:

    RVEC_ORDER_TYPE determines if an array is (non)strictly ascending/descending.

  Modified:

    14 September 2000

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of entries of the array.

    Input, float X[N], the array to be checked.

    Output, int *ORDER, order indicator:
    -1, no discernable order;
    0, all entries are equal;
    1, ascending order;
    2, strictly ascending order;
    3, descending order;
    4, strictly descending order.
*/
  int i;

/* Search for the first value not equal to X(0). */

  i = 0;

  for (;;) {

    i = i + 1;
    if ( i > n-1 ) {
      *order = 0;
      return;
    }

    if ( x[i] > x[0] ) {
      if ( i == 1 ) {
        *order = 2;
        break;
      }
      else {
        *order = 1;
        break;
      }
    }
    else if ( x[i] < x[0] ) {
      if ( i == 1 ) {
        *order = 4;
        break;
      }
      else {
        *order = 3;
        break;
      }
    }
  }
/* 
  Now we have a "direction".  Examine subsequent entries. 
*/
  for (;;) {

    i = i + 1;
    if ( i > n - 1 ) {
      return;
    }

    if ( *order == 1 ) {

      if ( x[i] < x[i-1] ) {
        *order = -1;
        return;
      }
    }
    else if ( *order == 2 ) {

      if ( x[i] < x[i-1] ) {
        *order = -1;
        return;
      }
      else if ( x[i] == x[i-1] ) {
        *order = 1;
      }
    }
    else if ( *order == 3 ) {

      if ( x[i] > x[i-1] ) {
        *order = -1;
        return;
      }
    }
    else if ( *order == 4 ) {

      if ( x[i] > x[i-1] ) {
        *order = -1;
        return;
      }
      else if ( x[i] == x[i-1] ) {
        *order = 3;
      }

    }

  }

}
/**********************************************************************/

int s3_fs ( float a1[], float a2[], float a3[], int n, float b[], float x[] ) {

/**********************************************************************/

/*
  Purpose:

    S3_FS factors and solves a tridiagonal linear system.

  Note:

    This algorithm requires that each diagonal entry be nonzero.

  Modified:

    03 February 1999

  Author:

    John Burkardt

  Parameters:

    Input/output, float A1[N], A2[N], A3[N].
    On input, the nonzero diagonals of the linear system.
    On output, the data in these vectors has been overwritten
    by factorization information.  Entries A1[0] and A3[N-1] are not used.

    Input, int N, the order of the linear system.

    Input/output, float B[N].
    On input, B contains the right hand side of the linear system.
    On output, B has been overwritten by factorization information.

    Output, float X[N], the solution of the linear system.

    Output, int S3_FS, error indicator:
    0, no error,
    1, a diagonal entry was zero.
*/
  int i;
  float xmult;
/*
  Check.
*/
  for ( i = 0; i < n; i++ ) {
    if ( a2[i] == 0.0 ) {
      printf ( "\n" );
      printf ( "S3_FS - Fatal error!\n" );
      printf ( "  A2[%d] = 0.\n", i );
      return 1;
    }
  }

  for ( i = 1; i < n-1; i++ ) {

    xmult = a1[i] / a2[i-1];
    a2[i] = a2[i] - xmult * a3[i-1];

    b[i] = b[i] - xmult * b[i-1];

  }

  xmult = a1[n-1] / a2[n-2];
  a2[n-1] = a2[n-1] - xmult * a3[n-2];

  x[n-1] = ( b[n-1] - xmult * b[n-2] ) / a2[n-1];
  for ( i = n-2; i >= 0; i-- ) {
    x[i] = ( b[i] - a3[i] * x[i+1] ) / a2[i];
  }

  return 0;
}
/**********************************************************************/

void setdif ( float diftab[], int ntab, float xtab[], float ytab[] ) {

/**********************************************************************/

/*

  Purpose:

    SETDIF sets up a divided difference table from raw data.

  Discussion:

    Space can be saved by using a single array for both the DIFTAB and
    YTAB dummy parameters.  In that case, the difference table will
    overwrite the Y data without interfering with the computation.

  Modified:

    12 April 1999

  Author:

    John Burkardt

  Parameters:

    Output, real DIFTAB[NTAB], the divided difference coefficients 
    corresponding to the input (XTAB,YTAB).

    Input, integer NTAB, the number of pairs of points
    (XTAB[I],YTAB[I]) which are to be used as data.  

    Input, real XTAB[NTAB], the X values at which data was taken.
    These values must be distinct.

    Input, real YTAB[NTAB], the corresponding Y values.
*/
  int i;
  int j;
/*
  Copy the data values into DIFTAB.
*/
  for ( i = 0; i < ntab; i++ ) {
    diftab[i] = ytab[i];
  }
/*
  Make sure the abscissas are distinct.
*/
  for ( i = 0; i < ntab; i++ ) {
    for ( j = i+1; j < ntab; j++ ) {

      if ( xtab[i] - xtab[j] == 0.0 ) {
        printf ( "\n" );
        printf ( "SETDIF - Fatal error!\n" );
        printf ( "  Two entries of XTAB are equal!\n" );
        printf ( "  XTAB[%d] = %f\n", i, xtab[i] );
        printf ( "  XTAB[%d] = %f\n", j, xtab[j] );
        exit ( 1 );
      }
    }
  }
/*
  Compute the divided differences.
*/
  for ( i = 1; i <= ntab-1; i++ ) {
    for ( j = ntab-1; j >= i; j-- ) {
      diftab[j] = ( diftab[j] - diftab[j-1] ) / ( xtab[j] - xtab[j-i] );

    }
  }
 
  return;
}
/**********************************************************************/

void spline_b_val ( int ndata, float tdata[], float ydata[], float tval, 
  float *yval ) {

/**********************************************************************/

/*
  Purpose:

    SPLINE_B_VAL evaluates a cubic B spline approximant.

  Discussion:

    The cubic B spline will approximate the data, but is not
    designed to interpolate it.

    In effect, two "phantom" data values are appended to the data,
    so that the spline will interpolate the first and last data values.

  Modified:

    12 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, int NDATA, the number of data values.

    Input, float TDATA[NDATA], the abscissas of the data.

    Input, float YDATA[NDATA], the data values.

    Input, float TVAL, a point at which the spline is to be evaluated.

    Output, float *YVAL, the value of the function at TVAL.
*/
  float bval;
  int left;
  int right;
  float u;
/*
  Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL.
*/
  rvec_bracket ( ndata, tdata, tval, &left, &right );
/*
  Evaluate the 5 nonzero B spline basis functions in the interval,
  weighted by their corresponding data values.
*/
  u = ( tval - tdata[left] ) / ( tdata[right] - tdata[left] );
  *yval = 0.0;
/*
  B function associated with node LEFT - 1, (or "phantom node"),
  evaluated in its 4th interval.
*/
  bval = ( 1.0 + u * ( - 3.0 + u * ( + 3.0 - u ) ) ) / 6.0;

  if ( left-1 >= 0 ) {
    *yval = *yval + ydata[left-1] * bval;
  }
  else {
    *yval = *yval + ( 2.0 * ydata[0] - ydata[1] ) * bval;
  }
/*
  B function associated with node LEFT,
  evaluated in its third interval.
*/
  bval = ( 4.0 + u * u * ( - 6.0 + 3.0 * u ) ) / 6.0;
  *yval = *yval + ydata[left] * bval;
/*
  B function associated with node RIGHT,
  evaluated in its second interval.
*/
  bval = ( 1.0 + u * ( 3.0 + u * ( 3.0 - 3.0 * u ) ) ) / 6.0;
  *yval = *yval + ydata[right] * bval;
/*
  B function associated with node RIGHT+1, (or "phantom node"),
  evaluated in its first interval.
*/
  bval = u * u * u / 6.0;
  if ( right+1 < ndata ) {
    *yval = *yval + ydata[right+1] * bval;
  }
  else {
    *yval = *yval + ( 2.0 * ydata[ndata-1] - ydata[ndata-2] ) * bval;
  }

  return;
}
/**********************************************************************/

void spline_beta_val ( float beta1, float beta2, int ndata, float tdata[], 
  float ydata[], float tval, float *yval ) {

/**********************************************************************/

/*
  Purpose:

    SPLINE_BETA_VAL evaluates a cubic beta spline approximant.

  Discussion:

    The cubic beta spline will approximate the data, but is not
    designed to interpolate it.

    If BETA1 = 1 and BETA2 = 0, the cubic beta spline will be the
    same as the cubic B spline approximant.

    With BETA1 = 1 and BETA2 large, the beta spline becomes more like
    a linear spline.

    In effect, two "phantom" data values are appended to the data,
    so that the spline will interpolate the first and last data values.

  Modified:

    12 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, float BETA1, the skew or bias parameter.
    BETA1 = 1 for no skew or bias.

    Input, float BETA2, the tension parameter.
    BETA2 = 0 for no tension.

    Input, int NDATA, the number of data values.

    Input, float TDATA[NDATA], the abscissas of the data.

    Input, float YDATA[NDATA], the data values.

    Input, float TVAL, a point at which the spline is to be evaluated.

    Output, float YVAL, the value of the function at TVAL.
*/
  float a;
  float b;
  float bval;
  float c;
  float d;
  float delta;
  int left;
  int right;
  float u;
/*
  Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL.
*/
  rvec_bracket ( ndata, tdata, tval, &left, &right );
/*
  Evaluate the 5 nonzero beta spline basis functions in the interval,
  weighted by their corresponding data values.
*/
  u = ( tval - tdata[left] ) / ( tdata[right] - tdata[left] );

  delta = 2.0 + beta2 + beta1 * ( 4.0 + beta1 * ( 4.0 + 2.0 * beta1 ) );

  *yval = 0.0;
/*
  Beta function associated with node LEFT - 1, (or "phantom node"),
  evaluated in its 4th interval.
*/
  bval = ( 2.0 * beta1 * beta1 * beta1 
    * ( 1.0 - u ) * ( 1.0 - u ) * ( 1.0 - u ) ) / delta;

  if ( left-1 >= 0 ) {
    *yval = *yval + ydata[left-1] * bval;
  }
  else {
    *yval = *yval + ( 2.0 * ydata[0] - ydata[1] ) * bval;
  }
/*
  Beta function associated with node LEFT,
  evaluated in its third interval.
*/
  a = beta2 + 4.0 * beta1 + 4.0 * beta1 * beta1;

  b = - 6.0 * beta1 * ( 1.0 + beta1 ) * ( 1.0 - beta1 );

  c = - 3.0 * ( beta2 + 2.0 * beta1 * beta1 * ( 1.0 + beta1 ) );

  d = 2.0 * ( beta2 + beta1 * ( 1.0 + beta1 * ( 1.0 + beta1 ) ) );

  bval = ( a + u * ( b + u * ( c + u * d ) )  ) / delta;

  *yval = *yval + ydata[left] * bval;
/*
  Beta function associated with node RIGHT,
  evaluated in its second interval.
*/
  a = 2.0;

  b = 6.0 * beta1;

  c = 3.0 * beta2 + 6.0 * beta1 * beta1;

  d = - 2.0 * ( 1.0 + beta2 + beta1 * ( 1.0 + beta1 ) );

  bval = ( a + u * ( b + u * ( c + u * d ) )  ) / delta;

  *yval = *yval + ydata[right] * bval;
/*
  Beta function associated with node RIGHT+1, (or "phantom node"),
  evaluated in its first interval.
*/
  bval = 2.0 * u * u * u / delta;

  if ( right+1 <= ndata-1 ) {
    *yval = *yval + ydata[right+1] * bval;
  }
  else {
    *yval = *yval + ( 2.0 * ydata[ndata-1] - ydata[ndata-2] ) * bval;
  }

  return;
}
/**********************************************************************/

int spline_cubic_set ( float diag[], int ibcbeg, int ibcend, int n, 
  float sub[], float sup[], float t[], float y[], float ybcbeg, 
  float ybcend, float ypp[] ) {

/**********************************************************************/

/*
  Purpose:

    SPLINE_CUBIC_SET computes the second derivatives of a cubic spline.

  Discussion:

    For data interpolation, the user must call SPLINE_SET to determine
    the second derivative data, passing in the data to be interpolated,
    and the desired boundary conditions.

    The data to be interpolated, plus the SPLINE_SET output, defines
    the spline.  The user may then call SPLINE_VAL to evaluate the
    spline at any point.

    The cubic spline is a piecewise cubic polynomial.  The intervals
    are determined by the "knots" or abscissas of the data to be
    interpolated.  The cubic spline has continous first and second
    derivatives over the entire interval of interpolation.

    For any point T in the interval T(IVAL), T(IVAL+1), the form of
    the spline is

      SPL(T) = A(IVAL)
             + B(IVAL) * ( T - T(IVAL) )
             + C(IVAL) * ( T - T(IVAL) )**2
             + D(IVAL) * ( T - T(IVAL) )**3

    If we assume that we know the values Y(*) and YPP(*), which represent
    the values and second derivatives of the spline at each knot, then
    the coefficients can be computed as:

      A(IVAL) = Y(IVAL)
      B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
      C(IVAL) = YPP(IVAL) / 2
      D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )

    Since the first derivative of the spline is

      SPL'(T) =     B(IVAL)
              + 2 * C(IVAL) * ( T - T(IVAL) )
              + 3 * D(IVAL) * ( T - T(IVAL) )**2,

    the requirement that the first derivative be continuous at interior
    knot I results in a total of N-2 equations, of the form:

      B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1))
      + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))**2 = B(IVAL)

    or, setting H(IVAL) = T(IVAL+1) - T(IVAL)

      ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
      - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6
      + YPP(IVAL-1) * H(IVAL-1)
      + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2
      =
      ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
      - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6

    or

      YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) )
      + YPP(IVAL) * H(IVAL)
      =
      6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
      - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)

    Boundary conditions must be applied at the first and last knots.  
    The resulting tridiagonal system can be solved for the YPP values.

  Modified:

    04 February 1999

  Author:

    John Burkardt

  Parameters:

    Workspace, float DIAG[N].

    Input, int IBCBEG, left boundary condition flag:

      0: the cubic spline should be a quadratic over the first interval;
      1: the first derivative at the left endpoint should be YBCBEG;
      2: the second derivative at the left endpoint should be YBCBEG.

    Input, int IBCEND, right boundary condition flag:

      0: the cubic spline should be a quadratic over the last interval;
      1: the first derivative at the right endpoint should be YBCEND;
      2: the second derivative at the right endpoint should be YBCEND.

    Input, int N, the number of data points.  N must be at least 2.
    In the special case where N = 2 and IBCBEG = IBCEND = 0, the
    spline will actually be linear.

    Workspace, float SUB[N], SUP[N].

    Input, float T[N], the knot values, that is, the points were data is
    specified.  The knot values should be distinct, and increasing.

    Input, float Y[N], the data values to be interpolated.

    Input, float YBCBEG, YBCEND, the values to be used in the boundary
    conditions if IBCBEG or IBCEND is equal to 1 or 2.

    Output, float YPP[N], the second derivatives of the cubic spline.

    Output, int SPLSET, error flag.
    0, no error;
    1, N is illegal;
    2, the T's are not monotonically increasing;
    3, IBCBEG is illegal;
    4, IBCEND is illegal;
    5, the linear system could not be solved properly.
*/
  int i;
/*
  Check.
*/
  if ( n <= 1 ) {
    printf ( "\n" );
    printf ( "SPLINE_CUBIC_SET - Fatal error!\n" );
    printf ( "  The number of data points N must be at least 2.\n" );
    printf ( "  The input value is %d.\n", n );
    return 1;
  }

  for ( i = 0; i < n - 1; i++ ) {
    if ( t[i] >= t[i+1] ) {
      printf ( "\n" );
      printf ( "SPLINE_CUBIC_SET - Fatal error!\n" );
      printf ( "  The knots must be strictly increasing, but\n" );
      printf ( "  T(%d) = %f\n", i, t[i] );
      printf ( "  T(%d) = %f\n", i+1, t[i+1] );
      return 2;
    }
  }
/*
  Set up the first equation.
*/
  if ( ibcbeg == 0 ) {
    ypp[0] = 0.0;
    diag[0] = 1.0;
    sup[0] = -1.0;
  }
  else if ( ibcbeg == 1 ) {
    ypp[0] = ( y[1] - y[0] ) / ( t[1] - t[0] ) - ybcbeg;
    diag[0] = ( t[1] - t[0] ) / 3.0;
    sup[0] = ( t[1] - t[0] ) / 6.0;
  }
  else if ( ibcbeg == 2 ) {
    ypp[0] = ybcbeg;
    diag[0] = 1.0;
    sup[0] = 0.0;
  }
  else {
    printf ( "\n" );
    printf ( "SPLINE_CUBIC_SET - Fatal error!\n" );
    printf ( "  IBCBEG must be 0, 1 or 2.\n" );
    printf ( "  The input value is %d.\n", ibcbeg );
    return 3;
  }
/*
  Set up the intermediate equations.
*/
  for ( i = 1; i < n-1; i++ ) {
    ypp[i] = ( y[i+1] - y[i] ) / ( t[i+1] - t[i] )
      - ( y[i] - y[i-1] ) / ( t[i] - t[i-1] );
    sub[i] = ( t[i] - t[i-1] ) / 6.0;
    diag[i] = ( t[i+1] - t[i-1] ) / 3.0;
    sup[i] = ( t[i+1] - t[i] ) / 6.0;
  }
/*
  Set up the last equation.
*/
  if ( ibcend == 0 ) {
    ypp[n-1] = 0.0;
    sub[n-1] = -1.0;
    diag[n-1] = 1.0;
  }
  else if ( ibcend == 1 ) {
    ypp[n-1] = ybcend - ( y[n-1] - y[n-2] ) / ( t[n-1] - t[n-2] );
    sub[n-1] = ( t[n-1] - t[n-2] ) / 6.0;
    diag[n-1] = ( t[n-1] - t[n-2] ) / 3.0;
  }
  else if ( ibcend == 2 ) {
    ypp[n-1] = ybcend;
    sub[n-1] = 0.0;
    diag[n-1] = 1.0;
  }
  else {
    printf ( "\n" );
    printf ( "SPLINE_CUBIC_SET - Fatal error!\n" );
    printf ( "  IBCEND must be 0, 1 or 2.\n" );
    printf ( "  The input value is %d.\n", ibcend );
    return 4;
  }
/*
  Solve the linear system.
*/
  if ( n == 2 && ibcbeg == 0 && ibcend == 0 ) {

    ypp[0] = 0.0;
    ypp[1] = 0.0;

  }
  else {

    i = s3_fs ( sub, diag, sup, n, ypp, ypp );

    if ( i != 0 ) {
      printf ( "\n" );
      printf ( "SPLINE_CUBIC_SET - Fatal error!\n" );
      printf ( "  The linear system could not be solved.\n" );
      return 5;
    }

  }

  return 0;
}
/**********************************************************************/

float spline_cubic_val ( int n, float t[], float tval, float y[], 
  float ypp[], float *ypval, float *yppval ) {

/**********************************************************************/

/*
  Purpose:

    SPLINE_CUBIC_VAL evaluates a cubic spline at a specific point.

  Discussion:

    SPLINE_CUBIC_SET must have already been called to define the values of YPP.

    For any point T in the interval T(IVAL), T(IVAL+1), the form of
    the spline is

      SPL(T) = A
             + B * ( T - T(IVAL) )
             + C * ( T - T(IVAL) )**2
             + D * ( T - T(IVAL) )**3

    Here:
      A = Y(IVAL)
      B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
      C = YPP(IVAL) / 2
      D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )

  Modified:

    04 February 1999

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of knots.

    Input, float Y[N], the data values at the knots.

    Input, float T[N], the knot values.

    Input, float TVAL, a point, typically between T[0] and T[N-1], at
    which the spline is to be evalulated.  If TVAL lies outside
    this range, extrapolation is used.

    Input, float Y[N], the data values at the knots.

    Input, float YPP[N], the second derivatives of the spline at
    the knots.

    Output, float *YPVAL, the derivative of the spline at TVAL.

    Output, float *YPPVAL, the second derivative of the spline at TVAL.

    Output, float SPLINE_VAL, the value of the spline at TVAL.
*/
  float dt;
  float h;
  int i;
  int ival;
  float yval;
/*
  Determine the interval [ T(I), T(I+1) ] that contains TVAL.
  Values below T[0] or above T[N-1] use extrapolation.
*/
  ival = n - 2;

  for ( i = 0; i < n-1; i++ ) {
    if ( tval < t[i+1] ) {
      ival = i;
      break;
    }
  }
/*
  In the interval I, the polynomial is in terms of a normalized
  coordinate between 0 and 1.
*/
  dt = tval - t[ival];
  h = t[ival+1] - t[ival];

  yval = y[ival]
    + dt * ( ( y[ival+1] - y[ival] ) / h
           - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
    + dt * ( 0.5 * ypp[ival]
    + dt * ( ( ypp[ival+1] - ypp[ival] ) / ( 6.0 * h ) ) ) );

  *ypval = ( y[ival+1] - y[ival] ) / h
    - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
    + dt * ( ypp[ival]
    + dt * ( 0.5 * ( ypp[ival+1] - ypp[ival] ) / h ) );

  *yppval = ypp[ival] + dt * ( ypp[ival+1] - ypp[ival] ) / h;

  return yval;
}
/**********************************************************************/

void spline_cubic_val2 ( int n, float t[], float tval, int *left, float y[], 
  float ypp[], float *yval, float *ypval, float *yppval ) {

/**********************************************************************/

/*
  Purpose:

   SPLINE_CUBIC_VAL2 evaluates a cubic spline at a specific point.

  Discussion:

    This routine is a modification of SPLINE_CUBIC_VAL; it allows the
    user to speed up the code by suggesting the appropriate T interval
    to search first.

    SPLINE_CUBIC_SET must have already been called to define the
    values of YPP.

    In the LEFT interval, let RIGHT = LEFT+1.  The form of the spline is

    SPL(T) =
      A
    + B * ( T - T[LEFT] )
    + C * ( T - T[LEFT] )**2
    + D * ( T - T[LEFT] )**3

    Here:
      A = Y[LEFT]
      B = ( Y[RIGHT] - Y[LEFT] ) / ( T[RIGHT] - T[LEFT] )
        - ( YPP[RIGHT] + 2 * YPP[LEFT] ) * ( T[RIGHT] - T[LEFT] ) / 6
      C = YPP[LEFT] / 2
      D = ( YPP[RIGHT] - YPP[LEFT] ) / ( 6 * ( T[RIGHT] - T[LEFT] ) )

  Modified:

    04 March 1999

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of knots.

    Input, float T[N], the knot values.

    Input, float TVAL, a point, typically between T[0] and T[N-1], at
    which the spline is to be evalulated.  If TVAL lies outside
    this range, extrapolation is used.

    Input/output, int *LEFT, the suggested T interval to search.
    LEFT should be between 0 and N-2.  If LEFT is not in this range,
    then its value will be ignored.  On output, LEFT is set to the
    actual interval in which TVAL lies.

    Input, float Y[N], the data values at the knots.

    Input, float YPP[N], the second derivatives of the spline at
    the knots.

    Output, float *YVAL, *YPVAL, *YPPVAL, the value of the spline, and
    its first two derivatives at TVAL.
*/
  float dt;
  float h;
  int right;

/* 
  Determine the interval [T[LEFT], T[RIGHT]] that contains TVAL.  
*/

/*  
  What you want from RVEC_BRACKET3 is that TVAL is to be computed
  by the data in interval [T[LEFT], T[RIGHT]].  
*/

  rvec_bracket3 ( n, t, tval, left );
 
  right = *left + 1;

/*
 In the interval LEFT, the polynomial is in terms of a normalized
 coordinate  ( DT / H ) between 0 and 1.
*/

  right = *left + 1;

  dt = tval - t[*left];
  h = t[right] - t[*left];

  *yval = y[*left]
     + dt * ( ( y[right] - y[*left] ) / h
        - ( ypp[right] / 6.0 + ypp[*left] / 3.0 ) * h
     + dt * ( 0.5 * ypp[*left]
     + dt * ( ( ypp[right] - ypp[*left] ) / ( 6.0 * h ) ) ) );

  *ypval = ( y[right] - y[*left] ) / h
     - ( ypp[right] / 6.0 + ypp[*left] / 3.0 ) * h
     + dt * ( ypp[*left]
     + dt * ( 0.5 * ( ypp[right] - ypp[*left] ) / h ) );

  *yppval = ypp[*left] + dt * ( ypp[right] - ypp[*left] ) / h;

  return;
}
/****************************************************************/

void spline_hermite_set ( int ndata, float tdata[], float c0[],
  float c1[], float c2[], float c3[] ) {

/****************************************************************/

/*

  Purpose:

    SPLINE_HERMITE_SET sets up a piecewise cubic Hermite interpolant.

  Reference:

    Conte and de Boor,
    Algorithm CALCCF,
    Elementary Numerical Analysis,
    1973, page 235.

  Modified:

    12 April 1999

  Parameters:

    Input, int NDATA, the number of data points.
    NDATA must be at least 2.

    Input, float TDATA[NDATA], the abscissas of the data points.
    The entries of TDATA are assumed to be strictly increasing.

    Input/output, float C0[NDATA], C1[NDATA}, C2[NDATA], C3[NDATA].

    On input, C0[I] and C1[I] should contain the value of the
    function and its derivative at TDATA[I], for I = 0 to NDATA-1.
    These values will not be changed by this routine.

    On output, C2[I] and C3[I] contain the quadratic
    and cubic coefficients of the Hermite polynomial
    in the interval (TDATA[I], TDATA[I+1]), for I = 0 to NDATA-2.
    C2[NDATA-1] and C3[NDATA-1] are set to 0.

    In the interval (TDATA[I], TDATA[I+1]), the interpolating Hermite
    polynomial is given by

    SVAL(TVAL) =             C0[I]
   + ( TVAL - TDATA(I) ) * ( C1[I]
   + ( TVAL - TDATA(I) ) * ( C2[I]
   + ( TVAL - TDATA(I) ) *   C3[I] ) )
*/
  float divdif1;
  float divdif3;
  float dt;
  int i;

  for ( i = 0; i < ndata-1; i++ ) {
    dt = tdata[i+1] - tdata[i];
    divdif1 = ( c0[i+1] - c0[i] ) / dt;
    divdif3 = c1[i] + c1[i+1] - 2.0 * divdif1;
    c2[i] = ( divdif1 - c1[i] - divdif3 ) / dt;
    c3[i] = divdif3 / ( dt * dt );
  }

  c2[ndata-1] = 0.0;
  c3[ndata-1] = 0.0;

  return;
}
/****************************************************************/

void spline_hermite_val ( int ndata, float tdata[], float c0[], 
  float c1[], float c2[], float c3[], float tval, float *sval ) {

/****************************************************************/

/*
  Purpose:

    SPLINE_HERMITE_VAL evaluates a piecewise cubic Hermite interpolant.

  Discussion:

    SPLINE_HERMITE_SET must be called first, to set up the
    spline data from the raw function and derivative data.

  Reference:

    Conte and de Boor,
    Algorithm PCUBIC,
    Elementary Numerical Analysis,
    1973, page 234.

  Modified:

    12 April 1999

  Parameters:

    Input, int NDATA, the number of data points.
    NDATA is assumed to be at least 2.

    Input, float TDATA(NDATA), the abscissas of the data points.
    The entries of TDATA are assumed to be strictly increasing.

    Input, float C0[NDATA], C1[NDATA], C2[NDATA], C3[NDATA], contains 
    the piecewise polynomial coefficent data computed by SPLINE_HERMITE_SET.

    Input, float TVAL, the point where the interpolant is to
    be evaluated.

    Output, float *SVAL, the value of the interpolant at TVAL.
*/
  float dt;
  int left;
  int right;
/*
  Find the interval [ TDATA[LEFT], TDATA[RIGHT] ] that contains
  or is nearest to TVAL.
*/
  rvec_bracket ( ndata, tdata, tval, &left, &right );
/*
  Evaluate the cubic polynomial.
*/
  dt = tval - tdata[left];

  *sval =         c0[left]
         + dt * ( c1[left]
         + dt * ( c2[left]
         + dt *   c3[left] ) );

  return;
}
/**********************************************************************/

void spline_linear_val ( int ndata, float tdata[], float ydata[], 
  float tval, float *yval, float *ypval ) {

/**********************************************************************/

/*
  Purpose:

    SPLINE_LINEAR_VAL evaluates a linear spline at a specific point.

  Discussion:

    Because of the extremely simple form of the linear spline,
    the raw data points ( TDATA(I), YDATA(I)) can be used directly to
    evaluate the spline at any point.  No processing of the data
    is required.

  Modified:

    12 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, int NDATA, the number of data points defining the spline.

    Input, float TDATA[NDATA], YDATA[NDATA], the values of the independent
    and dependent variables at the data points.  The values of TDATA should
    be distinct and increasing.

    Input, float TVAL, the point at which the spline is to be evaluated.

    Output, float *YVAL, *YPVAL, the value of the spline and its first
    derivative dYdT at TVAL.  YPVAL is not reliable if TVAL is exactly
    equal to TDATA(I) for some I.
*/
  int left;
  int right;
/*
  Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] that contains, or is
  nearest to, TVAL.
*/
  rvec_bracket ( ndata, tdata, tval, &left, &right );
/*
  Now evaluate the piecewise linear function.
*/
  *ypval = ( ydata[right] - ydata[left] ) / ( tdata[right] - tdata[left] );

  *yval = ydata[left] +  ( tval - tdata[left] ) * (*ypval);

  return;
}
/**********************************************************************/

void spline_overhauser_val ( int ndata, float tdata[], float ydata[], 
  float tval, float *yval ) {

/**********************************************************************/

/*
  Purpose:

    SPLINE_OVERHAUSER_VAL evaluates an Overhauser spline.

  Discussion:

    Over the first and last intervals, the Overhauser spline is a 
    quadratic.  In the intermediate intervals, it is a piecewise cubic.
    The Overhauser spline is also known as the Catmull-Rom spline.

  Reference:

    H Brewer and D Anderson,
    Visual Interaction with Overhauser Curves and Surfaces,
    SIGGRAPH 77, pages 132-137.

    E Catmull and R Rom,
    A Class of Local Interpolating Splines,
    in Computer Aided Geometric Design,
    edited by R Barnhill and R Reisenfeld,
    Academic Press, 1974, pages 317-326.

    David Rogers and Alan Adams,
    Mathematical Elements of Computer Graphics,
    McGraw Hill, 1990, Second Edition, pages 278-289.

  Modified:

    28 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, int NDATA, the number of data points.
    NDATA must be at least 3.

    Input, float TDATA[NDATA], the abscissas of the data points.  The
    values in TDATA must be in strictly ascending order.

    Input, float YDATA[NDATA], the data points corresponding to
    the abscissas.

    Input, float TVAL, the abscissa value at which the spline
    is to be evaluated.  Normally, TDATA[0] <= TVAL <= T[NDATA-1], and 
    the data will be interpolated.  For TVAL outside this range, 
    extrapolation will be used.

    Output, float *YVAL, the value of the spline at TVAL.
*/

  int left;
  int order;
  int right;
  float yl;
  float yr;

/* Check. */

  rvec_order_type ( ndata, tdata, &order );

  if ( order != 2 ) {
    printf ( "\n" );
    printf ( "SPLINE_OVERHAUSER_VAL - Fatal error!\n" );
    printf ( "  The data abscissas are not strictly ascending.\n" );
    exit ( EXIT_FAILURE );
  }

  if ( ndata < 3 ) {
    printf ( "\n" );
    printf ( "SPLINE_OVERHAUSER_VAL - Fatal error!\n" );
    printf ( "  NDATA < 3.\n" );
    exit ( EXIT_FAILURE );
  }

/* 
  Locate the abscissa interval T[LEFT], T[LEFT+1] nearest to or 
  containing TVAL. 
*/

  rvec_bracket ( ndata, tdata, tval, &left, &right );

/* 
  Evaluate the "left hand" quadratic defined at 
  T[LEFT-1], T[LEFT], T[RIGHT]. 
*/

  if ( left > 0 ) {
    parabola_val2 ( ndata, tdata, ydata, left-1, tval, &yl );
  }

/* 
  Evaluate the "right hand" quadratic defined at 
  T[LEFT], T[RIGHT], T[RIGHT+1]. 
*/

  if ( right < ndata - 1 ) {
    parabola_val2 ( ndata, tdata, ydata, left, tval, &yr );
  }

/* 
  Blend the quadratics. 
*/

  if ( left == 0 ) {
    *yval = yr;
  }
  else if ( right == ndata - 1 ) {
    *yval = yl;
  }
  else {
    *yval = (
        ( tdata[right] - tval ) * yl
        + ( tval - tdata[left] ) * yr 
     ) / ( tdata[right] - tdata[left] );
  }

  return;
}
/**********************************************************************/

void valdif ( float diftab[], int ntab, float xtab[], float xval, 
  float *yval ) {

/**********************************************************************/

/*
  Purpose:

    VALDIF evaluates a divided difference polynomial at a point.

  Discussion:

    SETDIF must be called first to set up the divided difference table.

  Modified:

    12 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, float DIFTAB[NTAB], the divided difference table.

    Input, integer NTAB, the number of divided difference
    coefficients in DIFTAB, and the number of points XTAB.

    Input, float XTAB[NTAB], the X values upon which the
    divided difference polynomial is based.

    Input, float XVAL, a value of X at which the polynomial
    is to be evaluated.

    Output, float *YVAL, the value of the polynomial at XVAL.
*/
  int i;

  *yval = diftab[ntab-1];
  for ( i = 2; i <= ntab; i++ ) {
    *yval = diftab[ntab-i] + ( xval - xtab[ntab-i] ) * (*yval);
  }
 
  return;
}