1 files changed, 103 insertions, 0 deletions
diff --git a/pkg/utilities/nttools/trebin/tuhunt.f b/pkg/utilities/nttools/trebin/tuhunt.f
new file mode 100644
index 00000000..3d8df648
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tuhunt.f
@@ -0,0 +1,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