; ----------------------------------------------------------------------------- ; ; Copyright (C) 2010-2011 ; Jacques Delabrouille ; Guillaume Castex ; ; This file is part of the Planck Sky Model (PSM). ; ; The Planck Sky Model is free software: you can redistribute it and/or modify ; it under the terms of the GNU General Public License as published by ; the Free Software Foundation, either version 3 of the License, or ; (at your option) any later version. ; ; The Planck Sky Model is distributed in the hope that it will be useful, ; but WITHOUT ANY WARRANTY; without even the implied warranty of ; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ; GNU General Public License for more details. ; ; You should have received a copy of the GNU General Public License ; along with the Planck Sky Model. If not, see . ; ; For more information about the PSM see: ; http://www.apc.univ-paris7.fr/~delabrou/PSM/psm.html ; or contact Jacques Delabrouille (delabrouille@apc.univ-paris7.fr) ; ; ----------------------------------------------------------------------------- ;+ ;- PRO PSM_WMAP_SOURCES, data_cat, model_cat, polarised=polarised COMMON psm, param gp = param.global runp = param.run names = param.names wmappscat = FIX_SEPARATOR(GETPSMDIR('psmdata')+'ancillary/observations/WMAP/PS-catalog')+'wmap_cat.txt' ok = FILE_TEST(wmappscat) IF NOT ok THEN GET_PSM_DATA, datafile=wmappscat READCOL,wmappscat,ra,dec,lon_w,lat_w,flux_k, flux_ka, flux_q, flux_v, flux_w,err_k,err_ka,err_q,err_v,err_w, /silent dist=11*!dpi/(180*60.) nsrc = N_ELEMENTS(lon_w) ;;create wmap catalog emlaw_list = ['POWERLAW','POWERLAW','POWERLAW','POWERLAW','POWERLAW','POWERLAW','POWERLAW'] wmap_cat = CREATE_PS_MODEL_CAT(nsrc, emlaw_list, name='strong'+names.wmapps_basename, id=id, polarised=polarised) *wmap_cat.e1.nuref.value = (*(*model_cat).e1.nuref.value) *wmap_cat.e1.numin.value = (*(*model_cat).e1.numin.value) *wmap_cat.e1.numax.value = 4.85e9 (*wmap_cat.e1.specind.value)=dblarr(nsrc) (*wmap_cat.e1.ampl.value)=dblarr(nsrc) (wmap_cat.e1.ampl.unit)='Jy' *wmap_cat.e2.nuref.value = 4.85e9 *wmap_cat.e2.numin.value = 4.85e9 *wmap_cat.e2.numax.value = 23e9 (*wmap_cat.e2.specind.value)=dblarr(nsrc) (*wmap_cat.e2.ampl.value)=dblarr(nsrc) (wmap_cat.e2.ampl.unit)='Jy' *wmap_cat.e3.nuref.value = 23e9 *wmap_cat.e3.numin.value = 23e9 *wmap_cat.e3.numax.value = 33e9 (*wmap_cat.e3.specind.value)=dblarr(nsrc) (*wmap_cat.e3.ampl.value)=dblarr(nsrc) (wmap_cat.e3.ampl.unit)='Jy' *wmap_cat.e4.nuref.value = 33e9 *wmap_cat.e4.numin.value = 33e9 *wmap_cat.e4.numax.value = 41e9 (*wmap_cat.e4.specind.value)=dblarr(nsrc) (*wmap_cat.e4.ampl.value)=dblarr(nsrc) (wmap_cat.e4.ampl.unit)='Jy' *wmap_cat.e5.nuref.value = 41e9 *wmap_cat.e5.numin.value = 41e9 *wmap_cat.e5.numax.value = 61e9 (*wmap_cat.e5.specind.value)=dblarr(nsrc) (*wmap_cat.e5.ampl.value)=dblarr(nsrc) (wmap_cat.e5.ampl.unit)='Jy' *wmap_cat.e6.nuref.value = 61e9 *wmap_cat.e6.numin.value = 61e9 *wmap_cat.e6.numax.value = 94e9 (*wmap_cat.e6.specind.value)=dblarr(nsrc) (*wmap_cat.e6.ampl.value)=dblarr(nsrc) (wmap_cat.e6.ampl.unit)='Jy' *wmap_cat.e7.nuref.value = 94e9 *wmap_cat.e7.numin.value = 94e9 *wmap_cat.e7.numax.value = (*(*model_cat).e4.numax.value) (*wmap_cat.e7.specind.value)=dblarr(nsrc) (*wmap_cat.e7.ampl.value)=dblarr(nsrc) (wmap_cat.e7.ampl.unit)='Jy' candidates=lonarr(nsrc) candidates[*]=-1 no_match=[-1] FOR i=0,nsrc-1 DO BEGIN lon=lon_w[i]*!dpi/180 lat=lat_w[i]*!dpi/180 wh_lat=WHERE(abs((*(*data_cat).lat.value)-lat) lt dist) if wh_lat[0] ne -1 then begin gcirc, 0, lon, lat,(*(*data_cat).lon.value)[wh_lat],(*(*data_cat).lat.value)[wh_lat],dist_pos wh_lon=WHERE(dist_pos le dist) if wh_lon[0] ne -1 then begin wh_can=wh_lat[wh_lon] ncan=n_elements(wh_can) IF ncan GT 1 THEN BEGIN good_fluxes=where((finite((*(*data_cat).obs3.iflux)[wh_can]) eq 1) and ((*(*data_cat).obs3.iflux)[wh_can] ne 0)) IF good_fluxes[0] eq -1 THEN BEGIN good_fluxes=where((finite((*(*data_cat).obs2.iflux)[wh_can]) eq 1) and ((*(*data_cat).obs2.iflux)[wh_can] ne 0)) IF good_fluxes[0] eq -1 THEN BEGIN good_fluxes=where((finite((*(*data_cat).obs1.iflux)[wh_can]) eq 1) and ((*(*data_cat).obs1.iflux)[wh_can] ne 0)) IF good_fluxes[0] eq -1 THEN BEGIN print, 'WTF?' stop ENDIF ELSE BEGIN best_can=where((*(*data_cat).obs1.iflux)[wh_can[good_fluxes]] eq max(((*(*data_cat).obs1.iflux)[wh_can[good_fluxes]]))) best_can=wh_can[good_fluxes[best_can]] ENDELSE ENDIF ELSE BEGIN best_can=where((*(*data_cat).obs2.iflux)[wh_can[good_fluxes]] eq max(((*(*data_cat).obs2.iflux)[wh_can[good_fluxes]]))) best_can=wh_can[good_fluxes[best_can]] ENDELSE ENDIF ELSE BEGIN best_can=where((*(*data_cat).obs3.iflux)[wh_can[good_fluxes]] eq max(((*(*data_cat).obs3.iflux)[wh_can[good_fluxes]]))) best_can=wh_can[good_fluxes[best_can]] ENDELSE ENDIF ELSE best_can=wh_can[0] candidates[i]=best_can ;;power law for 5-23GHz if finite((*(*data_cat).obs3.iflux)[best_can]) then begin data_flux=(*(*data_cat).obs3.iflux)[best_can] / 1e3 ; mJy to Jy data_freq=(*(*data_cat).obs3.band).nu_c endif else begin if finite((*(*data_cat).obs2.iflux)[best_can]) then begin data_flux=(*(*data_cat).obs2.iflux)[best_can] / 1e3 ; mJy to Jy data_freq=(*(*data_cat).obs2.band).nu_c endif else begin if finite((*(*data_cat).obs1.iflux)[best_can]) then begin data_flux=(*(*data_cat).obs1.iflux)[best_can] / 1e3 ; mJy to Jy data_freq=(*(*data_cat).obs1.band).nu_c endif else begin print, 'WTF?' stop endelse endelse endelse if flux_k[i] lt 0 then begin if flux_ka[i] lt 0 then begin if flux_q[i] lt 0 then begin if flux_v[i] lt 0 then begin if flux_w[i] lt 0 then begin print, 'WTF?' stop endif else begin specind=alog10(flux_w[i]/data_flux)/alog10(94.*1e9/data_freq) flux_k[i]=10^(specind*alog10(23.*1e9/data_freq))*data_flux flux_ka[i]=10^(specind*alog10(33.*1e9/data_freq))*data_flux flux_q[i]=10^(specind*alog10(41.*1e9/data_freq))*data_flux flux_v[i]=10^(specind*alog10(61.*1e9/data_freq))*data_flux endelse endif else begin specind=alog10(flux_v[i]/data_flux)/alog10(61.*1e9/data_freq) flux_k[i]=10^(specind*alog10(23.*1e9/data_freq))*data_flux flux_ka[i]=10^(specind*alog10(33.*1e9/data_freq))*data_flux flux_q[i]=10^(specind*alog10(41.*1e9/data_freq))*data_flux endelse endif else begin specind=alog10(flux_q[i]/data_flux)/alog10(41.*1e9/data_freq) flux_k[i]=10^(specind*alog10(23.*1e9/(*(*data_cat).obs3.band).nu_c))*data_flux flux_ka[i]=10^(specind*alog10(33.*1e9/(*(*data_cat).obs3.band).nu_c))*data_flux endelse endif else begin specind=alog10(flux_ka[i]/data_flux)/alog10(33.*1e9/data_freq) flux_k[i]=10^(specind*alog10(23.*1e9/data_freq))*data_flux endelse endif else begin specind=alog10(flux_k[i]/data_flux)/alog10(23.*1e9/data_freq) endelse if data_freq ne (*(*data_cat).obs3.band).nu_c then begin data_flux_485=10^(specind*alog10((*(*data_cat).obs3.band).nu_c/data_freq))*data_flux endif else begin data_flux_485=data_flux endelse (*wmap_cat.e2.specind.value)[i]=specind (*wmap_cat.e2.ampl.value)[i]=data_flux_485 specind_back = specind ;;power law under 4.85GHz if finite((*(*data_cat).obs3.iflux)[best_can]) then begin if finite((*(*data_cat).obs2.iflux)[best_can]) then begin specind=alog10((*(*data_cat).obs3.iflux)[best_can]/(*(*data_cat).obs2.iflux)[best_can])/alog10((*(*data_cat).obs3.band).nu_c/(*(*data_cat).obs2.band).nu_c) endif else begin if finite((*(*data_cat).obs1.iflux)[best_can]) then begin specind=alog10((*(*data_cat).obs3.iflux)[best_can]/(*(*data_cat).obs1.iflux)[best_can])/alog10((*(*data_cat).obs3.band).nu_c/(*(*data_cat).obs1.band).nu_c) endif endelse endif else begin if finite((*(*data_cat).obs2.iflux)[best_can]) then begin if finite((*(*data_cat).obs1.iflux)[best_can]) then begin specind=alog10((*(*data_cat).obs2.iflux)[best_can]/(*(*data_cat).obs1.iflux)[best_can])/alog10((*(*data_cat).obs2.band).nu_c/(*(*data_cat).obs1.band).nu_c) endif endif endelse if finite(specind) eq 0 then specind = specind_back (*wmap_cat.e1.specind.value)[i]=specind (*wmap_cat.e1.ampl.value)[i]=data_flux_485 ;;power law for 23-33GHz if flux_ka[i] lt 0 then begin if flux_q[i] lt 0 then begin if flux_v[i] lt 0 then begin if flux_w[i] lt 0 then begin specind=(*wmap_cat.e2.specind.value)[i] flux_ka[i]=10^(specind*alog10(33./23.))*flux_k[i] flux_q[i]=10^(specind*alog10(41./23.))*flux_k[i] flux_v[i]=10^(specind*alog10(61./23.))*flux_k[i] endif else begin specind=alog10(flux_w[i]/flux_k[i])/alog10(94./23.) flux_ka[i]=10^(specind*alog10(33./23.))*flux_k[i] flux_q[i]=10^(specind*alog10(41./23.))*flux_k[i] flux_v[i]=10^(specind*alog10(61./23.))*flux_k[i] endelse endif else begin specind=alog10(flux_v[i]/flux_k[i])/alog10(61./23.) flux_ka[i]=10^(specind*alog10(33./23.))*flux_k[i] flux_q[i]=10^(specind*alog10(41./23.))*flux_k[i] endelse endif else begin specind=alog10(flux_q[i]/flux_k[i])/alog10(41./23.) flux_ka[i]=10^(specind*alog10(33./23.))*flux_k[i] endelse endif else begin specind=alog10(flux_ka[i]/flux_k[i])/alog10(33./23.) endelse (*wmap_cat.e3.specind.value)[i]=specind (*wmap_cat.e3.ampl.value)[i]=flux_k[i] ;;power law for 33-41GHz if flux_q[i] lt 0 then begin if flux_v[i] lt 0 then begin if flux_w[i] lt 0 then begin specind=(*wmap_cat.e3.specind.value)[i] flux_q[i]=10^(specind*alog10(41./33.))*flux_ka[i] flux_v[i]=10^(specind*alog10(61./33.))*flux_ka[i] endif else begin specind=alog10(flux_w[i]/flux_ka[i])/alog10(94./33) flux_q[i]=10^(specind*alog10(41./33.))*flux_ka[i] flux_v[i]=10^(specind*alog10(61./33.))*flux_ka[i] endelse endif else begin specind=alog10(flux_v[i]/flux_ka[i])/alog10(61./33) flux_q[i]=10^(specind*alog10(41./33.))*flux_ka[i] endelse endif else begin specind=alog10(flux_q[i]/flux_ka[i])/alog10(41./33) endelse (*wmap_cat.e4.specind.value)[i]=specind (*wmap_cat.e4.ampl.value)[i]=flux_ka[i] ;;power law for 41-61GHz if flux_v[i] lt 0 then begin if flux_w[i] lt 0 then begin specind=(*wmap_cat.e4.specind.value)[i] flux_v[i]=10^(specind*alog10(61./41.))*flux_q[i] endif else begin specind=alog10(flux_w[i]/flux_q[i])/alog10(94./41) flux_v[i]=10^(specind*alog10(61./41.))*flux_q[i] endelse endif else begin specind=alog10(flux_v[i]/flux_q[i])/alog10(61./41.) endelse (*wmap_cat.e5.specind.value)[i]=specind (*wmap_cat.e5.ampl.value)[i]=flux_q[i] ;;power law for 61-94GHz if flux_w[i] lt 0 then begin specind=(*wmap_cat.e5.specind.value)[i] flux_w[i]=10^(specind*alog10(94./61.))*flux_v[i] endif else begin specind=alog10(flux_w[i]/flux_v[i])/alog10(94./61.) endelse (*wmap_cat.e6.specind.value)[i]=specind (*wmap_cat.e6.ampl.value)[i]=flux_v[i] if specind lt 0 then (*wmap_cat.e7.specind.value)[i]=specind else (*wmap_cat.e7.specind.value)[i]=0 (*wmap_cat.e7.ampl.value)[i]=flux_w[i] endif else begin no_match=[no_match,i] endelse endif else begin no_match=[no_match,i] endelse ENDFOR ;;copy of the caracteristics of the sources (*wmap_cat.objectname) = (*(*model_cat).objectname)[candidates] IF wmap_cat.polarised THEN BEGIN (*wmap_cat.pol_fraction) = (*(*model_cat).pol_fraction)[candidates] (*wmap_cat.pol_angle) = (*(*model_cat).pol_angle)[candidates] ENDIF (*wmap_cat.lon.value) = (*(*model_cat).lon.value)[candidates] (*wmap_cat.lat.value) = (*(*model_cat).lat.value)[candidates] indix = where(candidates ge 0) candidates = candidates[indix] remove_src_cat, candidates, model_cat wcat = PTR_NEW(wmap_cat) remove_src_cat, no_match[1:n_elements(no_match)-1], wcat filename = FIX_SEPARATOR(GETPSMDIR('ps')) +'strong'+ names.wmapps_basename + '.sav' WRITE_CATALOG, *wcat, filename, /sav, name='comp' DESTROY_DATA, wmap_cat DESTROY_DATA, wcat IF runp.verbosity GE 1 THEN PSM_INFO_MESSAGE, 'WMAP radiosource model catalogue written' END