;+ ; NAME: ; find_primary_field ; ; PURPOSE: ; Find the field information for the "primary" SDSS field at coordinates ; ; CALLING SEQUENCE: ; find_primary_field, ra, dec, [ ilist=, fout=, select_tags=, except_tags= ,$ ; /unified, iunified=, ibalkan= ] ; ; INPUTS: ; ra - Right ascension(s) [deg] ; dec - Declination(s) [deg] ; ; OPTIONAL INPUTS: ; select_tags- Tag names to include in output structure FOUT ; except_tags- Tag names to exclude in output structure FOUT ; unified - If set, then use the file window_unified.fits, which is ; typically faster ; ; OUTPUTS: ; ; OPTIONAL OUTPUTS: ; ilist - Field index number for primary field, line-by-line matched ; to the input RA,DEC; set to -1 for coordinates with no ; primary SDSS field ; fout - Field list information for each primary field, line-by-line ; matched to the input RA,DEC; zero values for coordinates ; with no primary SDSS field ; iunified - Balkan identifier for the unified balkan ID; set only ; if procedure called with /UNIFIED ; ibalkan - Balkan identifier for the full balkan ID; set only ; if procedure not called with /UNIFIED ; ; COMMENTS: ; This procedure uses the geometrical information from the SDSS "resolve" ; algorithm: ; https://www.sdss3.org/dr10/algorithms/resolve.php ; For each RA,DEC coordinate on the sky, find the appropriate SDSS field ; that is "primary" at that coordinate. ; ; The output structure FOUT has all of the fields found in the ; file "window_flist.fits", with the addition of AIRMASS and IMAGE_DEPTH. ; AIRMASS is computed for the timestamp of the r-band image. ; IMAGE_DEPTH is the 5-sigma point-source depth in linear units ; of nanomaggies/arcsec^2 (1.0 being equivalent to 22.5 mag/arcsec^2), ; computed according to Schlegel note on 7 Oct 2013 "A Question of Depth" ; ; On a 2013 Macbook, this takes about 100 sec to compute for 100k points ; that are not on stripe 82, longer if on stripe 82. ; Untested for more than 1 million points. ; ; EXAMPLES: ; Test which of 100k random points on the sky are in the BOSS footprint, ; and plot them. Limit to only unbinned imaging. ; IDL> n=100000L ; IDL> ra=randomu(1234,n)*360 & dec=asin(randomu(2345,n)*2-1)*180d0/!dpi ; IDL> find_primary_field,ra,dec,ilist=ilist,fout=fout ; IDL> i=where(ilist NE -1 AND fout.xbin EQ 1) ; IDL> plot,ra[i],sin(dec[i]/!radeg),xtitle='RA',ytitle='sin(Dec)',ps=3 ; ; BUGS: ; ; DATA FILES: ; $PHOTO_RESOLVE/window_flist.fits ; $PHOTO_RESOLVE/window_blist.fits ; $PHOTO_RESOLVE/window_bindx.fits ; $PHOTO_RESOLVE/window_bcaps.fits ; $PHOTO_RESOLVE/window_unified.fits ; ; PROCEDURES CALLED: ; ; REVISION HISTORY: ; 09-Feb-2014 Checked in older revision of photoop product find_primary_field.pro that was used to make dr12v4 catalogs (155612). Latest version does not run properly/it's too slow (?). Did not fully diagnose the problem. Could just be the /rescore keyword. ; ; 02-Apr-2014 Written by D. Schlegel, LBL ;- ;------------------------------------------------------------------------------ pro find_primary_field1, ra_list, dec_list, b_ra, b_dec, flist, ilist, $ iunified=iunified pfile = djs_filepath('window_unified.fits', $ root_dir=getenv('PHOTO_RESOLVE')) if (file_test(pfile) EQ 0) then $ message, 'Unified field list missing from $PHOTO_RESOLVE' read_fits_polygons, djs_filepath('window_unified.fits', $ root_dir=getenv('PHOTO_RESOLVE')), unif iunified = lonarr(n_elements(ra_list)) - 1L b_ra = flist[unif.ifield ].ra b_dec = flist[unif.ifield ].dec spherematch, ra_list, dec_list, b_ra, b_dec, 0.16, m1, m2, max=0 isort = sort(m2) iuniq = uniq(m2[isort]) ist = 0L for i=0L, n_elements(iuniq)-1L do begin ind = iuniq[i] icurr = isort[ist:ind] racurr = ra_list[m1[icurr]] deccurr = dec_list[m1[icurr]] ucurr = unif[m2[icurr[0]]] is_in = is_in_polygon(ra=racurr, dec=deccurr, ucurr) iin = where(is_in, nin) if (nin GT 0) then begin ilist[m1[icurr[iin]]] = ucurr.ifield iunified[m1[icurr[iin]]] = m2[icurr[0]] endif ist= ind+1L endfor return end ;------------------------------------------------------------------------------ pro find_primary_field2, ra_list, dec_list, b_ra, b_dec, blist, ilist, $ ibalkan=ibalkan window_read, bcaps=bcaps qexist = keyword_set(bcaps) if (qexist EQ 0) then $ message, 'Field list or balkan files missing from $PHOTO_RESOLVE' ibalkan = lonarr(n_elements(ra_list)) - 1L ; Find all balkans that might possibly overlap each input coordinate ; Hardwire 0.26 deg as the maximum distance from a coordinate to ; the primary field of a balkan. spherematch, b_ra, b_dec, ra_list, dec_list, 0.16, im1, im2, d12, maxmatch=0 nmatch = n_elements(im1) d12 = 0 ; free memory isort = sort(im2) im1 = im1[isort] im2 = im2[isort] j0 = 0L j1 = 0L nobj = n_elements(ra_list) for iobj=0L, nobj-1L do begin if ((iobj MOD 100) EQ 0) then print,iobj,nobj,string(13b),format='(i,i,a,$)' ; Get the list of balkans that might overlap this position while (im2[j0<(nmatch-1)] LT iobj AND j0 LT nmatch-1) do j0++ if (im2[j0] EQ iobj) then begin j1 = j0 while (im2[j1<(nmatch-1)] EQ iobj AND j1 LT nmatch) do j1++ j1-- ; The balkans that may overlap are indexed by im1[j0:j1] for ib=j0, j1 do begin ; Test if this object is in this particular balkan qinside = 1B For k=0, blist[im1[ib]].ncaps-1 do $ qinside *= is_in_cap(ra=ra_list[iobj], dec=dec_list[iobj], $ bcaps[blist[im1[ib]].icap+k]) if (qinside) then begin if (ilist[iobj] EQ -1) then begin ilist[iobj] = blist[im1[ib]].iprimary ibalkan[iobj] = im1[ib] endif else begin print & print, 'Warning: object in multiple balkans' endelse endif endfor endif j0 = j1 ; up the counters for the next object endfor blist = 0 ; free memory bcaps = 0 ; free memory return end ;------------------------------------------------------------------------------ pro find_primary_field, ra_list, dec_list, ilist=ilist, $ fout=fout, select_tags=select_tags, except_tags=except_tags, unified=unified, $ iunified=iunified, ibalkan=ibalkan nobj = n_elements(ra_list) if (nobj EQ 0 OR nobj NE n_elements(dec_list)) then $ message, 'Number of elements in RA,DEC must agree' print, 'Using PHOTO_RESOLVE='+getenv('PHOTO_RESOLVE') window_read, flist=flist, blist=blist, bindx=bindx qexist = keyword_set(flist) * keyword_set(blist) * keyword_set(bindx) if (qexist EQ 0) then $ message, 'Field list or balkan files missing from $PHOTO_RESOLVE' ; Add AIRMASS to the structure, computing for the r-band image timestamp, ; and IMAGE_DEPTH, the 5-sigma point source depth addtags = replicate(create_struct('AIRMASS', 0., $ 'IMAGE_DEPTH', fltarr(5)), n_elements(flist)) addtags.airmass = tai2airmass(flist.ra, flist.dec, tai=flist.tai[2]) anorm = [0.387, 0.218, 0.241, 0.297, 0.665] kterm = [0.49, 0.17, 0.10, 0.06, 0.06] ; From arxiv.org/abs/astro-ph/0703454) for i=0,4 do $ addtags.image_depth[i] = anorm[i] * flist.psf_fwhm[i] $ * sqrt(flist.skyflux[i]) $ * 10.^(0.4*kterm[i]*tai2airmass(flist.ra, flist.dec, tai=flist.tai[i])) flist = struct_addtags(flist, addtags) ; Output will be primary field number for each input coordinate, ; set to -1 if there is none ilist = lonarr(n_elements(ra_list)) - 1L ; Get the coordinates of the primary field for each balkan b_ra = flist[ bindx[blist.ibindx].ifield ].ra b_dec = flist[ bindx[blist.ibindx].ifield ].dec bindx = 0 ; free memory if (keyword_set(unified)) then begin find_primary_field1, ra_list, dec_list, b_ra, b_dec, flist, ilist, $ iunified=iunified endif else begin find_primary_field2, ra_list, dec_list, b_ra, b_dec, blist, ilist, $ ibalkan=ibalkan endelse blist = 0 ; free memory if (arg_present(fout)) then begin ; First trim the field list down in size to save memory if (keyword_set(select_tags) OR keyword_set(except_tags)) then $ flist = struct_selecttags(flist, select_tags=select_tags, $ except_tags=except_tags) fout = struct_selecttags(flist[ilist>0], select_tags=select_tags, $ except_tags=except_tags) ; Set output structure to all zeros when there is no primary field ibad = where(ilist EQ -1, nbad) if (nbad GT 0) then begin blankdat = fout[0] struct_assign, {junk:0}, blankdat fout[ibad] = blankdat endif endif print return end ;------------------------------------------------------------------------------