| Line |
Branch |
Exec |
Source |
| 1 |
|
|
! $Id$ |
| 2 |
|
|
module interpolation |
| 3 |
|
|
|
| 4 |
|
|
! From Press et al., 1996, version 2.10a |
| 5 |
|
|
! B3 Interpolation and Extrapolation |
| 6 |
|
|
|
| 7 |
|
|
IMPLICIT NONE |
| 8 |
|
|
|
| 9 |
|
|
contains |
| 10 |
|
|
|
| 11 |
|
✗ |
pure FUNCTION locate(xx,x) |
| 12 |
|
|
|
| 13 |
|
|
REAL, DIMENSION(:), INTENT(IN) :: xx |
| 14 |
|
|
REAL, INTENT(IN) :: x |
| 15 |
|
|
INTEGER locate |
| 16 |
|
|
|
| 17 |
|
|
! Given an array xx(1:N), and given a value x, returns a value j, |
| 18 |
|
|
! between 0 and N, such that x is between xx(j) and xx(j + 1). xx |
| 19 |
|
|
! must be monotonic, either increasing or decreasing. j = 0 or j = |
| 20 |
|
|
! N is returned to indicate that x is out of range. This |
| 21 |
|
|
! procedure should not be called with a zero-sized array argument. |
| 22 |
|
|
! See notes. |
| 23 |
|
|
|
| 24 |
|
|
INTEGER n,jl,jm,ju |
| 25 |
|
|
LOGICAL ascnd |
| 26 |
|
|
|
| 27 |
|
|
!---------------------------- |
| 28 |
|
|
|
| 29 |
|
✗ |
n=size(xx) |
| 30 |
|
✗ |
ascnd = (xx(n) >= xx(1)) |
| 31 |
|
|
! (True if ascending order of table, false otherwise.) |
| 32 |
|
|
! Initialize lower and upper limits: |
| 33 |
|
|
jl=0 |
| 34 |
|
✗ |
ju=n+1 |
| 35 |
|
✗ |
do while (ju-jl > 1) |
| 36 |
|
✗ |
jm=(ju+jl)/2 ! Compute a midpoint, |
| 37 |
|
✗ |
if (ascnd .eqv. (x >= xx(jm))) then |
| 38 |
|
|
jl=jm ! and replace either the lower limit |
| 39 |
|
|
else |
| 40 |
|
|
ju=jm ! or the upper limit, as appropriate. |
| 41 |
|
|
end if |
| 42 |
|
|
end do |
| 43 |
|
|
! {ju == jl + 1} |
| 44 |
|
|
|
| 45 |
|
|
! {(ascnd .and. xx(jl) <= x < xx(jl+1)) |
| 46 |
|
|
! .neqv. |
| 47 |
|
|
! (.not. ascnd .and. xx(jl+1) <= x < xx(jl))} |
| 48 |
|
|
|
| 49 |
|
|
! Then set the output, being careful with the endpoints: |
| 50 |
|
✗ |
if (x == xx(1)) then |
| 51 |
|
|
locate=1 |
| 52 |
|
✗ |
else if (x == xx(n)) then |
| 53 |
|
✗ |
locate=n-1 |
| 54 |
|
|
else |
| 55 |
|
|
locate=jl |
| 56 |
|
|
end if |
| 57 |
|
|
|
| 58 |
|
✗ |
END FUNCTION locate |
| 59 |
|
|
|
| 60 |
|
|
!*************************** |
| 61 |
|
|
|
| 62 |
|
✗ |
pure SUBROUTINE hunt(xx,x,jlo) |
| 63 |
|
|
|
| 64 |
|
|
! Given an array xx(1:N ), and given a value x, returns a value |
| 65 |
|
|
! jlo such that x is between xx(jlo) and xx(jlo+1). xx must be |
| 66 |
|
|
! monotonic, either increasing or decreasing. jlo = 0 or jlo = N is |
| 67 |
|
|
! returned to indicate that x is out of range. jlo on input is taken as |
| 68 |
|
|
! the initial guess for jlo on output. |
| 69 |
|
|
! Modified so that it uses the information "jlo = 0" on input. |
| 70 |
|
|
|
| 71 |
|
|
INTEGER, INTENT(INOUT) :: jlo |
| 72 |
|
|
REAL, INTENT(IN) :: x |
| 73 |
|
|
REAL, DIMENSION(:), INTENT(IN) :: xx |
| 74 |
|
|
INTEGER n,inc,jhi,jm |
| 75 |
|
|
LOGICAL ascnd, hunt_up |
| 76 |
|
|
|
| 77 |
|
|
!----------------------------------------------------- |
| 78 |
|
|
|
| 79 |
|
✗ |
n=size(xx) |
| 80 |
|
✗ |
ascnd = (xx(n) >= xx(1)) |
| 81 |
|
|
! (True if ascending order of table, false otherwise.) |
| 82 |
|
✗ |
if (jlo < 0 .or. jlo > n) then |
| 83 |
|
|
! Input guess not useful. Go immediately to bisection. |
| 84 |
|
✗ |
jlo=0 |
| 85 |
|
✗ |
jhi=n+1 |
| 86 |
|
|
else |
| 87 |
|
|
inc=1 ! Set the hunting increment. |
| 88 |
|
✗ |
if (jlo == 0) then |
| 89 |
|
|
hunt_up = .true. |
| 90 |
|
|
else |
| 91 |
|
✗ |
hunt_up = x >= xx(jlo) .eqv. ascnd |
| 92 |
|
|
end if |
| 93 |
|
✗ |
if (hunt_up) then ! Hunt up: |
| 94 |
|
✗ |
do |
| 95 |
|
✗ |
jhi=jlo+inc |
| 96 |
|
✗ |
if (jhi > n) then ! Done hunting, since off end of table. |
| 97 |
|
✗ |
jhi=n+1 |
| 98 |
|
✗ |
exit |
| 99 |
|
|
else |
| 100 |
|
✗ |
if (x < xx(jhi) .eqv. ascnd) exit |
| 101 |
|
✗ |
jlo=jhi ! Not done hunting, |
| 102 |
|
✗ |
inc=inc+inc ! so double the increment |
| 103 |
|
|
end if |
| 104 |
|
|
end do ! and try again. |
| 105 |
|
|
else ! Hunt down: |
| 106 |
|
|
jhi=jlo |
| 107 |
|
✗ |
do |
| 108 |
|
✗ |
jlo=jhi-inc |
| 109 |
|
✗ |
if (jlo < 1) then ! Done hunting, since off end of table. |
| 110 |
|
✗ |
jlo=0 |
| 111 |
|
✗ |
exit |
| 112 |
|
|
else |
| 113 |
|
✗ |
if (x >= xx(jlo) .eqv. ascnd) exit |
| 114 |
|
|
jhi=jlo ! Not done hunting, |
| 115 |
|
✗ |
inc=inc+inc ! so double the increment |
| 116 |
|
|
end if |
| 117 |
|
|
end do ! and try again. |
| 118 |
|
|
end if |
| 119 |
|
|
end if ! Done hunting, value bracketed. |
| 120 |
|
|
|
| 121 |
|
|
do ! Hunt is done, so begin the final bisection phase: |
| 122 |
|
✗ |
if (jhi-jlo <= 1) then |
| 123 |
|
✗ |
if (x == xx(n)) jlo=n-1 |
| 124 |
|
✗ |
if (x == xx(1)) jlo=1 |
| 125 |
|
|
exit |
| 126 |
|
|
else |
| 127 |
|
✗ |
jm=(jhi+jlo)/2 |
| 128 |
|
✗ |
if (x >= xx(jm) .eqv. ascnd) then |
| 129 |
|
✗ |
jlo=jm |
| 130 |
|
|
else |
| 131 |
|
|
jhi=jm |
| 132 |
|
|
end if |
| 133 |
|
|
end if |
| 134 |
|
|
end do |
| 135 |
|
|
|
| 136 |
|
✗ |
END SUBROUTINE hunt |
| 137 |
|
|
|
| 138 |
|
|
end module interpolation |
| 139 |
|
|
|