47 |
|
!! precomputation of spline parameters. |
48 |
|
!! |
49 |
|
!! @author Charles F. Vardeman II |
50 |
< |
!! @version $Id: interpolation.F90,v 1.1 2006-04-14 19:57:04 gezelter Exp $ |
50 |
> |
!! @version $Id: interpolation.F90,v 1.2 2006-04-14 20:04:31 gezelter Exp $ |
51 |
|
|
52 |
|
|
53 |
|
module INTERPOLATION |
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(4,:) :: c => null() |
67 |
> |
real (kind=dp), pointer,dimension(:,:) :: c => null() |
68 |
|
end type cubicSpline |
69 |
|
|
70 |
– |
interface splineLookup |
71 |
– |
module procedure multiSplint |
72 |
– |
module procedure splintd |
73 |
– |
module procedure splintd1 |
74 |
– |
module procedure splintd2 |
75 |
– |
end interface |
76 |
– |
|
70 |
|
interface newSpline |
71 |
|
module procedure newSplineWithoutDerivs |
72 |
|
module procedure newSplineWithDerivs |
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,n) = 1.0_DP |
185 |
> |
cs%c(4,cs%np) = 1.0_DP |
186 |
|
! |
187 |
|
! Set the off-diagonal coefficients. |
188 |
|
! |
216 |
|
cs%c(4,i) = divdif3 / ( dx * dx ) |
217 |
|
end do |
218 |
|
|
219 |
< |
cs%c(3,np) = 0.0_DP |
220 |
< |
cs%c(4,np) = 0.0_DP |
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%dxi = 1.0_DP / dx |
223 |
> |
cs%dx_i = 1.0_DP / dx |
224 |
|
return |
225 |
|
end subroutine newSplineWithoutDerivs |
226 |
|
|
279 |
|
do i = 2, cs%np - 1 |
280 |
|
cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) ) |
281 |
|
end do |
282 |
< |
cs%c(4,n) = 1.0_DP |
282 |
> |
cs%c(4,cs%np) = 1.0_DP |
283 |
|
! |
284 |
|
! Set the off-diagonal coefficients. |
285 |
|
! |
313 |
|
cs%c(4,i) = divdif3 / ( dx * dx ) |
314 |
|
end do |
315 |
|
|
316 |
< |
cs%c(3,np) = 0.0_DP |
317 |
< |
cs%c(4,np) = 0.0_DP |
316 |
> |
cs%c(3,cs%np) = 0.0_DP |
317 |
> |
cs%c(4,cs%np) = 0.0_DP |
318 |
|
|
319 |
|
cs%dx = dx |
320 |
< |
cs%dxi = 1.0_DP / dx |
320 |
> |
cs%dx_i = 1.0_DP / dx |
321 |
|
|
322 |
|
return |
323 |
< |
end subroutine newSplineWithoutDerivs |
323 |
> |
end subroutine newSplineWithDerivs |
324 |
|
|
325 |
|
subroutine deleteSpline(this) |
326 |
|
|
368 |
|
type (cubicSpline), intent(in) :: cs |
369 |
|
real( kind = DP ), intent(in) :: xval |
370 |
|
real( kind = DP ), intent(out) :: yval |
371 |
+ |
real( kind = DP ) :: dx |
372 |
|
integer :: i, j |
373 |
|
! |
374 |
|
! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains |
423 |
|
type (cubicSpline), intent(in) :: cs |
424 |
|
real( kind = DP ), intent(in) :: xval |
425 |
|
real( kind = DP ), intent(out) :: yval |
426 |
+ |
real( kind = DP ) :: dx |
427 |
|
integer :: i, j |
428 |
|
! |
429 |
|
! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains |
430 |
|
! or is nearest to xval. |
431 |
|
|
432 |
< |
j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dxi) + 1)) |
432 |
> |
j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1)) |
433 |
|
|
434 |
|
dx = xval - cs%x(j) |
435 |
|
|