Limber Approximation Library & 2D interpolation library

The cosmosis-standard-library/shear/limber directory contains a C library that can be incorporated by more than one module. It has facilities for:

• 2D interpolation onto a regular grid (adapted from the cosmocalc code by Matt Becker)
• the Limber approximation integral

2D interpolation

The 2D interpolation code used here is based on that used in Matt Becker's cosmocalc code. It uses GSL to interpolate first along the x axis and then along the y, using (by default) the akima spline interpolation scheme.

We have made notes on the accuracy of this integrator for the specific case of the matter power spectrum P(k,z).

To use it you should include the "interp2d.h" and "utils.h" headers. In the C/C++ code:

#include "limber.h"
#include "utils.h"

USER_CFLAGS = -I${COSMOSIS_SRC_DIR}/cosmosis-standard-library/shear/limber
USER_LDFLAGS = -L${COSMOSIS_SRC_DIR}/cosmosis-standard-library/shear/limber -llimber

The basic idea is that you first load the set up the interpolator, with two 1D arrays, (e.g. k and z), and one that can be either 1D or 2D (e.g. P(k,z))

Interpolator2D * init_interp_2d_akima(double *x1, double *x2, double *y, int N1, int N2);
Interpolator2D * init_interp_2d_akima_grid(double *x1, double *x2, double **y, int N1, int N2);

Once you have this Interpolator2D object you can use this function to do the actual interpolation:

double interp_2d(double x1, double x2, Interpolator2D * interp2d);

And then when finished deallocate the memory with:

void destroy_interp_2d(Interpolator2D * interp2d);

There are also some helper functions to get splines and interpolator functions directly from a datablock:

gsl_spline * spline_from_arrays(int n, double * x, double *y);

DATABLOCK_STATUS save_spline(c_datablock * block, const char * section, 
    const char * n_name, const char * x_name, const char * y_name,
    gsl_spline * s);

gsl_spline * load_spline(c_datablock * block, const char * section, 
    const char * x_name, const char * y_name);

void reverse(double * x, int n);


// This loads P(k,z) directly
Interpolator2D * 
load_interpolator(c_datablock * block,
    const char * section,
    const char * k_name, const char * z_name, const char * P_name);

// This loads P(k,z) and turns it into P(k,chi)
Interpolator2D * 
load_interpolator_chi(c_datablock * block, gsl_spline * chi_of_z_spline, 
    const char * section,
    const char * k_name, const char * z_name, const char * P_name);

##The Limber approximation The Limber approximation projects the power spectrum of a 3D field into 2D by integrating over radial distance using some weighting (or kernel) function that depends on the application. The form of the approximation implemented here is:

where the two W functions are the kernel weights and P is the 3D power spectrum, chi is the radial coordinate, and ell the angular frequency.

You can call the library function by including the limber.h and utils.h headers from the cosmosis-standard-library/shear/limber directory and linking to liblimber in the same directory.

In the C/C++ code:

#include "limber.h"
#include "utils.h"

In a Makefile:

USER_CFLAGS = -I${COSMOSIS_SRC_DIR}/cosmosis-standard-library/shear/limber
USER_LDFLAGS = -L${COSMOSIS_SRC_DIR}/cosmosis-standard-library/shear/limber -llimber

The library is used by the module spectra, which computes the shear-shear power spectra from P(k,z).

The main function to use is:

gsl_spline * limber_integral(limber_config * config, gsl_spline * WX, gsl_spline * WY, Interpolator2D * P);

You need to create two GSL spline objects to represent the WX and WY kernels, a limber_config structure and an Interpolator2D which can be loaded from a data block using the load_interpolator function in utils.h or created in a variety of other ways - see interp2d.h. The config structure is defined like this:

typedef struct limber_config{
    bool xlog;  // The output spline will be in terms of log(ell) not ell
    bool ylog;  // The output spline will return log(C_ell) not C_ell
    int n_ell;  // Number of ell values you want in the spline
    double * ell;  // The chosen ell values you want
    double prefactor; //Scaling prefactor
} limber_config;

and you need to fill in all the parts of it to tell the library what options you want.

Several full examples of this usage can be found in the cosmosis-standard-library/shear/spectra/interface.c file.