ccl_correlation.c File Reference

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_legendre.h>
#include "fftlog.h"
#include "ccl.h"
#include "ccl_params.h"

Include dependency graph for ccl_correlation.c:

Go to the source code of this file.

Classes
struct	corr_int_par

Functions
static int	taper_cl (int n_ell, double *ell, double *cl, double *ell_limits)
static void	ccl_tracer_corr_fftlog (ccl_cosmology *cosmo, int n_ell, double *ell, double *cls, int n_theta, double *theta, double *wtheta, int corr_type, int do_taper_cl, double *taper_cl_limits, int *status)
static double	corr_bessel_integrand (double l, void *params)
static void	ccl_tracer_corr_bessel (ccl_cosmology *cosmo, int n_ell, double *ell, double *cls, int n_theta, double *theta, double *wtheta, int corr_type, int *status)
static void	ccl_compute_legendre_polynomial (int corr_type, double theta, int ell_max, double *Pl_theta)
static void	ccl_tracer_corr_legendre (ccl_cosmology *cosmo, int n_ell, double *ell, double *cls, int n_theta, double *theta, double *wtheta, int corr_type, int do_taper_cl, double *taper_cl_limits, int *status)
void	ccl_correlation (ccl_cosmology *cosmo, int n_ell, double *ell, double *cls, int n_theta, double *theta, double *wtheta, int corr_type, int do_taper_cl, double *taper_cl_limits, int flag_method, int *status)
void	ccl_correlation_3d (ccl_cosmology *cosmo, double a, int n_r, double *r, double *xi, int do_taper_pk, double *taper_pk_limits, int *status)

Function Documentation

static void ccl_compute_legendre_polynomial ( int  corr_type, double  theta, int  ell_max, double *  Pl_theta  )

static

Definition at line 267 of file ccl_correlation.c.

References CCL_CORR_GG, CCL_CORR_GL, and M_PI.

Referenced by ccl_tracer_corr_legendre().

{
  int i,j;
  double k=0;
  double cth=cos(theta*M_PI/180);
  
  //Initialize Pl_theta
  for (j=0;j<=ell_max;j++)
      Pl_theta[j]=0.;

  if(corr_type==CCL_CORR_GG) {
    gsl_sf_legendre_Pl_array(ell_max,cth,Pl_theta);
    for (j=0;j<=ell_max;j++)
      Pl_theta[j]*=(2*j+1);
  }
  else if(corr_type==CCL_CORR_GL) {
    for (j=2;j<=ell_max;j++) {//https://arxiv.org/pdf/1007.4809.pdf
      Pl_theta[j]=gsl_sf_legendre_Plm(j,2,cth);
      Pl_theta[j]*=(2*j+1.)/((j+0.)*(j+1.));
    }
  }
}

void ccl_correlation	(	ccl_cosmology *	cosmo,
		int	n_ell,
		double *	ell,
		double *	cls,
		int	n_theta,
		double *	theta,
		double *	wtheta,
		int	corr_type,
		int	do_taper_cl,
		double *	taper_cl_limits,
		int	flag_method,
		int *	status
	)

Computes the correlation function (wrapper)

Parameters

cosmo	:Cosmological parameters
n_ell	: number of multipoles in the input power spectrum
ell	: multipoles at which the power spectrum is evaluated
cls	: input power spectrum
n_theta	: number of output values of the separation angle (theta)
theta	: values of the separation angle in degrees.
wtheta	: the values of the correlation function at the angles above will be returned in this array, which should be pre-allocated
do_taper_cl	:
taper_cl_limits
flag_method	: method to compute the correlation function. Choose between: CCL_CORR_FFTLOG : fast integration with FFTLog CCL_CORR_BESSEL : direct integration over the Bessel function CCL_CORR_LGNDRE : brute-force sum over legendre polynomials
corr_type	: type of correlation function. Choose between: CCL_CORR_GG : spin0-spin0 CCL_CORR_GL : spin0-spin2 CCL_CORR_LP : spin2-spin2 (xi+) CCL_CORR_LM : spin2-spin2 (xi-) Currently supported spin-0 fields are number counts and CMB lensing. The only spin-2 is currently shear.

Definition at line 391 of file ccl_correlation.c.

References ccl_check_status(), CCL_CORR_BESSEL, CCL_CORR_FFTLOG, CCL_CORR_LGNDRE, ccl_cosmology_set_status_message(), CCL_ERROR_INCONSISTENT, ccl_tracer_corr_bessel(), ccl_tracer_corr_fftlog(), and ccl_tracer_corr_legendre().

Referenced by compare_corr(), and main().

{
  switch(flag_method) {
  case CCL_CORR_FFTLOG :
    ccl_tracer_corr_fftlog(cosmo,n_ell,ell,cls,n_theta,theta,wtheta,corr_type,
               do_taper_cl,taper_cl_limits,status);
    break;
  case CCL_CORR_LGNDRE :
    ccl_tracer_corr_legendre(cosmo,n_ell,ell,cls,n_theta,theta,wtheta,corr_type,
                 do_taper_cl,taper_cl_limits,status);
    break;
  case CCL_CORR_BESSEL :
    ccl_tracer_corr_bessel(cosmo,n_ell,ell,cls,n_theta,theta,wtheta,corr_type,status);
    break;
  default :
    *status=CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation. Unknown algorithm\n");
  }

  ccl_check_status(cosmo,status);
}

void ccl_correlation_3d	(	ccl_cosmology *	cosmo,
		double	a,
		int	n_r,
		double *	r,
		double *	xi,
		int	do_taper_pk,
		double *	taper_pk_limits,
		int *	status
	)

Computes the 3dcorrelation function (wrapper)

Parameters

cosmo	:Cosmological parameters
a	: scale factor
n_r	: number of output values of distance r
r	: values of the distance in Mpc
xi	: the values of the correlation function at the distances above will be returned in this array, which should be pre-allocated
do_taper_pk	: key for tapering (using cosine tapering by default)
taper_pk_limits	limits of tapering

Definition at line 427 of file ccl_correlation.c.

References ccl_check_status(), ccl_cosmology_set_status_message(), CCL_ERROR_MEMORY, ccl_log_spacing(), ccl_nonlin_matter_power(), ccl_spline_eval(), ccl_spline_free(), ccl_spline_init(), ccl_splines, ccl_spline_params::K_MAX, ccl_spline_params::K_MIN, ccl_spline_params::N_K_3DCOR, pk2xi(), and taper_cl().

Referenced by compare_correlation_3d(), and main().

{
  int i,N_ARR;
  double *k_arr,*pk_arr,*r_arr,*xi_arr;

  //number of data points for k and pk array
  N_ARR=(int)(ccl_splines->N_K_3DCOR*log10(ccl_splines->K_MAX/ccl_splines->K_MIN));  

  k_arr=ccl_log_spacing(ccl_splines->K_MIN,ccl_splines->K_MAX,N_ARR);
  if(k_arr==NULL) {
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }

  pk_arr=malloc(N_ARR*sizeof(double));
  if(pk_arr==NULL) {
    free(k_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }  

  for (i=0; i<N_ARR; i++)
    pk_arr[i] = ccl_nonlin_matter_power(cosmo, k_arr[i], a, status);

  if (do_taper_pk)
    taper_cl(N_ARR,k_arr,pk_arr,taper_pk_limits);

  r_arr=malloc(sizeof(double)*N_ARR);
  if(r_arr==NULL) {
    free(k_arr);
    free(pk_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }
  xi_arr=malloc(sizeof(double)*N_ARR);
  if(xi_arr==NULL) {
    free(k_arr); free(pk_arr); free(r_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }

  for(i=0;i<N_ARR;i++)
    r_arr[i]=0;

  pk2xi(N_ARR,k_arr,pk_arr,r_arr,xi_arr);

  // Interpolate to output values of r
  SplPar *xi_spl=ccl_spline_init(N_ARR,r_arr,xi_arr,xi_arr[0],0);
  for(i=0;i<n_r;i++)
    xi[i]=ccl_spline_eval(r[i],xi_spl);
  ccl_spline_free(xi_spl);

  free(k_arr); free(pk_arr);
  free(r_arr); free(xi_arr);

  ccl_check_status(cosmo,status);

  return;
}

static void ccl_tracer_corr_bessel ( ccl_cosmology *  cosmo, int  n_ell, double *  ell, double *  cls, int  n_theta, double *  theta, double *  wtheta, int  corr_type, int *  status  )

static

Definition at line 184 of file ccl_correlation.c.

References CCL_CORR_GG, CCL_CORR_GL, CCL_CORR_LM, CCL_CORR_LP, ccl_cosmology_set_status_message(), CCL_ERROR_MEMORY, ccl_gsl, ccl_raise_gsl_warning(), ccl_spline_free(), ccl_spline_init(), ccl_splines, corr_int_par::cl0, corr_int_par::cl_spl, corr_int_par::clf, corr_bessel_integrand(), corr_int_par::ell0, ccl_spline_params::ELL_MAX_CORR, corr_int_par::ellf, corr_int_par::extrapol_0, corr_int_par::extrapol_f, corr_int_par::i_bessel, ccl_gsl_params::INTEGRATION_EPSREL, ccl_gsl_params::INTEGRATION_GAUSS_KRONROD_POINTS, M_PI, ccl_gsl_params::N_ITERATION, corr_int_par::nell, corr_int_par::th, corr_int_par::tilt0, corr_int_par::tiltf, and w.

Referenced by ccl_correlation().

{
  corr_int_par *cp=malloc(sizeof(corr_int_par));
  if(cp==NULL) {
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_bessel ran out of memory\n");
    return;
  }

  cp->nell=n_ell;
  cp->ell0=ell[0];
  cp->ellf=ell[n_ell-1];
  cp->cl0=cls[0];
  cp->clf=cls[n_ell-1];
  cp->cl_spl=ccl_spline_init(n_ell,ell,cls,cls[0],0);
  if(cp->cl_spl==NULL) {
    free(cp);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_bessel ran out of memory\n");
    return;
  }
  switch(corr_type) {
  case CCL_CORR_GG :
    cp->i_bessel=0;
    break;
  case CCL_CORR_GL :
    cp->i_bessel=2;
    break;
  case CCL_CORR_LP :
    cp->i_bessel=0;
    break;
  case CCL_CORR_LM :
    cp->i_bessel=4;
    break;
  }

  if(cls[0]*cls[1]<=0)
    cp->extrapol_0=0;
  else {
    cp->extrapol_0=1;
    cp->tilt0=log10(cls[1]/cls[0])/log10(ell[1]/ell[0]);
  }

  if(cls[n_ell-2]*cls[n_ell-1]<=0)
    cp->extrapol_f=0;
  else {
    cp->extrapol_f=1;
    cp->tiltf=log10(cls[n_ell-1]/cls[n_ell-2])/log10(ell[n_ell-1]/ell[n_ell-2]);
  }

  int ith, gslstatus;
  double result,eresult;
  gsl_function F;
  gsl_integration_workspace *w=gsl_integration_workspace_alloc(ccl_gsl->N_ITERATION);
  for(ith=0;ith<n_theta;ith++) {
    cp->th=theta[ith]*M_PI/180;
    F.function=&corr_bessel_integrand;
    F.params=cp;
    //TODO: Split into intervals between first bessel zeros before integrating
    //This will help both speed and accuracy of the integral.
    gslstatus = gsl_integration_qag(&F, 0, ccl_splines->ELL_MAX_CORR, 0,
                                    ccl_gsl->INTEGRATION_EPSREL, ccl_gsl->N_ITERATION,
                                    ccl_gsl->INTEGRATION_GAUSS_KRONROD_POINTS,
                                    w, &result, &eresult);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_correlation.c: ccl_tracer_corr_bessel():");
      *status |= gslstatus;
    }
    wtheta[ith]=result/(2*M_PI);
  }
  gsl_integration_workspace_free(w);
  ccl_spline_free(cp->cl_spl);
  free(cp);
}

static void ccl_tracer_corr_fftlog ( ccl_cosmology *  cosmo, int  n_ell, double *  ell, double *  cls, int  n_theta, double *  theta, double *  wtheta, int  corr_type, int  do_taper_cl, double *  taper_cl_limits, int *  status  )

static

Definition at line 50 of file ccl_correlation.c.

References CCL_CORR_GG, CCL_CORR_GL, CCL_CORR_LM, CCL_CORR_LP, ccl_cosmology_set_status_message(), CCL_ERROR_LINSPACE, CCL_ERROR_MEMORY, ccl_log_spacing(), ccl_spline_eval(), ccl_spline_free(), ccl_spline_init(), ccl_splines, ccl_spline_params::ELL_MAX_CORR, ccl_spline_params::ELL_MIN_CORR, fftlog_ComputeXi2D(), M_PI, ccl_spline_params::N_ELL_CORR, pow(), and taper_cl().

Referenced by ccl_correlation().

{
  int i;
  double *l_arr,*cl_arr,*th_arr,*wth_arr;

  l_arr=ccl_log_spacing(ccl_splines->ELL_MIN_CORR,ccl_splines->ELL_MAX_CORR,ccl_splines->N_ELL_CORR);
  if(l_arr==NULL) {
    *status=CCL_ERROR_LINSPACE;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }
  cl_arr=malloc(ccl_splines->N_ELL_CORR*sizeof(double));
  if(cl_arr==NULL) {
    free(l_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }

  //Interpolate input Cl into array needed for FFTLog
  SplPar *cl_spl=ccl_spline_init(n_ell,ell,cls,cls[0],0);
  if(cl_spl==NULL) {
    free(l_arr);
    free(cl_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }

  double cl_tilt,l_edge,cl_edge;
  l_edge=ell[n_ell-1];
  if((cls[n_ell-1]*cls[n_ell-2]<0) || (cls[n_ell-2]==0)) {
    cl_tilt=0;
    cl_edge=0;
  }
  else {
    cl_tilt=log(cls[n_ell-1]/cls[n_ell-2])/log(ell[n_ell-1]/ell[n_ell-2]);
    cl_edge=cls[n_ell-1];
  }
  for(i=0;i<ccl_splines->N_ELL_CORR;i++) {
    if(l_arr[i]>=l_edge)
      cl_arr[i]=cl_edge*pow(l_arr[i]/l_edge,cl_tilt);
    else
      cl_arr[i]=ccl_spline_eval(l_arr[i],cl_spl);
  }
  ccl_spline_free(cl_spl);

  if (do_taper_cl)
    taper_cl(ccl_splines->N_ELL_CORR,l_arr,cl_arr,taper_cl_limits);

  th_arr=malloc(sizeof(double)*ccl_splines->N_ELL_CORR);
  if(th_arr==NULL) {
    free(l_arr);
    free(cl_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }
  wth_arr=(double *)malloc(sizeof(double)*ccl_splines->N_ELL_CORR);
  if(wth_arr==NULL) {
    free(l_arr); free(cl_arr); free(th_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }

  for(i=0;i<ccl_splines->N_ELL_CORR;i++)
    th_arr[i]=0;
  //Although set here to 0, theta is modified by FFTlog to obtain the correlation at ~1/l

  int i_bessel=0;
  if(corr_type==CCL_CORR_GG) i_bessel=0;
  if(corr_type==CCL_CORR_GL) i_bessel=2;
  if(corr_type==CCL_CORR_LP) i_bessel=0;
  if(corr_type==CCL_CORR_LM) i_bessel=4;
  fftlog_ComputeXi2D(i_bessel,ccl_splines->N_ELL_CORR,l_arr,cl_arr,th_arr,wth_arr);

  // Interpolate to output values of theta
  SplPar *wth_spl=ccl_spline_init(ccl_splines->N_ELL_CORR,th_arr,wth_arr,wth_arr[0],0);
  for(i=0;i<n_theta;i++)
    wtheta[i]=ccl_spline_eval(theta[i]*M_PI/180.,wth_spl);
  ccl_spline_free(wth_spl);

  free(l_arr); free(cl_arr);
  free(th_arr); free(wth_arr);

  return;
}

static void ccl_tracer_corr_legendre ( ccl_cosmology *  cosmo, int  n_ell, double *  ell, double *  cls, int  n_theta, double *  theta, double *  wtheta, int  corr_type, int  do_taper_cl, double *  taper_cl_limits, int *  status  )

static

Definition at line 295 of file ccl_correlation.c.

References ccl_compute_legendre_polynomial(), CCL_CORR_LM, CCL_CORR_LP, ccl_cosmology_set_status_message(), CCL_ERROR_MEMORY, CCL_ERROR_NOT_IMPLEMENTED, ccl_spline_eval(), ccl_spline_free(), ccl_spline_init(), ccl_splines, ccl_spline_params::ELL_MAX_CORR, cl_cmbl_bm::l, M_PI, pow(), and taper_cl().

Referenced by ccl_correlation().

{
  int i;
  double *l_arr,*cl_arr,*Pl_theta;
  SplPar *cl_spl;

  if(corr_type==CCL_CORR_LM || corr_type==CCL_CORR_LP){
    *status=CCL_ERROR_NOT_IMPLEMENTED;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: CCL does not support full-sky xi+- calcuations.\nhttps://arxiv.org/abs/1702.05301 indicates flat-sky to be sufficient.\n");
  }
  
  if(*status==0) {
    l_arr=malloc(((int)(ccl_splines->ELL_MAX_CORR)+1)*sizeof(double));
    if(l_arr==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_legendre ran out of memory\n");
    }
  }
  
  if(*status==0) {
    cl_arr=malloc(((int)(ccl_splines->ELL_MAX_CORR)+1)*sizeof(double));
    if(cl_arr==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_legendre ran out of memory\n");
    }
  }

  if(*status==0) {
    //Interpolate input Cl into 
    cl_spl=ccl_spline_init(n_ell,ell,cls,cls[0],0);
    if(cl_spl==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_legendre ran out of memory\n");
    }
  }

  if(*status==0) {
    double cl_tilt,l_edge,cl_edge;
    l_edge=ell[n_ell-1];
    if((cls[n_ell-1]*cls[n_ell-2]<0) || (cls[n_ell-2]==0)) {
      cl_tilt=0;
      cl_edge=0;
    }
    else {
      cl_tilt=log(cls[n_ell-1]/cls[n_ell-2])/log(ell[n_ell-1]/ell[n_ell-2]);
      cl_edge=cls[n_ell-1];
    }
    for(i=0;i<=(int)(ccl_splines->ELL_MAX_CORR);i++) {
      double l=(double)i;
      l_arr[i]=l;
      if(l>=l_edge)
    cl_arr[i]=cl_edge*pow(l/l_edge,cl_tilt);
      else
    cl_arr[i]=ccl_spline_eval(l,cl_spl);
    }
    ccl_spline_free(cl_spl);

    if (do_taper_cl)
      *status=taper_cl((int)(ccl_splines->ELL_MAX_CORR)+1,l_arr,cl_arr,taper_cl_limits);
  }

  if(*status==0) {
    Pl_theta=malloc(sizeof(double)*((int)(ccl_splines->ELL_MAX_CORR)+1));
    if(Pl_theta==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_legendre ran out of memory\n");
    }
  }
  
  if(*status==0) {
    for (int i=0;i<n_theta;i++) {
      wtheta[i]=0;
      ccl_compute_legendre_polynomial(corr_type,theta[i],(int)(ccl_splines->ELL_MAX_CORR),Pl_theta);
      for(int i_L=1;i_L<(int)(ccl_splines->ELL_MAX_CORR);i_L+=1)
    wtheta[i]+=cl_arr[i_L]*Pl_theta[i_L];
      wtheta[i]/=(M_PI*4);
    }
  }

  free(Pl_theta);
  free(l_arr);
  free(cl_arr);
}

static double corr_bessel_integrand ( double  l, void *  params  )

static

Definition at line 158 of file ccl_correlation.c.

References ccl_spline_eval(), corr_int_par::cl0, corr_int_par::cl_spl, corr_int_par::clf, corr_int_par::ell0, corr_int_par::ellf, corr_int_par::extrapol_0, corr_int_par::extrapol_f, corr_int_par::i_bessel, p, pow(), corr_int_par::th, corr_int_par::tilt0, corr_int_par::tiltf, and x.

Referenced by ccl_tracer_corr_bessel().

{
  double cl,jbes;
  corr_int_par *p=(corr_int_par *)params;
  double x=l*p->th;

  if(l<p->ell0) {
    if(p->extrapol_0)
      cl=p->cl0*pow(l/p->ell0,p->tilt0);
    else
      cl=0;
  }
  else if(l>p->ellf) {
    if(p->extrapol_f)
      cl=p->clf*pow(l/p->ellf,p->tiltf);
    else
      cl=0;
  }
  else
    cl=ccl_spline_eval(l,p->cl_spl);

  jbes=gsl_sf_bessel_Jn(p->i_bessel,x);

  return l*jbes*cl;
}

static int taper_cl ( int  n_ell, double *  ell, double *  cl, double *  ell_limits  )

static

Definition at line 23 of file ccl_correlation.c.

References M_PI.

Referenced by ccl_correlation_3d(), ccl_tracer_corr_fftlog(), and ccl_tracer_corr_legendre().

{

  for(int i=0;i<n_ell;i++) {
    if(ell[i]<ell_limits[0] || ell[i]>ell_limits[3]) {
      cl[i]=0;//ell outside desirable range
      continue;
    }
    if(ell[i]>=ell_limits[1] && ell[i]<=ell_limits[2])
      continue;//ell within good ell range
    
    if(ell[i]<ell_limits[1])//tapering low ell
      cl[i]*=cos((ell[i]-ell_limits[1])/(ell_limits[1]-ell_limits[0])*M_PI/2.);
    
    if(ell[i]>ell_limits[2])//tapering high ell
      cl[i]*=cos((ell[i]-ell_limits[2])/(ell_limits[3]-ell_limits[2])*M_PI/2.);
  }

  return 0;
}