;+
; NAME:
;   VALUE_LOCATE
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;
;   Locate one or more values in a reference array (IDL LE 5.2 compatibility)
;
; CALLING SEQUENCE:
;
;   INDICES = VALUE_LOCATE(REF, VALUES)
;
; DESCRIPTION: 
;
;   VALUE_LOCATE locates the positions of given values within a
;   reference array.  The reference array need not be regularly
;   spaced.  This is useful for various searching, sorting and
;   interpolation algorithms.
;
;   The reference array should be a monotonically increasing or
;   decreasing list of values which partition the real numbers.  A
;   reference array of NBINS numbers partitions the real number line
;   into NBINS+1 regions, like so:
;
;
; REF:           X[0]         X[1]   X[2] X[3]     X[NBINS-1]
;      <----------|-------------|------|---|----...---|--------------->
; INDICES:  -1           0          1    2       3        NBINS-1
;
;
;   VALUE_LOCATE returns which partition each of the VALUES falls
;   into, according to the figure above.  For example, a value between
;   X[1] and X[2] would return a value of 1.  Values below X[0] return
;   -1, and above X[NBINS-1] return NBINS-1.  Thus, besides the value
;   of -1, the returned INDICES refer to the nearest reference value
;   to the left of the requested value.
;
;   If the reference array is monotonically decreasing then the
;   partitions are numbered starting at -1 from the right instead (and
;   the returned INDICES refer to the nearest reference value to the
;   *right* of the requested value).  If the reference array is
;   neither monotonically increasing or decreasing the results of
;   VALUE_LOCATE are undefined.
;
;   VALUE_LOCATE appears as a built-in funcion in IDL v5.3 and later.
;   This version of VALUE_LOCATE should work under IDL v4 and later,
;   and is intended to provide a portable solution for users who do
;   not have the latest version of IDL.  The algrorithm in this file
;   is slower but not terribly so, than the built-in version.
;
;   Users should be able to place this file in their IDL path safely:
;   under IDL 5.3 and later, the built-in function will take
;   precedence; under IDL 5.2 and earlier, this function will be used.
;
; INPUTS:
;
;   REF - the reference array of monotonically increasing or
;         decreasing values.
;
;   VALUES - a scalar value or array of values to be located in the
;            reference array.
;
;
; KEYWORDS:
;
;   L64 - (ignored) for compatibility with built-in version. 
;
;   NO_CROP - if set, and VALUES is outside of the region between X[0]
;             and X[NBINS-1], then the returned indices may be *less
;             than* -1 or *greater than* NBINS-1.  The user is the
;             responsible for cropping these values appropriately.
;
; RETURNS:
;
;   An array of indices between -1L and NBINS-1.  If VALUES is an
;   array then the returned array will have the same dimensions.
;
;
; EXAMPLE:
;
;   Cast random values into a histogram with bins from 1-10, 10-100,
;   100-1000, and 1000-10,000.
;
;     ;; Make bin edges - this is the ref. array
;     xbins = 10D^dindgen(5)  
;
;     ;; Make some random data that ranges from 1 to 10,000
;     x     = 10D^(randomu(seed,1000)*4)
;
;     ;; Find the bin number of each random value
;     ii    = value_locate(xbins, x)
;
;     ;; Histogram the data
;     hh    = histogram(ii)
;
;
; SEE ALSO:
;
;   VALUE_LOCATE (IDL 5.3 and later), HISTOGRAM, CMHISTOGRAM
;
;
; MODIFICATION HISTORY:
;   Written and documented, 21 Jan 2001
;   Case of XBINS having only one element, CM, 29 Apr 2001
;   Handle case of VALUES exactly hitting REF points, CM, 13 Oct 2001
; 
;  $Id: value_locate.pro,v 1.5 2001/10/13 17:59:34 craigm Exp $
;
;-
; Copyright (C) 2001, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
function value_locate, xbins, x, l64=l64, no_crop=nocrop, _EXTRA=extra

  on_error, 2
  nbins = n_elements(xbins)
  sz = size(xbins)

  ;; Error checking
  if nbins EQ 0 then message, 'ERROR: XBINS must have at least one element'
  if nbins EQ 1 then return, (x GE xbins(0)) - 1L
  
  ;; The values are computed by spline interpolation.  Here is the "y"
  ;; value of the spline, which is just the bin position.
  tp = sz(sz(0)+1)
  if tp EQ 1 OR tp EQ 2 OR tp EQ 12 then begin
      yy = findgen(nbins) - 0.5
      eps = (machar()).eps
  endif else begin
      yy = dindgen(nbins) - 0.5D
      eps = (machar(/double)).eps
  endelse      

  ;; Check if we are reversing.
  if xbins(nbins-1) GT xbins(0) then rev = 0 else rev = 1

  ;; Compute the spline interpolation.  Note here that we set the 2nd
  ;; derivative value to zero since the derivative computed by
  ;; SPL_INIT seems to be royally screwed up.  Also we do separate
  ;; computations for the increasing and non-increasing cases, since
  ;; SPL_INTERP seems to choke on the later.
  if rev EQ 0 then $
    ii = round(spl_interp(xbins, yy, yy*0, x) + eps) $
  else $
    ii = round(spl_interp(reverse(xbins), yy, yy*0, x) + eps)

  ;; Crop the end values appropriately
  if NOT keyword_set(nocrop) then begin
      ii = ii > (-1L) < (nbins-1)
  endif
      
  ;; Reverse the array
  if rev EQ 0 then $
    ret = temporary(ii) $
  else $
    ret = (nbins-2)-temporary(ii)

  ;; Reform the array to the correct dimensions (ie, add trailing
  ;; dimensions)
  sz = size(x)
  if sz(0) GT 0 then $
    ret = reform(ret, sz(1:sz(0)), /overwrite)

  return, ret
end