source: trunk/SRC/Interpolation/inrecgrid.pro

;+
;
; @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
; Without loop
;
; @param x1d {in}{required}{type=1d array}
; the x position on the points
;
; @param y1d {in}{required}{type=1d array}
; the y position on the points
;
; @param left {in}{required}{type=1d monotonically increasing array}
; the position of the "left" border of each cell.
;
; @param bottom {in}{required}{type=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 according to the 1d array defined by
; left and res[1,*] the y index according 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 according 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