1
-
!! Provides spline functions
4
-
!! Cubic spline interpolation functions and procedures
11
10
SUBROUTINE CSPLINE (TAU, C, N, IBCBEG, IBCEND, NDIM)
13
-
! TAKEN FROM "A PRACTICAL GUIDE TO SPLINES", BY CARL DE BOOR. 1978.
14
-
! SPRINGER-VERLAG. THE INPUT PARAMETER "NDIM" HAS BEEN ADDED TO
15
-
! ALLOW FOR MULTIPLE CALLS WITH DIFFERENT VALUES OF N. - DENNIS DUNDORE
17
-
! SUBSTANTIAL MODIFICATIONS MADE BY STEVE WALES, APRIL 1983,
18
-
! PRINCIPALLY TO HANDLE COMPLEX NUMBERS (C) & UPDATE THE FORTRAN.
20
-
! ***************************** I N P U T ****************************
22
-
! N = NUMBER OF DATA POINTS. ASSUMED TO BE .GE. 2.
24
-
! (TAU(I), C(1,I),I=1,...,N) = ABSCISSAE AND ORDINATES OF THE DATA
25
-
! POINTS. TAU IS ASSUMED TO BE STRICTLY INCREASING.
27
-
! IBCBEG, IBCEND = BOUNDARY CONDITION INDICATORS, AND
28
-
! C(2,1), C(2,N) = BOUNDARY CONDITION INFORMATION. SPECIFICALLY,
29
-
! IBCBEG = 0 IMPLIES NO BOUNDARY CONDITION AT TAU(1) IS GIVEN.
30
-
! IN THIS CASE, THE "NOT-A-KNOT" CONDITION IS USED, IE THE
31
-
! JUMP IN THE 3-RD DERIVATIVE ACROSS TAU(2) IS FORCED TO 0.,
32
-
! THUS THE 1-ST AND 2-ND POLYNOMIAL PIECES ARE MADE TO
34
-
! IBCBEG = 1 IMPLIES THAT THE SLOPE AT TAU(1) IS SET EQUAL TO C(2,1)
36
-
! IBCBEG = 2 IMPLIES THAT THE 2-ND DERIVATIVE AT TAU(1) IS SET EQUAL
37
-
! TO C(2,1), SUPPLIED BY INPUT.
38
-
! IBCEND = 0, 1, OR 2 HAS ANALOGOUS MEANING CONCERNING THE BOUNDARY
39
-
! CONDITION AT TAU(N), WITH INFORMATION SUPPLIED BY C(2,N).
41
-
! NDIM = ROW DIMENSION OF C MATRIX: C(4,NDIM)
43
-
! ************************** O U T P U T ****************************
45
-
! C(J,I), J=1,...,4; I=1,...,L=N-1 = THE POLY COEFFS OF THE CUBIC
46
-
! SPLINE WITH INTERIOR KNOTS TAU(2),...,TAU(N-1). PRECISELY, IN THE
47
-
! INTERVAL (TAU(I), TAU(I+1)), THE SPLINE F IS GIVEN BY
49
-
! F(X) = C(1,I) + H*(C(2,I) + H*(C(3,I)/2. + H*C(4,I)/6.))
51
-
! WHERE H = X - TAU(I).
53
-
! THE COEFFICIENTS CALCULATED ARE, 1) THE VALUE, 2) THE SLOPE, AND
54
-
! 3) THE CURVATURE AT EACH OF THE KNOTS 1 TO N-1, AND 4) THE RATE OF
55
-
! CHANGE OF THE CURVATURE OVER THE FOLLOWING INTERVAL. IN ADDITION,
56
-
! WE HAVE THE VALUE AND THE SLOPE AT THE LAST POINT. THE LAST TWO
57
-
! POSTIONS AT THE LAST POINT ARE THEN SET TO THE CURVATURE AT THAT
58
-
! POINT (IN C(3,N)) AND THE MEAN VALUE OVER THE ENTIRE INTERVAL,
59
-
! CALCULATED AS THE INTEGRAL OVER THE INTERVAL DIVIDED BY THE LENGTH
60
-
! OF THE INTERVAL (IN C(4,N)).
62
-
! **********************************************************************
64
-
INTEGER, INTENT(IN) :: N, IBCBEG, IBCEND, NDIM
65
-
REAL (KIND=8), INTENT(IN) :: TAU(N)
66
-
COMPLEX (KIND=8), INTENT(INOUT) :: C(4,NDIM)
69
-
INTEGER :: L, M, I, J
70
-
COMPLEX (KIND=8) :: G, DTAU, DIVDF1, DIVDF3
75
13*
C(3,M) = TAU(M) - TAU(M-1)
76
23*
C(4,M) = (C(1,M) - C(1,M-1)) / C(3,M)
79
-
! * BEGINNING BOUNDARY CONDITION SECTION *
81
10
IF (IBCBEG==0) THEN ! IBCBEG = 0
82
10
IF (N>2) THEN ! N > 2
84
1*
C(3,1) = C(3,2) + C(3,3)
85
2*
C(2,1) = ((C(3,2) + 2.0*C(3,1))*C(4,2)*C(3,3) + &
86
2*
C(3,2)**2 * C(4,3)) / C(3,1)
90
9*
C(2,1) = 2.0 * C(4,2)
93
#####
ELSE IF (IBCBEG==1) THEN ! IBCBEG = 1
94
#####
C(4,1) = (1.0,0.0)
95
#####
C(3,1) = (0.0,0.0)
97
#####
ELSE IF (IBCBEG==2) THEN ! IBCBEG = 2
98
#####
C(4,1) = (2.0,0.0)
99
#####
C(3,1) = (1.0,0.0)
100
#####
C(2,1) = 3.0*C(4,2) - C(2,1)*C(3,2)/2.0
103
-
! * RUNNING CALCULATIONS TO N-1 - LOOP IS NOT EXECUTED IF N = 2 *
106
3*
G = -C(3,M+1) / C(4,M-1)
107
3*
C(2,M) = G*C(2,M-1) + 3.0*(C(3,M)*C(4,M+1) + C(3,M+1)*C(4,M))
108
13*
C(4,M) = G*C(3,M-1) + 2.0*(C(3,M) + C(3,M+1))
111
-
! * ENDING BOUNDARY CONDITION SECTION *
113
10
IF (IBCEND /= 1) THEN
114
10
IF (IBCEND==0) THEN
115
10
IF (N==2 .AND. IBCBEG==0) THEN
117
1
ELSE IF ((N==3 .AND. IBCBEG==0) .OR. N==2) THEN
118
#####
C(2,N) = 2.0 * C(4,N)
119
#####
C(4,N) = (1.,0.)
120
#####
G = -1.0 / C(4,N-1)
122
1*
G = C(3,N-1) + C(3,N)
123
4*
C(2,N) = ((C(3,N) + 2.0*G) * C(4,N)*C(3,N-1) + &
124
4*
C(3,N)**2 * (C(1,N-1)-C(1,N-2)) / C(3,N-1)) / G
129
#####
ELSE IF (IBCEND==2) THEN
130
#####
C(2,N) = 3.0 * C(4,N) + C(2,N)*C(3,N)/2.0
131
#####
C(4,N) = (2.0,0.0)
132
#####
G = -1.0 / C(4,N-1)
135
10
IF ( IBCBEG>0 .OR. N>2) THEN
136
1*
C(4,N) = G*C(3,N-1) + C(4,N)
137
1*
C(2,N) = (G*C(2,N-1) + C(2,N)) / C(4,N)
141
-
! * RUN THE ENDING BOUNDARY EFFECT BACK THROUGH THE EQUATIONS *
144
23*
C(2,J) = (C(2,J) - C(3,J)*C(2,J+1)) / C(4,J)
147
-
! * FINAL CALCULATIONS *
151
13*
DIVDF1 = (C(1,I)-C(1,I-1)) / DTAU
152
13*
DIVDF3 = C(2,I-1) + C(2,I) - 2.0*DIVDF1
153
13*
C(3,I-1) = 2.0 * (DIVDF1-C(2,I-1)-DIVDF3) / DTAU
154
23*
C(4,I-1) = (DIVDF3/DTAU) * (6.0/DTAU)
157
-
! * ADD THE CURVATURE AT THE LAST POINT IN THE THIRD POSITION OF THE
160
10*
C(3,N) = C(3,L) + (TAU(N)-TAU(L)) * C(4,L)
163
-
! * ADD THE MEAN VALUE OF THE ENTIRE INTERVAL IN THE FOURTH POSITION OF
164
-
! THE LAST NODE. MEAN VALUE IS CALCULATED AS THE INTEGRAL OVER THE
165
-
! INTERVAL DIVIDED BY THE LENGTH OF THE INTERVAL. *
167
10*
C(4,N) = (0.0,0.0)
168
23*
DO I = 1,L ! INTEGRATE OVER THE INTERVAL
169
13*
DTAU = TAU(I+1) - TAU(I)
170
52*
C(4,N) = C(4,N) + DTAU*(C(1,I) + DTAU*(C(2,I)/2.0 + &
171
75*
DTAU*(C(3,I)/6.0 + DTAU*C(4,I)/24.0)))
173
10*
C(4,N) = C(4,N) / (TAU(N)-TAU(1)) ! DIVIDE BY LENGTH OF INTERVAL
176
-
END SUBROUTINE CSPLINE
180
#####
SUBROUTINE VSPLINE (TAU, C, M, MDIM, F, N)
182
-
! VSPLINE CALCULATES THE CUBIC SPLINE VALUES FOR A SET OF N POINTS
183
-
! IN F FROM THE M-POINT CUBI! SPLINE FIT IN ! AND THE NODES IN TAU.
184
-
! THE POINTS ARE RETURNED IN F. ALL OF THE POINTS IN F MUST LIE
185
-
! BETWEEN TAU(1) AND TAU(M).
187
-
! * * * * * * * * * * * * * WARNINGS * * * * * * * * * * * * *
189
-
! POINTS OUTSIDE OF THE SPLINE FIT REGION ARE EXTRAPOLATED FROM THE END
190
-
! INTERVALS. THIS CAN RESULT IN WILD VALUES IF EXTRAPOLATED TOO FAR.
191
-
! ALSO THE POINTS MUST BE IN STRICTLY ASCENDING ORDER, IF NOT THE
192
-
! POINTS WHICH ARE OUT OF ORDER WILL BE EXTRAPOLATED FROM THE CURRENT
193
-
! INTERVAL AGAIN RESULTING IN WILD VALUES.
195
-
INTEGER, INTENT(IN) :: M, N, MDIM
196
-
REAL (KIND=8), INTENT(IN) :: TAU(M)
197
-
COMPLEX (KIND=8), INTENT(IN) :: C(4,MDIM)
198
-
COMPLEX (KIND=8), INTENT(INOUT) :: F(N)
201
-
INTEGER :: I, J, J1
207
#####
IF (TAU(J1) < REAL(F(I)) .AND. J1 < M) THEN ! CHECK TO MAKE SURE
208
#####
J = J + 1 ! THIS POINT IS NOT
209
#####
GO TO 10 ! IN THE NEXT INTERVAL.
211
#####
H = DBLE (F(I)) - TAU(J) ! DISTANCE FROM START OF INTERVAL
212
#####
F(I) = SPLINE (C(1,J), H)
215
-
END SUBROUTINE VSPLINE
218
-
! **********************************************************************C
219
#####
FUNCTION SPLINE ( C, H )
221
-
! THIS FUNCTION EVALUATES THE SPLINE AT THE POINT H
223
-
COMPLEX (KIND=8), INTENT(IN) :: C(4)
224
-
REAL (KIND=8), INTENT(IN) :: H
225
-
COMPLEX (KIND=8) :: SPLINE
227
#####
SPLINE = C(1) + H * ( C(2) + H * ( C(3) / 2.0 + H * C(4) / 6.0 ) )
229
-
END FUNCTION SPLINE
231
#####
FUNCTION SPLINEX ( C, H )
233
-
! THIS FUNCTION EVALUATES THE SPLINE DERIVATIVE AT THE POINT H
235
-
COMPLEX (KIND=8), INTENT(IN) :: C(4)
236
-
REAL (KIND=8), INTENT(IN) :: H
237
-
COMPLEX (KIND=8) :: SPLINEX
239
#####
SPLINEX = C(2) + H * ( C(3) + H * C(4) / 2.0 )
241
-
END FUNCTION SPLINEX
243
#####
FUNCTION SPLINEXX ( C, H )
245
-
! THIS FUNCTION EVALUATES THE SPLINE 2ND DERIVATIVE AT THE POINT H
247
-
COMPLEX (KIND=8), INTENT(IN) :: C(4)
248
-
REAL (KIND=8), INTENT(IN) :: H
249
-
COMPLEX (KIND=8) :: SPLINEXX
251
#####
SPLINEXX = C(3) + H * C(4)
253
-
END FUNCTION SPLINEXX
254
-
! **********************************************************************C
255
34117873
SUBROUTINE SPLINEALL ( C, H, F, FX, FXX )
257
-
! THIS ROUTINE EVALUATES THE
259
-
! SPLINE DERIVATIVE, AND
260
-
! SPLINE 2ND DERIVATIVE AT THE POINT H
262
-
COMPLEX (KIND=8), INTENT(IN) :: C(4)
263
-
REAL (KIND=8), INTENT(IN) :: H
264
-
COMPLEX (KIND=8), INTENT(OUT) :: F, FX, FXX
266
-
REAL (KIND=8), PARAMETER :: HALF = 0.5, SIXTH = 1.0 / 6.0
268
34117873
F = C(1) + H * ( C(2) + H * ( HALF * C(3) + SIXTH * H * C(4) ) )
269
34117873
FX = C(2) + H * ( C(3) + H * HALF * C(4) )
270
34117873
FXX = C(3) + H * C(4)
273
-
END SUBROUTINE SPLINEALL