Commit 2a867cda authored by Julien Lesgourgues's avatar Julien Lesgourgues
Browse files

added write-in-file for thermo, fixed minor bugs and improved minor stylistic details

* devel:
  added version nunber 2.1.1
  Improved testings wiht nose
  fixed memory leak that appeared only for isocurvature modes and no Pk
  Added feature: write thermodynamics file
  improved details of testing input values of w_fld
  improved style of test_loops
  fixed misplaced free in transfer_free
  Macro-ised the write background part.
  background_functions is now an internal function, never called eslewhere
parents b619af39 7db712fc
......@@ -397,35 +397,39 @@ bias = 1.
s_bias = 0.
7) file name root 'root' for all output files (default: set 'root' to 'output/') (if Cl requested, written to '<root>cl.dat'; if P(k) requested, written to '<root>pk.dat'; plus similar files for scalars, tensors, pairs of initial conditions, etc.; if file with input parameters requested, written to '<root>parameters.ini')
7a) file name root 'root' for all output files (default: set 'root' to 'output/') (if Cl requested, written to '<root>cl.dat'; if P(k) requested, written to '<root>pk.dat'; plus similar files for scalars, tensors, pairs of initial conditions, etc.; if file with input parameters requested, written to '<root>parameters.ini')
root = output/test_
8) do you want headers at the beginning of each output file (giving precisions on the output units/ format) ? If 'headers' set to something containing the letter 'y' or 'Y', headers written, otherwise not written (default: written)
7b) do you want headers at the beginning of each output file (giving precisions on the output units/ format) ? If 'headers' set to something containing the letter 'y' or 'Y', headers written, otherwise not written (default: written)
headers = yes
9) in all output files, do you want columns to be normalized and ordered with the default CLASS definitions or with the CAMB definitions (often idential to the CMBFAST one) ? Set 'format' to either 'class', 'CLASS', 'camb' or 'CAMB' (default: 'class')
7c) in all output files, do you want columns to be normalized and ordered with the default CLASS definitions or with the CAMB definitions (often idential to the CMBFAST one) ? Set 'format' to either 'class', 'CLASS', 'camb' or 'CAMB' (default: 'class')
format = class
10) Do you want to write a table of background quantitites in a file? This will include H, densities, Omegas, various cosmological distances, sound horizon, etc., as a function of conformal time, proper time, scale factor. File created if 'write background' set to something containing the letter 'y' or 'Y', file written, otherwise not written (default: not written)
7d) Do you want to write a table of background quantitites in a file? This will include H, densities, Omegas, various cosmological distances, sound horizon, etc., as a function of conformal time, proper time, scale factor. File created if 'write background' set to something containing the letter 'y' or 'Y', file written, otherwise not written (default: not written)
write background = no
11) Do you want to write a table of perturbations to files for certain wavenumbers k? Dimension of k is 1/Mpc. The actual wave numbers are chosen such that they are as close as possible to the requested k-values.
7e) Do you want to write a table of thermodynamics quantitites in a file? File is created if 'write thermodynamics' is set to something containing the letter 'y' or 'Y'. (default: not written)
write thermodynamics = no
7f) Do you want to write a table of perturbations to files for certain wavenumbers k? Dimension of k is 1/Mpc. The actual wave numbers are chosen such that they are as close as possible to the requested k-values.
k_output_values = #0.01, 0.1, 0.0001
12) Do you want to write the primordial scalar(/tensor) spectrum in a file, with columns k [1/Mpc], P_s(k) [dimensionless], ( P_t(k) [dimensionless])? File created if 'write primordial' set to something containing the letter 'y' or 'Y', file written, otherwise not written (default: not written)
7g) Do you want to write the primordial scalar(/tensor) spectrum in a file, with columns k [1/Mpc], P_s(k) [dimensionless], ( P_t(k) [dimensionless])? File created if 'write primordial' set to something containing the letter 'y' or 'Y', file written, otherwise not written (default: not written)
write primordial = no
13) Do you want to have all input/precision parameters which have been read written in file '<root>parameters.ini', and those not written in file '<root>unused_parameters' ? If 'write parameters' set to something containing the letter 'y' or 'Y', file written, otherwise not written (default: not written)
7h) Do you want to have all input/precision parameters which have been read written in file '<root>parameters.ini', and those not written in file '<root>unused_parameters' ? If 'write parameters' set to something containing the letter 'y' or 'Y', file written, otherwise not written (default: not written)
write parameters = yeap
14) Do you want a warning written in the standard output when an input parameter or value could not be interpreted ? If 'write warnings' set to something containing the letter 'y' or 'Y', warnings written, otherwise not written (default: not written)
7i) Do you want a warning written in the standard output when an input parameter or value could not be interpreted ? If 'write warnings' set to something containing the letter 'y' or 'Y', warnings written, otherwise not written (default: not written)
write warnings =
......
......@@ -320,13 +320,6 @@ extern "C" {
double * pvecback
);
int background_functions(
struct background *pba,
double a,
short return_format,
double * pvecback
);
int background_tau_of_z(
struct background *pba,
double z,
......@@ -346,6 +339,13 @@ extern "C" {
struct background *pba
);
int background_functions(
struct background *pba,
double a,
short return_format,
double * pvecback
);
int background_ncdm_distribution(
void *pba,
double q,
......
......@@ -15,7 +15,7 @@
#ifndef __COMMON__
#define __COMMON__
#define _VERSION_ "v2.1.0"
#define _VERSION_ "v2.1.1"
#define _TRUE_ 1 /**< integer associated to true statement */
#define _FALSE_ 0 /**< integer associated to false statement */
......@@ -49,9 +49,9 @@ typedef char FileName[_FILENAMESIZE_];
#define _HUGE_ 1.e99
#define _OUTPUTPRECISION_ 10 /**< Number of significant digits in some output files */
#define _OUTPUTPRECISION_ 12 /**< Number of significant digits in some output files */
#define _COLUMNWIDTH_ 20 /**< Must be at least _OUTPUTPRECISION_+8 for guaranteed fixed width columns */
#define _COLUMNWIDTH_ 24 /**< Must be at least _OUTPUTPRECISION_+8 for guaranteed fixed width columns */
#ifndef __CLASSDIR__
#define __CLASSDIR__ "." /**< The directory of CLASS. This is set to the absolute path to the CLASS directory so this is just a failsafe. */
......@@ -237,7 +237,9 @@ void* class_protect_memcpy(void* dest, void* from, size_t sz);
title, \
condition){ \
if (condition == _TRUE_) \
fprintf(file,"%*s%-*s ",MAX(0,_COLUMNWIDTH_-_OUTPUTPRECISION_-6),"",_OUTPUTPRECISION_+6,title); \
fprintf(file,"%*s%-*s ",\
MAX(0,MIN(_COLUMNWIDTH_-_OUTPUTPRECISION_-6,_COLUMNWIDTH_-((int) strlen(title)))), \
"",_OUTPUTPRECISION_+6,title); \
}
/** parameters related to the precision of the code and to the method of calculation */
......
......@@ -52,6 +52,7 @@ struct output {
enum file_format output_format;
short write_background;
short write_thermodynamics;
short write_primordial;
//@}
......@@ -130,6 +131,12 @@ extern "C" {
struct output * pop
);
int output_thermodynamics(
struct background * pba,
struct thermo * pth,
struct output * pop
);
int output_primordial(
struct perturbs * ppt,
struct primordial * ppm,
......@@ -213,6 +220,21 @@ extern "C" {
double * pvecback
);
int output_open_thermodynamics_file(
struct thermo * pth,
struct output * pop,
FILE ** thermofile,
FileName filename
);
int output_one_line_of_thermodynamics(
struct thermo * pth,
FILE * thermofile,
double tau,
double z,
double * pvecthermo
);
int output_open_primordial_file(
struct perturbs * ppt,
struct primordial * ppm,
......
......@@ -39,6 +39,7 @@ struct spectra {
//@{
int md_size; /**< number of modes (scalar, tensor, ...) included in computation */
int index_md_scalars; /**< index for scalar modes */
int * ic_size; /**< for a given mode, ic_size[index_md] = number of initial conditions included in computation */
int * ic_ic_size; /**< for a given mode, ic_ic_size[index_md] = number of pairs of (index_ic1, index_ic2) with index_ic2 >= index_ic1; this number is just N(N+1)/2 where N = ic_size[index_md] */
......@@ -142,8 +143,6 @@ struct spectra {
//@{
int index_md_scalars; /**< index for scalar modes (the matter power spectrum refers by construction to scalar modes) */
int ln_k_size; /**< number ln(k) values */
double * ln_k; /**< list of ln(k) values ln_k[index_k] */
......
from classy import Class
from classy import CosmoSevereError
import itertools
import sys
import numpy as np
import unittest
from nose_parameterized import parameterized
......@@ -42,19 +43,15 @@ class TestClass(unittest.TestCase):
self.cosmo.cleanup()
del self.scenario
#'Mnu',
#'Positive_Omega_k',
#'Negative_Omega_k'),
#@parameterized.expand([
@parameterized.expand(
itertools.product(
('LCDM',
'Mnu',
'Positive_Omega_k',
'Negative_Omega_k',
'Isocurvature_modes'),
({}, {'output': 'mPk'}, {'output': 'tCl'},
{'output': 'tCl pCl lCl'}, {'output': 'mPk tCl lCl','P_k_max_h/Mpc':10},
'Isocurvature_modes', ),
({'output': ''}, {'output': 'mPk'}, {'output': 'tCl'},
{'output': 'tCl pCl lCl'}, {'output': 'mPk tCl lCl', 'P_k_max_h/Mpc':10},
{'output': 'nCl sCl'}, {'output': 'tCl pCl lCl nCl sCl'}),
({'gauge': 'newtonian'}, {'gauge': 'sync'}),
({}, {'non linear': 'halofit'})))
......@@ -67,19 +64,19 @@ class TestClass(unittest.TestCase):
elif name == 'Negative_Omega_k':
self.scenario.update({'Omega_k': -0.01})
elif name == 'Isocurvature_modes':
self.scenario.update({'ic':'ad,cdi,nid','c_ad_cdi':-0.5})
self.scenario.update({'ic': 'ad,nid,cdi', 'c_ad_cdi': -0.5})
self.scenario.update(scenario)
if scenario != {}:
self.scenario.update(gauge)
self.scenario.update(nonlinear)
print '\n\n--------------------------'
print '| Test case %s |' % name
print '--------------------------'
sys.stderr.write('\n\n---------------------------------\n')
sys.stderr.write('| Test case %s |\n' % name)
sys.stderr.write('---------------------------------\n')
for key, value in self.scenario.iteritems():
print key, '=', value
print
sys.stderr.write("%s = %s\n" % (key, value))
sys.stderr.write("\n")
setting = self.cosmo.set(
dict(self.verbose.items()+self.scenario.items()))
......@@ -87,7 +84,22 @@ class TestClass(unittest.TestCase):
cl_list = ['tCl', 'lCl', 'pCl', 'nCl', 'sCl']
self.cosmo.compute()
# Depending on the cases, the compute should fail or not
should_fail = True
output = self.scenario['output'].split()
for elem in output:
if elem in ['tCl', 'pCl']:
for elem2 in output:
if elem2 == 'lCl':
should_fail = False
break
if not should_fail:
self.cosmo.compute()
else:
self.assertRaises(CosmoSevereError, self.cosmo.compute)
return
self.assertTrue(
self.cosmo.state,
"Class failed to go through all __init__ methods")
......
......@@ -54,10 +54,9 @@
* any other qunatity, we can call background_at_tau(), which
* interpolates in the table and returns all values. Finally it can be
* useful to get 'tau' for a given redshift 'z': this can be done with
* background_tau_of_z(). So, in principle, if we know z, calling
* background_tau_of_z() and then background_at_tau() will provide
* roughly the same information as calling directly
* background_functions() with a=a_0/(1+z).
* background_tau_of_z(). So if we are somewhere in the code, knowing
* z and willing to get background quantitites, we should call first
* background_tau_of_z() and then background_at_tau().
*
*
* In order to save time, background_at_tau() can be called in three
......@@ -75,7 +74,7 @@
* In summary, the following functions can be called from other modules:
*
* -# background_init() at the beginning
* -# background_at_tau(), background_functions(), background_tau_of_z() at any later time
* -# background_at_tau(), background_tau_of_z() at any later time
* -# background_free() at the end, when no more calls to the previous functions are needed
*/
......@@ -225,181 +224,6 @@ int background_tau_of_z(
return _SUCCESS_;
}
/**
* Background quantities at given a.
*
* Function evaluating all background quantities which can be computed
* analytically as a function of parameters like the scale factor 'a'
* (see discussion at the beginnign of this file). In extended
* comsological models, we would include other input parameters than
* just 'a', e.g. (phi, phidot) for quintessence, some temperature of
* exotic relics, etc...
*
* @param pba Input: pointer to background structure
* @param a Input: value of scale factor
* @param return_format Input: format of output vector
* @param pvecback Output: vector of background quantities (assmued to be already allocated)
* @return the error status
*/
int background_functions(
struct background *pba,
double a, /* in extended models there could be more than one argument: phi, phidot of quintessence; temperature of some particles; etc. */
short return_format,
double * pvecback /* vector with argument pvecback[index_bg] (must be already allocated with a size comptible with return_format) */
) {
/** Summary: */
/** - define local variables */
/* total density */
double rho_tot;
/* total pressure */
double p_tot;
/* total relativistic density */
double rho_r;
/* total non-relativistic density */
double rho_m;
/* scale factor relative to scale factor today */
double a_rel;
/* background ncdm quantities */
double rho_ncdm,p_ncdm,pseudo_p_ncdm;
/* index for n_ncdm species */
int n_ncdm;
/** - initialize local variables */
rho_tot = 0.;
p_tot = 0.;
rho_r=0.;
rho_m=0.;
a_rel = a / pba->a_today;
class_test(a_rel <= 0.,
pba->error_message,
"a = %e instead of strictly positive",a_rel);
/** - pass value of a to output */
pvecback[pba->index_bg_a] = a;
/** - compute each component's density and pressure */
/* photons */
pvecback[pba->index_bg_rho_g] = pba->Omega0_g * pow(pba->H0,2) / pow(a_rel,4);
rho_tot += pvecback[pba->index_bg_rho_g];
p_tot += (1./3.) * pvecback[pba->index_bg_rho_g];
rho_r += pvecback[pba->index_bg_rho_g];
/* baryons */
pvecback[pba->index_bg_rho_b] = pba->Omega0_b * pow(pba->H0,2) / pow(a_rel,3);
rho_tot += pvecback[pba->index_bg_rho_b];
p_tot += 0;
rho_m += pvecback[pba->index_bg_rho_b];
/* cdm */
if (pba->has_cdm == _TRUE_) {
pvecback[pba->index_bg_rho_cdm] = pba->Omega0_cdm * pow(pba->H0,2) / pow(a_rel,3);
rho_tot += pvecback[pba->index_bg_rho_cdm];
p_tot += 0.;
rho_m += pvecback[pba->index_bg_rho_cdm];
}
/* ncdm */
if (pba->has_ncdm == _TRUE_) {
/* Loop over species: */
for(n_ncdm=0; n_ncdm<pba->N_ncdm; n_ncdm++){
/* function returning background ncdm[n_ncdm] quantities (only
those for which non-NULL pointers are passed) */
class_call(background_ncdm_momenta(
pba->q_ncdm_bg[n_ncdm],
pba->w_ncdm_bg[n_ncdm],
pba->q_size_ncdm_bg[n_ncdm],
pba->M_ncdm[n_ncdm],
pba->factor_ncdm[n_ncdm],
1./a_rel-1.,
NULL,
&rho_ncdm,
&p_ncdm,
NULL,
&pseudo_p_ncdm),
pba->error_message,
pba->error_message);
pvecback[pba->index_bg_rho_ncdm1+n_ncdm] = rho_ncdm;
rho_tot += rho_ncdm;
pvecback[pba->index_bg_p_ncdm1+n_ncdm] = p_ncdm;
p_tot += p_ncdm;
pvecback[pba->index_bg_pseudo_p_ncdm1+n_ncdm] = pseudo_p_ncdm;
/* (3 p_ncdm1) is the "relativistic" contrinution to rho_ncdm1 */
rho_r += 3.* p_ncdm;
/* (rho_ncdm1 - 3 p_ncdm1) is the "non-relativistic" contribution
to rho_ncdm1 */
rho_m += rho_ncdm - 3.* p_ncdm;
}
}
/* Lambda */
if (pba->has_lambda == _TRUE_) {
pvecback[pba->index_bg_rho_lambda] = pba->Omega0_lambda * pow(pba->H0,2);
rho_tot += pvecback[pba->index_bg_rho_lambda];
p_tot -= pvecback[pba->index_bg_rho_lambda];
}
/* fluid with w=w0+wa(1-a/a0) and constant cs2 */
if (pba->has_fld == _TRUE_) {
pvecback[pba->index_bg_rho_fld] = pba->Omega0_fld * pow(pba->H0,2)
/ pow(a_rel,3.*(1.+pba->w0_fld+pba->wa_fld))
* exp(3.*pba->wa_fld*(a_rel-1.));
rho_tot += pvecback[pba->index_bg_rho_fld];
p_tot += (pba->w0_fld+pba->wa_fld*(1.-a_rel)) * pvecback[pba->index_bg_rho_fld];
}
/* relativistic neutrinos (and all relativistic relics) */
if (pba->has_ur == _TRUE_) {
pvecback[pba->index_bg_rho_ur] = pba->Omega0_ur * pow(pba->H0,2) / pow(a_rel,4);
rho_tot += pvecback[pba->index_bg_rho_ur];
p_tot += (1./3.) * pvecback[pba->index_bg_rho_ur];
rho_r += pvecback[pba->index_bg_rho_ur];
}
/** - compute expansion rate H from Friedmann equation: this is the
unique place where the Friedmann equation is assumed. Remember
that densities are all expressed in units of [3c^2/8piG], ie
rho_class = [8 pi G rho_physical / 3 c^2] */
pvecback[pba->index_bg_H] = sqrt(rho_tot-pba->K/a/a);
/** - compute derivative of H with respect to conformal time */
pvecback[pba->index_bg_H_prime] = - (3./2.) * (rho_tot + p_tot) * a + pba->K/a;
/** - compute relativistic density to total density ratio */
pvecback[pba->index_bg_Omega_r] = rho_r / rho_tot;
/** - compute other quantities in the exhaustive, redundent format: */
if (return_format == pba->long_info) {
/** - compute critical density */
pvecback[pba->index_bg_rho_crit] = rho_tot-pba->K/a/a;
class_test(pvecback[pba->index_bg_rho_crit] <= 0.,
pba->error_message,
"rho_crit = %e instead of strictly positive",pvecback[pba->index_bg_rho_crit]);
/** - compute Omega_m */
pvecback[pba->index_bg_Omega_m] = rho_m / rho_tot;
/* one can put other variables here */
/* */
/* */
}
return _SUCCESS_;
}
/**
* Initialize the background structure, and in particular the
* background interpolation table.
......@@ -491,6 +315,13 @@ int background_init(
pba->error_message,
"T_cmb=%g out of bounds (%g<T_cmb<%g)",pba->T_cmb,_TCMB_SMALL_,_TCMB_BIG_);
if (pba->has_fld == _TRUE_) {
class_test(pba->w0_fld+pba->wa_fld>=1./3.,
pba->error_message,
"Your choice for w0_fld+wa_fld=%g is suspicious, there would not be radiation domination at early times\n",
pba->w0_fld+pba->wa_fld);
}
/* in verbose mode, inform the user about the value of the ncdm
masses in eV and about the ratio [m/omega_ncdm] in eV (the usual
93 point something)*/
......@@ -763,6 +594,181 @@ int background_indices(
}
/**
* Background quantities at given a.
*
* Function evaluating all background quantities which can be computed
* analytically as a function of parameters like the scale factor 'a'
* (see discussion at the beginnign of this file). In extended
* comsological models, we would include other input parameters than
* just 'a', e.g. (phi, phidot) for quintessence, some temperature of
* exotic relics, etc...
*
* @param pba Input: pointer to background structure
* @param a Input: value of scale factor
* @param return_format Input: format of output vector
* @param pvecback Output: vector of background quantities (assmued to be already allocated)
* @return the error status
*/
int background_functions(
struct background *pba,
double a, /* in extended models there could be more than one argument: phi, phidot of quintessence; temperature of some particles; etc. */
short return_format,
double * pvecback /* vector with argument pvecback[index_bg] (must be already allocated with a size comptible with return_format) */
) {
/** Summary: */
/** - define local variables */
/* total density */
double rho_tot;
/* total pressure */
double p_tot;
/* total relativistic density */
double rho_r;
/* total non-relativistic density */
double rho_m;
/* scale factor relative to scale factor today */
double a_rel;
/* background ncdm quantities */
double rho_ncdm,p_ncdm,pseudo_p_ncdm;
/* index for n_ncdm species */
int n_ncdm;
/** - initialize local variables */
rho_tot = 0.;
p_tot = 0.;
rho_r=0.;
rho_m=0.;
a_rel = a / pba->a_today;
class_test(a_rel <= 0.,
pba->error_message,
"a = %e instead of strictly positive",a_rel);
/** - pass value of a to output */
pvecback[pba->index_bg_a] = a;
/** - compute each component's density and pressure */
/* photons */
pvecback[pba->index_bg_rho_g] = pba->Omega0_g * pow(pba->H0,2) / pow(a_rel,4);
rho_tot += pvecback[pba->index_bg_rho_g];
p_tot += (1./3.) * pvecback[pba->index_bg_rho_g];
rho_r += pvecback[pba->index_bg_rho_g];
/* baryons */
pvecback[pba->index_bg_rho_b] = pba->Omega0_b * pow(pba->H0,2) / pow(a_rel,3);
rho_tot += pvecback[pba->index_bg_rho_b];
p_tot += 0;
rho_m += pvecback[pba->index_bg_rho_b];
/* cdm */
if (pba->has_cdm == _TRUE_) {
pvecback[pba->index_bg_rho_cdm] = pba->Omega0_cdm * pow(pba->H0,2) / pow(a_rel,3);
rho_tot += pvecback[pba->index_bg_rho_cdm];
p_tot += 0.;
rho_m += pvecback[pba->index_bg_rho_cdm];
}
/* ncdm */
if (pba->has_ncdm == _TRUE_) {
/* Loop over species: */
for(n_ncdm=0; n_ncdm<pba->N_ncdm; n_ncdm++){
/* function returning background ncdm[n_ncdm] quantities (only
those for which non-NULL pointers are passed) */
class_call(background_ncdm_momenta(
pba->q_ncdm_bg[n_ncdm],
pba->w_ncdm_bg[n_ncdm],
pba->q_size_ncdm_bg[n_ncdm],
pba->M_ncdm[n_ncdm],
pba->factor_ncdm[n_ncdm],
1./a_rel-1.,
NULL,
&rho_ncdm,
&p_ncdm,
NULL,
&pseudo_p_ncdm),
pba->error_message,
pba->error_message);
pvecback[pba->index_bg_rho_ncdm1+n_ncdm] = rho_ncdm;
rho_tot += rho_ncdm;
pvecback[pba->index_bg_p_ncdm1+n_ncdm] = p_ncdm;
p_tot += p_ncdm;
pvecback[pba->index_bg_pseudo_p_ncdm1+n_ncdm] = pseudo_p_ncdm;
/* (3 p_ncdm1) is the "relativistic" contrinution to rho_ncdm1 */
rho_r += 3.* p_ncdm;
/* (rho_ncdm1 - 3 p_ncdm1) is the "non-relativistic" contribution
to rho_ncdm1 */
rho_m += rho_ncdm - 3.* p_ncdm;
}
}
/* Lambda */
if (pba->has_lambda == _TRUE_) {
pvecback[pba->index_bg_rho_lambda] = pba->Omega0_lambda * pow(pba->H0,2);
rho_tot += pvecback[pba->index_bg_rho_lambda];
p_tot -= pvecback[pba->index_bg_rho_lambda];
}
/* fluid with w=w0+wa(1-a/a0) and constant cs2 */
if (pba->has_fld == _TRUE_) {
pvecback[pba->index_bg_rho_fld] = pba->Omega0_fld * pow(pba->H0,2)
/ pow(a_rel,3.*(1.+pba->w0_fld+pba->wa_fld))
* exp(3.*pba->wa_fld*(a_rel-1.));
rho_tot += pvecback[pba->index_bg_rho_fld];
p_tot += (pba->w0_fld+pba->wa_fld*(1.-a_rel)) * pvecback[pba->index_bg_rho_fld];
}
/* relativistic neutrinos (and all relativistic relics) */
if (pba->has_ur == _TRUE_) {
pvecback[pba->index_bg_rho_ur] = pba->Omega0_ur * pow(pba->H0,2) / pow(a_rel,4);
rho_tot += pvecback[pba->index_bg_rho_ur];
p_tot += (1./3.) * pvecback[pba->index_bg_rho_ur];
rho_r += pvecback[pba->index_bg_rho_ur];
}