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)'