locate.f

37 lines · 1 definitions · on GitHub (7be4133)

Python version: locate.py, 52 lines, on the right (switch to Python above the code).

Purpose

LOCATE finds, by bisection, the index $j$ such that $x$ lies between XX(j) and XX(j+1) in a monotonic array (ascending or descending), as in Numerical Recipes. $j = 0$ or $N$ means out of range.

How it fits

Called by the Livermore-table routine GTPURZ (GTLAOL and GTLAOL2 also called it; removed since rp_1) and by ALSURFP (Allard atmospheres). No COMMON blocks, no calls.

The method

Standard bisection. Then two tolerances are added: if $j = 0$ but $x > 0.99\,\texttt{XX(1)}$, it returns $j = 1$; if $j = N$ but $x < 1.01\,\texttt{XX(N)}$, it returns $N-1$.

Definitions

subroutineLOCATE

lines 9–37 · GitHub

Bisection search with 1 % end tolerances (above).

Called by ALSURFP ×2, GTPURZ

Change log

The changes since rp_1

Changed since rp_1. Declarations only: IMPLICIT NONE with explicit declarations. Two of its callers (GTLAOL, GTLAOL2) are removed since rp_1.

1C 2CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 3C LOCATE 4CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 5C GIVEN AN ARRAY XX OF LENGTH N, AND GIVEN A VALUE X, RETURNS 6C A VALUE J SUCH THAT X IS BETWEEN XX(J) AND XX(J+1). X=0 OR N 7C THEN OUT OF RANGE. 8C 9 SUBROUTINE LOCATE(XX, N, X, J) 10C 10/20/26 MHP via Claude 11C code changed to replace implicit typing with IMPLICIT NONE and 12C explicit declarations. 13 IMPLICIT NONE 14 REAL*8 X,XX 15 INTEGER J,JL,JM,JU,N 16 DIMENSION XX(N) 17 SAVE 18 JL = 0 19 JU = N+1 20 10 10 IF (JU-JL.GT.1) THEN 21 JM = (JU+JL)/2 22 IF((XX(N).GT.XX(1)).EQV.(X.GT.XX(JM))) THEN 23 JL = JM 24 ELSE 25 JU = JM 26 END IF 27 GO TO 10 28 END IF 29 J = JL 30 IF((J .EQ. 0) .AND. (X .GT. 0.99D0*XX(1)))THEN 31 J = 1 32 END IF 33 IF((J .EQ. N) .AND. (X .LT. 1.01D0*XX(N)))THEN 34 J = N-1 35 END IF 36 RETURN 37 END

1"""Python version of locate.f (YREC development, 7be4133).""" 2 3import numpy as np 4 5 6def locate(xx: np.ndarray, n: int, x: float, j: int) -> int: 7 """Bisection search in a monotonic table, with 1 % end tolerances. 8 9 Finds ``j`` such that ``x`` lies between ``xx[j]`` and ``xx[j+1]`` 10 (ascending or descending table); ``j = 0`` or ``n`` means out of 11 range, except that values within 1 % of the end values are pulled 12 back in (``j = 1`` or ``n - 1``). 13 14 Translates: locate.f#LOCATE (7be4133) 15 16 Parameters 17 ---------- 18 xx : ndarray, shape (n + 1,) 19 The table, ``xx[1..n]``. 20 n : int 21 Its length. 22 x : float 23 The value to locate. 24 j : int 25 Output only (passed in with any value). 26 27 Returns 28 ------- 29 j : int 30 The lower index of the bracketing interval. 31 32 """ 33 34 # Bisection (Fortran lines 18–29). 35 jl = 0 36 ju = n + 1 37 while ju - jl > 1: # label 10 38 jm = (ju + jl) // 2 39 if (xx[n] > xx[1]) == (x > xx[jm]): 40 jl = jm 41 else: 42 ju = jm 43 j = jl 44 45 # The 1 % end tolerances (Fortran lines 30–35). 46 # Question (see the page): the tolerances are relative to the end 47 # values, so they act only for positive ends of an ascending table. 48 if j == 0 and x > 0.99 * xx[1]: 49 j = 1 50 if j == n and x < 1.01 * xx[n]: 51 j = n - 1 52 return j