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 933 by chrisfen, Fri Apr 14 21:06:55 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.3 2006-04-14 21:06:55 chrisfen 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 newSpline
72 <  end interface
73 <
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 newSpline(cs, x, y, yp1, ypn)
80 <
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 103 | Line 104 | contains
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 +    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 <
135 <    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 158 | Line 162 | contains
162      cs%c(2,1) = yp1
163      cs%c(2,np) = ypn
164      
161 –
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      !
# Line 211 | Line 216 | contains
216      cs%c(3,cs%np) = 0.0_DP
217      cs%c(4,cs%np) = 0.0_DP
218  
214 –    cs%dx = dx
219      cs%dx_i = 1.0_DP / dx
220 +
221      return
222 <  end subroutine newSplineWithoutDerivs
222 >  end subroutine newSpline
223  
224    subroutine deleteSpline(this)
225  
# Line 233 | 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 286 | 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 317 | 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 327 | 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