;+ ; STACK_IN_REDSHIFT_SLICES, cmap, hd, cube, fhwm,$ ; cnoise=cnoise,mask=mask,beam_area=beam_area,$ ; err_ss=err_ss, quiet=quiet ; ; PURPOSE: Takes input map with header and cube containing layers and regresses ; with MPFITFUN to find the best stacked solution. Noise and masking ; are both optional. ; ; INPUTS: ; cmap : cropped map in Jy/beam. If in MJy/sr, must supply beam_area ; hd : header file, must contain CD2_2 or CDELT2 to determine pixel size ; cube : data cube, N sources per list, M lists, RA and DEC. Lists will have ; different number of sources, so N = longest list, and others have zero ; in place of values. ; fwhm : resolution of map. This is important to get right and a potential source ; of systematic uncertainty (i.e., consider simulating to measure potential error) ; ; OPTIONAL OUTPUTS: ; mask : Optional places where you don't want to stack ; beam_area : To convert from MJy/sr (like spitzer or BLAST) to Jy/beam (like SPIRE) ; err_ss: Errors from mpfit ; ; OUTPUTS: ; cov_ss: M stacked flux densities ; ; USAGE: ; stacked_fluxes = STACK_IN_REDSHIFT_SLICES(cmap,hd,cube,fwhm,cnoise=cnoise,mask=mask,quiet=1) ; ; CALLS: ; NAN_FINDER ; GAUSS_KERN ; CIRCLE_MASK ; MPFITFUN ; SIMULTANEOUS_STACK_ARRAY ; ; ; HISTORY: ; CREATED BY marco.viero@caltech.edu (2013-03-03) ; ADDED ROUNDING OF FLOATING POING AFTER ADXY (2014-04-07) ; ADDED MEAN OF ind_mean SUBSET OF tmap/cmap (2014-04-07) ; ADDED Second condition to if statement (nt0 gt 0) on line 99 (2014-30-06) ;- FUNCTION STACK_IN_REDSHIFT_SLICES, $ cmap, $ ; CROPPED MAP hd, $ ; NEED HEADER FOR ADXY AND PIXSIZE cube, $ ; N x M x 2 ARRAY, DEFINED BELOW fwhm, $ ; IN ARCSEC. TRY TO USE ACTUAL, NOT JUST NOMINAL cnoise=cnoise, $ ; OPTIONAL mask=mask, $ ; OPTIONAL beam_area=beam_area, $ ; INCLUDE IF YOU WANT TO CONVERT MAPS TO JY/beam err_ss=err_ss, $ ; Errors from MPFIT quiet=quiet ; QUIET=1 makes function silent ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;CHECK MAP FOR ZEROS OR NANS ind_map_zero0=where(nan_finder(cmap) eq 0,nimz) zeromask=cmap*0.+1. if keyword_set(nimz) then zeromask[ind_map_zero0]=0. if keyword_set(mask) then zeromask=zeromask*mask if not(keyword_set(quiet)) then quiet=0 ind_map_zero=where(zeromask eq 0,nzero) ;FIND SIZES OF MAP AND LISTS ms=size(cmap,/dim) size_cube=size(cube,/dim) nsrcmax=size_cube[0] nlists=size_cube[1] if keyword_set(hd) and keyword_set(pix) then begin print, 'cant have pix and header!' ;& RETURN stop endif if not(keyword_set(cnoise)) then cnoise=cmap*0.+1.0 pix=sxpar(hd, 'CD2_2')*3600. if pix eq 0 then pix=sxpar(hd, 'CDELT2')*3600. ;[STEP 0] - Calibrate maps (i.e., convert from MJy/sr to Jy/beam) if keyword_set(beam_area) then begin cmap=cmap*beam_area*1d6 cnoise=cnoise*beam_area*1d6 endif ; ; STEP 1 - Make Layers Cube layers=dblarr(nlists,ms[0],ms[1]) for s=0,nlists-1 do begin ind_src=where(cube[*,s,0] ne 0, nsrc) ra=cube[ind_src,s,0] dec=cube[ind_src,s,1] ; CONVERT FROM RA/DEC to X/Y adxy,hd,ra,dec,tx,ty ; CHECK FOR SOURCES THAT FALL OUTSIDE MAP ind_keep=where(tx ge 0 and tx lt ms[0] and ty ge 0 and ty lt ms[1],nt0) real_x=round(tx[ind_keep]) real_y=round(ty[ind_keep]) ; CHECK FOR SOURCES THAT FALL ON ZEROS MAP if keyword_set(nzero) and nt0 gt 0 then begin tally=dblarr(nt0) for d=0l,nt0-1 do $ if cmap[real_x[d],real_y[d]] ne 0 then tally[d]=1. ind_nz=where(tally eq 1,nt) real_x=real_x[ind_nz] real_y=real_y[ind_nz] endif else nt=nt0 if keyword_set(verbose) then $ print, 'of '+ strcompress(string(nt,format='(i10)'),/remove_all)+$ ' sources in list '+$ strcompress(string(nt0-nt,format='(i10)'),/remove_all)+' fall outside' for ni=0,nt-1 do $ layers[s, real_x[ni],real_y[ni]]+=1 endfor ; STEP 2 - Convolve Layers and put pixels to regress into vectors radius=1.1 kern=gauss_kern(fwhm,ms[0]