V3FIT
spline_ez1.f
1 C******************** START FILE SPLINE.FOR ; GROUP TRKRLIB ******************
2 CC
3 CCCCCCCCCCCCCCCC
4 C THE CODES (SPLINE & SEVAL) ARE TAKEN FROM:
5 C FORSYTHE,MALCOLM AND MOLER, "COMPUTER METHODS FOR
6 C MATHEMATICAL COMPUTATIONS",PRENTICE-HALL, 1977.
7 C
8 C THE CODES (SPLEEN,SPLAAN & SPEVAL) ARE ADAPTATIONS
9 C BY R.M. WIELAND FOR SPECIAL CASES ... SEE COMMENTS
10 
11 
12  SUBROUTINE spline_ez1 (N, X, Y, B, C, D)
13  INTEGER N
14  REAL X(N), Y(N), B(N), C(N), D(N)
15 C
16 C THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED
17 C FOR A CUBIC INTERPOLATING SPLINE
18 C
19 C S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3
20 C
21 C FOR X(I) .LE. X .LE. X(I+1)
22 C
23 C INPUT..
24 C
25 C N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2)
26 C X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER
27 C Y = THE ORDINATES OF THE KNOTS
28 C
29 C OUTPUT..
30 C
31 C B, C, D = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE.
32 C
33 C USING P TO DENOTE DIFFERENTIATION,
34 C
35 C Y(I) = S(X(I))
36 C B(I) = SP(X(I))
37 C C(I) = SPP(X(I))/2
38 C D(I) = SPPP(X(I))/6 (DERIVATIVE FROM THE RIGHT)
39 C
40 CCCCCCCCCCCCCCC
41 C THE ACCOMPANYING FUNCTION SUBPROGRAM SEVAL CAN BE USED
42 C TO EVALUATE THE SPLINE.
43 C
44 C
45  INTEGER NM1, IB, I
46  REAL T
47 C
48  nm1 = n-1
49  IF ( n .LT. 2 ) RETURN
50 C
51 C
52 C DMC - CHECK ORDINATES
53  DO 5 i=1,nm1
54  ip1=i+1
55  IF(x(ip1).LE.x(i)) then
56  write(6,9000)
57  write(6,9001) (x(ii),ii=1,n)
58  zbomb=sqrt(-1.0/(x(1)*x(1)))
59  write(6,*) zbomb
60  write(6,*) ' spline.for: code will stop.'
61  call exit(1)
62  9000 FORMAT(/
63  >' ? UNORDERED ORDINATE ARRAY PASSED TO SPLINE SUBROUTINE:')
64  9001 FORMAT(1x,5(1x,1pe12.5))
65  endif
66  5 CONTINUE
67 C
68  IF ( n .LT. 3 ) GO TO 50
69 C
70 C SET UP TRIDIAGONAL SYSTEM
71 C
72 C B = DIAGONAL, D = OFFDIAGONAL, C = RIGHT HAND SIDE.
73 C
74  d(1) = x(2) - x(1)
75  c(2) = (y(2) - y(1))/d(1)
76  DO 10 i = 2, nm1
77  d(i) = x(i+1) - x(i)
78  b(i) = 2.*(d(i-1) + d(i))
79  c(i+1) = (y(i+1) - y(i))/d(i)
80  c(i) = c(i+1) - c(i)
81  10 CONTINUE
82 C
83 C END CONDITIONS. THIRD DERIVATIVES AT X(1) AND X(N)
84 C OBTAINED FROM DIVIDED DIFFERENCES
85 C
86  b(1) = -d(1)
87  b(n) = -d(n-1)
88  c(1) = 0.
89  c(n) = 0.
90  IF ( n .EQ. 3 ) GO TO 15
91  c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1))
92  c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3))
93  c(1) = c(1)*d(1)**2/(x(4)-x(1))
94  c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3))
95 C
96 C FORWARD ELIMINATION
97 C
98  15 DO 20 i = 2, n
99  t = d(i-1)/b(i-1)
100  b(i) = b(i) - t*d(i-1)
101  c(i) = c(i) - t*c(i-1)
102  20 CONTINUE
103 C
104 C BACK SUBSTITUTION
105 C
106  c(n) = c(n)/b(n)
107  DO 30 ib = 1, nm1
108  i = n-ib
109  c(i) = (c(i) - d(i)*c(i+1))/b(i)
110  30 CONTINUE
111 C
112 C C(I) IS NOW THE SIGMA(I) OF THE TEXT
113 C
114 C COMPUTE POLYNOMIAL COEFFICIENTS
115 C
116  b(n) = (y(n) - y(nm1))/d(nm1) + d(nm1)*(c(nm1) + 2.*c(n))
117  DO 40 i = 1, nm1
118  b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.*c(i))
119  d(i) = (c(i+1) - c(i))/d(i)
120  c(i) = 3.*c(i)
121  40 CONTINUE
122  c(n) = 3.*c(n)
123  d(n) = d(n-1)
124  RETURN
125 C
126  50 b(1) = (y(2)-y(1))/(x(2)-x(1))
127  c(1) = 0.
128  d(1) = 0.
129  b(2) = b(1)
130  c(2) = 0.
131  d(2) = 0.
132  RETURN
133  END
134 C******************** END FILE SPLINE.FOR ; GROUP TRKRLIB ******************