GCC Code Coverage Report


Directory: ./
File: misc/interpolation.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 42 0.0%
Branches: 0 36 0.0%

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