ccl_utils.c File Reference

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include "ccl.h"

Include dependency graph for ccl_utils.c:

Go to the source code of this file.

Macros
#define	CCL_GAMMA1   2.6789385347077476336556
#define	CCL_GAMMA2   1.3541179394264004169452
#define	CCL_ROOTPI12   21.269446210866192327578

Functions
double *	ccl_linear_spacing (double xmin, double xmax, int N)
double *	ccl_linlog_spacing (double xminlog, double xmin, double xmax, int Nlog, int Nlin)
double *	ccl_log_spacing (double xmin, double xmax, int N)
SplPar *	ccl_spline_init (int n, double *x, double *y, double y0, double yf)
double	ccl_spline_eval (double x, SplPar *spl)
double	ccl_j_bessel (int l, double x)
void	ccl_spline_free (SplPar *spl)

Macro Definition Documentation

#define CCL_GAMMA1   2.6789385347077476336556

Definition at line 187 of file ccl_utils.c.

Referenced by ccl_j_bessel().

#define CCL_GAMMA2   1.3541179394264004169452

Definition at line 188 of file ccl_utils.c.

Referenced by ccl_j_bessel().

#define CCL_ROOTPI12   21.269446210866192327578

Definition at line 189 of file ccl_utils.c.

Referenced by ccl_j_bessel().

Function Documentation

double ccl_j_bessel	(	int	l,
		double	x
	)

Definition at line 190 of file ccl_utils.c.

References beta, beta2, CCL_GAMMA1, CCL_GAMMA2, CCL_ROOTPI12, M_PI, pow(), sqrt(), and x.

Referenced by __ctest_spherical_bessel_tests_compare_gsl_run().

{
  double jl;
  double ax=fabs(x);
  double ax2=x*x;
  if(l<0) {
    fprintf(stderr,"CosmoMas: l>0 for Bessel function");
    exit(1);
  }

  if(l<7) {
    if(l==0) {
      if(ax<0.1) jl=1-ax2*(1-ax2/20.)/6.;
      else jl=sin(x)/x;
    }
    else if(l==1) {
      if(ax<0.2) jl=ax*(1-ax2*(1-ax2/28)/10)/3;
      else jl=(sin(x)/ax-cos(x))/ax;
    }
    else if(l==2) {
      if(ax<0.3) jl=ax2*(1-ax2*(1-ax2/36)/14)/15;
      else jl=(-3*cos(x)/ax-sin(x)*(1-3/ax2))/ax;
    }
    else if(l==3) {
      if(ax<0.4)
    jl=ax*ax2*(1-ax2*(1-ax2/44)/18)/105;
      else
    jl=(cos(x)*(1-15/ax2)-sin(x)*(6-15/ax2)/ax)/ax;
    }
    else if(l==4) {
      if(ax<0.6)
    jl=ax2*ax2*(1-ax2*(1-ax2/52)/22)/945;
      else
    jl=(sin(x)*(1-(45-105/ax2)/ax2)+cos(x)*(10-105/ax2)/ax)/ax;
    }
    else if(l==5) {
      if(ax<1.0)
    jl=ax2*ax2*ax*(1-ax2*(1-ax2/60)/26)/10395;
      else {
    jl=(sin(x)*(15-(420-945/ax2)/ax2)/ax-
        cos(x)*(1-(105-945/ax2)/ax2))/ax;
      }
    }
    else {
      if(ax<1.0)
    jl=ax2*ax2*ax2*(1-ax2*(1-ax2/68)/30)/135135;
      else {
    jl=(sin(x)*(-1+(210-(4725-10395/ax2)/ax2)/ax2)+
        cos(x)*(-21+(1260-10395/ax2)/ax2)/ax)/ax;
      }
    }
  }
  else {
    double nu=l+0.5;
    double nu2=nu*nu;

    if(ax<1.0E-40) jl=0;
    else if((ax2/l)<0.5) {
      jl=(exp(l*log(ax/nu)-M_LN2+nu*(1-M_LN2)-(1-(1-3.5/nu2)/(30*nu2))/(12*nu))/nu)*
    (1-ax2/(4*nu+4)*(1-ax2/(8*nu+16)*(1-ax2/(12*nu+36))));
    }
    else if((l*l/ax)<0.5) {
      double beta=ax-0.5*M_PI*(l+1);
      jl=(cos(beta)*(1-(nu2-0.25)*(nu2-2.25)/(8*ax2)*(1-(nu2-6.25)*(nu2-12.25)/(48*ax2)))-
      sin(beta)*(nu2-0.25)/(2*ax)*(1-(nu2-2.25)*(nu2-6.25)/(24*ax2)*
                       (1-(nu2-12.25)*(nu2-20.25)/(80*ax2))))/ax;
    }
    else {
      double l3=pow(nu,0.325);
      if(ax<nu-1.31*l3) {
    double cosb=nu/ax;
    double sx=sqrt(nu2-ax2);
    double cotb=nu/sx;
    double secb=ax/nu;
    double beta=log(cosb+sx/ax);
    double cot3b=cotb*cotb*cotb;
    double cot6b=cot3b*cot3b;
    double sec2b=secb*secb;
    double expterm=((2+3*sec2b)*cot3b/24
            -((4+sec2b)*sec2b*cot6b/16
              +((16-(1512+(3654+375*sec2b)*sec2b)*sec2b)*cot3b/5760
                +(32+(288+(232+13*sec2b)*sec2b)*sec2b)*sec2b*cot6b/(128*nu))*
              cot6b/nu)/nu)/nu;
    jl=sqrt(cotb*cosb)/(2*nu)*exp(-nu*beta+nu/cotb-expterm);
      }
      else if(ax>nu+1.48*l3) {
    double cosb=nu/ax;
    double sx=sqrt(ax2-nu2);
    double cotb=nu/sx;
    double secb=ax/nu;
    double beta=acos(cosb);
    double cot3b=cotb*cotb*cotb;
    double cot6b=cot3b*cot3b;
    double sec2b=secb*secb;
    double trigarg=nu/cotb-nu*beta-0.25*M_PI-
      ((2+3*sec2b)*cot3b/24+(16-(1512+(3654+375*sec2b)*sec2b)*sec2b)*
       cot3b*cot6b/(5760*nu2))/nu;
    double expterm=((4+sec2b)*sec2b*cot6b/16-
            (32+(288+(232+13*sec2b)*sec2b)*sec2b)*
            sec2b*cot6b*cot6b/(128*nu2))/nu2;
    jl=sqrt(cotb*cosb)/nu*exp(-expterm)*cos(trigarg);
      }
      else {
    double beta=ax-nu;
    double beta2=beta*beta;
    double sx=6/ax;
    double sx2=sx*sx;
    double secb=pow(sx,0.3333333333333333333333);
    double sec2b=secb*secb;

    jl=(CCL_GAMMA1*secb+beta*CCL_GAMMA2*sec2b
        -(beta2/18-1.0/45.0)*beta*sx*secb*CCL_GAMMA1
        -((beta2-1)*beta2/36+1.0/420.0)*sx*sec2b*CCL_GAMMA2
        +(((beta2/1620-7.0/3240.0)*beta2+1.0/648.0)*beta2-1.0/8100.0)*sx2*secb*CCL_GAMMA1
        +(((beta2/4536-1.0/810.0)*beta2+19.0/11340.0)*beta2-13.0/28350.0)*beta*sx2*sec2b*CCL_GAMMA2
        -((((beta2/349920-1.0/29160.0)*beta2+71.0/583200.0)*beta2-121.0/874800.0)*
          beta2+7939.0/224532000.0)*beta*sx2*sx*secb*CCL_GAMMA1)*sqrt(sx)/CCL_ROOTPI12;
      }
    }
  }
  if((x<0)&&(l%2!=0)) jl=-jl;

  return jl;
}

double* ccl_linear_spacing	(	double	xmin,
		double	xmax,
		int	N
	)

Compute bin edges of N-1 linearly spaced bins on the interval [xmin,xmax]

Parameters

xmin	minimum value of spacing
xmax	maximum value of spacing
N	number of bins plus one (number of bin edges)

Returns

x, bin edges in range [xmin, xmax]

Definition at line 14 of file ccl_utils.c.

References CCL_ERROR_MEMORY, ccl_raise_warning(), x, xmax, and xmin.

Referenced by __ctest_spacing_tests_linear_spacing_simple_run(), calculate_nu_phasespace_spline(), ccl_cosmology_compute_distances(), ccl_cosmology_compute_power_emu(), ccl_cosmology_compute_sigma(), clt_init_wL(), and clt_init_wM().

{
  double dx = (xmax-xmin)/(N -1.);

  double * x = malloc(sizeof(double)*N);
  if (x==NULL) {
    ccl_raise_warning(
      CCL_ERROR_MEMORY,
      "ERROR: Could not allocate memory for linear-spaced array (N=%d)\n", N);
    return x;
  }

  for (int i=0; i<N; i++) {
    x[i] = xmin + dx*i;
  }
  x[0]=xmin; //Make sure roundoff errors don't spoil edges
  x[N-1]=xmax; //Make sure roundoff errors don't spoil edges

  return x;
}

double* ccl_linlog_spacing	(	double	xminlog,
		double	xmin,
		double	xmax,
		int	Nlin,
		int	Nlog
	)

Compute bin edges of N-1 logarithmically and then linearly spaced bins on the interval [xmin,xmax]

Parameters

xminlog	minimum value of spacing
xmin	value when logarithmical ends and linear spacing begins
xmax	maximum value of spacing
Nlin	number of linear bins plus one (number of bin edges)
Nlog	number of logarithmic bins plus one (number of bin edges)

Returns

x, bin edges in range [xminlog, xmax]

Definition at line 43 of file ccl_utils.c.

References CCL_ERROR_LINLOGSPACE, CCL_ERROR_MEMORY, ccl_raise_warning(), x, xmax, and xmin.

Referenced by __ctest_spacing_tests_linlog_spacing_simple_run(), ccl_cosmology_compute_distances(), ccl_cosmology_compute_growth(), ccl_cosmology_compute_power_bbks(), ccl_cosmology_compute_power_class(), ccl_cosmology_compute_power_eh(), and ccl_cosmology_compute_power_emu().

{
  if (Nlog<2) {
    ccl_raise_warning(
      CCL_ERROR_LINLOGSPACE,
      "ERROR: Cannot make log-spaced array with %d points - need at least 2\n", Nlog);
    return NULL;
  }

  if (!(xminlog>0 && xmin>0)) {
    ccl_raise_warning(
      CCL_ERROR_LINLOGSPACE,
      "ERROR: Cannot make log-spaced array xminlog or xmin non-positive (had %le, %le)\n", xminlog, xmin);
    return NULL;
  }

  if (xminlog>xmin){
    ccl_raise_warning(CCL_ERROR_LINLOGSPACE, "ERROR: xminlog must be smaller as xmin");
    return NULL;
  }

  if (xmin>xmax){
    ccl_raise_warning(CCL_ERROR_LINLOGSPACE, "ERROR: xmin must be smaller as xmax");
    return NULL;
  }

  double * x = malloc(sizeof(double)*(Nlin+Nlog-1));
  if (x==NULL) {
    ccl_raise_warning(
      CCL_ERROR_MEMORY,
      "ERROR: Could not allocate memory for array of size (Nlin+Nlog-1)=%d)\n", (Nlin+Nlog-1));
    return x;
  }

  double dx = (xmax-xmin)/(Nlin -1.);
  double log_xchange = log(xmin);
  double log_xmin = log(xminlog);
  double dlog_x = (log_xchange - log_xmin) /  (Nlog-1.);

  for (int i=0; i<Nlin+Nlog-1; i++) {
    if (i<Nlog)
        x[i] = exp(log_xmin + dlog_x*i);
    if (i>=Nlog)
        x[i] = xmin + dx*(i-Nlog+1);
  }

  x[0]=xminlog; //Make sure roundoff errors don't spoil edges
  x[Nlog-1]=xmin; //Make sure roundoff errors don't spoil edges
  x[Nlin+Nlog-2]=xmax; //Make sure roundoff errors don't spoil edges

  return x;
}

double* ccl_log_spacing	(	double	xmin,
		double	xmax,
		int	N
	)

Compute bin edges of N-1 logarithmically spaced bins on the interval [xmin,xmax]

Parameters

xmin	minimum value of spacing
xmax	maximum value of spacing
N	number of bins plus one (number of bin edges)

Returns

x, bin edges in range [xmin, xmax]

Definition at line 102 of file ccl_utils.c.

References CCL_ERROR_LOGSPACE, CCL_ERROR_MEMORY, ccl_raise_warning(), x, xmax, and xmin.

Referenced by __ctest_spacing_tests_log_spacing_simple_run(), ccl_correlation_3d(), ccl_cosmology_compute_power_bbks(), ccl_cosmology_compute_power_class(), ccl_cosmology_compute_power_eh(), ccl_cosmology_compute_power_emu(), ccl_tracer_corr_fftlog(), and main().

{
  if (N<2) {
    ccl_raise_warning(
      CCL_ERROR_LOGSPACE,
      "ERROR: Cannot make log-spaced array with %d points - need at least 2\n", N);
    return NULL;
  }

  if (!(xmin>0 && xmax>0)) {
    ccl_raise_warning(
      CCL_ERROR_LOGSPACE,
      "ERROR: Cannot make log-spaced array xmax or xmax non-positive (had %le, %le)\n", xmin, xmax);
    return NULL;
  }

  double log_xmax = log(xmax);
  double log_xmin = log(xmin);
  double dlog_x = (log_xmax - log_xmin) /  (N-1.);

  double * x = malloc(sizeof(double)*N);
  if (x==NULL) {
    ccl_raise_warning(
      CCL_ERROR_MEMORY,
      "ERROR: Could not allocate memory for log-spaced array (N=%d)\n", N);
    return x;
  }

  double xratio = exp(dlog_x);
  x[0] = xmin; //Make sure roundoff errors don't spoil edges
  for (int i=1; i<N-1; i++) {
    x[i] = x[i-1] * xratio;
  }
  x[N-1]=xmax; //Make sure roundoff errors don't spoil edges

  return x;
}

double ccl_spline_eval	(	double	x,
		SplPar *	spl
	)

Definition at line 170 of file ccl_utils.c.

References ccl_raise_gsl_warning(), SplPar::intacc, SplPar::spline, SplPar::xf, SplPar::y0, and SplPar::yf.

Referenced by ccl_angular_cls(), ccl_correlation_3d(), ccl_get_tracer_fa(), ccl_get_tracer_fas(), ccl_tracer_corr_fftlog(), ccl_tracer_corr_legendre(), corr_bessel_integrand(), f_dens(), f_IA_NLA(), f_lensing(), f_mag(), f_rsd(), integrand_mag(), integrand_wl(), and speval_bis().

{
  if(x<=spl->x0)
    return spl->y0;
  else if(x>=spl->xf)
    return spl->yf;
  else {
    double y;
    int stat=gsl_spline_eval_e(spl->spline,x,spl->intacc,&y);
    if (stat!=GSL_SUCCESS) {
      ccl_raise_gsl_warning(stat, "ccl_utils.c: ccl_splin_eval():");
      return NAN;
    }
    return y;
  }
}

void ccl_spline_free

(

SplPar *

spl

)

Definition at line 316 of file ccl_utils.c.

References SplPar::intacc, and SplPar::spline.

Referenced by ccl_angular_cls(), ccl_cl_tracer_free(), ccl_correlation_3d(), ccl_tracer_corr_bessel(), ccl_tracer_corr_fftlog(), ccl_tracer_corr_legendre(), and clt_init_nz().

{
  gsl_spline_free(spl->spline);
  gsl_interp_accel_free(spl->intacc);
  free(spl);
}

SplPar* ccl_spline_init	(	int	n,
		double *	x,
		double *	y,
		double	y0,
		double	yf
	)

Definition at line 146 of file ccl_utils.c.

References SplPar::intacc, SplPar::spline, SplPar::x0, SplPar::xf, SplPar::y0, and SplPar::yf.

Referenced by ccl_angular_cls(), ccl_correlation_3d(), ccl_tracer_corr_bessel(), ccl_tracer_corr_fftlog(), ccl_tracer_corr_legendre(), clt_init_ba(), clt_init_bz(), clt_init_nz(), clt_init_rf(), clt_init_wL(), and clt_init_wM().

{
  SplPar *spl=malloc(sizeof(SplPar));
  if(spl==NULL)
    return NULL;

  spl->intacc=gsl_interp_accel_alloc();
  spl->spline=gsl_spline_alloc(gsl_interp_cspline,n);
  int parstatus=gsl_spline_init(spl->spline,x,y,n);
  if(parstatus) {
    gsl_interp_accel_free(spl->intacc);
    gsl_spline_free(spl->spline);
    return NULL;
  }

  spl->x0=x[0];
  spl->xf=x[n-1];
  spl->y0=y0;
  spl->yf=yf;

  return spl;
}