#
# pointsource_utilities.PY
# Utilties functions for data analysis of point sources.
#
# date: 2017-12-21
# author: GIUSEPPE PUGLISI
#
# Copyright (C) 2017 Giuseppe Puglisi giuspugl@sissa.it
#
import numpy as np
import scipy.optimize as so
import astropy
from astropy import units as u, constants as C
import scipy as sp
from scipy import integrate, interpolate
from IO import *
from Stats import *
from healpy import read_map, query_strip, query_disc, get_nside, mollview
[docs]def compute_variance_ratios(fsky, dir_ps ='./data/'):
"""
Given a fraction of sky :math:`f_{sky}` if computes the ratio of the variance of the synchrotron and dust maps
within a patch with size ``fsky`` and the one within ``fsky0`` which is the value from Planck collaboration et al. 2015.X .
All the considered patches exclude the Galactic plane.
**Parameters**
- `fsky`: {float}
- `dir_ps`:{string}
path to the PS4C data folder
**Returns**
-`fs`: {float}
ratio of synchrotron variances
-`fd`: {float}
ratio of thermal dust variances
"""
sizeradpatch=np.sqrt(fsky).value
fskyval= fsky.value /(4.*np.pi)
synchro=hp.read_map(dir_ps+'synchro_30GHz_256.fits', field=[1,2],verbose=False)
Ps=np.sqrt(synchro[1] **2 +synchro[0]**2)
dust=hp.read_map(dir_ps+'thermaldust_353GHz_256.fits', field=[1,2],verbose=False)
Pd=np.sqrt(dust[1] **2 +dust[0]**2)
nside=get_nside(synchro)
listpix0=hp.query_strip(nside, (105.*u.deg).to(u.rad).value, (75.*u.deg).to(u.rad).value )
galstrip=hp.query_strip(nside, (75.*u.deg).to(u.rad).value, (105.*u.deg).to(u.rad).value )
center=(180.-30.)*u.deg
fsky2patch={fskyval<0.02: 'small', 0.02<=fskyval<0.1:'large', 0.1<=fskyval<=0.5:'verylarge',0.5<fskyval<=1:'full'}
listp_dict= {'small': query_disc(nside, hp.ang2vec(center.to(u.rad).value,(35*u.deg).to(u.rad)),sizeradpatch),
'large': query_disc(nside, hp.ang2vec(center.to(u.rad).value,(35*u.deg).to(u.rad)),sizeradpatch),
'verylarge': query_disc(nside, hp.ang2vec(center.to(u.rad).value,(35*u.deg).to(u.rad)),sizeradpatch),
'full': listpix0 }
patch=fsky2patch[True]
listpix=listp_dict[patch ]
mask=np.zeros(hp.nside2npix(nside))
mask[listpix]=1.
mask[galstrip]=0.
select_pix=np.ma.masked_equal(mask,1).mask
fs= (np.var(Ps[select_pix])/np.var (Ps[listpix0]))
fd= ( np.var(Pd[select_pix])/np.var(Pd[listpix0]))
return fs,fd
[docs]def forecast_Galaxy(ell,nu ,fsky, **kwargs ):
"""
Compute the Galactic contribution in multipole ``ell`` , frequency ``nu`` and
fraction of sky ``fsky`` given eq.(22) of Planck 2015 X data are taken from table 4, table 5, table 11.
**Parameters**
- `ell` : {float or array}
multipole orders
- `nu` :{float}
- `fsky`:{float}
"""
Asynch=20 *u.K #data in RJ
nuscale=0.408*u.GHz
nuscale2=30*u.GHz
beta_s=3.12
s_sync=lambda nu : Asynch *(nu/nuscale ) **(-( beta_s +2) )
qs=1.3 *(u.uK)**2
alpha_s=-0.31
qd=208 *(u.uK)**2
alpha_d=-0.59
gamma=(C.h.cgs/C.k_B.cgs /(20 *u.K)).to(1./u.GHz)
nudust=353 *u.GHz
beta_d=1.55
Adust= 15
fsky0=0.73 *4.*np.pi *u.sr
s_dust=lambda nu: Adust *(nu /nudust)**(beta_d +1)* (np.exp(gamma *nudust ) -1 ) / (np.exp(gamma *nu ) -1 )
a=1
b=200
c=140
r=.10
correctionfunction=lambda x: (a*np.exp(c*r) +b*np.exp(r*x)) /(np.exp(c*r) + np.exp(r*x))
factor=correctionfunction(nu.value)/correctionfunction(b)
fs, fd= compute_variance_ratios(fsky,**kwargs)
return ( qs*(ell/80.)**alpha_s * s_sync(nu)/s_sync(nuscale2) *fs
+ fd*qd *(ell/80.)**alpha_d * s_dust(nu)/s_dust(nudust))*(fsky/fsky0)
[docs]def fit_numbercounts_w_powerlaw(S, n , Smax, **kwargs):
"""
Fit the differential number counts as a single power law :math:`dN/dS \propto K S^{-\gamma}` up to a maximum flux density.
**Parameters**
- `S`: {array}
array of fluxes
- `n`: {function}
differential number counts computed from interpolating number counts from model
- `Smax`: {float}
the flux limit you may want to estimate the power law fit
**Returns**
- `K`:{float}
- `Gamma`: {float}
"""
Smin=min(S)*2
if not 'verbose' in kwargs: kwargs['verbose']=False
f=lambda x,A, gamma:A*x**(-gamma)
idval=np.argmin(abs(Smax.value - S.value ) )
popt,cov= so.curve_fit(f , S[:idval], n(S[:idval]),p0=[20,2])
K, Gamma=popt[0], popt[1]
if kwargs['verbose']: print "Fitting up to S< %g Jy, Gamma= %g, K=%g"%(S[idval].value, Gamma, K)
return K, Gamma
[docs]def compute_confusion_limit(K, Gamma, omega_b, q=5):
"""
Compute the confusion limit :math:`S_c=5\sigma_c` by means of the definition in Condon, 1974.
.. math::
\sigma_c = ( {q ^{3- \gamma}} k \Omega_e/({3- \gamma}) )^{1/ (\gamma-1)}
where we set :math:`q=5` for the :math:`5\sigma` *truncation limit*.
**Parameters**
- `K`, `Gamma` as computed from :func:`fit_numbercounts_w_powerlaw`.
"""
K*= np.power((u.Jy ), (Gamma-1)) /u.sr
Omega_e= omega_b /(Gamma-1)
sigmac= (q**(3-Gamma)*K*Omega_e/(3-Gamma))**(1/(Gamma-1))
return sigmac
[docs]def beamsolidangle(theta):
"""
Given a FWHM resolution `theta` in radians, it computes the solid angle in steradians
"""
return (np.pi* theta**2/4./np.log(2)).to(u.steradian)
[docs]def Krj2Kcmb(nu, Trj):
"""
Convert antenna temperature ( Rayleigh-Jeans) into the physical one
"""
nu0 =56.8 *u.GHz
x=nu/nu0
return Trj / (x**2*np.exp(x)/(np.exp(x) -1))
[docs]def Kcmb2Krj(nu, Tcmb):
"""
Convert physical temperature into the antenna one.
"""
nu0 =56.8 *u.GHz
x=nu/nu0
return (x**2*np.exp(x)/(np.exp(x) -1) )*Tcmb
[docs]def Krj2brightness(temp, freq, theta_beam):
"""
Convert brightness temperature to brightness measured in [Jy /beam]
"""
# uK_RJ----------> mJy /beam
# Astropy here is misleading, it requires Jy quantity but physically should be Jy/sr
#sigmaT [uK arcmin]
Omega_beam = beamsolidangle(theta_beam)
equiv = u.brightness_temperature(Omega_beam, freq)
return temp.to(u.milliJansky , equivalencies=equiv)
[docs]def brightness2Krj(sigmaS, freq, theta_beam):
"""
Convert brightness to brightness temperature
"""
# mJy /beam --------> K_RJ
# Astropy here is misleading, it requires Jy quantity but physically should be Jy/sr
#sigmaS [mJy beam-1]
Omega_beam = beamsolidangle(theta_beam)
equiv = u.brightness_temperature(Omega_beam, freq)
return sigmaS.to(u.Kelvin, equivalencies=equiv)
[docs]def brightness2Kcmb (nu,Ib ):
"""
Returns conversion factor to pass from ``Ib`` brightness units to physical temperature at a given frequency ``nu``
"""
nu0=56.8 *u.gigahertz
x= nu/nu0
return Ib* 4.e-2*(np.exp(x) -1)**2/( x**4 *np.exp(x)) *u.uK /(u.Jy/ u.sr )
[docs]def estimate_power_spectrum_contribution(dNdS, nu,omega_b, Smin, Smax ):
"""
Given differential number counts :math:`dN/dS`, estimated at certain flux densities :math:`S\in [S_{min}, S_{max}]`,
it estimates the integral
.. math::
C= \int _{S_{min}} ^{ S_{max}} dS {n(S)} S^2
**Parameters**
- `dNdS`: {func}
differential number counts interpolated by means of :func:`numpy.interpolate.interp1d`
- ` nu`:{float}
frequency
- `omega_b`:{float}
solid angle of the beam
- `Smin` and `Smax`:{float}
range of integration
"""
# function to integrate
integrand= lambda s: dNdS(s) * s**2
#integral
Integral = sp.integrate.quad( integrand, Smin,Smax,limit=1000 ,epsrel=1.e-3) [0] *u.Jy**2/u.sr
flux2Kcmb= brightness2Kcmb(nu, 1.)
#flux2Kcmb = brightness2Krj(1*u.Jy, nu,fwhm).to(u.uK)
return flux2Kcmb**2 *Integral
[docs]def linfunc(x,m,q):
return m*x +q
[docs]def sigmaf(P,I, sigmaP, sigmaI):
"""
Given polarized and intensity flux , P and I and their uncertainties sigmaP and sigmaI, propagates
the uncertainties to polarization fraction :math:`\Pi= P/I`.
"""
return np.sqrt(sigmaP**2 / I**2 + (P/I**2 * sigmaI)**2 )
[docs]def get_spectral_index(v1,v2, S1,S2):
"""
Computes the spectral index :math:`alpha_{v1} ^{v2}` of fluxes S1, S2 estimated at two different frequencies, v1 and v2.
"""
return np.log(S1/S2)/ np.log(v1/v2)
[docs]def compute_polfrac_scaling( dir_ps,exclude_HFI_data=True, include_steep=False, **kwargs ):
"""
Considering observations from several catalogues, it computes the observed polarization fractions
and fit a lognormal distribution in order to get from the best fit parameters the average
values of :math:`< \Pi> ,< \Pi^2>, \Pi_{med}` and their errors. Once all the datasets have been fitted
it computes the scaling of :math:`< \Pi^2>^{1/2}` by means of a linear function.
**Parameters**
- `dir_ps`: {string}
path to the PS4C data folder
- `exclude_HFI_data`:{bool}
if set to `True` it exclude the data from HFI channels of Planck.
- `include_steep`:{bool}
if set to `True` it includes in the computation even Steep Spectrum Radio Sources.
We recommend you to use this option **if** you are running at low radio frequencies, i.e. :math:`<20` GHz.
**Returns**
- `A, sigmaA`:{float}
the slope and the error on the slope as computed from the linear fit
- `B, sigmaB`:{float}
the constant term fitted from the linear function
"""
fname=dir_ps+'spass_nvss_catalog.tsv'
cat,dict_cat=read_spass_nvss_catalog(fname, **kwargs)
nfreq=18
if include_steep :
flat=np.ma.masked_less (cat['alpha'], 0)
nfreq-=1
else:
flat= np.ma.masked_greater(cat['alpha'], -0.5 )
upmaj= np.ma.masked_equal(cat['f_majaxis'],-1)
upmin= np.ma.masked_equal(cat['f_minaxis'],-1)
unres = np.logical_and(upmaj.mask, upmin.mask)
flat_unres=np.logical_and(flat.mask, unres )
notneg=np.ma.masked_greater(cat['PSP'],0)
mask=np.logical_and(flat_unres, notneg.mask)
if exclude_HFI_data:
nfreq-=3
bins={30:15,44:10,70:15}
bonaveravals=[4.03, 4.64,4.48]
bonaveraerrs=[0.65, 1.58,0.49]
else:
bins={30:15,44:10,70:15, 100:7,143:6,217:15}
bonaveravals=[4.03, 4.64,4.48,5.75,4.91,4.54]
bonaveraerrs=[0.65, 1.58,0.49,0.47,0.51,0.41]
pimedian=np.zeros(nfreq)
pisquare=pimedian*0.
deltapisq=pimedian*0.
deltapim =pimedian*0.
freqs=pimedian*0.
deltapi=np.zeros_like(pisquare)
pimean=np.zeros_like(pisquare)
idx=0
freqs[idx]=1.4
pimedian[idx],pisquare[idx],deltapisq[idx],deltapim[idx], pimean[idx], deltapi[idx]=fitting_lognormal_from_fluxes(cat['INV'][mask], cat['PNV'][mask],
cat['e_INV'][mask],cat['e_PNV'][mask],
'NVSS',workdir=dir_ps,fig=0, **kwargs)
idx+=1
freqs[idx]=2.3
pimedian[idx ],pisquare[idx ],deltapisq[idx ],deltapim[idx ], pimean[idx], deltapi[idx]=fitting_lognormal_from_fluxes(cat['ISP'][mask], cat['PSP'][mask],
cat['e_ISP'][mask],cat['e_PSP'][mask],
'SPASS', fig=1, workdir=dir_ps, **kwargs)
for f in [18,100]:
flag,fpol,errfpol= read_galluzzi_fpol_dat(dir_ps+'m'+str(f)+'_resampling.dat', **kwargs)
if f==100:
resampbool=True
nbins=10
string='_ALMA'
else:
nbins=15
string=''
resampbool=True
idx+=1
freqs[idx]=f
pimedian[idx],pisquare[idx ], deltapisq[idx ],deltapim[idx ], pimean[idx], deltapi[idx]=fitting_lognormal_from_fpol(fpol,errfpol,
flag,str(f)+'GHz'+string,workdir=dir_ps,
fig=f, nbins=nbins, resampling=resampbool, **kwargs)
data= read_jvas_catalog(dir_ps+'jvas_fsrq_ned.dat', **kwargs)
flagjvas = np.ma.masked_less_equal( value=0,x= data['efpol']).mask
idx+=1
freqs[idx]=8.4
pimedian[idx],pisquare[idx], deltapisq[idx], deltapim[idx], pimean[idx], deltapi[idx]=fitting_lognormal_from_fpol(data['fpol']/100., data['efpol']/100., flagjvas,
'JVAS', fig=8, nbins=100,resampling=True,workdir=dir_ps, **kwargs )
planck_catalogs= np.array(glob.glob(dir_ps+'COM_PCCS_*'))
channels =[catalog.split('_')[-2] for catalog in planck_catalogs]
idx_sort=np.argsort((channels),axis=None,kind='quicksort')
for j,catalog in np.ndenumerate(planck_catalogs[idx_sort]):
#READ CATALOGUES
ch =catalog.split('_')[-2]
freq=int(ch)
if exclude_HFI_data and freq>=100:
continue
Iflux,Polflux, eI, eP=extract_flux_from_planck_catalog(catalog,freq, **kwargs)
idx+=1
freqs[idx ]=int(ch)
pimedian[idx],pisquare[idx],deltapisq[idx],deltapim[idx], pimean[idx], deltapi[idx]=fitting_lognormal_from_fluxes(np.array(Iflux),np.array(Polflux),
eI, eP,ch+'GHz',workdir=dir_ps,
fig=freq, nbins=bins[freq], **kwargs)
cat,dict_cat= read_atca_catalog(dir_ps+'at20_catalog.tsv', **kwargs)
alpha=get_spectral_index(4.8,8.6, cat['S5'], cat['S8'])
# define FSRQs at low freqs and consider them at 20 GHz
if include_steep :
maskflat=np.ma.masked_less((alpha) ,0).mask
else:
maskflat=np.ma.masked_inside((alpha) ,-0.5 ,0).mask
# mask nan AND SSRQ in fpol
mask1 = np.logical_and( np.logical_not( np.ma.masked_invalid(cat['m5\n']).mask) , maskflat )
mask2 = np.logical_and(np.logical_not( np.ma.masked_invalid(cat['m8']).mask) , maskflat)
# flagging the upper limits
flag1=( np.ma.masked_invalid(cat['e_P5'][mask1]).mask )
flag2=( np.ma.masked_invalid(cat['e_P8'][mask2]).mask )
fpol1=np.array(cat['m5\n'][mask1] )/100.
fpol2=np.array(cat['m8'][mask2] ) /100.
erfpol1=sigmaf( np.array(cat['P5'][mask1] ),np.array(cat['S5'][mask1] ),np.array(cat['e_P5'][mask1] ),np.array(cat['e_S5'][mask1] ))
erfpol2=sigmaf( np.array(cat['P8'][mask2] ),np.array(cat['S8'][mask2] ),np.array(cat['e_P8'][mask2] ),np.array(cat['e_S8'][mask2] ))
idx +=1
freqs[idx]=4.8
pimedian[idx],pisquare[idx],deltapisq[idx],deltapim[idx], pimean[idx], deltapi[idx]=fitting_lognormal_from_fpol(fpol1,
erfpol1, flag1, '5GHz', fig=5,workdir=dir_ps,
nbins=20, resampling=True, **kwargs)
idx +=1
freqs[idx]=9.
pimedian[idx],pisquare[idx],deltapisq[idx],deltapim[idx], pimean[idx], deltapi[idx]=fitting_lognormal_from_fpol(fpol2,
erfpol2, flag2, '8GHz', fig=9,workdir=dir_ps,
nbins=20, resampling=True , **kwargs)
idx+=1
cat, dict_catalog= read_IRAM_catalog(dir_ps+'iram_catalog.tsv')
msk =np.ma.masked_not_equal( cat['mL'] , -1 ).mask
fpol= cat['mL'][msk]
errfpol= cat['e_mL'][msk]
flag= np.ma.masked_equal(errfpol , -1 ).mask
freqs[idx]=89
pimedian[idx],pisquare[idx], deltapisq[idx], deltapim[idx], pimean[idx], deltapi[idx]=fitting_lognormal_from_fpol(fpol/100., errfpol/100.,
flag,'89GHz', nbins=8,workdir=dir_ps,
**kwargs)
cat, dict_catalog= read_VLA_catalog(dir_ps+'vla_catalog.tsv')
alpha1=get_spectral_index(S1=cat['SC'], S2=cat['SX'],v1=4.86 , v2=8.46 )
if include_steep :
flats=np.ma.masked_less(alpha1, 0 ).mask
vlabins={'C':4,'X':8, 'K':8, 'Q':6}
bands={'C':5.6,'X':8.46, 'K':22.46}
else :
vlabins={'C':4,'X':8, 'K':12, 'Q':12}
flats=np.ma.masked_inside(alpha1, -0.5,0 ).mask
bands={'C':5.6,'X':8.46, 'K':22.46, 'Q':43.34}
for b,v in bands.items():
msk =np.ma.masked_not_equal( cat['p'+b] , -1 ).mask
mask1 = np.logical_and(flats, msk )
fpol= cat['p'+b][mask1]
errfpol= cat['e_p'+b][mask1]
flag= np.ma.masked_equal(errfpol , -1 ).mask
idx+=1
freqs[idx]=v
pimedian[idx],pisquare[idx], deltapisq[idx], deltapim[idx], pimean[idx], deltapi[idx]= fitting_lognormal_from_fpol(fpol/100, errfpol/100.,
flag, 'VLA_'+b+'_band',workdir=dir_ps,
nbins=vlabins[b] , **kwargs)
sqrtp, deltasqrtp=np.sqrt(pisquare), (deltapisq)/2./np.sqrt(pisquare)
fs= np.sort(bins.keys())
fitfreq=( np.concatenate((freqs, fs+3)) )
data= np.concatenate((sqrtp, np.array(bonaveravals)))
errors= np.concatenate((deltasqrtp,bonaveraerrs))
pop,cov=so.curve_fit(linfunc, fitfreq ,data ,sigma=errors)
A , B= pop[0]/u.gigahertz, pop[1]
sigmaA, sigmaB=np.sqrt(cov[0,0])/u.gigahertz,np.sqrt(cov[1,1])
return A, B, sigmaA, sigmaB
[docs]def compute_dNdP(A,mu, sigma, dNdS, S, Np_quad=32):
"""
Compute the convolution of differential number counts with the lognormal distribution function, i.e.
.. math::
n(P) = \int_{P=S_0} ^{\infty} {dS}/ {S} \mathcal{P} (\Pi) n(S)
**Parameters**
- `A, mu, sigma`: {floats}
lognormal best fit values computed from :func:`Stats.fitting_lognormal_from_fpol` and :func:`Stats.fitting_lognormal_from_fluxes`.
- `dNdS`:{func}
interpolated differential number counts
- `S`:{array}
flux densities
- `Np_quad`:{int}
number of quadrature points.
**Returns**
- `Pbin`: {array}
polarized fluxes
- `interpolated_dNdP`:{func}
interpolated polarized differential number counts.
"""
logf= lambda x: lognormal_distribution((x)*100., A, mu,sigma)
Normalization= integrate.quad(logf,a=1e-7, b=1 )[0]
dNdP=np.zeros(Np_quad)
Pbin =np.zeros(Np_quad)
for i,P in np.ndenumerate(np.logspace(-3,1,Np_quad)):
logfunc= lambda x: lognormal_distribution((P/x)*100., A, mu,sigma)/Normalization
integrand=lambda x: logfunc(x) *dNdS(x)/x
try:
dNdP[i],err= integrate.quad(integrand,a=P, b=S.max().value )
except ValueError:
dNdP[i],err= integrate.quad(integrand,a=S.min().value, b=S.max().value )
Pbin[i]= P
interpolated_dNdP= interpolate.interp1d(Pbin , dNdP )
return Pbin, interpolated_dNdP