function qmatch, ra1, dec1, ra2, dec2, dmatch ; inputs ra1, dec1, ra2, dec2 in degrees; dmatch in arcseconds ; this catalog-matching scheme is designed to be quick ; and should not be used as a general neighbor-finding routine. ; hardwired slowdown at dec=85 deg (currently inaccurate; need to fix) ; set up match1 array and offset1 array (init to NaN) nobj1 = n_elements(ra1) match1 = lonarr(nobj1) - 1l offset1 = fltarr(nobj1) -1. & offset1[*] = !VALUES.F_NAN ; calculate ra offsets from dec1 scaled by cos dec1 ddtor = !DPI / 180d dramatch = dmatch / (3600. * (cos(dec1 * ddtor) > 0.1) ) ; bin each dec array into dmatch-sized chunks dmatchdeg = dmatch / 3600. hdec1 = histogram(dec1, min=-90d, max=90d, binsize=dmatchdeg, rev=rev1) hdec2 = histogram(dec2, min=-90d, max=90d, binsize=dmatchdeg, rev=rev2) ; find where hdec1 contains objects and hdec2 contains objects ; within 3-chunk box centered on hdec1,2 bins ; this step assures detection of matches to companions at bin edges checklist = where(smooth(float(hdec2),3,/edge_truncate) * hdec1 ne 0., count) if count eq 0 then return, {index:match1, offset:offset1, nmatch:0L} ; step through checklist and search only objects within dec bin ncheck = n_elements(checklist) ncheck1 = ncheck - 1 for i = 0l, ncheck1 do begin icheck = checklist[i] ; who's here from list 1 j1 = rev1[rev1[icheck] : rev1[icheck+1] - 1] n1 = n_elements(j1) ; who's here from list 2 (again check one up to one down) icheckup = icheck + 1 icheckup2 = icheck + 2 icheckdown = icheck - 1 if rev2[icheck] ne rev2[icheckup] then begin j2 = rev2[rev2[icheck] : rev2[icheckup] - 1] endif if icheck ne ncheck1 then begin if rev2[icheckup] ne rev2[icheckup2] then begin j2 = (size(j2))[0] ? [j2, rev2[rev2[icheckup] : rev2[icheckup2] - 1]] \$ : [rev2[rev2[icheckup] : rev2[icheckup2] - 1]] endif endif if icheck ne 0 then begin if rev2[icheckdown] ne rev2[icheck] then begin j2 = (size(j2))[0] ? [j2, rev2[rev2[icheckdown] : rev2[icheck] - 1]] \$ : [rev2[rev2[icheckdown] : rev2[icheck] - 1]] endif endif n2 = n_elements(j2) ; if n2 ne 1 then print, icheck, n2 ; check for candidates that have dra < dramatch, get true angular separation ; in arcseconds for those that do. for i1 = 0L, n1 - 1 do begin dra = ra2[j2] - ra1[j1[i1]] candidates = where(abs(dra) le dramatch[j1[i1]], count) ; print, count if count ne 0 then begin ang = djs_diff_angle(ra2[j2[candidates]], dec2[j2[candidates]], \$ ra1[j1[i1]], dec1[j1[i1]], units=degrees) * 3600d ; if ang[0] ne 0d then print, ra2[j2[0]], dec2[j2[0]], \$ ; ra1[j1[i1]], dec1[j1[i1]], ang[0], format='(5f16.11)' ; keep the closest one that falls within dmatch dmin = min(ang, match2) if dmin le dmatch then begin match1[j1[i1]] = j2[candidates[match2]] offset1[j1[i1]] = ang[match2] endif endif endfor endfor matched_indices = where(match1 ne -1, count) return, {index:match1, offset:offset1, nmatch: count} end