https://github.com/EthanAnderes/BayesianCmbLensing
Tip revision: 4f7a07119712138b0974e7564dbcdf4fd0f25092 authored by Ethan Anderes on 21 October 2015, 19:14:35 UTC
replaced v0.3+ with v0.3 since the code isn't currently supported for v0.4
replaced v0.3+ with v0.3 since the code isn't currently supported for v0.4
Tip revision: 4f7a071
cmb.jl
using PyCall, Interp
@pyimport scipy.interpolate as scii
## pixel grid and Fourier conjugate grid in radians
immutable Grid_xandk
pixel_size_arcmin::Float64
n::Float64
deltx::Float64
deltk::Float64
period::Float64
nyq::Float64
kside::Array{Float64,1}
x::Array{Float64,2}
y::Array{Float64,2}
k1::Array{Float64,2}
k2::Array{Float64,2}
r::Array{Float64,2}
end
function Grid_xandk(pixel_size_arcmin, n)
deltx = pixel_size_arcmin * pi / (180 * 60) #this is in radians
x,y = meshgrid(1:n, 1:n)
x .*= deltx
x .-= minimum(x)
y .*= deltx
y .-= minimum(y)
kside = (2 * pi / deltx) .* (mod(1/2 .+ [0:(n-1)]./n, 1) .- 1/2)
k1,k2 = meshgrid(kside, kside)
period = deltx * n
deltk = 2 * pi / period
nyq = 2 * pi / (2 * deltx)
r = sqrt(k1.^2 .+ k2.^2)
return Grid_xandk(pixel_size_arcmin, n, deltx, deltk, period, nyq, kside, x, y, k1, k2, r)
end
immutable SpectrumGrids # this contains all the signal + noise spectral information
grd::Grid_xandk # everything is defined in reference to a grid
ell::Array{Float64,1}
ellLen::Array{Float64,1}
#------------------------
cMaskBool::BitArray{2} # this gives l_max for cmb
pMaskBool::BitArray{2} # this gives l_max for phi
cMaskInf::Array{Float64,2}
#---------- noise
nugget_at_each_pixel::Float64
beamFWHM::Float64
beam::Array{Float64,2}
# radial profile
CTell2d::Array{Float64,1}
CTell2dLen::Array{Float64,1}
CPell2d::Array{Float64,1}
# matrices
cTT::Array{Float64,2}
cTTLen::Array{Float64,2}
cPP::Array{Float64,2}
cNT::Array{Float64,2}
end
function SpectrumGrids(
ell,
ellLen,
CPell2d,
CTell2d,
CTell2dLen,
pixel_size_arcmin,
n,
beamFWHM,
nugget_at_each_pixel,
maskupC,
maskupP
)
d = 2.0
grd = Grid_xandk(pixel_size_arcmin, n)
sig=(beamFWHM * (pi / (180 * 60))) / ( 2 * sqrt(2 * log(2)))
beam = exp(-(sig^2) .* (grd.k1.^2 .+ grd.k2.^2) ./ 2)
beamSQ =abs2(beam)
cNT = nugget_at_each_pixel .* (grd.deltx^d) ./ beamSQ
cMaskBool = !((grd.r.<= maskupC) & (grd.r.>= 1.0)) # this is true when the masking is active Mat[cMaskBool] = 0 or Inf
pMaskBool = !((grd.r.<= maskupP) & (grd.r.>= 1.0)) # this is true when the masking is active Mat[cMaskBool] = 0 or Inf
cMaskInf = 1.0./((grd.r.<=maskupC) & (grd.r.>= 1.0)) # this equals Inf when the masking is active and 1.0 otherwise...
index=ceil(grd.r)
index[find(index.==0)]=1
logCPP = linear_interp1(ell,log(CPell2d), index)
logCPP[find(logCPP .== 0)] =-Inf
logCPP[find(isnan(logCPP))]=-Inf
cPP = exp(logCPP)
logCTT = linear_interp1(ell,log(CTell2d), index)
logCTT[find(logCTT .== 0)] =-Inf
logCTT[find(isnan(logCTT))]=-Inf
cTT = exp(logCTT)
logCTTLen = linear_interp1(ellLen,log(CTell2dLen), index)
logCTTLen[find(logCTTLen .== 0)] =-Inf
logCTTLen[find(isnan(logCTTLen))]=-Inf
logCTTLen[grd.r .> 8000]=-Inf
cTTLen = exp(logCTTLen)
inputs = {grd,
ell,
ellLen,
cMaskBool,
pMaskBool,
cMaskInf,
nugget_at_each_pixel,
beamFWHM,
beam,
CTell2d,
CTell2dLen,
CPell2d,
cTT,
cTTLen,
cPP,
cNT }
SpectrumGrids(inputs...)
end
# get spectra and an instance of SpectrumGrids
function setpar(pixel_size_arcmin, n, beamFWHM, nugget_at_each_pixel, maskupC, maskupP)
Sfile = open("src/camb/test_scalCls.dat")
lines = readlines(Sfile)
close(Sfile)
CLs = Array(Any,(length(lines),6))
for i in 1:length(lines)
CLs[i,:] = lines[i] |> strip |> chomp |> split |> transpose
end
CLs = map(parsefloat, CLs)
ell=[1;CLs[:,1]]
CTell2d=[0;CLs[:,2].*(2*pi)./(CLs[:,1].*(CLs[:,1] .+ 1))]
CPell2d=[0;CLs[:,5]./(7.4311e12*CLs[:,1].^4)]
LenSfile = open("src/camb/test_lensedCls.dat")
Lenlines = readlines(LenSfile)
close(LenSfile)
LenCLs = Array(Any,(length(Lenlines),5))
for i in 1:length(Lenlines)
LenCLs[i,:] = Lenlines[i] |> strip |> chomp |> split |> transpose
end
LenCLs = map(parsefloat, LenCLs) #[ell, CTT, CEE, CBB, CTE]
ellLen = [1;LenCLs[:,1]]
CTell2dLen =[0;LenCLs[:,2].*(2*pi)./(LenCLs[:,1].*(LenCLs[:,1].+1))]
inputs = {ell,
ellLen,
CPell2d,
CTell2d,
CTell2dLen,
pixel_size_arcmin,
n,
beamFWHM,
nugget_at_each_pixel,
maskupC,
maskupP}
SpectrumGrids(inputs...)
end
# if you need to specify where src is...used in graphics
function setpar(pixel_size_arcmin, n, beamFWHM, nugget_at_each_pixel, maskupC, maskupP, pathtosrc::String)
Sfile = open(joinpath(pathtosrc,"camb/test_scalCls.dat"))
lines = readlines(Sfile)
close(Sfile)
CLs = Array(Any,(length(lines),6))
for i in 1:length(lines)
CLs[i,:] = lines[i] |> strip |> chomp |> split |> transpose
end
CLs = map(parsefloat, CLs)
ell=[1;CLs[:,1]]
CTell2d=[0;CLs[:,2].*(2*pi)./(CLs[:,1].*(CLs[:,1] .+ 1))]
CPell2d=[0;CLs[:,5]./(7.4311e12*CLs[:,1].^4)]
LenSfile = open(joinpath(pathtosrc,"camb/test_lensedCls.dat"))
Lenlines = readlines(LenSfile)
close(LenSfile)
LenCLs = Array(Any,(length(Lenlines),5))
for i in 1:length(Lenlines)
LenCLs[i,:] = Lenlines[i] |> strip |> chomp |> split |> transpose
end
LenCLs = map(parsefloat, LenCLs) #[ell, CTT, CEE, CBB, CTE]
ellLen = [1;LenCLs[:,1]]
CTell2dLen =[0;LenCLs[:,2].*(2*pi)./(LenCLs[:,1].*(LenCLs[:,1].+1))]
inputs = {ell,
ellLen,
CPell2d,
CTell2d,
CTell2dLen,
pixel_size_arcmin,
n,
beamFWHM,
nugget_at_each_pixel,
maskupC,
maskupP}
SpectrumGrids(inputs...)
end