source: trunk/ToBeReviewed/MATRICE/inrecgrid.pro.save @ 67

;+
; NAME: inrecgrid
;
; PURPOSE: 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.
;
; CATEGORY: no DO loop
;
; CALLING SEQUENCE:res = inrecgrid(xin, yin, wb1d, eb1d, sb1d, nb1d)
;
; INPUTS: 
;    x1d: a 1d array, the x position on the points
;    y1d: a 1d array, the y position on the points
;    wb1d: a 1d sorted array, the position of the "west" border of each cell.
;    eb1d: a 1d sorted array, the position of the "east" border of each cell.
;    sb1d: a 1d sorted array, the position of the "south" border of each cell.
;    nb1d: a 1d sorted array, the position of the "north" border of each cell.
;    + see RESTRICTIONS part    
;
; OPTIONAL INPUTS:
;
; KEYWORD PARAMETERS:
;    /output2d: to the output as a 2d array (2,n_elements(x1d)), with
;    res[0,*] the x index accoring to the 1d array defined by wb1d or
;    eb1d and res[1,*] the y index accoring to the 1d array defined by
;    sb1d or nb1d.
;
;    /check: to check if the inputs parameters verify all the test
;    described in section RESTRICTIONS
;
; OUTPUTS:the index on the cell accoring to the 2d array defined by
; wb1d, eb1d, sb1d and nb1d
;
; OPTIONAL OUTPUTS:
;
; COMMON BLOCKS: no
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;    x1d and y1d must have the same number of elements
;    wb1d and eb1d must have the same number of elements
;    sb1d and nb1d must have the same number of elements
;    the cells form a 2d array with the size:
;      n_elements(wb1d)*n_elements(sb1d)
;    we must have:
;      total(eb1d le wb1d) eq 0
;      total(nb1d le sb1d) eq 0
;      min(wb1d) le min(x1d)
;      max(eb1d) ge max(x1d)
;      min(sb1d) le min(y1d)
;      max(nb1d) ge max(y1d)
;      total(sort(wb1d) ne lindgen(n_elements(wb1d))) eq 0
;      total(sort(eb1d) ne lindgen(n_elements(eb1d))) eq 0
;      total(sort(sb1d) ne lindgen(n_elements(sb1d))) eq 0
;      total(sort(nb1d) ne lindgen(n_elements(nb1d))) eq 0
;
; PROCEDURE:
;
; EXAMPLE:
;
;  IDL> a=indgen(5)
;  IDL> b=indgen(7)
;  IDL> r=inrecgrid([0.25,3.25,2],[4.25,2.8,1.4],a,a+1,b,b+1)
;  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
;  
; MODIFICATION HISTORY:
;            S. Masson (smasson@lodyc.jussieu.fr)
;                      July 3rd, 2002
;-

FUNCTION inrecgrid, x1d, y1d, wb1d, eb1d, sb1d, nb1d, output2d = output2d,  check = check
; test
  IF keyword_set(check) THEN BEGIN
    IF total(eb1d le wb1d) ne 0 then stop
    IF total(nb1d le sb1d) ne 0 then stop
    IF min(wb1d) gt min(x1d) then stop
    IF max(eb1d) lt max(x1d) then stop
    IF min(sb1d) gt min(y1d) then stop
    IF max(nb1d) lt max(y1d) then stop
    IF total(sort(wb1d) ne lindgen(n_elements(wb1d))) ne 0 then stop
    IF total(sort(eb1d) ne lindgen(n_elements(eb1d))) ne 0 then stop
    IF total(sort(sb1d) ne lindgen(n_elements(sb1d))) ne 0 then stop
    IF total(sort(nb1d) ne lindgen(n_elements(nb1d))) ne 0 then stop
  ENDIF
;
  ncellx = n_elements(wb1d)
  ncelly = n_elements(sb1d)
  nx = n_elements(x1d)
  ny = n_elements(y1d)
; x dimension
  wb = wb1d[*]#replicate(1, nx)
  eb = eb1d[*]#replicate(1, nx)
  x = replicate(1, ncellx)#x1d[*]
  testx = (temporary(wb) LE x)*(temporary(eb) GT x)
  x = -1 ; free memory
  testx = where(temporary(testx) EQ 1)
  IF n_elements(testx) NE nx THEN BEGIN
    print, 'One point is out of the grid'
    stop
  ENDIF
  testx = temporary(testx) MOD ncellx
; y dimension
  sb = sb1d[*]#replicate(1, ny)
  nb = nb1d[*]#replicate(1, ny)
  y = replicate(1, ncelly)#y1d[*]
  testy = (temporary(sb) LE y)*(temporary(nb) GT y)
  y = -1 ; free memory
  testy = where(temporary(testy) EQ 1)
  IF n_elements(testy) NE ny THEN BEGIN
    print, 'One point is out of the grid'
    stop
  ENDIF
  testy = temporary(testy) MOD ncelly

  IF keyword_set(output2d) THEN return, [transpose(testx), transpose(testy)] $
  ELSE RETURN, long(testx)+ncellx*long(testy)

END