1 |
!! |
2 |
!! Copyright (c) 2006 The University of Notre Dame. All Rights Reserved. |
3 |
!! |
4 |
!! The University of Notre Dame grants you ("Licensee") a |
5 |
!! non-exclusive, royalty free, license to use, modify and |
6 |
!! redistribute this software in source and binary code form, provided |
7 |
!! that the following conditions are met: |
8 |
!! |
9 |
!! 1. Acknowledgement of the program authors must be made in any |
10 |
!! publication of scientific results based in part on use of the |
11 |
!! program. An acceptable form of acknowledgement is citation of |
12 |
!! the article in which the program was described (Matthew |
13 |
!! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher |
14 |
!! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented |
15 |
!! Parallel Simulation Engine for Molecular Dynamics," |
16 |
!! J. Comput. Chem. 26, pp. 252-271 (2005)) |
17 |
!! |
18 |
!! 2. Redistributions of source code must retain the above copyright |
19 |
!! notice, this list of conditions and the following disclaimer. |
20 |
!! |
21 |
!! 3. Redistributions in binary form must reproduce the above copyright |
22 |
!! notice, this list of conditions and the following disclaimer in the |
23 |
!! documentation and/or other materials provided with the |
24 |
!! distribution. |
25 |
!! |
26 |
!! This software is provided "AS IS," without a warranty of any |
27 |
!! kind. All express or implied conditions, representations and |
28 |
!! warranties, including any implied warranty of merchantability, |
29 |
!! fitness for a particular purpose or non-infringement, are hereby |
30 |
!! excluded. The University of Notre Dame and its licensors shall not |
31 |
!! be liable for any damages suffered by licensee as a result of |
32 |
!! using, modifying or distributing the software or its |
33 |
!! derivatives. In no event will the University of Notre Dame or its |
34 |
!! licensors be liable for any lost revenue, profit or data, or for |
35 |
!! direct, indirect, special, consequential, incidental or punitive |
36 |
!! damages, however caused and regardless of the theory of liability, |
37 |
!! arising out of the use of or inability to use software, even if the |
38 |
!! University of Notre Dame has been advised of the possibility of |
39 |
!! such damages. |
40 |
!! |
41 |
!! |
42 |
!! interpolation.F90 |
43 |
!! |
44 |
!! Created by Charles F. Vardeman II on 03 Apr 2006. |
45 |
!! |
46 |
!! PURPOSE: Generic Spline interpolation routines. These routines |
47 |
!! assume that we are on a uniform grid for precomputation of |
48 |
!! spline parameters. |
49 |
!! |
50 |
!! @author Charles F. Vardeman II |
51 |
!! @version $Id: interpolation.F90,v 1.4 2006-04-14 21:49:54 gezelter Exp $ |
52 |
|
53 |
|
54 |
module INTERPOLATION |
55 |
use definitions |
56 |
use status |
57 |
implicit none |
58 |
PRIVATE |
59 |
|
60 |
character(len = statusMsgSize) :: errMSG |
61 |
|
62 |
type, public :: cubicSpline |
63 |
private |
64 |
logical :: isUniform = .false. |
65 |
integer :: np = 0 |
66 |
real(kind=dp) :: dx_i |
67 |
real (kind=dp), pointer,dimension(:) :: x => null() |
68 |
real (kind=dp), pointer,dimension(:,:) :: c => null() |
69 |
end type cubicSpline |
70 |
|
71 |
public :: newSpline |
72 |
public :: deleteSpline |
73 |
public :: lookup_spline |
74 |
public :: lookup_uniform_spline |
75 |
public :: lookup_nonuniform_spline |
76 |
|
77 |
contains |
78 |
|
79 |
|
80 |
subroutine newSpline(cs, x, y, yp1, ypn, isUniform) |
81 |
|
82 |
!************************************************************************ |
83 |
! |
84 |
! newSpline solves for slopes defining a cubic spline. |
85 |
! |
86 |
! Discussion: |
87 |
! |
88 |
! A tridiagonal linear system for the unknown slopes S(I) of |
89 |
! F at x(I), I=1,..., N, is generated and then solved by Gauss |
90 |
! elimination, with S(I) ending up in cs%C(2,I), for all I. |
91 |
! |
92 |
! Reference: |
93 |
! |
94 |
! Carl DeBoor, |
95 |
! A Practical Guide to Splines, |
96 |
! Springer Verlag. |
97 |
! |
98 |
! Parameters: |
99 |
! |
100 |
! Input, real x(N), the abscissas or X values of |
101 |
! the data points. The entries of x are assumed to be |
102 |
! strictly increasing. |
103 |
! |
104 |
! Input, real y(I), contains the function value at x(I) for |
105 |
! I = 1, N. |
106 |
! |
107 |
! Input, real yp1 contains the slope at x(1) |
108 |
! Input, real ypn contains the slope at x(N) |
109 |
! |
110 |
! On output, the slopes at x(I) have been stored in |
111 |
! cs%C(2,I), for I = 1 to N. |
112 |
|
113 |
implicit none |
114 |
|
115 |
type (cubicSpline), intent(inout) :: cs |
116 |
real( kind = DP ), intent(in) :: x(:), y(:) |
117 |
real( kind = DP ), intent(in) :: yp1, ypn |
118 |
logical, intent(in) :: isUniform |
119 |
real( kind = DP ) :: g, divdif1, divdif3, dx |
120 |
integer :: i, alloc_error, np |
121 |
|
122 |
alloc_error = 0 |
123 |
|
124 |
if (cs%np .ne. 0) then |
125 |
call handleWarning("interpolation::newSpline", & |
126 |
"cubicSpline struct was already created") |
127 |
call deleteSpline(cs) |
128 |
end if |
129 |
|
130 |
! make sure the sizes match |
131 |
|
132 |
np = size(x) |
133 |
|
134 |
if ( size(y) .ne. np ) then |
135 |
call handleError("interpolation::newSpline", & |
136 |
"Array size mismatch") |
137 |
end if |
138 |
|
139 |
cs%np = np |
140 |
cs%isUniform = isUniform |
141 |
|
142 |
allocate(cs%x(np), stat=alloc_error) |
143 |
if(alloc_error .ne. 0) then |
144 |
call handleError("interpolation::newSpline", & |
145 |
"Error in allocating storage for x") |
146 |
endif |
147 |
|
148 |
allocate(cs%c(4,np), stat=alloc_error) |
149 |
if(alloc_error .ne. 0) then |
150 |
call handleError("interpolation::newSpline", & |
151 |
"Error in allocating storage for c") |
152 |
endif |
153 |
|
154 |
do i = 1, np |
155 |
cs%x(i) = x(i) |
156 |
cs%c(1,i) = y(i) |
157 |
enddo |
158 |
|
159 |
! Set the first derivative of the function to the second coefficient of |
160 |
! each of the endpoints |
161 |
|
162 |
cs%c(2,1) = yp1 |
163 |
cs%c(2,np) = ypn |
164 |
|
165 |
! |
166 |
! Set up the right hand side of the linear system. |
167 |
! |
168 |
|
169 |
do i = 2, cs%np - 1 |
170 |
cs%c(2,i) = 3.0_DP * ( & |
171 |
(x(i) - x(i-1)) * (cs%c(1,i+1) - cs%c(1,i)) / (x(i+1) - x(i)) + & |
172 |
(x(i+1) - x(i)) * (cs%c(1,i) - cs%c(1,i-1)) / (x(i) - x(i-1))) |
173 |
end do |
174 |
|
175 |
! |
176 |
! Set the diagonal coefficients. |
177 |
! |
178 |
cs%c(4,1) = 1.0_DP |
179 |
do i = 2, cs%np - 1 |
180 |
cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) ) |
181 |
end do |
182 |
cs%c(4,cs%np) = 1.0_DP |
183 |
! |
184 |
! Set the off-diagonal coefficients. |
185 |
! |
186 |
cs%c(3,1) = 0.0_DP |
187 |
do i = 2, cs%np |
188 |
cs%c(3,i) = x(i) - x(i-1) |
189 |
end do |
190 |
! |
191 |
! Forward elimination. |
192 |
! |
193 |
do i = 2, cs%np - 1 |
194 |
g = -cs%c(3,i+1) / cs%c(4,i-1) |
195 |
cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1) |
196 |
cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1) |
197 |
end do |
198 |
! |
199 |
! Back substitution for the interior slopes. |
200 |
! |
201 |
do i = cs%np - 1, 2, -1 |
202 |
cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i) |
203 |
end do |
204 |
! |
205 |
! Now compute the quadratic and cubic coefficients used in the |
206 |
! piecewise polynomial representation. |
207 |
! |
208 |
do i = 1, cs%np - 1 |
209 |
dx = x(i+1) - x(i) |
210 |
divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx |
211 |
divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1 |
212 |
cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx |
213 |
cs%c(4,i) = divdif3 / ( dx * dx ) |
214 |
end do |
215 |
|
216 |
cs%c(3,cs%np) = 0.0_DP |
217 |
cs%c(4,cs%np) = 0.0_DP |
218 |
|
219 |
cs%dx_i = 1.0_DP / dx |
220 |
|
221 |
return |
222 |
end subroutine newSplineWithoutDerivs |
223 |
|
224 |
subroutine deleteSpline(this) |
225 |
|
226 |
type(cubicSpline) :: this |
227 |
|
228 |
if(associated(this%x)) then |
229 |
deallocate(this%x) |
230 |
this%x => null() |
231 |
end if |
232 |
if(associated(this%c)) then |
233 |
deallocate(this%c) |
234 |
this%c => null() |
235 |
end if |
236 |
|
237 |
this%np = 0 |
238 |
|
239 |
end subroutine deleteSpline |
240 |
|
241 |
subroutine lookup_nonuniform_spline(cs, xval, yval) |
242 |
|
243 |
!************************************************************************* |
244 |
! |
245 |
! lookup_nonuniform_spline evaluates a piecewise cubic Hermite interpolant. |
246 |
! |
247 |
! Discussion: |
248 |
! |
249 |
! newSpline must be called first, to set up the |
250 |
! spline data from the raw function and derivative data. |
251 |
! |
252 |
! Modified: |
253 |
! |
254 |
! 06 April 1999 |
255 |
! |
256 |
! Reference: |
257 |
! |
258 |
! Conte and de Boor, |
259 |
! Algorithm PCUBIC, |
260 |
! Elementary Numerical Analysis, |
261 |
! 1973, page 234. |
262 |
! |
263 |
! Parameters: |
264 |
! |
265 |
implicit none |
266 |
|
267 |
type (cubicSpline), intent(in) :: cs |
268 |
real( kind = DP ), intent(in) :: xval |
269 |
real( kind = DP ), intent(out) :: yval |
270 |
real( kind = DP ) :: dx |
271 |
integer :: i, j |
272 |
! |
273 |
! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains |
274 |
! or is nearest to xval. |
275 |
! |
276 |
j = cs%np - 1 |
277 |
|
278 |
do i = 1, cs%np - 2 |
279 |
|
280 |
if ( xval < cs%x(i+1) ) then |
281 |
j = i |
282 |
exit |
283 |
end if |
284 |
|
285 |
end do |
286 |
! |
287 |
! Evaluate the cubic polynomial. |
288 |
! |
289 |
dx = xval - cs%x(j) |
290 |
|
291 |
yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) ) |
292 |
|
293 |
return |
294 |
end subroutine lookup_nonuniform_spline |
295 |
|
296 |
subroutine lookup_uniform_spline(cs, xval, yval) |
297 |
|
298 |
!************************************************************************* |
299 |
! |
300 |
! lookup_uniform_spline evaluates a piecewise cubic Hermite interpolant. |
301 |
! |
302 |
! Discussion: |
303 |
! |
304 |
! newSpline must be called first, to set up the |
305 |
! spline data from the raw function and derivative data. |
306 |
! |
307 |
! Modified: |
308 |
! |
309 |
! 06 April 1999 |
310 |
! |
311 |
! Reference: |
312 |
! |
313 |
! Conte and de Boor, |
314 |
! Algorithm PCUBIC, |
315 |
! Elementary Numerical Analysis, |
316 |
! 1973, page 234. |
317 |
! |
318 |
! Parameters: |
319 |
! |
320 |
implicit none |
321 |
|
322 |
type (cubicSpline), intent(in) :: cs |
323 |
real( kind = DP ), intent(in) :: xval |
324 |
real( kind = DP ), intent(out) :: yval |
325 |
real( kind = DP ) :: dx |
326 |
integer :: i, j |
327 |
! |
328 |
! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains |
329 |
! or is nearest to xval. |
330 |
|
331 |
j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1)) |
332 |
|
333 |
dx = xval - cs%x(j) |
334 |
|
335 |
yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) ) |
336 |
|
337 |
return |
338 |
end subroutine lookup_uniform_spline |
339 |
|
340 |
subroutine lookup_spline(cs, xval, yval) |
341 |
|
342 |
type (cubicSpline), intent(in) :: cs |
343 |
real( kind = DP ), intent(inout) :: xval |
344 |
real( kind = DP ), intent(inout) :: yval |
345 |
|
346 |
if (cs%isUniform) then |
347 |
call lookup_uniform_spline(cs, xval, yval) |
348 |
else |
349 |
call lookup_nonuniform_spline(cs, xval, yval) |
350 |
endif |
351 |
|
352 |
return |
353 |
end subroutine lookup_spline |
354 |
|
355 |
end module INTERPOLATION |