See wiki page: https://wiki.cosmos.esa.int/planckpla2015/index.php/Main_Page
Links for masks: http://irsa.ipac.caltech.edu/data/Planck/release_2/ancillary-data/previews/HFI_Mask_PointSrc_2048_R2.00/header.txt http://irsa.ipac.caltech.edu/data/Planck/release_2/ancillary-data/previews/HFI_Mask_GalPlane-apo5_2048_R2.00/header.txt
SZ Source Catalog: https://wiki.cosmos.esa.int/planckpla2015/index.php/Catalogues#.282015.29_Second_SZ_Catalogue
Format: Healpix format, with Nside of 1024 (LFI 30, 44 and 70) and 2048 (LFI 70 and HFI), in Galactic coordinates, and Nested ordering.
In [1]:
path = "/Users/inchani/Desktop/UC\ Davis/My\ Courses/STA\ 250\ (AstroStatistics)/Project/";
include("$path""lib.jl")
map_name = "HFI_SkyMap_100_2048_R2.02_full.fits"; # Polarization map at 100 GHz
using PyCall, PyPlot, StatsBase
@pyimport healpy as hp
@pyimport healpy.fitsfunc as fitsfunc;
@pyimport healpy.pixelfunc as pixelfunc;
@pyimport numpy as np;
dθ = dϕ = 0.0007669903939429012; # 5 arcmin of resolution to radian
res_highest = 3. / 60. * pi / 180; # 3 arcmin (actual highest resolution = 2.637 arcmin at equator)
NSIDE = 2048;
Nested = false;
I_STOKES = hp.read_map("$path$map_name", field = 0, memmap = true);
dim = length(I_STOKES);
const dtheta_deg = 15.
const dtheta = dtheta_deg / 180 * pi; # 15 degrees in radian
const FWHM = 30. / 60. * pi / 180; # 30. arcmin to radian unit (See Planck 2013 ISW paper)
const σ = FWHM / 2.355 # 2.355 ~ 2√(2log(2))
const σnorm2 = 2.*σ^2.
const σlim2 = (3σ)^2.
const angsize = copy(dtheta_deg*2) # width and height in degree
const XYsize = angsize * pi / 180; # in radian
#const res = 6. / 60. * pi / 180;
const res = 3. / 60. * pi / 180;
# 6 arcmin of pixel size in radian unit (See Planck 2013 ISW paper)
const Nsize = Int64(round(XYsize/res)); # 600
const Tmin = -593.5015506111085; # mu K
const Tmax = 709.0113358572125; # mu K
# min and max temperature if we take out 40 % of sky as foreground
Out[1]:
In [3]:
function PlotTempProfile(Rmax::Float64,dR::Float64,Temperature::Array{Float64, 2};vmin=1.,vmax=0.)
x1dang = linspace(-dtheta_deg,dtheta_deg,Nsize); y1dang = linspace(-dtheta_deg,dtheta_deg,Nsize);# dR = 0.1;
R0 = dR;# Rmax = 15.
dim = Int((Rmax - R0) / dR) + 1
Rarr = linspace(R0,Rmax,dim) - dR/2
numarr = zeros(dim)
Tarr = zeros(dim)
for i=1:Nsize
for j=1:Nsize
n = Int(floor( √(x1dang[j]^2 +y1dang[i]^2) / dR) ) + 1
if n <= dim
Tarr[n] += Temperature[i,j]
numarr[n] += 1.
end
end
end
for i=1:dim
Tarr[i] = Tarr[i] / numarr[i]
end
PyPlot.plot(Rarr,Tarr)
if vmin > vmax
vmin = np.min(Tarr)
vmax = np.max(Tarr)
end
PyPlot.xlabel("Distance from center (\$\\degree\$\)")
PyPlot.ylabel("Average Temperature (\$\\mu\$\K\$\_{CMB}\$\)")
PyPlot.xlim(0,Rmax)
PyPlot.ylim(vmin,vmax)
end
function PlotTmap(Temperature::Array{Float64, 2},vmin::Float64,vmax::Float64,Rmax::Float64)
PyPlot.imshow(Temperature, origin='l',vmin = vmin, vmax = vmax, extent= [-15,15,-15,15])
PyPlot.xlabel("degree (\$\\degree\$\)")
PyPlot.ylabel("degree (\$\\degree\$\)")
PyPlot.colorbar(label= " \$\\mu\$\K\$\_{CMB}\$\ ");
end
Out[3]:
In [ ]:
Q_STOKES = hp.read_map("$path$map_name", field = 1, memmap = true);
U_STOKES = hp.read_map("$path$map_name", field = 2, memmap = true);
Each fits file contains:
Field 0 = I_STOKES, the intensity in each specific band
Field 1 = Q_STOKES, the polarized brightness
Field 2 = U_STOKES,
Field 3 = HITS , the number of observations
Field 4 = II_COV , the variance in the corresponding Stokes parameter
Field 5 = IQ_COV
Field 6 = IU_COV , the covariance inbetween the corresponding Stokes parameter
Field 7 = QQ_COV
Field 8 = QU_COV
Field 9 = UU_COV
In [2]:
GalMapFile = "HFI_Mask_GalPlane-apo5_2048_R2.00.fits"
PtMapFile = "HFI_Mask_PointSrc_2048_R2.00.fits"
GalMap = hp.read_map("$path$GalMapFile", field = 2, memmap = true); # 40% of sky is ruled out
PtMap = hp.read_map("$path$PtMapFile", field = 0, memmap = true);
(Ordering converted to RING)
In [3]:
hp.mollview(I_STOKES*1e6, xsize = 800, title = "Intensity Map (\$\\mu\$\K\$\_{CMB}\$\)", min = -300, max = 300)
savefig("$path/allskymap.png")
masking galactic foregrounds and point sources
In [3]:
planck_map_mask = copy(I_STOKES)*1e6
cnt = 0
for i = 1:max(length(GalMap),length(PtMap))
if (GalMap[i] == 0.) | (PtMap[i] == 0.)
planck_map_mask[i] = hp.UNSEEN
cnt += 1
end
end
In [4]:
hp.mollview(planck_map_mask, xsize = 800, title = "Intensity Map (\$\\mu\$\K\$\_{CMB}\$\)", min = -300, max = 300)
savefig("$path/allskymap_mask.png")
How to use pixelfunc for change index to coord, or vice versa. @pyimport healpy.pixelfunc as pixelfunc ind = 0 l,b = pixelfunc.pix2ang(2048,ind,nest=false) # pixel index > Spherical Coord, [θ, ϕ] [radian unit] ind_rtn = pixelfunc.ang2pix(2048,l,b,nest=false) # Spherical Coord, [θ, ϕ]> pixel index Note: ϕ = 0 (θ=90 deg) is toward galactic center > need to convert coordinates of SZ sources increasing ϕ = moving East, increasing θ = moving South
In [5]:
CoordSZ = np.load("$path""Coord_SZ.npy");
CoordGR08SC = np.load("$path""Coord_GR08_sc.npy");
CoordGR08Void = np.load("$path""Coord_GR08_void.npy");
In [6]:
Coord_obj = copy(CoordSZ);
Nsample = 10;
planck_map = copy(I_STOKES)*1e6;
i_pixel = Array(Int64, 0) # a list of pixel indices of a supercluster
ids = Array(Int64, 0) # survived index of Coord_obj after sorting out bad samples
N_obj = Array(Int64, 0) # number of pixels in a region
println("---| Start clipping regions of superclusters: 30 deg x 30 deg (15 deg = $dtheta radian)")
i = 0; ntry = 1;
if Nsample > length(Coord_obj[:,1])
Nsample = length(Coord_obj[:,1])
end
while (i < Nsample) & (ntry < length(Coord_obj))
l, b = Coord_obj[ntry,:]
θ = l * 180 / pi
ϕ = b * 180 / pi
#if (( θ > 15. ) & (θ < 75.)) | (( θ > 115. ) & (θ < 165.))
N, ind = SpheCoord2Index(Coord_obj[ntry,1]-dtheta,Coord_obj[ntry,1]+dtheta,
Coord_obj[ntry,2]-dtheta,Coord_obj[ntry,2]+dtheta;discrepancy=0.7)
if N > 0
i +=1
println(" |>> No. $i out of $ntry with angular position, (θ,ϕ) = ($θ, $ϕ)")
planck_map[ind] = hp.UNSEEN
ids = vcat(ids, i)
N_obj = vcat(N_obj, N)
i_pixel = vcat(i_pixel, ind)
println(" |>> total N = $N")
end
#end
ntry +=1
end
In [7]:
hp.mollview(planck_map, xsize = 800, title = "Intensity Map (\$\\mu\$\K\$\_{CMB}\$\)", min = -300, max = 300)
In [8]:
Coord_obj = copy(CoordGR08SC);
Nsample = 20;
planck_map = copy(I_STOKES)*1e6;
i_pixel = Array(Int64, 0) # a list of pixel indices of a supercluster
ids = Array(Int64, 0) # survived index of Coord_obj after sorting out bad samples
N_obj = Array(Int64, 0) # number of pixels in a region
println("---| Start clipping regions of superclusters: 30 deg x 30 deg (15 deg = $dtheta radian)")
i = 0; ntry = 1;
if Nsample > length(Coord_obj[:,1])
Nsample = length(Coord_obj[:,1])
end
while (i < Nsample) & (ntry < length(Coord_obj))
l, b = Coord_obj[ntry,:]
θ = l * 180 / pi
ϕ = b * 180 / pi
#if (( θ > 15. ) & (θ < 75.)) | (( θ > 115. ) & (θ < 165.))
N, ind = SpheCoord2Index(Coord_obj[ntry,1]-dtheta,Coord_obj[ntry,1]+dtheta,
Coord_obj[ntry,2]-dtheta,Coord_obj[ntry,2]+dtheta;discrepancy=0.7)
if N > 0
i +=1
println(" |>> No. $i out of $ntry with angular position, (θ,ϕ) = ($θ, $ϕ)")
planck_map[ind] = hp.UNSEEN
ids = vcat(ids, i)
N_obj = vcat(N_obj, N)
i_pixel = vcat(i_pixel, ind)
println(" |>> total N = $N")
end
#end
ntry +=1
end
In [9]:
hp.mollview(planck_map, xsize = 800, title = "Intensity Map (\$\\mu\$\K\$\_{CMB}\$\)", min = -300, max = 300)
In [13]:
Nrs = 50;
Ntry = 1; i = 0;
CoordRand = Array(Float64, Nrs, 2)
N_rs = Array(Int64, 0)
i_rs_pixel = Array(Int64, 0)
println("constructing random coordinates")
while i < Nrs
θ = (rand() * 178 + 1.) * pi / 180.
ϕ = rand() * 2pi
#if (θ < 60.) | (θ > 120.) # choose θ where 1 < θ < 75 or 105 < θ < 179 degree
if planck_map_mask[pixelfunc.ang2pix(2048,θ,ϕ)] > hp.UNSEEN
N, ind = SpheCoord2Index(θ-dtheta,θ+dtheta,ϕ-dtheta,ϕ+dtheta)
if N > 0
i +=1
println("No. $i out of $Ntry tries with angular position, (θ,ϕ) = (", θ*180/pi, ", ", ϕ*180/pi,")")
N_rs = vcat(N_rs, N)
i_rs_pixel = vcat(i_rs_pixel, ind)
CoordRand[i,:] = [θ, ϕ]
end
end
Ntry += 1
end
np.save("$path""CoordRand""$Nrs",CoordRand)
In [11]:
CoordRand = np.load("$path""CoordRand10.npy");
planck_map = copy(I_STOKES)*1e6
for i = 1:length(CoordRand[:,1])
l, b = CoordRand[i,:]
N, ind = SpheCoord2Index(CoordRand[i,1]-dtheta,CoordRand[i,1]+dtheta,
CoordRand[i,2]-dtheta,CoordRand[i,2]+dtheta)
planck_map[ind] = hp.UNSEEN
end
In [14]:
hp.mollview(planck_map, xsize = 800, title = "Intensity Map (\$\\mu\$\K\$\_{CMB}\$\)", min = -300, max = 300)
In [16]:
function GKernel(x::Float64,y::Float64,x₀::Float64,y₀::Float64)
r2 = (x-x₀)^2. + (y-y₀)^2.
if r2 < σlim2
return exp( -r2 / σnorm2 )
else
return 0.
end
end
Out[16]:
In [29]:
x1d = linspace(-XYsize*0.5,XYsize*0.5,Nsize)
y1d = linspace(XYsize*0.5,-XYsize*0.5,Nsize)
Tmap = Float64[ GKernel(xi,yi,0.,0.) for xi in x1d, yi in y1d];
Umap = Float64[ GKernel(xi,yi,0.,0.)>0.? 1.:0. for xi in x1d, yi in y1d];
TmapNorm = sum(Tmap);
Tmap /= TmapNorm;
In [326]:
PyPlot.figure(figsize=(5,5))
PyPlot.imshow(Tmap, origin='l', extent= [-dtheta_deg,dtheta_deg,-dtheta_deg,dtheta_deg]);
PyPlot.title("A Gaussian kernel with FWHM = 30'",fontsize=10)
PyPlot.xlabel("degree (\$\\degree\$\)")
PyPlot.ylabel("degree (\$\\degree\$\)")
PyPlot.colorbar(label= " \$\\mu\$\K\$\_{CMB}\$\ ");
In [346]:
#x = rand(30) * XYsize - 0.5XYsize
#y = rand(30) * XYsize - 0.5XYsize
x = [linspace(-XYsize/2,XYsize/2, 10); linspace(-XYsize/2,XYsize/2, 10);zeros(2)]
y = [linspace(-XYsize/2,XYsize/2, 10); zeros(10);zeros(2)]
θc, ϕc = 0., 0.;
TestMap = zeros(Nsize,Nsize)
for i =1:length(x)
θ, ϕ = x[i], y[i]
row_shift = Int64(round( (θ - θc) /res ))
col_shift = Int64(round( (ϕ - ϕc) /res ))
TestMap += ShiftArray(Tmap, row_shift, col_shift)
end
PyPlot.figure(figsize=(5,5))
PyPlot.imshow(TestMap, origin='l', extent= [-dtheta_deg,dtheta_deg,-dtheta_deg,dtheta_deg])
PyPlot.title("Testing superposition of kernels",fontsize=12)
PyPlot.xlabel("degree (\$\\degree\$\)")
PyPlot.ylabel("degree (\$\\degree\$\)")
PyPlot.colorbar(label= " \$\\mu\$\K\$\_{CMB}\$\ ");
In [17]:
planck_map = copy(I_STOKES) * 1e6 # Converting it to mu K scale;
Tmin = -593.5015506111085 # mu K
Tmax = 137299.59726333618 # mu K
for i = 1:length(PtMap)
if PtMap[i] == 0.
planck_map[i] = hp.UNSEEN
end
end
In [32]:
img = zeros(Nsize,Nsize)
cnt = 0
Nmin = 1; Nmax = 200;
for i = Nmin:Nmax
if cnt == 0
try
img = np.load("$path""img/targetimg_highres_Coord_SZ_$i""_HFI.npy");
cnt += 1
Nsize = size(img)[1]
end
else
try
img += np.rot90(np.load("$path""img/targetimg_highres_Coord_SZ_$i""_HFI.npy"),Int(floor(rand()*4.)))
cnt += 1
catch
end
end
end
img /= cnt
vmin = np.min(img); vmax = np.max(img)
vmin = -30.; vmax = 100.
#vmax = 0.5*(abs(vmin)+abs(vmax))
#vmin = -vmax
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
PlotTmap(img,vmin,vmax,dtheta_deg)
PyPlot.title("Stacked SZ sources",fontsize=15)
μ = mean(img); noise = std(img)
Ntot = copy(cnt);
println("Stacked $cnt images: mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,img;vmin=0.,vmax=60.)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [33]:
imgrand = zeros(Nsize,Nsize)
Nmin = 1; Nmax = 100; cnt =0;
ind_use = Array(Int64, 0)
add = 0
for i = 1:Nmax
try
np.load("$path""img/randimg_highres_CoordRand100_$i""_HFI.npy");
ind_use = vcat(ind_use, i)
catch
end
end
#ind_use = pick_random_ind(ind_use, Ntot)
for i = 1:Nmax
n = ind_use[i]
imgrand += np.rot90(np.load("$path""img/randimg_highres_CoordRand100_$n""_HFI.npy"),Int(floor(rand()*4.)))
#imgrand += np.load("$path""img/randimg_highres_CoordRand100_$n""_HFI.npy");
end
imgrand /= Nmax
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
PlotTmap(imgrand,vmin,vmax,dtheta_deg)
PyPlot.title("Stacked random positions",fontsize=15)
μ = mean(imgrand); noise = std(imgrand)
println("Stacked $Nmax images: mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,imgrand;vmin=-60,vmax=60)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [35]:
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
vmin = -50.; vmax = 50.;
PlotTmap(img-imgrand,vmin,vmax,dtheta_deg)
PyPlot.title("Residual image",fontsize=15)
μ = mean(img - imgrand); noise = std(img - imgrand)
println("mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,img-imgrand;vmin=-30.,vmax=30.)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [36]:
img = zeros(Nsize,Nsize)
cnt = 0
Nmin = 1; Nmax = 100;
for i = Nmin:Nmax
try
img += np.rot90(np.load("$path""img/targetimg_highres_Coord_GR08_sc_$i""_HFI.npy"),Int(floor(rand()*4.)))
cnt += 1
catch
end
end
img /= cnt
vmin = np.min(img); vmax = np.max(img)
vmin = -30.; vmax = 100.;
#vmax = 0.5*(abs(vmin)+abs(vmax))
#vmin = -vmax
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
PlotTmap(img,vmin,vmax,dtheta_deg)
PyPlot.title("Stacked GR08 sources",fontsize=15)
μ = mean(img); noise = std(img)
Ntot = copy(cnt);
println("Stacked $cnt images: mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,img;vmin=-20.,vmax=80)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [38]:
Nmin = 1; Nmax = 100; cnt =0;
ind_use = Array(Int64, 0)
add = 0
for i = 1:Nmax
try
np.load("$path""img/randimg_highres_CoordRand100_$i""_HFI.npy");
ind_use = vcat(ind_use, i)
catch
end
end
#imgrand = zeros(Nsize,Nsize)
#ind_use = pick_random_ind(ind_use, Ntot)
#for i = 1:Nmax
# n = ind_use[i]
# #imgrand += np.load("$path""img/randimg_highres_CoordRand100_$n""_HFI.npy");
# imgrand += np.rot90(np.load("$path""img/randimg_highres_CoordRand100_$n""_HFI.npy"),Int(floor(rand()*4.)))
#end
#imgrand /= Nmax
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
vmin = -50.; vmax = 50.;
PlotTmap(img-imgrand,vmin,vmax,dtheta_deg)
PyPlot.title("Residual image",fontsize=15)
μ = mean(img - imgrand); noise = std(img - imgrand)
println("mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,img - imgrand;vmin=-30.,vmax=30.)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [43]:
img = zeros(Nsize,Nsize)
cnt = 0
Nmin = 1; Nmax = 50;
for i = Nmin:Nmax
try
img += np.rot90(np.load("$path""img/targetimg_highres_Coord_GR08_void_$i""_HFI.npy"),Int(floor(rand()*4.)))
cnt += 1
catch
end
end
img /= cnt
vmin = np.min(img); vmax = np.max(img)
vmin = -30.; vmax = 100.;
#vmax = 0.5*(abs(vmin)+abs(vmax))
#vmin = -vmax
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
PlotTmap(img,vmin,vmax,dtheta_deg)
PyPlot.title("Stacked GR08 VOID sources",fontsize=15)
μ = mean(img); noise = std(img)
Ntot = copy(cnt);
println("Stacked $cnt images: mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,img;vmin=-20.,vmax=50.)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [44]:
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
vmin = -50.; vmax = 50.;
PlotTmap(img-imgrand,vmin,vmax,dtheta_deg)
PyPlot.title("Residual image",fontsize=15)
μ = mean(img - imgrand); noise = std(img - imgrand)
println("mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,img - imgrand;vmin=-50.,vmax=30.)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [48]:
SMICA="COM_CMB_IQU-smica-field-Int_2048_R2.01_full.fits"
I_SMICA = hp.read_map("$path$SMICA", field = 0, memmap = true);
In [298]:
hp.mollview(I_SMICA*1e6, xsize = 800, title = "(SMICA) Intensity Map (\$\\mu\$\K\$\_{CMB}\$\)", min = -300, max = 300)
In [48]:
img = zeros(Nsize,Nsize)
cnt = 0
Nmin = 1; Nmax = 200;
for i = Nmin:Nmax
try
img += np.rot90(np.load("$path""img/targetimg_highres_Coord_SZ_$i""_SMICA.npy"),Int(floor(rand()*4.)));
cnt += 1
catch
end
end
img /= cnt
vmin = np.min(img); vmax = np.max(img)
vmin = -80.; vmax = 100.
#vmax = 0.5*(abs(vmin)+abs(vmax))
#vmin = -vmax
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
PlotTmap(img,vmin,vmax,dtheta_deg)
PyPlot.title("Stacked SZ sources",fontsize=15)
μ = mean(img); noise = std(img)
Ntot = copy(cnt);
println("Stacked $cnt images: mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,img;vmin=-20.,vmax=20.)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [53]:
imgrand = zeros(Nsize,Nsize)
Nmin = 1; Nmax = 100; cnt =0;
ind_use = Array(Int64, 0)
add = 0
for i = 1:Nmax
try
np.load("$path""img/randimg_highres_CoordRand100_$i""_SMICA.npy");
ind_use = vcat(ind_use, i)
catch
end
end
#ind_use = pick_random_ind(ind_use, Ntot)
for i = 1:length(ind_use)
n = ind_use[i]
imgrand += np.rot90(np.load("$path""img/randimg_highres_CoordRand100_$n""_SMICA.npy"),Int(floor(rand()*4.)))
end
imgrand /= length(ind_use)
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
PlotTmap(imgrand,vmin,vmax,dtheta_deg)
PyPlot.title("Stacked random positions",fontsize=15)
cnt = length(ind_use)
μ = mean(imgrand); noise = std(imgrand)
println("Stacked $cnt images: mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,imgrand;vmin=-20.,vmax=20.)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [51]:
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
PlotTmap(img-imgrand,vmin,vmax,dtheta_deg)
PyPlot.title("Residual image",fontsize=15)
μ = mean(img - imgrand); noise = std(img - imgrand)
println("mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,img-imgrand;vmin=-30.,vmax=20.)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [13]:
img = zeros(Nsize,Nsize)
cnt = 0
Nmin = 1; Nmax = 100;
for i = Nmin:Nmax
try
img += np.rot90(np.load("$path""img/targetimg_highres_Coord_GR08_sc_$i""_SMICA.npy"),Int(floor(rand()*4.)))
cnt += 1
catch
end
end
img /= cnt
vmin = np.min(img); vmax = np.max(img)
vmin = -50.; vmax = 50.
#vmin = max(abs(vmin),abs(vmax))
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
PlotTmap(img,vmin,vmax,dtheta_deg)
PyPlot.title("Stacked GR08 sources",fontsize=15)
μ = mean(img); noise = std(img)
Ntot = copy(cnt);
println("Stacked $cnt images: mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,img;vmin=-20.,vmax=20.)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [14]:
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
PlotTmap(img-imgrand,vmin,vmax,dtheta_deg)
PyPlot.title("Residual image",fontsize=15)
μ = mean(img - imgrand); noise = std(img - imgrand)
println("mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,img-imgrand;vmin=-20.,vmax=20.)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [28]:
img = zeros(Nsize,Nsize)
cnt = 0
Nmin = 1; Nmax = 50;
for i = Nmin:Nmax
try
img += np.rot90(np.load("$path""img/targetimg_highres_Coord_GR08_void_$i""_SMICA.npy"),Int(floor(rand()*4.)))
cnt += 1
catch
end
end
img /= cnt
vmin = np.min(img); vmax = np.max(img)
vmin = -50.; vmax = 50.;
#vmax = 0.5*(abs(vmin)+abs(vmax))
#vmin = -vmax
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
PlotTmap(img,vmin,vmax,dtheta_deg)
PyPlot.title("Stacked GR08 VOID sources",fontsize=15)
μ = mean(img); noise = std(img)
Ntot = copy(cnt);
println("Stacked $cnt images: mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,img;vmin=-50.,vmax=20.)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
In [24]:
PyPlot.figure(figsize=(12,4))
PyPlot.subplot(1, 2, 1)
PlotTmap(img-imgrand,vmin,vmax,dtheta_deg)
PyPlot.title("Residual image",fontsize=15)
μ = mean(img - imgrand); noise = std(img - imgrand)
println("mean = $μ , std = $noise")
PyPlot.subplot(1, 2, 2)
PlotTempProfile(15.,0.2,img-imgrand;vmin=-50.,vmax=20.)
PyPlot.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
(in progress ... )
To check ISW effect on the CMB, I have used two different CMB maps: HFI map at 100 GHz and SMICA image which went through SMICA pipline. I stopped my analysis at this step because I could not see a clear difference between stacked images of superclusters/voids and random positions. The amount of excess or deficit in temperature is much smaller than GR08 paper or Planck results by two orders of magnitude. Possible reasons for this would be from:
1. coordinate transformation
2. resolution of stacked image
coordinate transformation : the coordinates of superclusters are given by Galactic coordinate system,(l, b) while Healpix uses spherical coordinates (θ, ϕ). I changed galactic coordinates to spherical ones by: ϕ [rad] = pi / 180 b [deg] θ [rad] = pi / 2 - pi / 180 l [deg] If the transformation is wrong, I stack wrong positions
resolution : Since the resolution of the stacked image is smaller than the Planck CMB map, I have taken the average of temperatures in the same pixel position. I think this would lead to small fluctuations of temperature as seen above. I want to see any changes if I raise the resolution as the same as the Planck CMB map (even after this quarter)
Although I upload my current progress to Github, I will check and fix the issues mentioned above.