source: trunk/SRC/Interpolation/inrecgrid.pro @ 134

;+
;
; @file_comments 
; given - a list of points, (x,y) position
;       - the x and y limits of a rectangular grid
; find in which cell is located each given point.
;
; @categories no DO loop, use the wonderfull value_locate function!
;
; @param x1d {in}{required}
; a 1d array, the x position on the points
;
; @param y1d {in}{required}
; a 1d array, the y position on the points
;
; @param left {in}{required}
; a 1d, monotonically increasing array,
; the position of the "left" border of each cell.
;
; @param bottom {in}{required}
; a 1d, monotonically increasing array,
; the position of the "bottom" border of each cell.
;
; @keyword OUTPUT2D
; to get the output as a 2d array (2,n_elements(x1d)),
;    with res[0,*] the x index accoring to the 1d array defined by
;    left and res[1,*] the y index accoring to the 1d array defined by
;    bottom.
;
; @keyword CHECKOUT
; = [rbgrid,ubgrid] specify the right and upper boundaries of
;    the grid and check if some points are out.
;
; @returns
; the index on the cell accoring to the 2d array defined by left and bottom.
;
; @examples
;
; IDL> a=indgen(5)
; IDL> b=indgen(7)
; IDL> r=inrecgrid([0.25,3.25,2],[4.25,2.8,1.4],a,b)
; IDL> print, r
;            20          13           7
; IDL> r=inrecgrid([0.25,3.25,2],[4.25,2.8,1.4],a,a+1,b,b+1,/output2d)
; IDL> print, r
;        0.00000      4.00000
;        3.00000      2.00000
;        2.00000      1.00000
;
; @history
; S. Masson (smasson\@lodyc.jussieu.fr)
; July 3rd, 2002
; October 3rd, 2003: use value_locate
;
; @version $Id$
;
;-
FUNCTION inrecgrid, x1d, y1d, left, bottom, OUTPUT2D = output2d, CHECKOUT = checkout
;
  compile_opt idl2, strictarrsubs
;
  ncellx = n_elements(left)
  ncelly = n_elements(bottom)
;
  xpos = value_locate(left[*], x1d[*])
  ypos = value_locate(bottom[*], y1d[*])
;
  IF n_elements(checkout) EQ 2 THEN BEGIN
    out = where(x1d GT checkout[0])
    IF out[0] NE -1 THEN xpos[out] = -1
    out = where(y1d GT checkout[1])
    IF out[0] NE -1 THEN ypos[out] = -1
  ENDIF
;
  IF keyword_set(output2d) THEN return, [transpose(xpos), transpose(ypos)]
;
  IF NOT keyword_set(checkout) THEN RETURN, xpos+ncellx*ypos
;
  res = xpos+ncellx*ypos
  out = where(xpos EQ -1 OR ypos EQ -1)
  IF out[0] NE -1 THEN res[out] = -1
;
  RETURN, res

END