ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/utils/interpolation.F90
(Generate patch)

Comparing trunk/src/utils/interpolation.F90 (file contents):
Revision 932 by gezelter, Fri Apr 14 20:04:31 2006 UTC vs.
Revision 938 by gezelter, Mon Apr 17 21:49:12 2006 UTC

# Line 43 | Line 43
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.
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.2 2006-04-14 20:04:31 gezelter Exp $
51 > !! @version $Id: interpolation.F90,v 1.6 2006-04-17 21:49:12 gezelter Exp $
52  
53  
54 < module  INTERPOLATION
54 > module interpolation
55    use definitions
56    use status
57    implicit none
# Line 59 | Line 60 | module  INTERPOLATION
60    character(len = statusMsgSize) :: errMSG
61  
62    type, public :: cubicSpline
63 <     private
63 >     logical :: isUniform = .false.
64       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 newSplineWithoutDerivs
72 <     module procedure newSplineWithDerivs
73 <  end interface
74 <
70 >  public :: newSpline
71    public :: deleteSpline
72 <
72 >  public :: lookupSpline
73 >  public :: lookupUniformSpline
74 >  public :: lookupNonuniformSpline
75 >  public :: lookupUniformSpline1d
76 >  
77   contains
78 +  
79  
80 <
81 <  subroutine newSplineWithoutDerivs(cs, x, y, yp1, ypn, boundary)
81 <
80 >  subroutine newSpline(cs, x, y, yp1, ypn, isUniform)
81 >    
82      !************************************************************************
83      !
84 <    ! newSplineWithoutDerivs solves for slopes defining a cubic spline.
84 >    ! newSpline solves for slopes defining a cubic spline.
85      !
86      !  Discussion:
87      !
# Line 98 | Line 98 | contains
98      !  Parameters:
99      !
100      !    Input, real x(N), the abscissas or X values of
101 <    !    the data points.  The entries of TAU are assumed to be
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 <    !    yp1 contains the slope at x(1) and ypn contains
108 <    !    the slope at x(N).
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 intermediate slopes at x(I) have been
111 <    !    stored in cs%C(2,I), for I = 2 to N-1.
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 <    character(len=*), intent(in) :: boundary
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::newSplineWithoutDerivs", &
126 <            "Type was already created")
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 <    if (size(x) .ne. size(y)) then
133 <       call handleError("interpolation::newSplineWithoutDerivs", &
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 <
137 <    np = size(x)
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::newSplineWithoutDerivs", &
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::newSplineWithoutDerivs", &
150 >       call handleError("interpolation::newSpline", &
151              "Error in allocating storage for c")
152      endif
153        
# Line 154 | Line 156 | contains
156         cs%c(1,i) = y(i)      
157      enddo
158  
159 <    if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
160 <         (boundary.eq.'b').or.(boundary.eq.'B')) then
159 <       cs%c(2,1) = yp1
160 <    else
161 <       cs%c(2,1) = 0.0_DP
162 <    endif
163 <    if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
164 <         (boundary.eq.'b').or.(boundary.eq.'B')) then
165 <       cs%c(2,1) = ypn
166 <    else
167 <       cs%c(2,1) = 0.0_DP
168 <    endif
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
178 –    !
179 –    !  Set the diagonal coefficients.
180 –    !
181 –    cs%c(4,1) = 1.0_DP
182 –    do i = 2, cs%np - 1
183 –       cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) )
184 –    end do
185 –    cs%c(4,cs%np) = 1.0_DP
186 –    !
187 –    !  Set the off-diagonal coefficients.
188 –    !
189 –    cs%c(3,1) = 0.0_DP
190 –    do i = 2, cs%np
191 –       cs%c(3,i) = x(i) - x(i-1)
192 –    end do
193 –    !
194 –    !  Forward elimination.
195 –    !
196 –    do i = 2, cs%np - 1
197 –       g = -cs%c(3,i+1) / cs%c(4,i-1)
198 –       cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1)
199 –       cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1)
200 –    end do
201 –    !
202 –    !  Back substitution for the interior slopes.
203 –    !
204 –    do i = cs%np - 1, 2, -1
205 –       cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i)
206 –    end do
207 –    !
208 –    !  Now compute the quadratic and cubic coefficients used in the
209 –    !  piecewise polynomial representation.
210 –    !
211 –    do i = 1, cs%np - 1
212 –       dx = x(i+1) - x(i)
213 –       divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx
214 –       divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1
215 –       cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx
216 –       cs%c(4,i) = divdif3 / ( dx * dx )
217 –    end do
174  
219 –    cs%c(3,cs%np) = 0.0_DP
220 –    cs%c(4,cs%np) = 0.0_DP
221 –
222 –    cs%dx = dx
223 –    cs%dx_i = 1.0_DP / dx
224 –    return
225 –  end subroutine newSplineWithoutDerivs
226 –
227 –  subroutine newSplineWithDerivs(cs, x, y, yp)
228 –
229 –    !************************************************************************
175      !
231 –    ! newSplineWithDerivs
232 –
233 –    implicit none
234 –
235 –    type (cubicSpline), intent(inout) :: cs
236 –    real( kind = DP ), intent(in) :: x(:), y(:), yp(:)
237 –    real( kind = DP ) :: g, divdif1, divdif3, dx
238 –    integer :: i, alloc_error, np
239 –
240 –    alloc_error = 0
241 –
242 –    if (cs%np .ne. 0) then
243 –       call handleWarning("interpolation::newSplineWithDerivs", &
244 –            "Type was already created")
245 –       call deleteSpline(cs)
246 –    end if
247 –
248 –    ! make sure the sizes match
249 –
250 –    if ((size(x) .ne. size(y)).or.(size(x) .ne. size(yp))) then
251 –       call handleError("interpolation::newSplineWithDerivs", &
252 –            "Array size mismatch")
253 –    end if
254 –    
255 –    np = size(x)
256 –    cs%np = np
257 –
258 –    allocate(cs%x(np), stat=alloc_error)
259 –    if(alloc_error .ne. 0) then
260 –       call handleError("interpolation::newSplineWithDerivs", &
261 –            "Error in allocating storage for x")
262 –    endif
263 –    
264 –    allocate(cs%c(4,np), stat=alloc_error)
265 –    if(alloc_error .ne. 0) then
266 –       call handleError("interpolation::newSplineWithDerivs", &
267 –            "Error in allocating storage for c")
268 –    endif
269 –    
270 –    do i = 1, np
271 –       cs%x(i) = x(i)
272 –       cs%c(1,i) = y(i)      
273 –       cs%c(2,i) = yp(i)
274 –    enddo
275 –    !
176      !  Set the diagonal coefficients.
177      !
178      cs%c(4,1) = 1.0_DP
# Line 316 | Line 216 | contains
216      cs%c(3,cs%np) = 0.0_DP
217      cs%c(4,cs%np) = 0.0_DP
218  
319 –    cs%dx = dx
219      cs%dx_i = 1.0_DP / dx
220  
221      return
222 <  end subroutine newSplineWithDerivs
222 >  end subroutine newSpline
223  
224    subroutine deleteSpline(this)
225  
# Line 339 | Line 238 | contains
238      
239    end subroutine deleteSpline
240  
241 <  subroutine lookup_nonuniform_spline(cs, xval, yval)
241 >  subroutine lookupNonuniformSpline(cs, xval, yval)
242      
243      !*************************************************************************
244      !
245 <    ! lookup_nonuniform_spline evaluates a piecewise cubic Hermite interpolant.
245 >    ! lookupNonuniformSpline evaluates a piecewise cubic Hermite interpolant.
246      !
247      !  Discussion:
248      !
# Line 392 | Line 291 | contains
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
294 >  end subroutine lookupNonuniformSpline
295  
296 <  subroutine lookup_uniform_spline(cs, xval, yval)
296 >  subroutine lookupUniformSpline(cs, xval, yval)
297      
298      !*************************************************************************
299      !
300 <    ! lookup_uniform_spline evaluates a piecewise cubic Hermite interpolant.
300 >    ! lookupUniformSpline evaluates a piecewise cubic Hermite interpolant.
301      !
302      !  Discussion:
303      !
# Line 423 | Line 322 | contains
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
325 >    real( kind = DP ) :: a, b, c, d, dx
326      integer :: i, j
327      !
328      !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
# Line 433 | Line 332 | contains
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) ) )
335 >    a = cs%c(1,j)
336 >    b = cs%c(2,j)
337 >    c = cs%c(3,j)
338 >    d = cs%c(4,j)
339 >
340 >    yval = c + dx * d
341 >    yval = b + dx * yval  
342 >    yval = a + dx * yval
343      
344      return
345 <  end subroutine lookup_uniform_spline
345 >  end subroutine lookupUniformSpline
346 >
347 >  subroutine lookupUniformSpline1d(cs, xval, yval, dydx)
348 >    
349 >    implicit none
350 >
351 >    type (cubicSpline), intent(in) :: cs
352 >    real( kind = DP ), intent(in)  :: xval
353 >    real( kind = DP ), intent(out) :: yval, dydx
354 >    real( kind = DP ) :: a, b, c, d, dx
355 >    integer :: i, j
356 >    
357 >    !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
358 >    !  or is nearest to xval.
359 >
360 >    j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))
361 >
362 >    dx = xval - cs%x(j)
363 >
364 >    a = cs%c(1,j)
365 >    b = cs%c(2,j)
366 >    c = cs%c(3,j)
367 >    d = cs%c(4,j)
368 >
369 >    yval = c + dx * d
370 >    yval = b + dx * yval  
371 >    yval = a + dx * yval
372 >
373 >    dydx = 2.0d0 * c + 3.0d0 * d * dx
374 >    dydx = b + dx * dydx
375 >      
376 >    return
377 >  end subroutine lookupUniformSpline1d
378 >
379 >  subroutine lookupSpline(cs, xval, yval)
380 >
381 >    type (cubicSpline), intent(in) :: cs
382 >    real( kind = DP ), intent(inout) :: xval
383 >    real( kind = DP ), intent(inout) :: yval
384 >    
385 >    if (cs%isUniform) then
386 >       call lookupUniformSpline(cs, xval, yval)
387 >    else
388 >       call lookupNonuniformSpline(cs, xval, yval)
389 >    endif
390 >
391 >    return
392 >  end subroutine lookupSpline
393    
394 < end module INTERPOLATION
394 > end module interpolation

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines