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 960 by gezelter, Wed May 17 15:37:15 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.
47   !!
48   !! @author Charles F. Vardeman II
49 < !! @version $Id: interpolation.F90,v 1.3 2006-04-14 21:06:55 chrisfen Exp $
49 > !! @version $Id: interpolation.F90,v 1.8 2006-05-17 15:37:15 gezelter Exp $
50  
51  
52 < module  INTERPOLATION
52 > module interpolation
53    use definitions
54    use status
55    implicit none
56    PRIVATE
57  
59 –  character(len = statusMsgSize) :: errMSG
60 –
58    type, public :: cubicSpline
59 <     private
60 <     integer :: np = 0
64 <     real(kind=dp) :: dx
59 >     logical :: isUniform = .false.
60 >     integer :: n = 0
61       real(kind=dp) :: dx_i
62       real (kind=dp), pointer,dimension(:)   :: x => null()
63 <     real (kind=dp), pointer,dimension(:,:) :: c => null()
63 >     real (kind=dp), pointer,dimension(:)   :: y => null()
64 >     real (kind=dp), pointer,dimension(:)   :: b => null()
65 >     real (kind=dp), pointer,dimension(:)   :: c => null()
66 >     real (kind=dp), pointer,dimension(:)   :: d => null()
67    end type cubicSpline
68  
69 <  interface newSpline
71 <     module procedure newSpline
72 <  end interface
73 <
69 >  public :: newSpline
70    public :: deleteSpline
71 <
71 >  public :: lookupSpline
72 >  public :: lookupUniformSpline
73 >  public :: lookupNonuniformSpline
74 >  public :: lookupUniformSpline1d
75 >  
76   contains
77 +  
78  
79 <
80 <  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 <
79 >  subroutine newSpline(cs, x, y, isUniform)
80 >    
81      implicit none
82  
83      type (cubicSpline), intent(inout) :: cs
84      real( kind = DP ), intent(in) :: x(:), y(:)
85 <    real( kind = DP ), intent(in) :: yp1, ypn
86 <    real( kind = DP ) :: g, divdif1, divdif3, dx
118 <    integer :: i, alloc_error, np
85 >    real( kind = DP ) :: fp1, fpn, p
86 >    REAL( KIND = DP), DIMENSION(size(x)-1) :: diff_y, H
87  
88 +    logical, intent(in) :: isUniform
89 +    integer :: i, alloc_error, n, k
90 +
91      alloc_error = 0
92  
93 <    if (cs%np .ne. 0) then
94 <       call handleWarning("interpolation::newSplineWithoutDerivs", &
95 <            "Type was already created")
93 >    if (cs%n .ne. 0) then
94 >       call handleWarning("interpolation::newSpline", &
95 >            "cubicSpline struct was already created")
96         call deleteSpline(cs)
97      end if
98  
99      ! make sure the sizes match
100  
101 <    if (size(x) .ne. size(y)) then
102 <       call handleError("interpolation::newSplineWithoutDerivs", &
101 >    n = size(x)
102 >    
103 >    if ( size(y) .ne. size(x) ) then
104 >       call handleError("interpolation::newSpline", &
105              "Array size mismatch")
106      end if
107 +    
108 +    cs%n = n
109 +    cs%isUniform = isUniform
110  
111 <    np = size(x)
136 <    cs%np = np
137 <
138 <    allocate(cs%x(np), stat=alloc_error)
111 >    allocate(cs%x(n), stat=alloc_error)
112      if(alloc_error .ne. 0) then
113 <       call handleError("interpolation::newSplineWithoutDerivs", &
113 >       call handleError("interpolation::newSpline", &
114              "Error in allocating storage for x")
115      endif
116  
117 <    allocate(cs%c(4,np), stat=alloc_error)
117 >    allocate(cs%y(n), stat=alloc_error)
118      if(alloc_error .ne. 0) then
119 <       call handleError("interpolation::newSplineWithoutDerivs", &
119 >       call handleError("interpolation::newSpline", &
120 >            "Error in allocating storage for y")
121 >    endif
122 >
123 >    allocate(cs%b(n), stat=alloc_error)
124 >    if(alloc_error .ne. 0) then
125 >       call handleError("interpolation::newSpline", &
126 >            "Error in allocating storage for b")
127 >    endif
128 >
129 >    allocate(cs%c(n), stat=alloc_error)
130 >    if(alloc_error .ne. 0) then
131 >       call handleError("interpolation::newSpline", &
132              "Error in allocating storage for c")
133      endif
134 <      
135 <    do i = 1, np
134 >
135 >    allocate(cs%d(n), stat=alloc_error)
136 >    if(alloc_error .ne. 0) then
137 >       call handleError("interpolation::newSpline", &
138 >            "Error in allocating storage for d")
139 >    endif
140 >
141 >    ! make sure we are monotinically increasing in x:
142 >
143 >    h = diff(x)
144 >    if (any(h <= 0)) then
145 >       call handleError("interpolation::newSpline", &
146 >            "Negative dx interval found")
147 >    end if
148 >
149 >    ! load x and y values into the cubicSpline structure:
150 >
151 >    do i = 1, n
152         cs%x(i) = x(i)
153 <       cs%c(1,i) = y(i)      
154 <    enddo
153 >       cs%y(i) = y(i)
154 >    end do
155  
156 <    ! Set the first derivative of the function to the second coefficient of
157 <    ! each of the endpoints
156 >    ! Calculate coefficients for the tridiagonal system: store
157 >    ! sub-diagonal in B, diagonal in D, difference quotient in C.
158  
159 <    cs%c(2,1) = yp1
160 <    cs%c(2,np) = ypn
159 >    cs%b(1:n-1) = h
160 >    diff_y = diff(y)
161 >    cs%c(1:n-1) = diff_y / h
162 >
163 >    if (n == 2) then
164 >       ! Assume the derivatives at both endpoints are zero
165 >       ! another assumption could be made to have a linear interpolant
166 >       ! between the two points.  In that case, the b coefficients
167 >       ! below would be diff_y(1)/h(1) and the c and d coefficients would
168 >       ! both be zero.
169 >       cs%b(1) = 0.0_dp
170 >       cs%c(1) = -3.0_dp * (diff_y(1)/h(1))**2
171 >       cs%d(1) = -2.0_dp * (diff_y(1)/h(1))**3
172 >       cs%b(2) = cs%b(1)
173 >       cs%c(2) = 0.0_dp
174 >       cs%d(2) = 0.0_dp
175 >       cs%dx_i = 1.0_dp / h(1)
176 >      return
177 >    end if
178 >
179 >    cs%d(1) = 2.0_dp * cs%b(1)
180 >    do i = 2, n-1
181 >      cs%d(i) = 2.0_dp * (cs%b(i) + cs%b(i-1))
182 >    end do
183 >    cs%d(n) = 2.0_dp * cs%b(n-1)
184 >
185 >    ! Calculate estimates for the end slopes using polynomials
186 >    ! that interpolate the data nearest the end.
187      
188 +    fp1 = cs%c(1) - cs%b(1)*(cs%c(2) - cs%c(1))/(cs%b(1) + cs%b(2))
189 +    if (n > 3) then
190 +      fp1 = fp1 + cs%b(1)*((cs%b(1) + cs%b(2))*(cs%c(3) - cs%c(2))/ &
191 +           (cs%b(2) + cs%b(3)) - cs%c(2) + cs%c(1))/(x(4) - x(1))
192 +    end if
193 +          
194 +    fpn = cs%c(n-1) + cs%b(n-1)*(cs%c(n-1) - cs%c(n-2))/(cs%b(n-2) + cs%b(n-1))
195 +    if (n > 3) then
196 +      fpn = fpn + cs%b(n-1)*(cs%c(n-1) - cs%c(n-2) - (cs%b(n-2) + cs%b(n-1))* &
197 +           (cs%c(n-2) - cs%c(n-3))/(cs%b(n-2) + cs%b(n-3)))/(x(n) - x(n-3))
198 +    end if
199  
200 <    !
201 <    !  Set up the right hand side of the linear system.
202 <    !
203 <    do i = 2, cs%np - 1
204 <       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)))
200 >    ! Calculate the right hand side and store it in C.
201 >
202 >    cs%c(n) = 3.0_dp * (fpn - cs%c(n-1))
203 >    do i = n-1,2,-1
204 >      cs%c(i) = 3.0_dp * (cs%c(i) - cs%c(i-1))
205      end do
206 <    !
207 <    !  Set the diagonal coefficients.
208 <    !
209 <    cs%c(4,1) = 1.0_DP
210 <    do i = 2, cs%np - 1
211 <       cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) )
206 >    cs%c(1) = 3.0_dp * (cs%c(1) - fp1)
207 >
208 >    ! Solve the tridiagonal system.
209 >
210 >    do k = 2, n
211 >      p = cs%b(k-1) / cs%d(k-1)
212 >      cs%d(k) = cs%d(k) - p*cs%b(k-1)
213 >      cs%c(k) = cs%c(k) - p*cs%c(k-1)
214      end do
215 <    cs%c(4,cs%np) = 1.0_DP
216 <    !
217 <    !  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)
215 >    cs%c(n) = cs%c(n) / cs%d(n)
216 >    do k = n-1, 1, -1
217 >      cs%c(k) = (cs%c(k) - cs%b(k) * cs%c(k+1)) / cs%d(k)
218      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
219  
220 <    cs%c(3,cs%np) = 0.0_DP
212 <    cs%c(4,cs%np) = 0.0_DP
220 >    ! Calculate the coefficients defining the spline.
221  
222 <    cs%dx = dx
223 <    cs%dx_i = 1.0_DP / dx
224 <    return
217 <  end subroutine newSplineWithoutDerivs
222 >    cs%d(1:n-1) = diff(cs%c) / (3.0_dp * h)
223 >    cs%b(1:n-1) = diff_y / h - h * (cs%c(1:n-1) + h * cs%d(1:n-1))
224 >    cs%b(n) = cs%b(n-1) + h(n-1) * (2.0_dp*cs%c(n-1) + h(n-1)*3.0_dp*cs%d(n-1))
225  
226 +    if (isUniform) then
227 +       cs%dx_i = 1.0_dp / (x(2) - x(1))
228 +    endif
229 +
230 +    return
231 +    
232 +  contains
233 +    
234 +    function diff(v)
235 +      ! Auxiliary function to compute the forward difference
236 +      ! of data stored in a vector v.
237 +      
238 +      implicit none
239 +      real (kind = dp), dimension(:), intent(in) :: v
240 +      real (kind = dp), dimension(size(v)-1) :: diff
241 +      
242 +      integer :: n
243 +      
244 +      n = size(v)
245 +      diff = v(2:n) - v(1:n-1)
246 +      return
247 +    end function diff
248 +    
249 +  end subroutine newSpline
250 +      
251    subroutine deleteSpline(this)
252  
253      type(cubicSpline) :: this
# Line 229 | Line 261 | contains
261         this%c => null()
262      end if
263      
264 <    this%np = 0
264 >    this%n = 0
265      
266    end subroutine deleteSpline
267  
268 <  subroutine lookup_nonuniform_spline(cs, xval, yval)
268 >  subroutine lookupNonuniformSpline(cs, xval, yval)
269      
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 –    !
270      implicit none
271  
272      type (cubicSpline), intent(in) :: cs
273      real( kind = DP ), intent(in)  :: xval
274      real( kind = DP ), intent(out) :: yval
275 <    real( kind = DP ) :: dx
275 >    real( kind = DP ) ::  dx
276      integer :: i, j
277      !
278      !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
279      !  or is nearest to xval.
280      !
281 <    j = cs%np - 1
281 >    j = cs%n - 1
282  
283 <    do i = 1, cs%np - 2
283 >    do i = 0, cs%n - 2
284  
285         if ( xval < cs%x(i+1) ) then
286            j = i
# Line 282 | Line 292 | contains
292      !  Evaluate the cubic polynomial.
293      !
294      dx = xval - cs%x(j)
295 <
286 <    yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
295 >    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))
296      
297      return
298 <  end subroutine lookup_nonuniform_spline
298 >  end subroutine lookupNonuniformSpline
299  
300 <  subroutine lookup_uniform_spline(cs, xval, yval)
300 >  subroutine lookupUniformSpline(cs, xval, yval)
301      
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 –    !
302      implicit none
303  
304      type (cubicSpline), intent(in) :: cs
305      real( kind = DP ), intent(in)  :: xval
306      real( kind = DP ), intent(out) :: yval
307 <    real( kind = DP ) :: dx
307 >    real( kind = DP ) ::  dx
308      integer :: i, j
309      !
310      !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
311      !  or is nearest to xval.
312 +    
313 +    j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
314 +    
315 +    dx = xval - cs%x(j)
316 +    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))
317 +    
318 +    return
319 +  end subroutine lookupUniformSpline
320  
321 <    j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))
321 >  subroutine lookupUniformSpline1d(cs, xval, yval, dydx)
322 >    
323 >    implicit none
324  
325 +    type (cubicSpline), intent(in) :: cs
326 +    real( kind = DP ), intent(in)  :: xval
327 +    real( kind = DP ), intent(out) :: yval, dydx
328 +    real( kind = DP ) :: dx
329 +    integer :: i, j
330 +    
331 +    !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
332 +    !  or is nearest to xval.
333 +
334 +
335 +    j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
336 +    
337      dx = xval - cs%x(j)
338 +    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))
339  
340 <    yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
340 >    dydx = cs%b(j) + dx*(2.0_dp * cs%c(j) + 3.0_dp * dx * cs%d(j))
341 >      
342 >    return
343 >  end subroutine lookupUniformSpline1d
344 >
345 >  subroutine lookupSpline(cs, xval, yval)
346 >
347 >    type (cubicSpline), intent(in) :: cs
348 >    real( kind = DP ), intent(inout) :: xval
349 >    real( kind = DP ), intent(inout) :: yval
350      
351 +    if (cs%isUniform) then
352 +       call lookupUniformSpline(cs, xval, yval)
353 +    else
354 +       call lookupNonuniformSpline(cs, xval, yval)
355 +    endif
356 +
357      return
358 <  end subroutine lookup_uniform_spline
358 >  end subroutine lookupSpline
359    
360 < end module INTERPOLATION
360 > end module interpolation

Diff Legend

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