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 |