| 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 |
|
|
