]>
Commit | Line | Data |
---|---|---|
0795afa3 | 1 | #include "isajet/pilot.h" |
2 | SUBROUTINE SPLINE(X,C,N,IBCBEG,IBCEND) | |
3 | C********************************************************************** | |
4 | C* Computes the coefficient of a cubic interpolating spline. The X(i),* | |
5 | C* i=1,...,N, are the knots or x-values of data points; C(1,i) are * | |
6 | C* the corresponding y-values. N is the number of data points (>3!). * | |
7 | C* IBCBEG = 0 means that the slope at X(1) is unknown, in which case * | |
8 | C* it is determined from requiring a smooth 3rd derivative at x(2); * | |
9 | C* IBCBEG = 1 means that the slope is known, in which case it has to * | |
10 | C* be stored in C(1,2). IBCEND has the same meaning for the end of * | |
11 | C* x-region; if IBCEND = 1, the slope is to be stored in C(2,N). The * | |
12 | C* routine then computes the coefficients C(l,i) of the i-th spline, * | |
13 | C* written in the form * | |
14 | C* f_i(x) = C(1,i) + h_i C(2,i) + h_i^2 C(3,i) + h_i^3 C(4,i), where * | |
15 | C* h_i = x - X(I). * | |
16 | C Modified from contributed subroutine by M. Drees, 1/14/99 | |
17 | C********************************************************************** | |
18 | #if defined(CERNLIB_IMPNONE) | |
19 | IMPLICIT NONE | |
20 | #endif | |
21 | INTEGER N,IBCBEG,IBCEND,I,L,M,J | |
22 | REAL C(4,N),X(N),G,DTAU,DIVDF1,DIVDF3 | |
23 | C | |
24 | L = N - 1 | |
25 | DO 10 M = 2, N | |
26 | C(3,M) = X(M) - X(M-1) | |
27 | 10 C(4,M) = (C(1,M) - C(1,M-1))/C(3,M) | |
28 | C First slope unknown | |
29 | IF(IBCBEG.EQ.0) THEN | |
30 | C(4,1) = C(3,3) | |
31 | C(3,1) = C(3,2) + C(3,3) | |
32 | C(2,1) = ( (C(3,2)+2.0*C(3,1))*C(4,2)*C(3,3) + | |
33 | & C(3,2)**2*C(4,3) )/C(3,1) | |
34 | C First slope already known | |
35 | ELSE | |
36 | C(4,1) = 1.0 | |
37 | C(3,1) = 0.0 | |
38 | ENDIF | |
39 | C Forward pass of Gauss elimination | |
40 | DO 20 M = 2, L | |
41 | G = -C(3,M+1)/C(4,M-1) | |
42 | C(2,M) = G*C(2,M-1) + 3.0*(C(3,M)*C(4,M+1)+C(3,M+1)*C(4,M)) | |
43 | 20 C(4,M) = G*C(3,M-1) + 2.0*(C(3,M)+C(3,M+1)) | |
44 | ||
45 | IF(IBCEND.EQ.0) THEN | |
46 | G = C(3,N-1) + C(3,N) | |
47 | C(2,N) = ( (C(3,N)+2.0*G)*C(4,N)*C(3,N-1) + | |
48 | & C(3,N)**2*(C(1,N-1)-C(1,N-2))/C(3,N-1) )/G | |
49 | G = -G/C(4,N-1) | |
50 | C(4,N) = (G+1.0)*C(3,N-1) + C(4,N) | |
51 | C(2,N) = ( G*C(2,N-1) + C(2,N) )/C(4,N) | |
52 | ENDIF | |
53 | C Back substitution | |
54 | DO 30 J = L,1,-1 | |
55 | 30 C(2,J) = ( C(2,J) - C(3,J)*C(2,J+1) )/C(4,J) | |
56 | C Computation of coefficients | |
57 | DO 40 I = 2,N | |
58 | DTAU = C(3,I) | |
59 | DIVDF1 = (C(1,I)-C(1,I-1))/DTAU | |
60 | DIVDF3 = C(2,I-1) + C(2,I) - 2.0*DIVDF1 | |
61 | C(3,I-1) = ( DIVDF1 - C(2,I-1) - DIVDF3 ) / DTAU | |
62 | 40 C(4,I-1) = DIVDF3/DTAU/DTAU | |
63 | C | |
64 | RETURN | |
65 | END |