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

;+
;
; @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!
;
; @examples 
; res = inrecgrid(xin, yin, left, bottom)
;
;    @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
;    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 bondaries 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
;-

FUNCTION inrecgrid, x1d, y1d, left, bottom, output2d = output2d, checkout = checkout
;
  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