SuperFaB
SuperFaB is a code for cosmological spherical Fourier-Bessel (SFB) analysis. The details of the code are presented in 2102.10079.
Installation
To install SuperFaB, start the Julia REPL and type
]add SphericalFourierBesselDecompositions
The package makes use of some python packages (e.g.
healpy) that are only supported on MacOSX
and Linux. The above command should work if healpy is already installed, but
if problems occur when first using the package, see
PyCall. Specifically, try
ENV["PYTHON"]=""; ]build PyCall
.
Basic Usage
Load the package and create a shortcut
julia> using SphericalFourierBesselDecompositions
julia> SFB = SphericalFourierBesselDecompositions
To perform a SFB decomposition, we create a modes object amodes
that contains
the modes and basis functions, and for the pseudo-SFB power spectrum we create
cmodes
,
julia> amodes = SFB.AnlmModes(kmax, rmin, rmax)
julia> cmodes = SFB.ClnnModes(amodes, Δnmax=0)
Here, kmax
is the maximum k-mode to be calculated, rmin
and rmax
are the
radial boundaries. Δnmax
specifies how many off-diagonals k ≠ k'
to
include.
The window function is described by
julia> wmodes = SFB.ConfigurationSpaceModes(rmin, rmax, nr, amodes.nside)
julia> win = SFB.SeparableArray(phi, mask, name1=:phi, name2=:mask)
julia> win_rhat_ln = SFB.win_rhat_ln(win, wmodes, amodes)
julia> Veff = SFB.integrate_window(win, wmodes)
where nr
is the number of radial bins, phi
is an array of length nr
, and
mask
is a HEALPix mask. In general, win
can be a 2D-array, where the first
dimension is radial, and the second dimension is the HEALPix mask at each
radius. Using a SeparableArray
uses Julia's dispatch mechanism to call more
efficient specialized algorithms when the radial and angular window are
separable. SFB.win_rhat_ln()
performs the radial transform of the window,
SFB.integrate_window()
is a convenient way to calculate the effective volume
Veff
.
The SFB decomposition is now performed with
julia> anlm = SFB.cat2amln(rθϕ, amodes, nbar, win_rhat_ln)
julia> CNobs = SFB.amln2clnn(anlm, cmodes)
where rθϕ
is a 3 × Ngalaxies
array with the r
, θ
, and ϕ
coordinates
of each galaxy in the survey, and nbar = Ngalaxies / Veff
is the average
number density. The second line calculates the pseudo-SFB power spectrum.
Shot noise and pixel window are calculated with
julia> Nobs_th = SFB.win_lnn(win, wmodes, cmodes) ./ nbar
julia> pixwin = SFB.pixwin(cmodes)
Window deconvolution is performed with bandpower binning:
julia> w̃mat, vmat = SFB.bandpower_binning_weights(cmodes; Δℓ=Δℓ, Δn=Δn)
julia> bcmodes = SFB.ClnnBinnedModes(w̃mat, vmat, cmodes)
julia> bcmix = SFB.power_win_mix(win, w̃mat, vmat, wmodes, bcmodes)
julia> Cobs = @. (CNobs - Nobs_th) / pixwin ^ 2
julia> C = bcmix \ (w̃mat * Cobs)
The first line calculates binning matrices w̃
and v
for bin sizes Δℓ ~ 1/fsky
and Δn = 1
, the second line describes modes similar to cmodes
but
for bandpower binned modes. The coupling matrix is calculated in the third
line, the shot noise and pixel window are corrected in the fourth line, and the
last line does the binning and deconvolves the window function.
To compare with a theoretical prediction, we calculate the deconvolved binning
matrix wmat
,
julia> using LinearAlgebra
julia> wmat = bcmix * SFB.power_win_mix(win, w̃mat, I, wmodes, bcmodes)
The modes of the pseudo-SFB power spectrum are given by
julia> lkk = SFB.getlkk(bcmodes)
where for a given i
the element lkk[1,i]
is the ℓ-mode, lkk[2,i]
is the
n
-mode, lkk[3,i]
is the n'
-mode of the pseudo-SFB power spectrum element
C[i]
.
An unoptimized way to calculate the covariance matrix is
julia> VW = SFB.calc_covariance_exact_chain(C_th, nbar, win, wmodes, cmodes)
julia> V = inv(bcmix) * w̃mat * VW * w̃mat' * inv(bcmix)'