Commit 20aed25e authored by Julien Lesgourgues's avatar Julien Lesgourgues
Browse files

Little bug fixes and improvement from Tokyo workshop

* 2.4:
  version number updated to 2.4.2
  Moved inflation tests from primordial module to input module and updated test_class.py accordingly.
  minor stylistic modifs of previous commit
  k_output_values are now added to actual k-list
  Fixed hardcoding in fix of issue #19
  Fixed a wrong allocation which would lead to segfaults on some 32 bit machines (if sizeof(double*)<sizeof(double) )
  added growth factors D, in write background
  Fixed Github issue #19 which was a bug for perturbations output
  removed deprecated classy function cleanup
  Fixed bug where the C-compiler defined in the Makefile would not be used by Cython when compiling the wrapper. Also avoids using --user when prefix have been set in distutils.cfg
  Fixed CPU.py, created script now works
  fixed compatibility with CosmoHammer
  changed the naming of output files from CPU.py
  modified default behaviour: if no field is asked, plot them all
parents befd7ba2 56d364ca
......@@ -30,20 +30,16 @@ cosmo_mini_toolbox, available under GPLv3 at
https://github.com/JesusTorrado/cosmo_mini_toolbox
"""
from __future__ import unicode_literals
# System imports
import os
import sys
import argparse
import itertools
from math import floor
from matplotlib import scale as mscale
from matplotlib.transforms import Transform
from matplotlib.ticker import FixedLocator
# Numerics
import numpy as np
from numpy import ma
from scipy.interpolate import splrep, splev, InterpolatedUnivariateSpline
from scipy.interpolate import InterpolatedUnivariateSpline
from math import floor
# Plotting
......@@ -122,8 +118,7 @@ def plot_CLASS_output(files, x_axis, y_axis, ratio=False, printing='',
"""
# Define the python script name, and the pdf path
python_script_path = files[0]+'.py'
pdf_path = files[0]+'.pdf'
python_script_path = os.path.splitext(files[0])[0]+'.py'
# The variable text will contain all the lines to be printed in the end to
# the python script path, joined with newline characters. Beware of the
......@@ -146,7 +141,7 @@ def plot_CLASS_output(files, x_axis, y_axis, ratio=False, printing='',
'for data_file in files:',
' data.append(np.loadtxt(data_file))']
# Recover the base name of the files, everything before the .
# Recover the base name of the files, everything before the dot
roots = [elem.split(os.path.sep)[-1].split('.')[0] for elem in files]
text += ['roots = [%s]' % ', '.join(["'%s'" % root for root in roots])]
......@@ -163,6 +158,7 @@ def plot_CLASS_output(files, x_axis, y_axis, ratio=False, printing='',
# title.
num_columns, names, tex_names = extract_headers(files[index])
text += ['', 'index, curve = %i, data[%i]' % (index, index)]
# Check if everything is in order
if num_columns == 2:
y_axis = [names[1]]
......@@ -176,14 +172,13 @@ def plot_CLASS_output(files, x_axis, y_axis, ratio=False, printing='',
# Store the selected text and tex_names to the script
selected = []
for elem in y_axis:
selected.extend([name for name in names if name.find(elem) != -1 and
name not in selected])
selected.extend(
[name for name in names if name.find(elem) != -1 and
name not in selected])
if not y_axis:
selected = names[1:]
y_axis = selected
text += ['y_axis = %s' % selected]
text += ['tex_names = %s' % [elem for (elem, name) in
zip(tex_names, names) if name in selected]]
# Decide for the x_axis, by default the index will be set to zero
x_index = 0
if x_axis:
......@@ -191,15 +186,17 @@ def plot_CLASS_output(files, x_axis, y_axis, ratio=False, printing='',
if name.find(x_axis) != -1:
x_index = index_name
break
# Store to text
text += ['y_axis = %s' % selected]
text += ['tex_names = %s' % [elem for (elem, name) in
zip(tex_names, names) if name in selected]]
text += ["x_axis = '%s'" % tex_names[x_index]]
# Store the limits variable in the text
text += ["ylim = %s" % ylim]
text += ["xlim = %s" % xlim]
for selec in y_axis:
index_selec = names.index(selec)
plot_line = ' ax.'
plot_line = 'ax.'
if scale == 'lin':
plot_line += 'plot(curve[:, %i], curve[:, %i])' % (
x_index, index_selec)
......@@ -456,15 +453,12 @@ def main():
# spectrum
if not args.y_axis:
if args.files[0].rfind('cl') != -1:
y_axis = ['TT']
scale = 'loglog'
elif args.files[0].rfind('pk') != -1:
y_axis = ['P']
scale = 'loglog'
else:
raise TypeError(
"Please specify a field to plot")
args.y_axis = y_axis
scale = 'lin'
args.y_axis = []
else:
scale = ''
if not args.scale:
......
......@@ -17,8 +17,8 @@ vpath .base build
########################################################
# your C compiler:
CC = gcc -Wall #-g -pg -E#-ggdb
#CC = icc #-ggdb
CC = gcc
#CC = icc
#CC = pgcc
# your tool for creating static libraries:
......@@ -58,6 +58,16 @@ INCLUDES = -I../include
# automatically add external programs if needed. First, initialize to blank.
EXTERNAL =
# Try to automatically avoid an error 'error: can't combine user with ...'
# which sometimes happens with brewed Python on OSX:
CFGFILE=$(shell $(PYTHON) -c "import sys; print sys.prefix+'/lib/'+'python'+'.'.join(['%i' % e for e in sys.version_info[0:2]])+'/distutils/distutils.cfg'")
PYTHONPREFIX=$(shell grep -s "prefix" $(CFGFILE))
ifeq ($(PYTHONPREFIX),)
PYTHONFLAGS=--user
else
PYTHONFLAGS=
endif
# eventually update flags for including HyRec
ifneq ($(HYREC),)
vpath %.c $(HYREC)
......@@ -175,7 +185,7 @@ tar: $(C_ALL) $(C_TEST) $(H_ALL) $(PRE_ALL) $(INI_ALL) $(MISC_FILES) $(HYREC) $(
tar czvf class.tar.gz $(C_ALL) $(H_ALL) $(PRE_ALL) $(INI_ALL) $(MISC_FILES) $(HYREC) $(PYTHON_FILES)
classy: libclass.a python/classy.pyx python/cclassy.pxd
cd python; $(PYTHON) setup.py install --user
cd python; export CC=$(CC); $(PYTHON) setup.py install $(PYTHONFLAGS)
clean: .base
rm -rf $(WRKDIR);
......
......@@ -15,7 +15,7 @@
#ifndef __COMMON__
#define __COMMON__
#define _VERSION_ "v2.4.1"
#define _VERSION_ "v2.4.2"
#define _TRUE_ 1 /**< integer associated to true statement */
#define _FALSE_ 0 /**< integer associated to false statement */
......
......@@ -267,6 +267,11 @@ extern "C" {
int * aux_flag,
ErrorMsg error_message);
int compare_integers (const void * elem1, const void * elem2);
int compare_doubles(const void *a,const void *b);
#ifdef __cplusplus
}
#endif
......
......@@ -162,7 +162,7 @@ struct perturbs
int store_perturbations; /**< Do we want to store perturbations? */
int k_output_values_num; /**< Number of perturbation outputs (default=0) */
double k_output_values[_MAX_NUMBER_OF_K_FILES_]; /**< List of k values where perturbation output is requested. */
int index_k_output_values[_MAX_NUMBER_OF_K_FILES_]; /**< List of indices corresponding to k-values close to k_output_values */
int *index_k_output_values; /**< List of indices corresponding to k-values close to k_output_values for each mode. [index_md*k_output_values_num+ik]*/
char scalar_titles[_MAXTITLESTRINGLENGTH_]; /**< _DELIMITER_ separated string of titles for scalar perturbation output files. */
char vector_titles[_MAXTITLESTRINGLENGTH_]; /**< _DELIMITER_ separated string of titles for vector perturbation output files. */
char tensor_titles[_MAXTITLESTRINGLENGTH_]; /**< _DELIMITER_ separated string of titles for tensor perturbation output files. */
......
......@@ -133,14 +133,7 @@ cdef class Class:
def empty(self):
self._pars = {}
self.ready=False
def cleanup(self):
if self.ready==False:
return True
for i in range(len(self._pars)):
if self.fc.read[i]==0:
del(self._pars[self.fc.name[i]])
self.ready = False
# Create an equivalent of the parameter file. Non specified values will be
# taken at their default (in Class)
......@@ -1401,7 +1394,8 @@ cdef class Class:
# Compute the derived paramter value and store them
params = ctx.getData()
self.get_current_derived_parameters(data)
self.get_current_derived_parameters(
data.get_mcmc_parameters(['derived']))
for elem in data.get_mcmc_parameters(['derived']):
data.mcmc_parameters[elem]['current'] /= \
data.mcmc_parameters[elem]['scale']
......
......@@ -91,15 +91,15 @@ COMPARE_OUTPUT = True
# against. Indeed, when not specifying a field, CLASS takes the default input.
CLASS_INPUT = {}
CLASS_INPUT['Mnu'] = (
[{'N_eff': 0.0, 'N_ncdm': 1, 'm_ncdm': 0.06, 'deg_ncdm': 3.0},
{'N_eff': 1.5, 'N_ncdm': 1, 'm_ncdm': 0.03, 'deg_ncdm': 1.5}],
'normal')
#CLASS_INPUT['Mnu'] = (
# [{'N_eff': 0.0, 'N_ncdm': 1, 'm_ncdm': 0.06, 'deg_ncdm': 3.0},
# {'N_eff': 1.5, 'N_ncdm': 1, 'm_ncdm': 0.03, 'deg_ncdm': 1.5}],
# 'normal')
CLASS_INPUT['Curvature'] = (
[{'Omega_k': 0.01},
{'Omega_k': -0.01}],
'normal')
#CLASS_INPUT['Curvature'] = (
# [{'Omega_k': 0.01},
# {'Omega_k': -0.01}],
# 'normal')
CLASS_INPUT['Isocurvature_modes'] = (
[{'ic': 'ad,nid,cdi', 'c_ad_cdi': -0.5}],
......@@ -375,6 +375,30 @@ class TestClass(unittest.TestCase):
if 'modes' in self.scenario.keys() and self.scenario['modes'].find('s') == -1:
should_fail = True
# If we specify initial conditions (for scalar modes), we must have
# perturbations and scalar modes.
if 'ic' in self.scenario.keys():
if 'modes' in self.scenario.keys() and self.scenario['modes'].find('s') == -1:
should_fail = True
if 'output' not in self.scenario.keys():
should_fail = True
# If we use inflation module, we must have scalar modes,
# tensor modes, no vector modes and we should only have adiabatic IC:
if 'P_k_ini type' in self.scenario.keys() and self.scenario['P_k_ini type'].find('inflation') != -1:
if 'modes' not in self.scenario.keys():
should_fail = True
else:
if self.scenario['modes'].find('s') == -1:
should_fail = True
if self.scenario['modes'].find('v') != -1:
should_fail = True
if self.scenario['modes'].find('t') == -1:
should_fail = True
if 'ic' in self.scenario.keys() and self.scenario['ic'].find('i') != -1:
should_fail = True
return should_fail
def compare_output(self, reference, candidate):
......
......@@ -1904,6 +1904,9 @@ int background_output_titles(struct background * pba,
class_store_columntitle(titles,"V'_scf",pba->has_scf);
class_store_columntitle(titles,"V''_scf",pba->has_scf);
class_store_columntitle(titles,"gr.fac. D",_TRUE_);
class_store_columntitle(titles,"gr.fac. f",_TRUE_);
return _SUCCESS_;
}
......@@ -1949,6 +1952,9 @@ int background_output_data(
class_store_double(dataptr,pvecback[pba->index_bg_V_scf],pba->has_scf,storeidx);
class_store_double(dataptr,pvecback[pba->index_bg_dV_scf],pba->has_scf,storeidx);
class_store_double(dataptr,pvecback[pba->index_bg_ddV_scf],pba->has_scf,storeidx);
class_store_double(dataptr,pvecback[pba->index_bg_D],_TRUE_,storeidx);
class_store_double(dataptr,pvecback[pba->index_bg_f],_TRUE_,storeidx);
}
return _SUCCESS_;
......
......@@ -1905,6 +1905,26 @@ int input_read_parameters(
class_read_double("custom10",ppm->custom10);
}
/** Tests moved from primordial module: */
if ((ppm->primordial_spec_type == inflation_V) || (ppm->primordial_spec_type == inflation_H) || (ppm->primordial_spec_type == inflation_V_end)) {
class_test(ppt->has_scalars == _FALSE_,
errmsg,
"inflationary module cannot work if you do not ask for scalar modes");
class_test(ppt->has_vectors == _TRUE_,
errmsg,
"inflationary module cannot work if you ask for vector modes");
class_test(ppt->has_tensors == _FALSE_,
errmsg,
"inflationary module cannot work if you do not ask for tensor modes");
class_test(ppt->has_bi == _TRUE_ || ppt->has_cdi == _TRUE_ || ppt->has_nid == _TRUE_ || ppt->has_niv == _TRUE_,
errmsg,
"inflationary module cannot work if you ask for isocurvature modes");
}
/** (e) parameters for final spectra */
if (ppt->has_cls == _TRUE_) {
......@@ -2510,10 +2530,15 @@ int input_read_parameters(
"you want to write some output for %d different values of k, hence you should increase _MAX_NUMBER_OF_K_FILES_ in include/perturbations.h to at least this number",
int1);
ppt->k_output_values_num = int1;
for (i=0; i<int1; i++) {
ppt->k_output_values[i] = pointer1[i];
}
free(pointer1);
/** Sort the k_array using qsort */
qsort (ppt->k_output_values, ppt->k_output_values_num, sizeof(double), compare_doubles);
ppt->store_perturbations = _TRUE_;
pop->write_perturbations = _TRUE_;
}
......@@ -2709,6 +2734,7 @@ int input_default_params(
ppt->vector_perturbations_data[filenum] = NULL;
ppt->tensor_perturbations_data[filenum] = NULL;
}
ppt->index_k_output_values=NULL;
/** - primordial structure */
......@@ -3588,3 +3614,21 @@ int input_auxillary_target_conditions(struct file_content * pfc,
}
return _SUCCESS_;
}
int compare_integers (const void * elem1, const void * elem2) {
int f = *((int*)elem1);
int s = *((int*)elem2);
if (f > s) return 1;
if (f < s) return -1;
return 0;
}
int compare_doubles(const void *a,const void *b) {
double *x = (double *) a;
double *y = (double *) b;
if (*x < *y)
return -1;
else if
(*x > *y) return 1;
return 0;
}
......@@ -1264,8 +1264,8 @@ int output_perturbations(
for (index_ikout=0; index_ikout<ppt->k_output_values_num; index_ikout++){
if (ppt->has_scalars == _TRUE_){
index_md = 0;
k = ppt->k[index_md][ppt->index_k_output_values[index_ikout]];
index_md = ppt->index_md_scalars;
k = ppt->k[index_md][ppt->index_k_output_values[index_md*ppt->k_output_values_num+index_ikout]];
sprintf(file_name,"%s%s%d%s",pop->root,"perturbations_k",index_ikout,"_s.dat");
class_open(out, file_name, "w", ppt->error_message);
fprintf(out,"#scalar perturbations for mode k = %.*e Mpc^(-1)\n",_OUTPUTPRECISION_,k);
......@@ -1277,8 +1277,8 @@ int output_perturbations(
fclose(out);
}
if (ppt->has_vectors == _TRUE_){
index_md = 1;
k = ppt->k[index_md][ppt->index_k_output_values[index_ikout]];
index_md = ppt->index_md_vectors;
k = ppt->k[index_md][ppt->index_k_output_values[index_md*ppt->k_output_values_num+index_ikout]];
sprintf(file_name,"%s%s%d%s",pop->root,"perturbations_k",index_ikout,"_v.dat");
class_open(out, file_name, "w", ppt->error_message);
fprintf(out,"#vector perturbations for mode k = %.*e Mpc^(-1)\n",_OUTPUTPRECISION_,k);
......@@ -1290,8 +1290,8 @@ int output_perturbations(
fclose(out);
}
if (ppt->has_tensors == _TRUE_){
index_md = 2;
k = ppt->k[index_md][ppt->index_k_output_values[index_ikout]];
index_md = ppt->index_md_tensors;
k = ppt->k[index_md][ppt->index_k_output_values[index_md*ppt->k_output_values_num+index_ikout]];
sprintf(file_name,"%s%s%d%s",pop->root,"perturbations_k",index_ikout,"_t.dat");
class_open(out, file_name, "w", ppt->error_message);
fprintf(out,"#tensor perturbations for mode k = %.*e Mpc^(-1)\n",_OUTPUTPRECISION_,k);
......@@ -1327,7 +1327,7 @@ int output_primordial(
pop->error_message);
number_of_titles = get_number_of_titles(titles);
size_data = number_of_titles*ppm->lnk_size;
class_alloc(data,sizeof(double*)*size_data,pop->error_message);
class_alloc(data,sizeof(double)*size_data,pop->error_message);
class_call(primordial_output_data(ppt,
ppm,
number_of_titles,
......
......@@ -468,6 +468,9 @@ int perturb_free(
/** Stuff related to perturbations output: */
/** Free non-NULL pointers: */
if (ppt->index_k_output_values != NULL)
free(ppt->index_k_output_values);
for (filenum = 0; filenum<_MAX_NUMBER_OF_K_FILES_; filenum++){
if (ppt->scalar_perturbations_data[filenum] != NULL)
free(ppt->scalar_perturbations_data[filenum]);
......@@ -1180,13 +1183,14 @@ int perturb_get_k_list(
struct thermo * pth,
struct perturbs * ppt
) {
int index_k, index_k_output;
int index_k, index_k_output, index_mode;
double k,k_min=0.,k_rec,step,tau1;
double k_max_cmb=0.;
double k_max_cl=0.;
double * k_max_cmb;
double * k_max_cl;
double k_max=0.;
double scale2;
double k_target;
double *tmp_k_list;
int newk_size, index_newk, add_k_output_value;
class_test(ppr->k_step_transition == 0.,
ppt->error_message,
......@@ -1211,6 +1215,15 @@ int perturb_get_k_list(
ppt->md_size*sizeof(double*),
ppt->error_message);
class_calloc(k_max_cmb,
ppt->md_size,
sizeof(double),
ppt->error_message);
class_calloc(k_max_cl,
ppt->md_size,
sizeof(double),
ppt->error_message);
/** - scalar modes */
if (ppt->has_scalars == _TRUE_) {
......@@ -1232,35 +1245,35 @@ int perturb_get_k_list(
k_min = sqrt((8.-1.e-4)*pba->K);
}
/** - find k_max (as well as k_max_cmb, k_max_cl) */
/** - find k_max (as well as k_max_cmb[ppt->index_md_scalars], k_max_cl[ppt->index_md_scalars]) */
k_rec = 2. * _PI_ / pth->rs_rec; /* comoving scale corresponding to sound horizon at recombination */
k_max_cmb = k_min;
k_max_cl = k_min;
k_max_cmb[ppt->index_md_scalars] = k_min;
k_max_cl[ppt->index_md_scalars] = k_min;
k_max = k_min;
if (ppt->has_cls == _TRUE_) {
/* find k_max_cmb: */
/* find k_max_cmb[ppt->index_md_scalars] : */
/* choose a k_max_cmb corresponding to a wavelength on the last
/* choose a k_max_cmb[ppt->index_md_scalars] corresponding to a wavelength on the last
scattering surface seen today under an angle smaller than
pi/lmax: this is equivalent to
k_max_cl*[comvoving.ang.diameter.distance] > l_max */
k_max_cl[ppt->index_md_scalars]*[comvoving.ang.diameter.distance] > l_max */
k_max_cmb = ppr->k_max_tau0_over_l_max*ppt->l_scalar_max
k_max_cmb[ppt->index_md_scalars] = ppr->k_max_tau0_over_l_max*ppt->l_scalar_max
/pba->conformal_age/pth->angular_rescaling;
k_max_cl = k_max_cmb;
k_max = k_max_cmb;
k_max_cl[ppt->index_md_scalars] = k_max_cmb[ppt->index_md_scalars];
k_max = k_max_cmb[ppt->index_md_scalars];
/* find k_max_cl: */
/* find k_max_cl[ppt->index_md_scalars] : */
/* if we need density/lensing Cl's, we must impose a stronger condition,
such that the minimum wavelength on the shell corresponding
to the center of smallest redshift bin is seen under an
angle smaller than pi/lmax. So we must mutiply our previous
k_max_cl by the ratio tau0/(tau0-tau[center of smallest
k_max_cl[ppt->index_md_scalars] by the ratio tau0/(tau0-tau[center of smallest
redhsift bin]). Note that we could do the same with the
lensing potential if we needed a very precise C_l^phi-phi at
large l. We don't do it by default, because the lensed ClT,
......@@ -1274,8 +1287,8 @@ int perturb_get_k_list(
pba->error_message,
ppt->error_message);
k_max_cl = MAX(k_max_cl,ppr->k_max_tau0_over_l_max*ppt->l_lss_max/(pba->conformal_age-tau1)); // to be very accurate we should use angular diameter distance to given redhsift instead of comoving radius: would implement corrections dependning on curvature
k_max = k_max_cl;
k_max_cl[ppt->index_md_scalars] = MAX(k_max_cl[ppt->index_md_scalars],ppr->k_max_tau0_over_l_max*ppt->l_lss_max/(pba->conformal_age-tau1)); // to be very accurate we should use angular diameter distance to given redhsift instead of comoving radius: would implement corrections dependning on curvature
k_max = k_max_cl[ppt->index_md_scalars];
}
}
......@@ -1313,7 +1326,7 @@ int perturb_get_k_list(
/* allocate array with, for the moment, the largest possible size */
class_alloc(ppt->k[ppt->index_md_scalars],
((int)((k_max_cmb-k_min)/k_rec/MIN(ppr->k_step_super,ppr->k_step_sub))+
((int)((k_max_cmb[ppt->index_md_scalars]-k_min)/k_rec/MIN(ppr->k_step_super,ppr->k_step_sub))+
(int)(MAX(ppr->k_per_decade_for_pk,ppr->k_per_decade_for_bao)*log(k_max/k_min)/log(10.))+3)
*sizeof(double),ppt->error_message);
......@@ -1324,9 +1337,9 @@ int perturb_get_k_list(
ppt->k[ppt->index_md_scalars][index_k] = k;
index_k++;
/* values until k_max_cmb */
/* values until k_max_cmb[ppt->index_md_scalars] */
while (k < k_max_cmb) {
while (k < k_max_cmb[ppt->index_md_scalars]) {
/* the linear step is not constant, it has a step-like shape,
centered around the characteristic scale set by the sound
......@@ -1370,9 +1383,9 @@ int perturb_get_k_list(
ppt->k_size_cmb[ppt->index_md_scalars] = index_k;
/* values until k_max_cl */
/* values until k_max_cl[ppt->index_md_scalars] */
while (k < k_max_cl) {
while (k < k_max_cl[ppt->index_md_scalars]) {
k *= pow(10.,1./(ppr->k_per_decade_for_pk
+(ppr->k_per_decade_for_bao-ppr->k_per_decade_for_pk)
......@@ -1402,7 +1415,6 @@ int perturb_get_k_list(
ppt->k[ppt->index_md_scalars],
ppt->k_size[ppt->index_md_scalars]*sizeof(double),
ppt->error_message);
}
/** - vector modes */
......@@ -1426,12 +1438,12 @@ int perturb_get_k_list(
k_min = sqrt((7.-1.e-4)*pba->K);
}
/** - find k_max (as well as k_max_cmb, k_max_cl) */
/** - find k_max (as well as k_max_cmb[ppt->index_md_vectors], k_max_cl[ppt->index_md_vectors]) */
k_rec = 2. * _PI_ / pth->rs_rec; /* comoving scale corresponding to sound horizon at recombination */
k_max_cmb = k_min;
k_max_cl = k_min;
k_max_cmb[ppt->index_md_vectors] = k_min;
k_max_cl[ppt->index_md_vectors] = k_min;
k_max = k_min;
if (ppt->has_cls == _TRUE_) {
......@@ -1443,10 +1455,10 @@ int perturb_get_k_list(
pi/lmax: this is equivalent to
k_max_cl*[comvoving.ang.diameter.distance] > l_max */
k_max_cmb = ppr->k_max_tau0_over_l_max*ppt->l_vector_max
k_max_cmb[ppt->index_md_vectors] = ppr->k_max_tau0_over_l_max*ppt->l_vector_max
/pba->conformal_age/pth->angular_rescaling;
k_max_cl = k_max_cmb;
k_max = k_max_cmb;
k_max_cl[ppt->index_md_vectors] = k_max_cmb[ppt->index_md_vectors];
k_max = k_max_cmb[ppt->index_md_vectors];
}
/** - test that result for k_min, k_max make sense */
......@@ -1475,7 +1487,7 @@ int perturb_get_k_list(
/* allocate array with, for the moment, the largest possible size */
class_alloc(ppt->k[ppt->index_md_vectors],
((int)((k_max_cmb-k_min)/k_rec/MIN(ppr->k_step_super,ppr->k_step_sub))+1)
((int)((k_max_cmb[ppt->index_md_vectors]-k_min)/k_rec/MIN(ppr->k_step_super,ppr->k_step_sub))+1)
*sizeof(double),ppt->error_message);
/* first value */
......@@ -1485,9 +1497,9 @@ int perturb_get_k_list(
ppt->k[ppt->index_md_vectors][index_k] = k;
index_k++;
/* values until k_max_cmb */
/* values until k_max_cmb[ppt->index_md_vectors] */
while (k < k_max_cmb) {
while (k < k_max_cmb[ppt->index_md_vectors]) {
/* the linear step is not constant, it has a step-like shape,
centered around the characteristic scale set by the sound
......@@ -1537,7 +1549,6 @@ int perturb_get_k_list(
ppt->k[ppt->index_md_vectors],
ppt->k_size[ppt->index_md_vectors]*sizeof(double),
ppt->error_message);
}
/** - tensor modes */
......@@ -1561,27 +1572,27 @@ int perturb_get_k_list(
k_min = sqrt((6.-1.e-4)*pba->K);
}
/** - find k_max (as well as k_max_cmb, k_max_cl) */
/** - find k_max (as well as k_max_cmb[ppt->index_md_tensors], k_max_cl[ppt->index_md_tensors]) */
k_rec = 2. * _PI_ / pth->rs_rec; /* comoving scale corresponding to sound horizon at recombination */
k_max_cmb = k_min;
k_max_cl = k_min;
k_max_cmb[ppt->index_md_tensors] = k_min;
k_max_cl[ppt->index_md_tensors] = k_min;
k_max = k_min;
if (ppt->has_cls == _TRUE_) {
/* find k_max_cmb: */
/* find k_max_cmb[ppt->index_md_tensors]: */
/* choose a k_max_cmb corresponding to a wavelength on the last
/* choose a k_max_cmb[ppt->index_md_tensors] corresponding to a wavelength on the last
scattering surface seen today under an angle smaller than
pi/lmax: this is equivalent to
k_max_cl*[comvoving.ang.diameter.distance] > l_max */
k_max_cl[ppt->index_md_tensors]*[comvoving.ang.diameter.distance] > l_max */
k_max_cmb = ppr->k_max_tau0_over_l_max*ppt->l_tensor_max