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 931 by gezelter, Fri Apr 14 19:57:04 2006 UTC vs.
Revision 934 by gezelter, Fri Apr 14 21:49:54 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.1 2006-04-14 19:57:04 gezelter Exp $
51 > !! @version $Id: interpolation.F90,v 1.4 2006-04-14 21:49:54 gezelter Exp $
52  
53  
54   module  INTERPOLATION
# Line 60 | Line 61 | module  INTERPOLATION
61  
62    type, public :: cubicSpline
63       private
64 +     logical :: isUniform = .false.
65       integer :: np = 0
64 –     real(kind=dp) :: dx
66       real(kind=dp) :: dx_i
67       real (kind=dp), pointer,dimension(:)   :: x => null()
68 <     real (kind=dp), pointer,dimension(4,:) :: c => null()
68 >     real (kind=dp), pointer,dimension(:,:) :: c => null()
69    end type cubicSpline
70  
71 <  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 <
71 >  public :: newSpline
72    public :: deleteSpline
73 <
73 >  public :: lookup_spline
74 >  public :: lookup_uniform_spline
75 >  public :: lookup_nonuniform_spline
76 >  
77   contains
78 +  
79  
80 <
81 <  subroutine newSplineWithoutDerivs(cs, x, y, yp1, ypn, boundary)
88 <
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 105 | 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 <
144 <    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 161 | 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
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
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
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
174  
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 –    !************************************************************************
175      !
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 –    !
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,n) = 1.0_DP
182 >    cs%c(4,cs%np) = 1.0_DP
183      !
184      !  Set the off-diagonal coefficients.
185      !
# Line 320 | Line 213 | contains
213         cs%c(4,i) = divdif3 / ( dx * dx )
214      end do
215  
216 <    cs%c(3,np) = 0.0_DP
217 <    cs%c(4,np) = 0.0_DP
216 >    cs%c(3,cs%np) = 0.0_DP
217 >    cs%c(4,cs%np) = 0.0_DP
218  
219 <    cs%dx = dx
327 <    cs%dxi = 1.0_DP / dx
219 >    cs%dx_i = 1.0_DP / dx
220  
221      return
222    end subroutine newSplineWithoutDerivs
# Line 375 | Line 267 | contains
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
# Line 429 | 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
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%dxi) + 1))
331 >    j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))
332  
333      dx = xval - cs%x(j)
334  
# Line 442 | Line 336 | contains
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

Diff Legend

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