/*******************************************************
*  CUBIC SPLINE INTERPOLATION OF A DISCRETE FUNCTION   *
*  F(X), GIVEN BY N POINTS X(I), Y(I)                  *
* ---------------------------------------------------- *
* SAMPLE RUN:                                          *
*                                                      *
* (Interpolate F(X), defined by 11 points:             *
*      0.0           0.0                               *
*      1.0           1.01                              *
*      2.0           3.99                              *
*      3.0           8.85                              *
*      4.0          15.0                               *
*      5.0          25.1                               *
*      6.0          37.0                               *
*      7.0          50.0                               *
*      8.0          63.9                               *
*      9.0          81.2                               *
*     10.0         100.5  for X = 8.7 )                *
*                                                      *
*  CUBIC SPLINE INTERPOLATION:                         *
*                                                      *
*  For X=  8.70000000   Y= 75.69211348                 *
*                                                      *
* ---------------------------------------------------- *
* Ref.: From Numath Library By Tuan Dang Trong in      *
*       Fortran 77.                                    *
*                                                      *
*                   C++ Release By J-P Moreau, Paris.  *
*                          (www.jpmoreau.fr)           *
*******************************************************/
#include <stdio.h>
#include <math.h>

#define  NMAX  25

  int    N;                       //Number of given points
  double X[NMAX],Y[NMAX];         //Coordinates of N points
  double B[NMAX],C[NMAX],D[NMAX]; //Coefficients of cubic spline
  double U,V;                     //Coordinates of interpolated point

  void VecPrint(char *s, int N, double *V) {
    int i;
    printf("\n%s\n", s);
	for (i=1; i<=N; i++) {
      printf(" %12.6f", V[i]);
	  if ((i%5)==0) printf("\n");
    }
	printf("\n");
  }

  double SEVAL(int N, double U, double *X, double *Y,
	           double *B, double *C, double *D)  {
/*----------------------------------------------------------------------
!     EVALUATE A CUBIC SPLINE INTERPOLATION OF A DISCRETE FUNCTION F(X),
!     GIVEN IN N POINTS X(I), Y(I). THE B, C AND D COEFFICIENTS DEFINING
!     THE BEST CUBIC SPLINE FOR THE GIVEN POINTS, ARE CALCULATED BEFORE
!     BY THE SPLINE SUBROUTINE.
!
!     INPUTS:
!     N       NUMBER OF POINTS OF CURVE Y = F(X)
!     U       ABSCISSA OF POINT TO BE INTERPOLATED
!     X,Y     POINTERS TO TABLES OF DIMENSION N, STORING THE 
!             COORDINATES OF CURVE F(X)
!     B,C,D   POINTERS TO TABLES STORING THE COEFFICIENTS DEFINING
!             THE CUBIC SPLINE
!
!     OUTPUT:
!     INTERPOLATED VALUE:
!             = Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I)))
!             WITH DX = U-X(I), U BETWEEN X(I) AND X(I+1)
!
!     REFERENCE :
!     FORSYTHE,G.E. (1977) COMPUTER METHODS FOR MATHEMATICAL
!     COMPUTATIONS. PRENTICE-HALL,INC.
!----------------------------------------------------------------------*/
//  Labels: e10, e20, e30
    double DX; int I,J,K;

      I=1;

//    BINARY SEARCH

      if (I >= N)  I = 1;
      if (U < X[I]) goto e10;
      if( U <= X[I+1]) goto e30;
e10:  I = 1;
      J = N+1;
e20:  K = (I+J) / 2;
      if (U < X[K])  J = K;
      if (U >= X[K]) I = K;
      if (J > I+1) goto e20;

//    SPLINE EVALUATION

e30:  DX = U-X[I];
      return (Y[I]+DX*(B[I]+DX*(C[I]+DX*D[I])));
  }


    void SPLINE (int N, double *X, double *Y, double *B, 
		         double *C, double *D)  {
/*--------------------------------------------------------------------
!     THIS SUBROUTINE CALCULATES THE COEFFICIENTS B,C,D OF A CUBIC
!     SPLINE TO BEST APPROXIMATE A DISCRETE FUNCTION GIVEN BY N POINTS
!
!     INPUTS:
!     N       NUMBER OF GIVEN POINTS
!     X,Y     VECTORS OF DIMENSION N, STORING THE COORDINATES
!             OF FUNCTION F(X)
!
!     OUTPUTS:
!     A,B,C   VECTORS OF DIMENSION N, STORING THE COEFFICIENTS
!             OF THE CUBIC SPLINE
!
!     REFERENCE:
!     FORSYTHE,G.E. (1977) COMPUTER METHODS FOR MATHEMATICAL
!     COMPUTATIONS. PRENTICE-HALL,INC.
!-------------------------------------------------------------------*/
//  Labels:  e15, e50
      int I, L, NM1;
      double T;

      NM1 = N - 1;

      if (N < 2) return;;
      if (N < 3) goto e50;

//    BUILD THE TRIDIAGONAL SYSTEM
//    B (DIAGONAL), D (UPPERDIAGONAL) , C (SECOND MEMBER)

      D[1] = X[2]-X[1];
      C[2] = (Y[2]-Y[1])/D[1];
      for (I = 2; I<N; I++) {
        D[I] = X[I+1]-X[I];
        B[I] = 2.0*(D[I-1]+D[I]);
        C[I+1] = (Y[I+1]-Y[I])/D[I];
        C[I] = C[I+1]-C[I];
      }

//    CONDITIONS AT LIMITS
//    THIRD DERIVATIVES OBTAINED BY DIVIDED DIFFERENCES

      B[1] = -D[1];
      B[N] = -D[N-1];
      C[1] = 0.0;
      C[N] = 0.0;
      if (N == 3) goto e15;
      C[1] = C[3]/(X[4]-X[2])-C[2]/(X[3]-X[1]);
      C[N] = C[N-1]/(X[N]-X[N-2])-C[N-2]/(X[N-1]-X[N-3]);
      C[1] = C[1]*D[1]*D[1]/(X[4]-X[1]);
      C[N] = -C[N]*D[N-1]*D[N-1]/(X[N]-X[N-3]);

//    FORWARD ELIMINATION

e15:  for (I = 2; I<=N; I++) {
        T = D[I-1]/B[I-1];
        B[I] = B[I]-T*D[I-1];
        C[I] = C[I]-T*C[I-1];
      }

//    BACK SUBSTITUTION

      C[N] = C[N]/B[N];
      for  (L = 1; L<N; L++) {
        I = N-L;
        C[I] = (C[I]-D[I]*C[I+1])/B[I];
      }

//    COEFFICIENTS OF 3RD DEGREE POLYNOMIAL

      B[N] = (Y[N]-Y[NM1])/D[NM1]+D[NM1]*(C[NM1]+2.0*C[N]);
      for (I = 1; I<N; I++) {
        B[I] = (Y[I+1]-Y[I])/D[I]-D[I]*(C[I+1]+2.0*C[I]);
        D[I] = (C[I+1]-C[I])/D[I];
        C[I] = 3.0*C[I];
      } 
      C[N] = 3.0*C[N];
      D[N] = D[NM1];
      return;

//    CASE N = 2

e50:  B[1] = (Y[2]-Y[1])/(X[2]-X[1]);
      C[1] = 0.0;
      D[1] = 0.0;
      B[2] = B[1];
      C[2] = 0.0;
      D[2] = 0.0;
   }

void main()  {

  N=11;

  X[1]=0.0; X[2]=1.0; X[3]=2.0; X[4]=3.0; X[5]=4.0; X[6]=5.0;
  X[7]=6.0; X[8]=7.0; X[9]=8.0; X[10]=9.0; X[11]=10.0;

  Y[1]=0.01; Y[2]=1.01; Y[3]=3.99; Y[4]=8.85; Y[5]=15.0; Y[6]=25.1;
  Y[7]=37.0; Y[8]=50.0; Y[9]=63.9; Y[10]=81.2; Y[11]=100.5;

  SPLINE (N,X,Y,B,C,D);

  VecPrint("B Vector",N,B);
  VecPrint("C Vector",N,C);
  VecPrint("D Vector",N,D);

  U = 9.2;

  V = SEVAL(N,U,X,Y,B,C,D);

  printf("\n  CUBIC SPLINE INTERPOLATION:\n");
  printf("\n  For X=%12.8f   Y=%12.8f\n", U, V);
  printf("\n");

}

// end of file tseval.cpp