PRO SKY_IN_BANDS, bandlist=bandlist, banddirnames=banddirnames, bandnames=bandnames, bandunits=bandunits, polarised=polarised, double=double, $ group_gal=group_gal, group_strong_ps=group_strong_ps, group_faint_ps=group_faint_ps COMMON psm, param COMMON obs, obser IF NOT DEFINED(polarised) THEN polarised = SKY_POLARISED() IF NOT DEFINED(double) THEN double = PRECISION_ISDOUBLE() ;; Find existing components, select those to be observed ;;------------------------------------------------------- possible_comp = PSM_ALLCOMP() npc = N_ELEMENTS(possible_comp) existing_comp = EXISTING_COMPONENT_LIST(ncomp=nec) complist = existing_comp ncomp = N_ELEMENTS(complist) IF nec EQ 0 THEN BEGIN PSM_INFO_MESSAGE, 'No existing component in sky model' RETURN ENDIF tobe_observed = BYTARR(nec) & tobe_observed[*] = 1B PSM_DEBUG_MESSAGE, 'Fix list of observed and coadded components', routine='SKY_IN_BANDS', infolimit=3, stoplimit=3 ;; Check whether band integration has already been done ;;------------------------------------------------------ nbands = N_ELEMENTS(bandlist) tobedone = BYTARR(nbands,ncomp)+1B IF param.tasks.sky_task EQ 'restore' THEN BEGIN tobedone = BYTARR(nbands,ncomp) tocheck = complist IF KEYWORD_SET(group_gal) THEN BEGIN FOR i=0, ncomp-1 DO BEGIN IF IS_GAL_COMP(complist[i]) THEN tocheck[i]='galaxy' ENDFOR ENDIF IF KEYWORD_SET(group_faint_ps) THEN BEGIN FOR i=0, ncomp-1 DO BEGIN IF IS_FAINTPS_COMP(complist[i]) THEN tocheck[i]='faintps' ENDFOR ENDIF IF KEYWORD_SET(group_strong_ps) THEN BEGIN FOR i=0, ncomp-1 DO BEGIN IF IS_STRONGPS_COMP(complist[i]) THEN tocheck[i]='strongps' ENDFOR ENDIF FOR i=0, N_ELEMENTS(bandlist)-1 DO BEGIN newband = *(bandlist[i]) subdir = GETPSMDIR('skyinbands')+banddirnames[i] ;; Select polarisation options skypol = SKY_POLARISED() polar = polarised AND skypol ;; Units to use in the coaddition bandunit = bandunits[i] newinfo = {polar:polar,bandunit:bandunit,complist:complist,double:double,group_gal:group_gal,group_strong_ps:group_strong_ps,group_faint_ps:group_faint_ps} ;; Test whether the band directory exists, the band is the ;; same, and the component has already been 'bandintegrated' ok = FILE_TEST(FIX_SEPARATOR(subdir)+'band.sav') IF NOT ok THEN tobedone [i,*] = 1B ELSE BEGIN RESTORE, file=FIX_SEPARATOR(subdir)+'band.sav' IF NOT BANDS_EQUAL(newband,band) THEN tobedone[i,*] = 1B ELSE BEGIN ok = FILE_TEST(FIX_SEPARATOR(subdir)+'band.sav') IF NOT ok THEN tobedone[i,*] = 1B ELSE BEGIN RESTORE, file=FIX_SEPARATOR(subdir)+'info.sav' IF NOT STRUCT_EQUAL(newinfo, info) THEN BEGIN tobedone[i,*] = 1B ;; empty the directory FILE_DELETE, subdir, /quiet, /recursive ENDIF ELSE BEGIN ;; Loop on components FOR j=0, ncomp-1 DO BEGIN compname = complist[j] neglect = 0B PSM_DEBUG_MESSAGE, 'Revise neglecting of components', routine='SKY_IN_BANDS', infolimit=2 IF PSM_CAREFULNESS() LT 1 THEN neglect = COMPONENT_NEGLIGIBLE(compname,band) IF NOT neglect THEN BEGIN obscompmapname = FIX_SEPARATOR(subdir) + tocheck[j] + '_map_' + bandnames[i] + '.fits' ok = FILE_TEST(obscompmapname) IF IS_STRONGPS_COMP(compname) THEN BEGIN IF NOT ok THEN IF (obser.what.strong_sources_to_map EQ 'yes') THEN tobedone[i,j] = 1B $ ELSE BEGIN PSM_DEBUG_MESSAGE, 'Check that the component is in the file', routine='SKY_IN_BANDS', infolimit=2, stoplimit=2 ;; Check that the component is in the ;; map (using CMIX in header) ENDELSE obscompcatname = FIX_SEPARATOR(subdir) + tocheck[j] + '_cat_' + bandnames[i] + '.sav' ok = FILE_TEST(obscompcatname) IF NOT ok THEN IF (obser.what.strong_sources_to_cat EQ 'yes') THEN tobedone[i,j] = 1B $ ELSE BEGIN ;; Check that the component is in the ;; catalogue (using CMIX in header) ENDELSE ENDIF ELSE BEGIN IF NOT ok THEN tobedone[i,j] = 1B $ ELSE BEGIN ;; Check that the component is in the ;; map (using CMIX in header) ENDELSE ENDELSE ENDIF ENDFOR ENDELSE ENDELSE ENDELSE DESTROY_DATA, band ENDELSE ENDFOR ENDIF IF MAX(tobedone) EQ 0 THEN GOTO, closing ;;----------------------------------------------------------------------------------------------------------- ;; Create beams on sky for strong sources ;;----------------------------------------------------------------------------------------------------------- ;; This needs to be done only if maps of strong sources are produced ;;------------------------------------------------------------------- CASE obser.what.strong_sources_to_map OF 'yes': do_strong_sources = 1B 'no': do_strong_sources = 0B ENDCASE IF do_strong_sources THEN BEGIN skyfwhm = SKY_RESOLUTION() ;; Infrared ;;---------- strongirps_pointing_filename = PSM_TMPFILENAME(exten='.sav') whir = WHERE(complist EQ 'strongirps', nwhir) IF nwhir NE 0 THEN BEGIN RESTORE_COMP, comp, compname='strongirps' polar = polarised AND (comp.polarised) GAUSSIAN_BEAM_ON_SKY, skyfwhm, !dpi/2d0-(*comp.lat.value), (*comp.lon.value), listpix, shape, nside=SKY_NSIDE(), polar=polar SAVE, file=strongirps_pointing_filename, listpix, shape DESTROY_DATA, listpix DESTROY_DATA, shape DESTROY_DATA, comp ENDIF ;; Radio ;;------- strongradiops_pointing_filename = PSM_TMPFILENAME(exten='.sav') whrad = WHERE(complist EQ 'strongradiops', nwhrad) IF nwhrad NE 0 THEN BEGIN RESTORE_COMP, comp, compname='strongradiops' polar = polarised AND (comp.polarised) GAUSSIAN_BEAM_ON_SKY, skyfwhm, !dpi/2d0-(*comp.lat.value), (*comp.lon.value), listpix, shape, nside=SKY_NSIDE(), polar=polar SAVE, file=strongradiops_pointing_filename, listpix, shape DESTROY_DATA, listpix DESTROY_DATA, shape DESTROY_DATA, comp ENDIF ;; uchii ;;------- stronguchii_pointing_filename = PSM_TMPFILENAME(exten='.sav') whhii = WHERE(complist EQ 'stronguchii', nwhhii) IF nwhhii NE 0 THEN BEGIN RESTORE_COMP, comp, compname='stronguchii' polar = polarised AND (comp.polarised) GAUSSIAN_BEAM_ON_SKY, skyfwhm, !dpi/2d0-(*comp.lat.value), (*comp.lon.value), listpix, shape, nside=SKY_NSIDE(), polar=polar SAVE, file=stronguchii_pointing_filename, listpix, shape DESTROY_DATA, listpix DESTROY_DATA, shape DESTROY_DATA, comp ENDIF ;; wmap ;;------- strongwmapps_pointing_filename = PSM_TMPFILENAME(exten='.sav') whwmap = WHERE(complist EQ 'strongwmapps', nwhwmap) IF nwhwmap NE 0 THEN BEGIN RESTORE_COMP, comp, compname='strongwmapps' polar = polarised AND (comp.polarised) GAUSSIAN_BEAM_ON_SKY, skyfwhm, !dpi/2d0-(*comp.lat.value), (*comp.lon.value), listpix, shape, nside=SKY_NSIDE(), polar=polar SAVE, file=strongwmapps_pointing_filename, listpix, shape DESTROY_DATA, listpix DESTROY_DATA, shape DESTROY_DATA, comp ENDIF ;; ercsc ;;------- strongercscps_pointing_filename = PSM_TMPFILENAME(exten='.sav') whercsc = WHERE(complist EQ 'strongercscps', nwhercsc) IF nwhercsc NE 0 THEN BEGIN RESTORE_COMP, comp, compname='strongercscps' polar = polarised AND (comp.polarised) elliptic = (comp.elliptic) IF elliptic EQ 0 THEN BEGIN GAUSSIAN_BEAM_ON_SKY, skyfwhm, !dpi/2d0-(*comp.lat.value), (*comp.lon.value), listpix, shape, nside=SKY_NSIDE(), polar=polar ENDIF ELSE BEGIN ELLIPTIC_GAUSSIAN_BEAM_ON_SKY, ((*comp.el_fwhm_maj)+skyfwhm), (*comp.el_fwhm_min)+skyfwhm, $ !dpi/2d0-(*comp.lat.value), (*comp.lon.value), (*comp.el_theta), listpix, shape, nside=SKY_NSIDE(), polar=polar ENDELSE SAVE, file=strongercscps_pointing_filename, listpix, shape DESTROY_DATA, listpix DESTROY_DATA, shape DESTROY_DATA, comp ENDIF ENDIF ;;----------------------------------------------------------------------------------------------------------- ;; Loop on bands ;;----------------------------------------------------------------------------------------------------------- FOR i=0, N_ELEMENTS(bandlist)-1 DO BEGIN IF MAX(tobedone[i,*]) EQ 0 THEN GOTO, skipthisband band = *(bandlist[i]) ;; Create directory if it does not exist subdir = GETPSMDIR('skyinbands')+banddirnames[i] exists = FILE_TEST(subdir,/directory) IF NOT exists THEN FILE_MKDIR, subdir ;; Save the band IF PSM_VERBOSITY() GE 1 THEN SAVE, file=FIX_SEPARATOR(subdir)+'band.sav', band, /verbose $ ELSE SAVE, file=FIX_SEPARATOR(subdir)+'band.sav', band ;; Select polarisation options skypol = SKY_POLARISED() polar = polarised AND skypol ;; Units to use in the coaddition bandunit = bandunits[i] info = {polar:polar,bandunit:bandunit,complist:complist,double:double,group_gal:group_gal,group_strong_ps:group_strong_ps,group_faint_ps:group_faint_ps} IF PSM_VERBOSITY() GE 1 THEN SAVE, file=FIX_SEPARATOR(subdir)+'info.sav', info, /verbose $ ELSE SAVE, file=FIX_SEPARATOR(subdir)+'info.sav', info ;;--------------------- ;; Loop on components | ;;--------------------- FOR j=0, ncomp-1 DO BEGIN IF tobedone[i,j] EQ 0B THEN GOTO, skipthiscomp compname = complist[j] neglect = 0B PSM_DEBUG_MESSAGE, 'Revise neglecting of components', routine='SKY_IN_BANDS' IF PSM_CAREFULNESS() LT 1 THEN neglect = COMPONENT_NEGLIGIBLE(compname,band) IF NOT neglect THEN BEGIN IF COMP_IS_CAT(compname) THEN BEGIN compol = ISPOLARISED(compname) ispol = (compol AND polarised) RESTORE_COMP, comp, compname=compname npol = 1+2*ispol channelbeam = SKY_RESOLUTION() CASE STRMID(compname,0,5) OF ;; strong sources ------------------------------------------------------------------ ;; draw them directly on the map, or make a catalogue, or both 'stron': BEGIN CASE obser.what.strong_sources_to_map OF 'yes': dopsmap = 1B 'no': dopsmap = 0B END CASE obser.what.strong_sources_to_cat OF 'yes': dopscat = 1B 'no': dopscat = 0B END ;; Catalogue of strong sources ;;------------------------------ IF dopscat THEN BEGIN catfilename = FIX_SEPARATOR(subdir) + compname + '_cat_' + bandnames[i] + '.sav' IF ispol THEN datacat = MODEL_TO_DATA_CAT(comp, band, /polar) ELSE datacat = MODEL_TO_DATA_CAT(comp, band, /nopolar) IF ispol THEN BEGIN result = [[*datacat.obs1.iflux],[*datacat.obs1.qflux],[*datacat.obs1.uflux]] ENDIF ELSE BEGIN result = [*datacat.obs1.iflux] ENDELSE WRITE_CATALOG, datacat, catfilename, /sav, /name DESTROY_DATA, datacat ENDIF ELSE BEGIN IF ispol THEN OBSERVE_MODEL_CAT, comp, band, result ELSE OBSERVE_MODEL_CAT, comp, band, result, /nopolar ENDELSE ;; Map of strong sources ;;----------------------- IF dopsmap THEN BEGIN ok = EXECUTE('RESTORE, file='+compname+'_pointing_filename, /verbose') nwhstrong = N_ELEMENTS(shape) IF PRECISION_ISDOUBLE() THEN strongpsmap = DBLARR(SKY_NSIDE()^2*12L,npol) ELSE strongpsmap = FLTARR(SKY_NSIDE()^2*12L,npol) FOR isource=0L, nwhstrong-1 DO BEGIN lis = *(listpix[isource]) shpi = result[isource] * (*(shape[isource])).amplitude if finite((*(shape[isource])).amplitude[0]) eq 0 then stop strongpsmap[lis,0] = strongpsmap[lis,0] + shpi ;; check dim of result IF ispol THEN BEGIN ang = (*(shape[isource])).rotation + (*comp.pol_angle)[isource] shpi = shpi*(*comp.pol_fraction)[isource] shpq = shpi*COS(2*ang) shpu = shpi*SIN(2*ang) strongpsmap[lis,1] = strongpsmap[lis,1] + shpq strongpsmap[lis,2] = strongpsmap[lis,2] + shpu ENDIF ENDFOR DESTROY_DATA, listpix DESTROY_DATA, shape ;; convert map to channel units: map in mJy/sr, convert to band units fac = CONVERSION_FACTOR('mJy/sr',bandunit, band, double=double) strongpsmap = TEMPORARY(strongpsmap)*fac compdata = PTR_NEW(strongpsmap,/no_copy) ENDIF END ;; faint sources ------------------------------------------------------------------- ;; put them on the map and do harmonic space convolution 'faint': BEGIN IF ispol THEN OBSERVE_MODEL_CAT, comp, band, result, /approx ELSE OBSERVE_MODEL_CAT, comp, band, result, /approx, /nopolar nsources = N_ELEMENTS(result[*,0]) IF PRECISION_ISDOUBLE() THEN faintpsmap = DBLARR(SKY_NSIDE()^2*12L,npol) ELSE faintpsmap = FLTARR(SKY_NSIDE()^2*12L,npol) ANG2PIX_RING, SKY_NSIDE(), !dpi/2d0 - *comp.lat.value, *comp.lon.value, ipix IF KEYWORD_SET(double) NE 0 THEN fourpi = 4d0*!dpi ELSE fourpi=4.*!pi omegapix = fourpi / NSIDE2NPIX(SKY_NSIDE()) FOR is=0L, nsources-1L DO faintpsmap[ipix[is],*] = faintpsmap[ipix[is],*] + result[is,*] fac = CONVERSION_FACTOR('mJy/sr',bandunit,band, double=double) fac = fac/omegapix faintpsmap = TEMPORARY(faintpsmap)*fac CHANGE_MAP_RESOLUTION, TEMPORARY(faintpsmap), newfaintpsmap, beam_in=GET_BEAM_STRUCT('GAUSSIAN',beamsize=0.,sizeunits='arcmin'), $ beam_out=GET_BEAM_STRUCT('GAUSSIAN',beamsize=SKY_RESOLUTION(),sizeunits='arcmin'), /nohdr, lmax=SKY_LMAX() compdata = PTR_NEW(newfaintpsmap,/no_copy) END ENDCASE DESTROY_DATA, comp ENDIF ELSE BEGIN ispol = ISPOLARISED(compname) COMP_IN_BAND, compname, band, compdata, units=bandunit, dtshp='array', dtfm='map', $ /pointer, polarised=ispol, hpx_nside=nside, lmax=lmax, double=double ENDELSE IF PTR_VALID(compdata) THEN BEGIN dimsout = SIZE(*compdata, /dimensions) ndimout = SIZE(*compdata, /n_dimensions) ;; Write individual component observation ;;---------------------------------------- IF tobe_observed[j] EQ 1 THEN BEGIN block1 = PSM_BASE_HDR(datatype='OBS',dataform='MAP',/setval) block2 = PSM_BAND_HDR(band) block3 = PSM_BEAM_HDR(SKY_RESOLUTION()) block4 = PSM_MAP_HDR(pixtype='HEALPIX',lmax=SKY_LMAX()) block5 = PSM_OBS_HDR(compname, obsid=obsid) ;block6 = PSM_SKYMIX_HDR(compname) MK_PSMHDR, xhdr, block1, block2, block3, block4, block5 ;, block6 ;; write the map in appropriate format obscompmapname = FIX_SEPARATOR(subdir) + compname + '_map_' + bandnames[i] + '.fits' skypix = SKY_PIX() IF ispol THEN WRITE_HPXMAP, obscompmapname, tqustokes=compdata, xhdr=xhdr, units=bandunit, coordsys=skypix.coordsys $ ELSE WRITE_HPXMAP, obscompmapname, tstokes=compdata, xhdr=xhdr, units=bandunit, coordsys=skypix.coordsys ENDIF ;; Free the pointer PTR_FREE, compdata ENDIF ENDIF skipthiscomp: ENDFOR ; (end loop on components) IF param.tasks.sky_task NE 'restore' THEN BEGIN IF KEYWORD_SET(group_gal) NE 0 THEN GROUP_GAL_INBAND, subdir;, /no_delete IF KEYWORD_SET(group_strong_ps) NE 0 THEN GROUP_STRONG_PS_INBAND, subdir;, /no_delete IF KEYWORD_SET(group_faint_ps) NE 0 THEN GROUP_FAINT_PS_INBAND, subdir;, /no_delete ENDIF ELSE BEGIN IF KEYWORD_SET(group_gal) NE 0 THEN BEGIN ok = FILE_TEST(FIX_SEPARATOR(subdir)+'galaxy_map_'+bandnames[i]+'.fits') IF NOT ok THEN GROUP_GAL_INBAND, subdir ENDIF IF KEYWORD_SET(group_strong_ps) NE 0 THEN BEGIN ok = FILE_TEST(FIX_SEPARATOR(subdir)+'strongps_map_'+bandnames[i]+'.fits') IF NOT ok THEN GROUP_STRONG_PS_INBAND, subdir ENDIF IF KEYWORD_SET(group_faint_ps) NE 0 THEN BEGIN ok = FILE_TEST(FIX_SEPARATOR(subdir)+'faintps_map_'+bandnames[i]+'.fits') IF NOT ok THEN GROUP_FAINT_PS_INBAND, subdir ENDIF ENDELSE skipthisband: ;; Clean memory leaks ;;-------------------- PSM_DEBUG_MESSAGE, 'Find memory leaks', routine='SKY_IN_BANDS' HEAP_GC ;; Make figures ;;-------------- ENDFOR ; (end loop on bands) closing: ;; Clean-up temporary pointing files for point sources ;;----------------------------------------------------- IF DEFINED(strongradiops_pointing_filename) THEN FILE_DELETE, strongradiops_pointing_filename, /quiet, /allow_nonexistent IF DEFINED(strongirps_pointing_filename) THEN FILE_DELETE, strongirps_pointing_filename, /quiet, /allow_nonexistent IF DEFINED(stronguchii_pointing_filename) THEN FILE_DELETE, stronguchii_pointing_filename, /quiet, /allow_nonexistent IF DEFINED(strongwmapps_pointing_filename) THEN FILE_DELETE, strongwmapps_pointing_filename, /quiet, /allow_nonexistent IF DEFINED(strongercscps_pointing_filename) THEN FILE_DELETE, strongercscps_pointing_filename, /quiet, /allow_nonexistent END