1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
subroutine tuhunt (xa, n, x, klo)
C
C The array XA is searched for an element KLO such that
C xa(klo) <= x <= xa(klo+1)
C If X < XA(1) then KLO is set to zero; if X > XA(N) then KLO is set to N.
C That is, KLO = 0 or N is a flag indicating X is out of bounds.
C
C KLO must be set to an initial guess on input; it will then be replaced
C by the correct value on output.
C
C This routine was copied with some modifications from the HUNT
C subroutine in Numerical Recipes by Press, Flannery, Teukolsky and
C Vetterling.
C
C N i: number of elements in each array
C XA i: array of independent-variable values
C X i: the value to be bracketed by elements in XA
C KLO io: the lower index in XA that brackets X
C
CH Phil Hodge, 14-Apr-1988 Subroutine copied from Numerical Recipes HUNT.
CH Phil Hodge, 21-May-1996 Don't flag endpoints of XA as out of bounds.
C
integer n, klo
double precision xa(n), x
C--
integer inc, km, khi
logical ascnd
C Set ASCND, and check for X out of bounds.
if (xa(n).gt.xa(1)) then
ascnd = .true.
if (x .lt. xa(1)) then
klo = 0
return
else if (x .gt. xa(n)) then
klo = n
return
end if
else
ascnd = .false.
if (x .gt. xa(1)) then
klo = 0
return
else if (x .lt. xa(n)) then
klo = n
return
end if
end if
if ((klo .le. 0) .or. (klo .gt. n)) then
klo = 1
khi = n
go to 3
endif
inc = 1
if ((x.ge.xa(klo) .and. ascnd) .or.
+ (x.lt.xa(klo) .and. .not.ascnd)) then
1 khi = klo + inc
if (khi .gt. n) then
khi = n + 1
else if ((x.ge.xa(khi) .and. ascnd) .or.
+ (x.lt.xa(khi) .and. .not.ascnd)) then
klo = khi
inc = inc + inc
go to 1
endif
else
khi = klo
2 klo = khi - inc
if (klo .lt. 1) then
klo = 0
else if ((x.lt.xa(klo) .and. ascnd) .or.
+ (x.ge.xa(klo) .and. .not.ascnd)) then
khi = klo
inc = inc + inc
go to 2
endif
endif
3 continue
C Before we return, make sure we don't return a value of KLO that
C implies X is out of bounds. We know it isn't because we checked
C at the beginning.
if (khi-klo .eq. 1) then
klo = max (klo, 1)
klo = min (klo, n-1)
return
end if
km = (khi + klo) / 2
if ((x .gt. xa(km) .and. ascnd) .or.
+ (x .le. xa(km) .and. .not.ascnd)) then
klo = km
else
khi = km
endif
go to 3
end
|
