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 interplelation routines. These routines assume that we are on a uniform grid for |
47 |
!! precomputation of spline parameters. |
48 |
!! |
49 |
!! @author Charles F. Vardeman II |
50 |
!! @version $Id: interpolation.F90,v 1.1 2006-04-14 19:57:04 gezelter Exp $ |
51 |
|
52 |
|
53 |
module INTERPOLATION |
54 |
use definitions |
55 |
use status |
56 |
implicit none |
57 |
PRIVATE |
58 |
|
59 |
character(len = statusMsgSize) :: errMSG |
60 |
|
61 |
type, public :: cubicSpline |
62 |
private |
63 |
integer :: np = 0 |
64 |
real(kind=dp) :: dx |
65 |
real(kind=dp) :: dx_i |
66 |
real (kind=dp), pointer,dimension(:) :: x => null() |
67 |
real (kind=dp), pointer,dimension(4,:) :: c => null() |
68 |
end type cubicSpline |
69 |
|
70 |
interface splineLookup |
71 |
module procedure multiSplint |
72 |
module procedure splintd |
73 |
module procedure splintd1 |
74 |
module procedure splintd2 |
75 |
end interface |
76 |
|
77 |
interface newSpline |
78 |
module procedure newSplineWithoutDerivs |
79 |
module procedure newSplineWithDerivs |
80 |
end interface |
81 |
|
82 |
public :: deleteSpline |
83 |
|
84 |
contains |
85 |
|
86 |
|
87 |
subroutine newSplineWithoutDerivs(cs, x, y, yp1, ypn, boundary) |
88 |
|
89 |
!************************************************************************ |
90 |
! |
91 |
! newSplineWithoutDerivs solves for slopes defining a cubic spline. |
92 |
! |
93 |
! Discussion: |
94 |
! |
95 |
! A tridiagonal linear system for the unknown slopes S(I) of |
96 |
! F at x(I), I=1,..., N, is generated and then solved by Gauss |
97 |
! elimination, with S(I) ending up in cs%C(2,I), for all I. |
98 |
! |
99 |
! Reference: |
100 |
! |
101 |
! Carl DeBoor, |
102 |
! A Practical Guide to Splines, |
103 |
! Springer Verlag. |
104 |
! |
105 |
! Parameters: |
106 |
! |
107 |
! Input, real x(N), the abscissas or X values of |
108 |
! the data points. The entries of TAU are assumed to be |
109 |
! strictly increasing. |
110 |
! |
111 |
! Input, real y(I), contains the function value at x(I) for |
112 |
! I = 1, N. |
113 |
! |
114 |
! yp1 contains the slope at x(1) and ypn contains |
115 |
! the slope at x(N). |
116 |
! |
117 |
! On output, the intermediate slopes at x(I) have been |
118 |
! stored in cs%C(2,I), for I = 2 to N-1. |
119 |
|
120 |
implicit none |
121 |
|
122 |
type (cubicSpline), intent(inout) :: cs |
123 |
real( kind = DP ), intent(in) :: x(:), y(:) |
124 |
real( kind = DP ), intent(in) :: yp1, ypn |
125 |
character(len=*), intent(in) :: boundary |
126 |
real( kind = DP ) :: g, divdif1, divdif3, dx |
127 |
integer :: i, alloc_error, np |
128 |
|
129 |
alloc_error = 0 |
130 |
|
131 |
if (cs%np .ne. 0) then |
132 |
call handleWarning("interpolation::newSplineWithoutDerivs", & |
133 |
"Type was already created") |
134 |
call deleteSpline(cs) |
135 |
end if |
136 |
|
137 |
! make sure the sizes match |
138 |
|
139 |
if (size(x) .ne. size(y)) then |
140 |
call handleError("interpolation::newSplineWithoutDerivs", & |
141 |
"Array size mismatch") |
142 |
end if |
143 |
|
144 |
np = size(x) |
145 |
cs%np = np |
146 |
|
147 |
allocate(cs%x(np), stat=alloc_error) |
148 |
if(alloc_error .ne. 0) then |
149 |
call handleError("interpolation::newSplineWithoutDerivs", & |
150 |
"Error in allocating storage for x") |
151 |
endif |
152 |
|
153 |
allocate(cs%c(4,np), stat=alloc_error) |
154 |
if(alloc_error .ne. 0) then |
155 |
call handleError("interpolation::newSplineWithoutDerivs", & |
156 |
"Error in allocating storage for c") |
157 |
endif |
158 |
|
159 |
do i = 1, np |
160 |
cs%x(i) = x(i) |
161 |
cs%c(1,i) = y(i) |
162 |
enddo |
163 |
|
164 |
if ((boundary.eq.'l').or.(boundary.eq.'L').or. & |
165 |
(boundary.eq.'b').or.(boundary.eq.'B')) then |
166 |
cs%c(2,1) = yp1 |
167 |
else |
168 |
cs%c(2,1) = 0.0_DP |
169 |
endif |
170 |
if ((boundary.eq.'u').or.(boundary.eq.'U').or. & |
171 |
(boundary.eq.'b').or.(boundary.eq.'B')) then |
172 |
cs%c(2,1) = ypn |
173 |
else |
174 |
cs%c(2,1) = 0.0_DP |
175 |
endif |
176 |
|
177 |
! |
178 |
! Set up the right hand side of the linear system. |
179 |
! |
180 |
do i = 2, cs%np - 1 |
181 |
cs%c(2,i) = 3.0_DP * ( & |
182 |
(x(i) - x(i-1)) * (cs%c(1,i+1) - cs%c(1,i)) / (x(i+1) - x(i)) + & |
183 |
(x(i+1) - x(i)) * (cs%c(1,i) - cs%c(1,i-1)) / (x(i) - x(i-1))) |
184 |
end do |
185 |
! |
186 |
! Set the diagonal coefficients. |
187 |
! |
188 |
cs%c(4,1) = 1.0_DP |
189 |
do i = 2, cs%np - 1 |
190 |
cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) ) |
191 |
end do |
192 |
cs%c(4,n) = 1.0_DP |
193 |
! |
194 |
! Set the off-diagonal coefficients. |
195 |
! |
196 |
cs%c(3,1) = 0.0_DP |
197 |
do i = 2, cs%np |
198 |
cs%c(3,i) = x(i) - x(i-1) |
199 |
end do |
200 |
! |
201 |
! Forward elimination. |
202 |
! |
203 |
do i = 2, cs%np - 1 |
204 |
g = -cs%c(3,i+1) / cs%c(4,i-1) |
205 |
cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1) |
206 |
cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1) |
207 |
end do |
208 |
! |
209 |
! Back substitution for the interior slopes. |
210 |
! |
211 |
do i = cs%np - 1, 2, -1 |
212 |
cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i) |
213 |
end do |
214 |
! |
215 |
! Now compute the quadratic and cubic coefficients used in the |
216 |
! piecewise polynomial representation. |
217 |
! |
218 |
do i = 1, cs%np - 1 |
219 |
dx = x(i+1) - x(i) |
220 |
divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx |
221 |
divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1 |
222 |
cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx |
223 |
cs%c(4,i) = divdif3 / ( dx * dx ) |
224 |
end do |
225 |
|
226 |
cs%c(3,np) = 0.0_DP |
227 |
cs%c(4,np) = 0.0_DP |
228 |
|
229 |
cs%dx = dx |
230 |
cs%dxi = 1.0_DP / dx |
231 |
return |
232 |
end subroutine newSplineWithoutDerivs |
233 |
|
234 |
subroutine newSplineWithDerivs(cs, x, y, yp) |
235 |
|
236 |
!************************************************************************ |
237 |
! |
238 |
! newSplineWithDerivs |
239 |
|
240 |
implicit none |
241 |
|
242 |
type (cubicSpline), intent(inout) :: cs |
243 |
real( kind = DP ), intent(in) :: x(:), y(:), yp(:) |
244 |
real( kind = DP ) :: g, divdif1, divdif3, dx |
245 |
integer :: i, alloc_error, np |
246 |
|
247 |
alloc_error = 0 |
248 |
|
249 |
if (cs%np .ne. 0) then |
250 |
call handleWarning("interpolation::newSplineWithDerivs", & |
251 |
"Type was already created") |
252 |
call deleteSpline(cs) |
253 |
end if |
254 |
|
255 |
! make sure the sizes match |
256 |
|
257 |
if ((size(x) .ne. size(y)).or.(size(x) .ne. size(yp))) then |
258 |
call handleError("interpolation::newSplineWithDerivs", & |
259 |
"Array size mismatch") |
260 |
end if |
261 |
|
262 |
np = size(x) |
263 |
cs%np = np |
264 |
|
265 |
allocate(cs%x(np), stat=alloc_error) |
266 |
if(alloc_error .ne. 0) then |
267 |
call handleError("interpolation::newSplineWithDerivs", & |
268 |
"Error in allocating storage for x") |
269 |
endif |
270 |
|
271 |
allocate(cs%c(4,np), stat=alloc_error) |
272 |
if(alloc_error .ne. 0) then |
273 |
call handleError("interpolation::newSplineWithDerivs", & |
274 |
"Error in allocating storage for c") |
275 |
endif |
276 |
|
277 |
do i = 1, np |
278 |
cs%x(i) = x(i) |
279 |
cs%c(1,i) = y(i) |
280 |
cs%c(2,i) = yp(i) |
281 |
enddo |
282 |
! |
283 |
! Set the diagonal coefficients. |
284 |
! |
285 |
cs%c(4,1) = 1.0_DP |
286 |
do i = 2, cs%np - 1 |
287 |
cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) ) |
288 |
end do |
289 |
cs%c(4,n) = 1.0_DP |
290 |
! |
291 |
! Set the off-diagonal coefficients. |
292 |
! |
293 |
cs%c(3,1) = 0.0_DP |
294 |
do i = 2, cs%np |
295 |
cs%c(3,i) = x(i) - x(i-1) |
296 |
end do |
297 |
! |
298 |
! Forward elimination. |
299 |
! |
300 |
do i = 2, cs%np - 1 |
301 |
g = -cs%c(3,i+1) / cs%c(4,i-1) |
302 |
cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1) |
303 |
cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1) |
304 |
end do |
305 |
! |
306 |
! Back substitution for the interior slopes. |
307 |
! |
308 |
do i = cs%np - 1, 2, -1 |
309 |
cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i) |
310 |
end do |
311 |
! |
312 |
! Now compute the quadratic and cubic coefficients used in the |
313 |
! piecewise polynomial representation. |
314 |
! |
315 |
do i = 1, cs%np - 1 |
316 |
dx = x(i+1) - x(i) |
317 |
divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx |
318 |
divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1 |
319 |
cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx |
320 |
cs%c(4,i) = divdif3 / ( dx * dx ) |
321 |
end do |
322 |
|
323 |
cs%c(3,np) = 0.0_DP |
324 |
cs%c(4,np) = 0.0_DP |
325 |
|
326 |
cs%dx = dx |
327 |
cs%dxi = 1.0_DP / dx |
328 |
|
329 |
return |
330 |
end subroutine newSplineWithoutDerivs |
331 |
|
332 |
subroutine deleteSpline(this) |
333 |
|
334 |
type(cubicSpline) :: this |
335 |
|
336 |
if(associated(this%x)) then |
337 |
deallocate(this%x) |
338 |
this%x => null() |
339 |
end if |
340 |
if(associated(this%c)) then |
341 |
deallocate(this%c) |
342 |
this%c => null() |
343 |
end if |
344 |
|
345 |
this%np = 0 |
346 |
|
347 |
end subroutine deleteSpline |
348 |
|
349 |
subroutine lookup_nonuniform_spline(cs, xval, yval) |
350 |
|
351 |
!************************************************************************* |
352 |
! |
353 |
! lookup_nonuniform_spline evaluates a piecewise cubic Hermite interpolant. |
354 |
! |
355 |
! Discussion: |
356 |
! |
357 |
! newSpline must be called first, to set up the |
358 |
! spline data from the raw function and derivative data. |
359 |
! |
360 |
! Modified: |
361 |
! |
362 |
! 06 April 1999 |
363 |
! |
364 |
! Reference: |
365 |
! |
366 |
! Conte and de Boor, |
367 |
! Algorithm PCUBIC, |
368 |
! Elementary Numerical Analysis, |
369 |
! 1973, page 234. |
370 |
! |
371 |
! Parameters: |
372 |
! |
373 |
implicit none |
374 |
|
375 |
type (cubicSpline), intent(in) :: cs |
376 |
real( kind = DP ), intent(in) :: xval |
377 |
real( kind = DP ), intent(out) :: yval |
378 |
integer :: i, j |
379 |
! |
380 |
! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains |
381 |
! or is nearest to xval. |
382 |
! |
383 |
j = cs%np - 1 |
384 |
|
385 |
do i = 1, cs%np - 2 |
386 |
|
387 |
if ( xval < cs%x(i+1) ) then |
388 |
j = i |
389 |
exit |
390 |
end if |
391 |
|
392 |
end do |
393 |
! |
394 |
! Evaluate the cubic polynomial. |
395 |
! |
396 |
dx = xval - cs%x(j) |
397 |
|
398 |
yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) ) |
399 |
|
400 |
return |
401 |
end subroutine lookup_nonuniform_spline |
402 |
|
403 |
subroutine lookup_uniform_spline(cs, xval, yval) |
404 |
|
405 |
!************************************************************************* |
406 |
! |
407 |
! lookup_uniform_spline evaluates a piecewise cubic Hermite interpolant. |
408 |
! |
409 |
! Discussion: |
410 |
! |
411 |
! newSpline must be called first, to set up the |
412 |
! spline data from the raw function and derivative data. |
413 |
! |
414 |
! Modified: |
415 |
! |
416 |
! 06 April 1999 |
417 |
! |
418 |
! Reference: |
419 |
! |
420 |
! Conte and de Boor, |
421 |
! Algorithm PCUBIC, |
422 |
! Elementary Numerical Analysis, |
423 |
! 1973, page 234. |
424 |
! |
425 |
! Parameters: |
426 |
! |
427 |
implicit none |
428 |
|
429 |
type (cubicSpline), intent(in) :: cs |
430 |
real( kind = DP ), intent(in) :: xval |
431 |
real( kind = DP ), intent(out) :: yval |
432 |
integer :: i, j |
433 |
! |
434 |
! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains |
435 |
! or is nearest to xval. |
436 |
|
437 |
j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dxi) + 1)) |
438 |
|
439 |
dx = xval - cs%x(j) |
440 |
|
441 |
yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) ) |
442 |
|
443 |
return |
444 |
end subroutine lookup_uniform_spline |
445 |
|
446 |
end module INTERPOLATION |