1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
|
subroutine tucspl (xa, ya, n, work, y2)
C
C Compute the second derivative of YA at each point in the array.
C This is the initialization needed in preparation for interpolating
C using cubic splines by the subroutine TUISPL. Input and output
C are all double precision.
C
C This routine was copied with slight modifications from the SPLINE
C subroutine in Numerical Recipes by Press, Flannery, Teukolsky and
C Vetterling.
C
C N i: number of elements in each array
C XA i: array of independent-variable values
C YA i: array of dependent-variable values
C WORK io: scratch array used for work space
C Y2 o: second derivative of YA at each point
C
CH Phil Hodge, 14-Apr-1988 Subroutine copied from Numerical Recipes SPLINE.
C
integer n
double precision xa(n), ya(n), work(n), y2(n)
C--
integer i
double precision p, sig
C These values (and y2(n) = 0) are for a "natural" spline.
y2(1) = 0.
work(1) = 0.
C
C This is the decomposition loop of the tridiagonal algorithm.
C Y2 and WORK are used for temporary storage of the decomposed factors.
C
do 10 i = 2, n-1
sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
p = sig * y2(i-1) + 2.
y2(i) = (sig - 1.) / p
work(i) = (6. * ((ya(i+1) - ya(i)) / (xa(i+1) - xa(i))
+ - (ya(i) - ya(i-1)) / (xa(i) - xa(i-1))) /
+ (xa(i+1) - xa(i-1)) - sig * work(i-1)) / p
10 continue
C "natural" spline
y2(n) = 0.
C
C This is the backsubstitution loop of the tridiagonal algorithm.
C
do 20 i = n-1, 1, -1
y2(i) = y2(i) * y2(i+1) + work(i)
20 continue
return
end
|
