PRO cross_cats, ra1, dec1, ra2, dec2, radius, index, multi = multi, r_min = r_min

;
; PURPOSE:
;  	Identify overlap sources in two catalogs by searching source
;	pairs with separations withing a certain value.
;
; INPUTS:
;  	ra1, dec1 - double precision vectors that contain source ra
;    		and dec in the first catalog, in unit of degree.
;  	ra2, dec2 - ra, dec of the second catalog
;  	radius - searching radius, in arcsec
;  	multi - set to 1 to avoid two sources having the same counterpart.
;    		In this case, the counterpart will be assigned to the source
;    		that is closest to the counterpart.  Could be time consuming.
;
; OUTPUT:
;  	index - integer vector that has the same number of elements
;    		as the first catalog indicating which source in cat2
;    		is the counterpart.  If -1, the corresponding cat1
;    		source doesn't have a counterpart in cat2.  
;
; NOTE:  1. It doesn't work if your sources are close (<1 arcmin, maybe) 
;           to the poles.
;        2. The subroutine cross_dogs.pro is required.
;
;




If n_elements(multi) EQ 0 THEN multi = 0

n1 = n_elements(ra1)
IF n1 LE 2000 THEN cross_dogs, ra1, dec1, ra2, dec2, radius, index, multi = multi
IF n1 LE 2000 THEN goto, finish


index = lon64arr(n1)*0-1

; roughly 1000-2000 sources per block, if uniformly distributed
ncut = ceil(sqrt(ceil(n1/2000.0)))

ra_max = max(ra1)
ra_min = min(ra1)
dec_max = max(dec1)
dec_min = min(dec1)
dra = (ra_max - ra_min)/ncut
ddec = (dec_max - dec_min)/ncut

;order1 = sort(ra1)
;order2 = sort(dec1)
;dn = n1*1.0/ncut

FOR i=0,ncut-1 DO BEGIN
FOR j=0,ncut-1 DO BEGIN
   A = where(ra1 GE i*dra+ra_min AND ra1 LE (i+1)*dra+ra_min AND $
             dec1 GE j*ddec+dec_min AND dec1 LE (j+1)*ddec+dec_min)
   IF total(A) NE -1 THEN BEGIN
      cross_dogs, ra1[A], dec1[A], ra2, dec2, radius, index2, multi = multi
      index[A] = index2
   ENDIF

;   why plan B doesn't work well??
;   B = where(order1 GE i*dn AND order1 LE (i+1)*dn AND $
;             order2 GE j*dn AND order2 LE (j+1)*dn)
;   IF total(B) NE -1 THEN BEGIN
;      cross_dogs, ra1[order1[B]], dec1[order2[B]], ra2, dec2, radius, index2, multi = multi
;      index[B] = index2
;   ENDIF
ENDFOR
ENDFOR




finish:

END