;+ ; NAME: ; msrcor ; PURPOSE: (one line only) ; Spatial correlation souce positions found in multiple lists ; DESCRIPTION: ; This routine is a conceptual extension of the Astronomy Users Library ; routine, srcor.pro. This extension of this routine is that you ; can correlate more than two lists of positions. ; Each call of this routine will serve to add one list to the description ; of the union of all the prior lists. The lists are copied into the ; information collected about the set of lists so you don't need to ; maintain a copy of the original lists. ; CATEGORY: ; Photometry ; CALLING SEQUENCE: ; msrcor,set,x,y,dcr ; INPUTS: ; set - This is the description of the set of all lists seen thus far. ; There are three ways to signal the first call to this routine. ; 1) set is undefined ; 2) set is not a structure ; 3) set is a valid structure but set.nlists is equal to zero. ; Note that there is absolutely no performance advantage to maintaining ; a previous structure and then setting set.nlists=0. If you want ; to restart it is just as effective to simply say set=0 prior to ; an initilizing call. ; If the structure exists already and set.nlists!=0 then the information ; will be added to and the anonymous structure will be modified ; before returning. ; x,y - Array of x and y coordinates for a new list to add to the set. ; The order in which these array are successively presented to msrcor ; will determine the ordinal number of the list within the set. ; dcr - Critical radious outside which correlations are rejected. ; x,y,dcr can be in any units that make sense to your problem but ; they must all be the same units. ; OPTIONAL INPUT PARAMETERS: ; KEYWORD INPUT PARAMETERS: ; OUTPUTS: ; set - Information structure, the following tags are defined: ; x - X position of source (all lists are concatenated here) ; y - Y position of source (all lists are concatenated here) ; objid - vector that matches the length of x and y. These values ; are unique object id numbers, starting at 0. ; objcnt - vector that indicates how many times each object ; appears in the lists (<= nlists) ; lidx - vector (same length as x,y) that indicates which list ; the position appears in. ; nobj - scalar, number of unique objects found. ; nlists - scalar, total number of lists included. ; ; There may be other related information you might like to keep along with ; this information, such as magnitude or FWHM. If you concatentate ; everything in list order you will have a vector whose structure matches ; set.x and set.y and the indexing for all will be the same. ; KEYWORD OUTPUT PARAMETERS: ; COMMON BLOCKS: ; SIDE EFFECTS: ; RESTRICTIONS: ; PROCEDURE: ; MODIFICATION HISTORY: ; 2008/06/18, Written by Marc W. Buie with algorithmic input from Leslie Young, ; Southwest Research Institute. ; 2011/08/02, MWB, promoted all counters and pointers to long, this caused ; some odd problems when linking a lot of lists. ;- pro msrcor,set,in_x,in_y,dcr,DEBUG=debug self='MSRCOR: ' if badpar(set,[0,1,2,3,8],[0,1],caller=self+'(set) ',type=stype) then return if badpar(in_x,[2,3,4,5],1,caller=self+'(x) ',npts=nx) then return if badpar(in_y,[2,3,4,5],1,caller=self+'(y) ',npts=ny) then return if badpar(dcr,[2,3,4,5],0,caller=self+'(dcr) ') then return if badpar(debug,[0,1,2,3],0,caller=self+'(DEBUG) ',default=0) then return if nx ne ny then begin print,self,' Error! Input x and y vectors must be the same length. Aborting.' return endif if stype eq 8 then begin if set.nlists eq 0 then newstart=1 else newstart=0 endif else begin newstart=1 endelse if debug then begin print,'Newstart=',strn(newstart) endif if newstart then begin x = in_x y = in_y objid = lindgen(nx) objcnt = replicate(1L,nx) lidx = replicate(0L,nx) nobj = nx nlists = 1L endif else begin if debug then begin help,set endif ; Here's the real work in this program. ; Step through the lists already in the set. For each list, only ; consider the objects in the new set that have NOT been matched up ; yet. Of course, that means we always check the entirety of the ; new list against the very first list in the set. After that ; the amount checked against subsequent lists should dwindle (or ; disappear completely). Anything that isn't matched after all old ; lists are checked become new object ids that have no links. ; This array keeps track of the object id number as we proceed. ; The starting value indicates this object hasn't been linked up yet. objid = replicate(-1L,nx) for i=0,set.nlists-1 do begin ; Find the unmatched objects from the new list. Exit loop if none. znew = where(objid eq -1,countnew,/NULL) if countnew eq 0 then break ; list 1 (from new list) x1 = in_x[znew] y1 = in_y[znew] ; list 2 (from old lists) z = where(set.lidx eq i,count,/NULL) x2 = set.x[z] y2 = set.y[z] srcor,x1,y1,x2,y2,dcr,ind1,ind2,option=1 if ind1[0] ne -1 then $ objid[znew[ind1]] = set.objid[z[ind2]] endfor ; deal with the unmatched objects at the end objcnt = set.objcnt z=where(objid eq -1,count,/NULL) if count ne 0 then begin objid[z] = lindgen(count)+set.nobj objcnt = [objcnt,replicate(0,count)] endif if debug then begin print,strn(count),' new unmatched objects after checking old lists.' help,set print,'objcnt[0]',objcnt[0] print,'objid',objid endif objcnt[objid] += 1 x = [set.x,in_x] y = [set.y,in_y] objid = [set.objid,objid] lidx = [set.lidx,replicate(set.nlists,nx)] nobj = set.nobj + count nlists = set.nlists+1 endelse set = { $ x: x, $ y: y, $ objid: objid, $ objcnt: objcnt, $ lidx: lidx, $ nobj: nobj, $ nlists: nlists $ } end