V3FIT
r8nspline.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 c
11 c this is nspline -- natural cubic spline
12 c second derivative of interpolant set to zero at endpts
13 c
14 
15  SUBROUTINE r8nspline (N, X, Y, B, C, D)
16 !============
17 ! idecl: explicitize implicit INTEGER declarations:
18  IMPLICIT NONE
19  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
20  INTEGER ip1,ii
21 !============
22 ! idecl: explicitize implicit REAL declarations:
23  real*8 zbomb
24 !============
25  INTEGER N
26  real*8 x(n), y(n), b(n), c(n), d(n)
27 C
28 C THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED
29 C FOR A CUBIC INTERPOLATING SPLINE
30 C
31 C S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3
32 C
33 C FOR X(I) .LE. X .LE. X(I+1)
34 C
35 C INPUT..
36 C
37 C N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2)
38 C X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER
39 C Y = THE ORDINATES OF THE KNOTS
40 C
41 C OUTPUT..
42 C
43 C B, C, D = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE.
44 C
45 C USING P TO DENOTE DIFFERENTIATION,
46 C
47 C Y(I) = S(X(I))
48 C B(I) = SP(X(I))
49 C C(I) = SPP(X(I))/2
50 C D(I) = SPPP(X(I))/6 (DERIVATIVE FROM THE RIGHT)
51 C
52 CCCCCCCCCCCCCCC
53 C THE ACCOMPANYING FUNCTION SUBPROGRAM SEVAL CAN BE USED
54 C TO EVALUATE THE SPLINE.
55 C
56 C
57  INTEGER NM1, NM2, IB, I
58  real*8 t
59 C
60  nm1 = n-1
61  IF ( n .LT. 2 ) RETURN
62 C
63 C
64 C DMC - CHECK ORDINATES
65  DO 5 i=1,nm1
66  ip1=i+1
67  IF(x(ip1).LE.x(i)) then
68  write(6,9000)
69  write(6,9001) (x(ii),ii=1,n)
70  zbomb=sqrt(-1.0_r8/(x(1)*x(1)))
71  write(6,*) zbomb
72  write(6,*) ' nspline -- code stop.'
73  call exit(1)
74  9000 FORMAT(/
75  >' ? UNORDERED ORDINATE ARRAY PASSED TO SPLINE SUBROUTINE:')
76  9001 FORMAT(1x,5(1x,1pe12.5))
77  endif
78  5 CONTINUE
79 C
80  IF ( n .LT. 3 ) GO TO 50
81 C
82 C SET UP TRIDIAGONAL SYSTEM
83 C
84 C B = DIAGONAL, D = OFFDIAGONAL, C = RIGHT HAND SIDE.
85 C
86  d(1) = x(2) - x(1)
87  c(2) = (y(2) - y(1))/d(1)
88  DO 10 i = 2, nm1
89  d(i) = x(i+1) - x(i)
90  b(i) = 2._r8*(d(i-1) + d(i))
91  c(i+1) = (y(i+1) - y(i))/d(i)
92  c(i) = c(i+1) - c(i)
93  10 CONTINUE
94 C
95 C END CONDITIONS: 2nd deriv = 0 at endpts; solve reduced n-2 x n-2
96 C system.
97  c(1) = 0._r8
98  c(n) = 0._r8
99 C
100 C FORWARD ELIMINATION
101 C
102  15 DO 20 i = 3, nm1
103  t = d(i-1)/b(i-1)
104  b(i) = b(i) - t*d(i-1)
105  c(i) = c(i) - t*c(i-1)
106  20 CONTINUE
107 C
108 C BACK SUBSTITUTION
109 C
110  nm2=nm1-1
111  c(nm1) = c(nm1)/b(nm1)
112  DO 30 ib = 2, nm2
113  i = n-ib
114  c(i) = (c(i) - d(i)*c(i+1))/b(i)
115  30 CONTINUE
116 C
117 C C(I) IS NOW THE SIGMA(I) OF THE TEXT
118 C
119 C COMPUTE POLYNOMIAL COEFFICIENTS
120 C
121  b(n) = (y(n) - y(nm1))/d(nm1) + d(nm1)*(c(nm1) + 2._r8*c(n))
122  DO 40 i = 1, nm1
123  b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2._r8*c(i))
124  d(i) = (c(i+1) - c(i))/d(i)
125  c(i) = 3._r8*c(i)
126  40 CONTINUE
127  c(n) = 3._r8*c(n)
128  d(n) = d(n-1)
129  RETURN
130 C
131  50 b(1) = (y(2)-y(1))/(x(2)-x(1))
132  c(1) = 0._r8
133  d(1) = 0._r8
134  b(2) = b(1)
135  c(2) = 0._r8
136  d(2) = 0._r8
137  RETURN
138  END
139 C******************** END FILE SPLINE.FOR ; GROUP TRKRLIB ******************