# TwoFAST

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")
```