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