TwoFAST

Master Branch:

Next Branch:

The 2-FAST (2-point function from Fast and Accurate Spherical bessel Transform) algorithm is implemented here in the Julia programming language. It computes integrals over one and two spherical Bessel functions as they frequently occur in cosmology.

The algorithm is documented in the paper Fast and Accurate Computation of Projected Two-point functions, arXiv link, ADS link.

To install in Julia-1.0, press ] to enter package mode, and then

   pkg> add TwoFAST

Minimal example

Load the module:

    using TwoFAST

For both minimal examples we need a power spectrum. For example, we can use the one in the test/ subdirectory of this project:

    using Dierckx
    using DelimitedFiles
    path = homedir() * "/.julia/packages/TwoFAST/"
    path *= readdir(path)[1]
    path *= "/test/data/planck_base_plikHM_TTTEEE_lowTEB_lensing_post_BAO_H070p6_JLA_matterpower.dat"
    d = readdlm(path, comments=true)
    pk = Spline1D(d[:,1], d[:,2])

To calculate the real-space correlation function, use

    N = 1024    # number of points to use in the Fourier transform
    kmax = 1e3  # maximum k-value
    kmin = 1e-5 # minimum k-value
    r0 = 1e-3   # minimum r-value (should be ~1/kmax)

    print("ξ(r), ℓ=0, ν=0: ")
    r00, xi00 = xicalc(pk, 0, 0; N=N, kmin=kmin, kmax=kmax, r0=r0)

    print("ξ(r), ℓ=0, ν=-2:")
    r, xi0m2 = xicalc(pk, 0, -2; N=N, kmin=kmin, kmax=kmax, r0=r0)

To calculate the integrals over two spherical Bessel functions, we first calculate the Fourier kernels at the highest needed ℓ. This is done with the structure F21EllCache. Then, we generate the full Mll-cache for each ℓ. This will automatically store the result in the file out/MlCache/MlCache.bin, and all related info will be stored in the structure MlCache. Finally, to actually calculate the wljj-terms we call the function calcwljj(). However, to store the wljj-terms, we need to create the output arrays, and write a function, outfunc(), that will store them in the arrays. The function outfunc() will be called for each ℓ in the array ell. Here's an example:

    N = 4096
    chi0 = 1e-3
    kmin = 1e-5
    kmax = 1e3
    q = 1.1
    ell = [42]  # only ell=42 for this run
    RR = [0.6, 0.7, 0.8, 0.9, 1.0]

    # calculate M_ll at high ell, result gets saved to a file:
    f21cache = F21EllCache(maximum(ell), RR, N; q=q, kmin=kmin, kmax=kmax, χ0=chi0)
    write("out/F21EllCache", f21cache)

    # calculate all M_ll, result gets saved to a file:
    mlcache = MlCache(ell, "out/F21EllCache", "out/MlCache")
    write("out/MlCache", mlcache)

    # calculate wljj:
    w00 = Array{Float64}(undef, N, length(RR))
    w02 = Array{Float64}(undef, N, length(RR))
    function outfunc(wjj, ell, rr, RR)
        if ell == 42
            w00[:,:] = wjj[1]
            w02[:,:] = wjj[2]
        end
    end
    rr = calcwljj(pk, RR; ell=ell, kmin=kmin, kmax=kmax, N=N, r0=chi0, q=q, outfunc=outfunc, cachefile="out/MlCache/MlCache.bin")