3
! Spline interpolation, based on code in "Numerical Recipes"
6
integer, parameter :: dp = selected_real_kind(16,100)
7
integer, parameter :: sp = selected_real_kind(5,10)
10
public :: generate_spline
11
interface generate_spline
12
module procedure generate_spline_sp
13
module procedure generate_spline_dp
14
end interface generate_spline
16
public :: evaluate_spline
17
interface evaluate_spline
18
module procedure evaluate_spline_sp
19
module procedure evaluate_spline_dp
20
end interface evaluate_spline
23
private :: generate_spline_sp
24
private :: generate_spline_dp
25
private :: evaluate_spline_sp
26
private :: evaluate_spline_dp
28
CONTAINS !=================
30
SUBROUTINE generate_spline_sp(X,Y,N,Y2,YP1,YPN,stat)
32
! Given arrays x and y of size n, this routine computes
33
! an array y2 holding information about the second derivative
34
! of the interpolating spline. YP1 and YPN, if present, are the
35
! requested values of the first derivative of the spline at the
36
! end points. If not given, a "natural" spline with zero second
37
! derivative at the extremes is constructed.
38
! The array y2 is needed to evaluate the spline, but it only needs
39
! to be computed once.
41
integer, intent(in) :: n
42
real(kind=sp), dimension(:), intent(in) :: x, y
43
real(kind=sp), dimension(:), intent(out) :: y2
44
real(kind=sp), intent(in), optional :: yp1, ypn
45
integer, intent(out), optional :: stat
47
real(kind=sp) :: P,QN,SIG,UN
51
real(kind=sp), parameter :: zero = 0.0_sp, &
58
real(kind=sp), dimension(n) :: U ! Automatic array
60
!--------------------------------------------------------
61
! Check that the nodes are different
65
if (x(i) == x(i-1)) then
67
if (present(stat)) stat = flag
72
IF (present(YP1)) THEN
74
U(1) = (three/ (X(2)-X(1)))* ((Y(2)-Y(1))/ (X(2)-X(1))-YP1)
81
SIG = (X(I)-X(I-1))/ (X(I+1)-X(I-1))
84
U(I) = ( six * ((Y(I+1)-Y(I))/ (X(I+1)-X(I))- (Y(I)-Y(I-1))/ &
85
(X(I)-X(I-1)))/ (X(I+1)-X(I-1))-SIG*U(I-1))/P
87
IF (present(YPN)) THEN
89
UN = (three/(X(N)-X(N-1)))* (YPN - (Y(N)-Y(N-1))/ (X(N)-X(N-1)))
95
Y2(N) = (UN-QN*U(N-1))/ (QN*Y2(N-1) + one)
97
Y2(K) = Y2(K)*Y2(K+1) + U(K)
100
if (present(stat)) stat = flag
102
end subroutine generate_spline_sp
104
!..............................................
105
SUBROUTINE generate_spline_dp(X,Y,N,Y2,YP1,YPN,stat)
107
! Given arrays x and y of size n, this routine computes
108
! an array y2 holding information about the second derivative
109
! of the interpolating spline. YP1 and YPN, if present, are the
110
! requested values of the first derivative of the spline at the
111
! end points. If not given, a "natural" spline with zero second
112
! derivative at the extremes is constructed.
113
! The array y2 is needed to evaluate the spline, but it only needs
114
! to be computed once.
116
integer, intent(in) :: n
117
real(kind=dp), dimension(:), intent(in) :: x, y
118
real(kind=dp), dimension(:), intent(out) :: y2
119
real(kind=dp), intent(in), optional :: yp1, ypn
120
integer, intent(out), optional :: stat
122
real(kind=dp) :: P,QN,SIG,UN
126
real(kind=dp), parameter :: zero = 0.0_dp, &
133
real(kind=dp), dimension(n) :: U ! Automatic array
135
!--------------------------------------------------------
136
! Check that the nodes are different
140
if (x(i) == x(i-1)) then
142
if (present(stat)) stat = flag
147
IF (present(YP1)) THEN
149
U(1) = (three/ (X(2)-X(1)))* ((Y(2)-Y(1))/ (X(2)-X(1))-YP1)
156
SIG = (X(I)-X(I-1))/ (X(I+1)-X(I-1))
157
P = SIG*Y2(I-1) + two
159
U(I) = ( six * ((Y(I+1)-Y(I))/ (X(I+1)-X(I))- (Y(I)-Y(I-1))/ &
160
(X(I)-X(I-1)))/ (X(I+1)-X(I-1))-SIG*U(I-1))/P
162
IF (present(YPN)) THEN
164
UN = (three/(X(N)-X(N-1)))* (YPN - (Y(N)-Y(N-1))/ (X(N)-X(N-1)))
170
Y2(N) = (UN-QN*U(N-1))/ (QN*Y2(N-1) + one)
172
Y2(K) = Y2(K)*Y2(K+1) + U(K)
175
if (present(stat)) stat = flag
177
end subroutine generate_spline_dp
179
!..............................................
180
!------------------------------------------------
181
SUBROUTINE evaluate_spline_sp(XA,YA,Y2A,N,X,Y)
183
! Given arrays xa and ya of size n, and an array y2a computed
184
! previously by routine spline, this routine computes the value
185
! at x of the interpolating spline. The value is returned in y.
187
integer, intent(in) :: n
188
real(kind=sp), dimension(:), intent(in) :: xa, ya, y2a
189
real(kind=sp), intent(in) :: x
190
real(kind=sp), intent(out) :: y
193
real(kind=sp) :: A,B,H
199
IF (KHI-KLO <= 1) exit
208
H = XA(KHI) - XA(KLO)
212
Y = A*YA(KLO) + B*YA(KHI) + ((A**3-A)*Y2A(KLO)+ &
213
(B**3-B)*Y2A(KHI))* (H**2)/6.0_sp
215
end subroutine evaluate_spline_sp
216
!..............................................
217
SUBROUTINE evaluate_spline_dp(XA,YA,Y2A,N,X,Y)
219
! Given arrays xa and ya of size n, and an array y2a computed
220
! previously by routine spline, this routine computes the value
221
! at x of the interpolating spline. The value is returned in y.
223
integer, intent(in) :: n
224
real(kind=dp), dimension(:), intent(in) :: xa, ya, y2a
225
real(kind=dp), intent(in) :: x
226
real(kind=dp), intent(out) :: y
229
real(kind=dp) :: A,B,H
235
IF (KHI-KLO <= 1) exit
244
H = XA(KHI) - XA(KLO)
248
Y = A*YA(KLO) + B*YA(KHI) + ((A**3-A)*Y2A(KLO)+ &
249
(B**3-B)*Y2A(KHI))* (H**2)/6.0_dp
251
end subroutine evaluate_spline_dp
252
!..............................................
253
!------------------------------------------------
254
end module flib_spline
added
removed
Pseudo/atom/Util/Ps2Inp/flib_spline.f90
