/******************************************************* * 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 #include #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