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.3 2006-04-14 21:06:55 chrisfen 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(:,:) :: c => null() |
68 |
end type cubicSpline |
69 |
|
70 |
interface newSpline |
71 |
module procedure newSpline |
72 |
end interface |
73 |
|
74 |
public :: deleteSpline |
75 |
|
76 |
contains |
77 |
|
78 |
|
79 |
subroutine newSpline(cs, x, y, yp1, ypn) |
80 |
|
81 |
!************************************************************************ |
82 |
! |
83 |
! newSplineWithoutDerivs solves for slopes defining a cubic spline. |
84 |
! |
85 |
! Discussion: |
86 |
! |
87 |
! A tridiagonal linear system for the unknown slopes S(I) of |
88 |
! F at x(I), I=1,..., N, is generated and then solved by Gauss |
89 |
! elimination, with S(I) ending up in cs%C(2,I), for all I. |
90 |
! |
91 |
! Reference: |
92 |
! |
93 |
! Carl DeBoor, |
94 |
! A Practical Guide to Splines, |
95 |
! Springer Verlag. |
96 |
! |
97 |
! Parameters: |
98 |
! |
99 |
! Input, real x(N), the abscissas or X values of |
100 |
! the data points. The entries of x are assumed to be |
101 |
! strictly increasing. |
102 |
! |
103 |
! Input, real y(I), contains the function value at x(I) for |
104 |
! I = 1, N. |
105 |
! |
106 |
! yp1 contains the slope at x(1) and ypn contains |
107 |
! the slope at x(N). |
108 |
! |
109 |
! On output, the intermediate slopes at x(I) have been |
110 |
! stored in cs%C(2,I), for I = 2 to N-1. |
111 |
|
112 |
implicit none |
113 |
|
114 |
type (cubicSpline), intent(inout) :: cs |
115 |
real( kind = DP ), intent(in) :: x(:), y(:) |
116 |
real( kind = DP ), intent(in) :: yp1, ypn |
117 |
real( kind = DP ) :: g, divdif1, divdif3, dx |
118 |
integer :: i, alloc_error, np |
119 |
|
120 |
alloc_error = 0 |
121 |
|
122 |
if (cs%np .ne. 0) then |
123 |
call handleWarning("interpolation::newSplineWithoutDerivs", & |
124 |
"Type was already created") |
125 |
call deleteSpline(cs) |
126 |
end if |
127 |
|
128 |
! make sure the sizes match |
129 |
|
130 |
if (size(x) .ne. size(y)) then |
131 |
call handleError("interpolation::newSplineWithoutDerivs", & |
132 |
"Array size mismatch") |
133 |
end if |
134 |
|
135 |
np = size(x) |
136 |
cs%np = np |
137 |
|
138 |
allocate(cs%x(np), stat=alloc_error) |
139 |
if(alloc_error .ne. 0) then |
140 |
call handleError("interpolation::newSplineWithoutDerivs", & |
141 |
"Error in allocating storage for x") |
142 |
endif |
143 |
|
144 |
allocate(cs%c(4,np), stat=alloc_error) |
145 |
if(alloc_error .ne. 0) then |
146 |
call handleError("interpolation::newSplineWithoutDerivs", & |
147 |
"Error in allocating storage for c") |
148 |
endif |
149 |
|
150 |
do i = 1, np |
151 |
cs%x(i) = x(i) |
152 |
cs%c(1,i) = y(i) |
153 |
enddo |
154 |
|
155 |
! Set the first derivative of the function to the second coefficient of |
156 |
! each of the endpoints |
157 |
|
158 |
cs%c(2,1) = yp1 |
159 |
cs%c(2,np) = ypn |
160 |
|
161 |
|
162 |
! |
163 |
! Set up the right hand side of the linear system. |
164 |
! |
165 |
do i = 2, cs%np - 1 |
166 |
cs%c(2,i) = 3.0_DP * ( & |
167 |
(x(i) - x(i-1)) * (cs%c(1,i+1) - cs%c(1,i)) / (x(i+1) - x(i)) + & |
168 |
(x(i+1) - x(i)) * (cs%c(1,i) - cs%c(1,i-1)) / (x(i) - x(i-1))) |
169 |
end do |
170 |
! |
171 |
! Set the diagonal coefficients. |
172 |
! |
173 |
cs%c(4,1) = 1.0_DP |
174 |
do i = 2, cs%np - 1 |
175 |
cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) ) |
176 |
end do |
177 |
cs%c(4,cs%np) = 1.0_DP |
178 |
! |
179 |
! Set the off-diagonal coefficients. |
180 |
! |
181 |
cs%c(3,1) = 0.0_DP |
182 |
do i = 2, cs%np |
183 |
cs%c(3,i) = x(i) - x(i-1) |
184 |
end do |
185 |
! |
186 |
! Forward elimination. |
187 |
! |
188 |
do i = 2, cs%np - 1 |
189 |
g = -cs%c(3,i+1) / cs%c(4,i-1) |
190 |
cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1) |
191 |
cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1) |
192 |
end do |
193 |
! |
194 |
! Back substitution for the interior slopes. |
195 |
! |
196 |
do i = cs%np - 1, 2, -1 |
197 |
cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i) |
198 |
end do |
199 |
! |
200 |
! Now compute the quadratic and cubic coefficients used in the |
201 |
! piecewise polynomial representation. |
202 |
! |
203 |
do i = 1, cs%np - 1 |
204 |
dx = x(i+1) - x(i) |
205 |
divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx |
206 |
divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1 |
207 |
cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx |
208 |
cs%c(4,i) = divdif3 / ( dx * dx ) |
209 |
end do |
210 |
|
211 |
cs%c(3,cs%np) = 0.0_DP |
212 |
cs%c(4,cs%np) = 0.0_DP |
213 |
|
214 |
cs%dx = dx |
215 |
cs%dx_i = 1.0_DP / dx |
216 |
return |
217 |
end subroutine newSplineWithoutDerivs |
218 |
|
219 |
subroutine deleteSpline(this) |
220 |
|
221 |
type(cubicSpline) :: this |
222 |
|
223 |
if(associated(this%x)) then |
224 |
deallocate(this%x) |
225 |
this%x => null() |
226 |
end if |
227 |
if(associated(this%c)) then |
228 |
deallocate(this%c) |
229 |
this%c => null() |
230 |
end if |
231 |
|
232 |
this%np = 0 |
233 |
|
234 |
end subroutine deleteSpline |
235 |
|
236 |
subroutine lookup_nonuniform_spline(cs, xval, yval) |
237 |
|
238 |
!************************************************************************* |
239 |
! |
240 |
! lookup_nonuniform_spline evaluates a piecewise cubic Hermite interpolant. |
241 |
! |
242 |
! Discussion: |
243 |
! |
244 |
! newSpline must be called first, to set up the |
245 |
! spline data from the raw function and derivative data. |
246 |
! |
247 |
! Modified: |
248 |
! |
249 |
! 06 April 1999 |
250 |
! |
251 |
! Reference: |
252 |
! |
253 |
! Conte and de Boor, |
254 |
! Algorithm PCUBIC, |
255 |
! Elementary Numerical Analysis, |
256 |
! 1973, page 234. |
257 |
! |
258 |
! Parameters: |
259 |
! |
260 |
implicit none |
261 |
|
262 |
type (cubicSpline), intent(in) :: cs |
263 |
real( kind = DP ), intent(in) :: xval |
264 |
real( kind = DP ), intent(out) :: yval |
265 |
real( kind = DP ) :: dx |
266 |
integer :: i, j |
267 |
! |
268 |
! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains |
269 |
! or is nearest to xval. |
270 |
! |
271 |
j = cs%np - 1 |
272 |
|
273 |
do i = 1, cs%np - 2 |
274 |
|
275 |
if ( xval < cs%x(i+1) ) then |
276 |
j = i |
277 |
exit |
278 |
end if |
279 |
|
280 |
end do |
281 |
! |
282 |
! Evaluate the cubic polynomial. |
283 |
! |
284 |
dx = xval - cs%x(j) |
285 |
|
286 |
yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) ) |
287 |
|
288 |
return |
289 |
end subroutine lookup_nonuniform_spline |
290 |
|
291 |
subroutine lookup_uniform_spline(cs, xval, yval) |
292 |
|
293 |
!************************************************************************* |
294 |
! |
295 |
! lookup_uniform_spline evaluates a piecewise cubic Hermite interpolant. |
296 |
! |
297 |
! Discussion: |
298 |
! |
299 |
! newSpline must be called first, to set up the |
300 |
! spline data from the raw function and derivative data. |
301 |
! |
302 |
! Modified: |
303 |
! |
304 |
! 06 April 1999 |
305 |
! |
306 |
! Reference: |
307 |
! |
308 |
! Conte and de Boor, |
309 |
! Algorithm PCUBIC, |
310 |
! Elementary Numerical Analysis, |
311 |
! 1973, page 234. |
312 |
! |
313 |
! Parameters: |
314 |
! |
315 |
implicit none |
316 |
|
317 |
type (cubicSpline), intent(in) :: cs |
318 |
real( kind = DP ), intent(in) :: xval |
319 |
real( kind = DP ), intent(out) :: yval |
320 |
real( kind = DP ) :: dx |
321 |
integer :: i, j |
322 |
! |
323 |
! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains |
324 |
! or is nearest to xval. |
325 |
|
326 |
j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1)) |
327 |
|
328 |
dx = xval - cs%x(j) |
329 |
|
330 |
yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) ) |
331 |
|
332 |
return |
333 |
end subroutine lookup_uniform_spline |
334 |
|
335 |
end module INTERPOLATION |