V3FIT
splaan.f
1 C******************** START FILE SPLAAN.FOR ; GROUP TRKRLIB ******************
2 C-------------------------------------------
3 C derived from Forsyth/Malcolm/Moler SPLINE.FOR -- different BC:
4 C y'=0 at LHS of data interval (low index side)
5 C
6 CCCCCCCCCCCCCCC
7  SUBROUTINE splaan (N, X, Y, B, C, D)
8  INTEGER N
9  REAL X(N), Y(N), B(N), C(N), D(N)
10 C
11 C THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED
12 C FOR A CUBIC INTERPOLATING SPLINE FOR WHICH S'(X1)=0.
13 C
14 C S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3
15 C
16 C FOR X(I) .LE. X .LE. X(I+1)
17 C
18 C INPUT..
19 C
20 C N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2)
21 C X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER
22 C Y = THE ORDINATES OF THE KNOTS
23 C
24 C OUTPUT..
25 C
26 C B, C, D = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE.
27 C
28 C USING P TO DENOTE DIFFERENTIATION,
29 C
30 C Y(I) = S(X(I))
31 C B(I) = SP(X(I))
32 C C(I) = SPP(X(I))/2
33 C D(I) = SPPP(X(I))/6 (DERIVATIVE FROM THE RIGHT)
34 C
35 CCCCCCCCCCCCCCC
36 C THE ACCOMPANYING FUNCTION SUBPROGRAM SEVAL CAN BE USED
37 C TO EVALUATE THE SPLINE.
38 C
39 C
40  INTEGER NM1, IB, I
41  REAL T
42 C
43  nm1 = n-1
44  IF ( n .LT. 2 ) RETURN
45  IF ( n .LT. 3 ) GO TO 50
46 C
47 C SET UP TRIDIAGONAL SYSTEM
48 C
49 C B = DIAGONAL, D = OFFDIAGONAL, C = RIGHT HAND SIDE.
50 C
51  d(1) = x(2) - x(1)
52  c(2) = (y(2) - y(1))/d(1)
53  DO 10 i = 2, nm1
54  d(i) = x(i+1) - x(i)
55  b(i) = 2.*(d(i-1) + d(i))
56  c(i+1) = (y(i+1) - y(i))/d(i)
57  c(i) = c(i+1) - c(i)
58  10 CONTINUE
59 C
60 C END CONDITIONS. THIRD DERIVATIVES AT X(1) AND X(N)
61 C OBTAINED FROM DIVIDED DIFFERENCES
62 C
63  b(1)=2.*d(1)
64  b(n) = -d(n-1)
65  c(1) = 0.
66  c(n) = 0.
67  IF ( n .EQ. 3 ) GO TO 15
68  c(1) = (y(2)-y(1))/d(1)
69  c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3))
70  c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3))
71 C
72 C FORWARD ELIMINATION
73 C
74  15 DO 20 i = 2, n
75  t = d(i-1)/b(i-1)
76  b(i) = b(i) - t*d(i-1)
77  c(i) = c(i) - t*c(i-1)
78  20 CONTINUE
79 C
80 C BACK SUBSTITUTION
81 C
82  c(n) = c(n)/b(n)
83  DO 30 ib = 1, nm1
84  i = n-ib
85  c(i) = (c(i) - d(i)*c(i+1))/b(i)
86  30 CONTINUE
87 C
88 C C(I) IS NOW THE SIGMA(I) OF THE TEXT
89 C
90 C COMPUTE POLYNOMIAL COEFFICIENTS
91 C
92  b(n) = (y(n) - y(nm1))/d(nm1) + d(nm1)*(c(nm1) + 2.*c(n))
93  DO 40 i = 1, nm1
94  b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.*c(i))
95  d(i) = (c(i+1) - c(i))/d(i)
96  c(i) = 3.*c(i)
97  40 CONTINUE
98  c(n) = 3.*c(n)
99  d(n) = d(n-1)
100  RETURN
101 C
102  50 b(1) = (y(2)-y(1))/(x(2)-x(1))
103  c(1) = 0.
104  d(1) = 0.
105  b(2) = b(1)
106  c(2) = 0.
107  d(2) = 0.
108  RETURN
109  END
110 C******************** END FILE SPLAAN.FOR ; GROUP TRKRLIB ******************