Commit 779fc65d authored by Maria Archidiacono's avatar Maria Archidiacono
Browse files

sigma8 and sigma8_cb now computed in two separate functions to avoid confusion

parent 96391c20
......@@ -425,12 +425,20 @@ extern "C" {
struct background * pba,
struct primordial * ppm,
struct spectra * psp,
short compute_sigma8_cb,
double R,
double z,
double *sigma
);
int spectra_sigma_cb(
struct background * pba,
struct primordial * ppm,
struct spectra * psp,
double R,
double z,
double *sigma_cb
);
int spectra_matter_transfers(
struct background * pba,
struct perturbs * ppt,
......
......@@ -191,6 +191,7 @@ cdef extern from "class.h":
int index_md_scalars
double* ln_k
double sigma8
double sigma8_cb
double alpha_II_2_20
double alpha_RI_2_20
double alpha_RR_2_20
......@@ -339,11 +340,18 @@ cdef extern from "class.h":
void * pba,
void * ppm,
void * psp,
short compute_sigma8_cb,
double R,
double z,
double * sigma)
int spectra_sigma_cb(
void * pba,
void * ppm,
void * psp,
double R,
double z,
double * sigma_cb)
int spectra_fast_pk_at_kvec_and_zvec(
void * pba,
void * psp,
......
......@@ -911,12 +911,12 @@ cdef class Class:
"In order to get sigma(R,z) you must set 'P_k_max_h/Mpc' to 1 or bigger, in order to have k_max > 1 h/Mpc."
)
if spectra_sigma(&self.ba,&self.pm,&self.sp,_FALSE_,R,z,&sigma)==_FAILURE_:
if spectra_sigma(&self.ba,&self.pm,&self.sp,R,z,&sigma)==_FAILURE_:
raise CosmoSevereError(self.sp.error_message)
return sigma
# Gives sigma(R,z) for a given (R,z)
# Gives sigma_cb(R,z) for a given (R,z)
def sigma_cb(self,double R,double z):
"""
Gives the pk for a given R and z
......@@ -946,7 +946,7 @@ cdef class Class:
"No massive neutrinos. You must use sigma, rather than sigma_cb."
)
if spectra_sigma(&self.ba,&self.pm,&self.sp,_TRUE_,R,z,&sigma_cb)==_FAILURE_:
if spectra_sigma_cb(&self.ba,&self.pm,&self.sp,R,z,&sigma_cb)==_FAILURE_:
raise CosmoSevereError(self.sp.error_message)
return sigma_cb
......@@ -990,6 +990,10 @@ cdef class Class:
self.compute(["spectra"])
return self.sp.sigma8
def sigma8_cb(self):
self.compute(["spectra"])
return self.sp.sigma8_cb
def rs_drag(self):
self.compute(["thermodynamics"])
return self.th.rs_d
......@@ -1643,6 +1647,8 @@ cdef class Class:
value = self.sp.alpha_RR_2_2500
elif name == 'sigma8':
value = self.sp.sigma8
elif name == 'sigma8_cb':
value = self.sp.sigma8
else:
raise CosmoSevereError("%s was not recognized as a derived parameter" % name)
derived[name] = value
......
......@@ -2949,7 +2949,6 @@ int spectra_pk(
double source_ic1_cb;
double source_ic2_cb;
double pk_cb_tot=0.,ln_pk_cb_tot=0.;
short compute_sigma8_cb=_FALSE_;
/** - check the presence of scalar modes */
......@@ -3219,7 +3218,7 @@ int spectra_pk(
/* compute sigma8 (mean variance today in sphere of radius 8/h Mpc */
class_call(spectra_sigma(pba,ppm,psp,compute_sigma8_cb,8./pba->h,0.,&(psp->sigma8)),
class_call(spectra_sigma(pba,ppm,psp,8./pba->h,0.,&(psp->sigma8)),
psp->error_message,
psp->error_message);
......@@ -3229,8 +3228,7 @@ int spectra_pk(
exp(psp->ln_k[psp->ln_k_size-1])/pba->h);
if(pba->has_ncdm){
compute_sigma8_cb = _TRUE_;
class_call(spectra_sigma(pba,ppm,psp,compute_sigma8_cb,8./pba->h,0.,&(psp->sigma8_cb)),
class_call(spectra_sigma_cb(pba,ppm,psp,8./pba->h,0.,&(psp->sigma8_cb)),
psp->error_message,
psp->error_message);
if (psp->spectra_verbose>0){
......@@ -3238,7 +3236,6 @@ int spectra_pk(
psp->sigma8_cb,
exp(psp->ln_k[psp->ln_k_size-1])/pba->h);
}
compute_sigma8_cb = _FALSE_;
}
/**- if interpolation of \f$ P_{NL}(k,\tau)\f$ will be needed (as a function of tau),
......@@ -3298,7 +3295,6 @@ int spectra_sigma(
struct background * pba,
struct primordial * ppm,
struct spectra * psp,
short compute_sigma8_cb,
double R,
double z,
double * sigma
......@@ -3319,10 +3315,15 @@ int spectra_sigma(
double k,W,x;
if (psp->ic_ic_size[psp->index_md_scalars]>1)
if (psp->ic_ic_size[psp->index_md_scalars]>1){
class_alloc(pk_ic,
psp->ic_ic_size[psp->index_md_scalars]*sizeof(double),
psp->error_message);
if (pba->has_ncdm)
class_alloc(pk_cb_ic,
psp->ic_ic_size[psp->index_md_scalars]*sizeof(double),
psp->error_message);
}
i=0;
index_k=i;
......@@ -3346,12 +3347,7 @@ int spectra_sigma(
psp->error_message,
psp->error_message);
array_for_sigma[i*index_num+index_k]=k;
if(compute_sigma8_cb == _TRUE_){
array_for_sigma[i*index_num+index_y]=k*k*pk_cb*W*W;
}
else{
array_for_sigma[i*index_num+index_y]=k*k*pk*W*W;
}
array_for_sigma[i*index_num+index_y]=k*k*pk*W*W;
}
class_call(array_spline(array_for_sigma,
......@@ -3390,6 +3386,102 @@ int spectra_sigma(
}
//**/
int spectra_sigma_cb(
struct background * pba,
struct primordial * ppm,
struct spectra * psp,
double R,
double z,
double * sigma_cb
) {
double pk;
double * pk_ic = NULL;
double pk_cb;
double * pk_cb_ic = NULL;
double * array_for_sigma;
int index_num;
int index_k;
int index_y;
int index_ddy;
int i;
double k,W,x;
if (psp->ic_ic_size[psp->index_md_scalars]>1){
class_alloc(pk_ic,
psp->ic_ic_size[psp->index_md_scalars]*sizeof(double),
psp->error_message);
if (pba->has_ncdm)
class_alloc(pk_cb_ic,
psp->ic_ic_size[psp->index_md_scalars]*sizeof(double),
psp->error_message);
}
i=0;
index_k=i;
i++;
index_y=i;
i++;
index_ddy=i;
i++;
index_num=i;
class_alloc(array_for_sigma,
psp->ln_k_size*index_num*sizeof(double),
psp->error_message);
for (i=0;i<psp->ln_k_size;i++) {
k=exp(psp->ln_k[i]);
if (i == (psp->ln_k_size-1)) k *= 0.9999999; // to prevent rounding error leading to k being bigger than maximum value
x=k*R;
W=3./x/x/x*(sin(x)-x*cos(x));
class_call(spectra_pk_at_k_and_z(pba,ppm,psp,k,z,&pk,pk_ic,&pk_cb,pk_cb_ic),
psp->error_message,
psp->error_message);
array_for_sigma[i*index_num+index_k]=k;
array_for_sigma[i*index_num+index_y]=k*k*pk_cb*W*W;
}
class_call(array_spline(array_for_sigma,
index_num,
psp->ln_k_size,
index_k,
index_y,
index_ddy,
_SPLINE_EST_DERIV_,
psp->error_message),
psp->error_message,
psp->error_message);
class_call(array_integrate_all_spline(array_for_sigma,
index_num,
psp->ln_k_size,
index_k,
index_y,
index_ddy,
sigma_cb,
psp->error_message),
psp->error_message,
psp->error_message);
free(array_for_sigma);
if (psp->ic_ic_size[psp->index_md_scalars]>1){
free(pk_ic);
if (pba->has_ncdm)
free(pk_cb_ic);
}
*sigma_cb = sqrt(*sigma_cb/(2.*_PI_*_PI_));
return _SUCCESS_;
}
/**
* This routine computes a table of values for all matter power spectra P(k),
* given the source functions and primordial spectra.
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment