diff --git a/CMakeLists.txt b/CMakeLists.txt index b37b2c3..2d651c4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -179,6 +179,9 @@ elseif(CODE_PRECISION STREQUAL "LONGDOUBLE") target_link_libraries(${PRGNAME} PRIVATE FFTW3::FFTW3_LONGDOUBLE_SERIAL) endif() +target_include_directories(${PRGNAME} PRIVATE ${GSL_INCLUDE_DIR}) +target_link_libraries(${PRGNAME} PRIVATE GSL::gsl) + if(HDF5_FOUND) target_link_libraries(${PRGNAME} PRIVATE ${HDF5_LIBRARIES}) target_include_directories(${PRGNAME} PRIVATE ${HDF5_INCLUDE_DIRS}) @@ -192,10 +195,22 @@ if(TIRPC_FOUND) target_compile_options(${PRGNAME} PRIVATE "-DHAVE_TIRPC") endif(TIRPC_FOUND) +######################################################################################################################## +# include PANPHASIA if(ENABLE_PANPHASIA) target_compile_options(${PRGNAME} PRIVATE "-DHAVE_PANPHASIA") endif(ENABLE_PANPHASIA) +######################################################################################################################## -target_include_directories(${PRGNAME} PRIVATE ${GSL_INCLUDE_DIR}) -target_link_libraries(${PRGNAME} PRIVATE GSL::gsl) + +######################################################################################################################## +# include CLASS +option(ENABLE_CLASS "Enable CLASS support (as a submodule)." ON) + +if(ENABLE_CLASS) + include(${CMAKE_CURRENT_SOURCE_DIR}/ext/class.cmake) + target_link_libraries(${PRGNAME} PRIVATE class::libclass_cpp) + target_compile_definitions(${PRGNAME} PRIVATE "USE_CLASS") +endif() +######################################################################################################################## diff --git a/ext/class.cmake b/ext/class.cmake new file mode 100644 index 0000000..ea6a8f5 --- /dev/null +++ b/ext/class.cmake @@ -0,0 +1,17 @@ +cmake_minimum_required(VERSION 3.11) +include(FetchContent) +FetchContent_Declare( + class + GIT_REPOSITORY https://github.com/ohahn/class_public.git + GIT_TAG master + GIT_SHALLOW YES + GIT_PROGRESS TRUE + USES_TERMINAL_DOWNLOAD TRUE # <---- this is needed only for Ninja +) + +FetchContent_GetProperties(class) +if(NOT class_POPULATED) + set(FETCHCONTENT_QUIET OFF) + FetchContent_Populate(class) + add_subdirectory(${class_SOURCE_DIR} ${class_BINARY_DIR}) +endif() diff --git a/src/general.hh b/src/general.hh index af1a146..f7eb32f 100644 --- a/src/general.hh +++ b/src/general.hh @@ -80,70 +80,26 @@ inline T POW4( T a ){ //return a*a*a*a; } +#if defined(_OPENMP) +#include +inline double get_wtime() +{ + return omp_get_wtime(); +} +#else +#include +inline double get_wtime() +{ + return std::clock() / double(CLOCKS_PER_SEC); +} +#endif -//! structure for cosmological parameters -/* -typedef struct cosmology_old{ - double - Omega_m, //!< baryon+dark matter density - Omega_b, //!< baryon matter density - Omega_DE, //!< dark energy density (cosmological constant or parameterised) - Omega_r, //!< photon + relativistic particle density - Omega_k, //!< curvature density - H0, //!< Hubble constant in km/s/Mpc - nspect, //!< long-wave spectral index (scale free is nspect=1) - sigma8, //!< power spectrum normalization - w_0, //!< dark energy equation of state parameter 1: w = w0 + a * wa - w_a, //!< dark energy equation of state parameter 2: w = w0 + a * wa - - //Gamma, //!< shape parameter (of historical interest, as a free parameter) - //fnl, //!< non-gaussian contribution parameter - //w0, //!< dark energy equation of state parameter (not implemented, i.e. =1 at the moment) - //wa, //!< dark energy equation of state parameter (not implemented, i.e. =1 at the moment) - dplus, //!< linear perturbation growth factor - pnorm, //!< actual power spectrum normalisation factor - vfact, //!< velocity<->displacement conversion factor in Zel'dovich approx. - WDMmass, //!< Warm DM particle mass - WDMg_x, //!< Warm DM particle degrees of freedom - astart; //!< expansion factor a for which to generate initial conditions - - cosmology( config_file cf ) - { - double zstart = cf.get_value( "setup", "zstart" ); - - astart = 1.0/(1.0+zstart); - Omega_b = cf.get_value( "cosmology", "Omega_b" ); - Omega_m = cf.get_value( "cosmology", "Omega_m" ); - Omega_DE = cf.get_value( "cosmology", "Omega_L" ); - w_0 = cf.get_value_safe( "cosmology", "w0", -1.0 ); - w_a = cf.get_value_safe( "cosmology", "wa", 0.0 ); - - Omega_r = cf.get_value_safe( "cosmology", "Omega_r", 0.0 ); // no longer default to nonzero (8.3e-5) - Omega_k = 1.0 - Omega_m - Omega_DE - Omega_r; - - H0 = cf.get_value( "cosmology", "H0" ); - sigma8 = cf.get_value( "cosmology", "sigma_8" ); - nspect = cf.get_value( "cosmology", "nspec" ); - WDMg_x = cf.get_value_safe( "cosmology", "WDMg_x", 1.5 ); - WDMmass = cf.get_value_safe( "cosmology", "WDMmass", 0.0 ); - - dplus = 0.0; - pnorm = 0.0; - vfact = 0.0; - } - - cosmology( void ) - { - - } -}Cosmology;*/ - -//! basic box/grid/refinement structure parameters -typedef struct { - unsigned levelmin, levelmax; - double boxlength; - std::vector offx,offy,offz,llx,lly,llz; -}Parameters; +// //! basic box/grid/refinement structure parameters +// typedef struct { +// unsigned levelmin, levelmax; +// double boxlength; +// std::vector offx,offy,offz,llx,lly,llz; +// }Parameters; //! measure elapsed wallclock time inline double time_seconds( void ) diff --git a/src/plugins/random_panphasia.cc b/src/plugins/random_panphasia.cc index 37e8fef..84dfd82 100644 --- a/src/plugins/random_panphasia.cc +++ b/src/plugins/random_panphasia.cc @@ -261,16 +261,6 @@ void RNG_panphasia::backward_transform_field(real_t *field, int nx, int ny, int FFTW_API(destroy_plan(ipf)); } -#include -inline double get_wtime(void) -{ -#ifdef _OPENMP - return omp_get_wtime(); -#else - return (double)clock() / CLOCKS_PER_SEC; -#endif -} - void RNG_panphasia::fill_grid(int level, DensityGrid &R) { real_t *pr0, *pr1, *pr2, *pr3, *pr4; diff --git a/src/plugins/transfer_CLASS.cc b/src/plugins/transfer_CLASS.cc new file mode 100644 index 0000000..fbf5c40 --- /dev/null +++ b/src/plugins/transfer_CLASS.cc @@ -0,0 +1,336 @@ +// This file is part of monofonIC (MUSIC2) +// A software package to generate ICs for cosmological simulations +// Copyright (C) 2020 by Oliver Hahn +// +// monofonIC is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// monofonIC is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifdef USE_CLASS + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +// #include + +#include + +class transfer_CLASS_plugin : public transfer_function_plugin +{ +private: + + using transfer_function_plugin::cosmo_params_; + + interpolated_function_1d delta_c_, delta_b_, delta_n_, delta_m_, theta_c_, theta_b_, theta_n_, theta_m_; + interpolated_function_1d delta_c0_, delta_b0_, delta_n0_, delta_m0_, theta_c0_, theta_b0_, theta_n0_, theta_m0_; + + double zstart_, ztarget_, astart_, atarget_, kmax_, kmin_, h_, tnorm_; + + ClassParams pars_; + std::unique_ptr the_ClassEngine_; + std::ofstream ofs_class_input_; + + template + void add_class_parameter(std::string parameter_name, const T parameter_value) + { + pars_.add(parameter_name, parameter_value); + ofs_class_input_ << parameter_name << " = " << parameter_value << std::endl; + } + + //! Set up class parameters from MUSIC cosmological parameters + void init_ClassEngine(void) + { + //--- general parameters ------------------------------------------ + add_class_parameter("z_max_pk", std::max(std::max(zstart_, ztarget_),199.0)); // use 1.2 as safety + add_class_parameter("P_k_max_h/Mpc", std::max(2.0,kmax_)); + add_class_parameter("output", "dTk,vTk"); + add_class_parameter("extra metric transfer functions","yes"); + // add_class_parameter("lensing", "no"); + + //--- choose gauge ------------------------------------------------ + // add_class_parameter("extra metric transfer functions", "yes"); + add_class_parameter("gauge", "synchronous"); + + //--- cosmological parameters, densities -------------------------- + add_class_parameter("h", cosmo_params_.get("h")); + + add_class_parameter("Omega_b", cosmo_params_.get("Omega_b")); + add_class_parameter("Omega_cdm", cosmo_params_.get("Omega_c")); + add_class_parameter("Omega_k", cosmo_params_.get("Omega_k")); + add_class_parameter("Omega_fld", 0.0); + add_class_parameter("Omega_scf", 0.0); + + + // add_class_parameter("fluid_equation_of_state","CLP"); + // add_class_parameter("w0_fld", -1 ); + // add_class_parameter("wa_fld", 0. ); + // add_class_parameter("cs2_fld", 1); + + //--- massive neutrinos ------------------------------------------- +#if 0 + //default off + // add_class_parameter("Omega_ur",0.0); + add_class_parameter("N_ur", cosmo_params_.get("N_ur")); + add_class_parameter("N_ncdm", 0); + +#else + + add_class_parameter("N_ur", cosmo_params_.get("N_ur")); + add_class_parameter("N_ncdm", cosmo_params_.get("N_nu_massive")); + if( cosmo_params_.get("N_nu_massive") > 0 ){ + std::stringstream sstr; + if( cosmo_params_.get("m_nu1") > 1e-9 ) sstr << cosmo_params_.get("m_nu1"); + if( cosmo_params_.get("m_nu2") > 1e-9 ) sstr << ", " << cosmo_params_.get("m_nu2"); + if( cosmo_params_.get("m_nu3") > 1e-9 ) sstr << ", " << cosmo_params_.get("m_nu3"); + add_class_parameter("m_ncdm", sstr.str().c_str()); + } + + // change above to enable + //add_class_parameter("omega_ncdm", 0.0006451439); + //add_class_parameter("m_ncdm", "0.4"); + //add_class_parameter("T_ncdm", 0.71611); +#endif + + //--- cosmological parameters, primordial ------------------------- + add_class_parameter("P_k_ini type", "analytic_Pk"); + + if( cosmo_params_.get("A_s") > 0.0 ){ + add_class_parameter("A_s", cosmo_params_.get("A_s")); + }else{ + add_class_parameter("sigma8", cosmo_params_.get("sigma_8")); + } + add_class_parameter("n_s", cosmo_params_.get("n_s")); + add_class_parameter("alpha_s", 0.0); + add_class_parameter("T_cmb", cosmo_params_.get("Tcmb")); + add_class_parameter("YHe", cosmo_params_.get("YHe")); + + // additional parameters + add_class_parameter("reio_parametrization", "reio_none"); + + // precision parameters + add_class_parameter("k_per_decade_for_pk", 100); + add_class_parameter("k_per_decade_for_bao", 100); + add_class_parameter("compute damping scale", "yes"); + add_class_parameter("tol_perturb_integration", 1.e-8); + add_class_parameter("tol_background_integration", 1e-9); + + // high precision options from cl_permille.pre: + // precision file to be passed as input in order to achieve at least percent precision on scalar Cls + add_class_parameter("hyper_flat_approximation_nu", 7000.); + add_class_parameter("transfer_neglect_delta_k_S_t0", 0.17); + add_class_parameter("transfer_neglect_delta_k_S_t1", 0.05); + add_class_parameter("transfer_neglect_delta_k_S_t2", 0.17); + add_class_parameter("transfer_neglect_delta_k_S_e", 0.13); + add_class_parameter("delta_l_max", 1000); + + int class_verbosity = 0; + + add_class_parameter("background_verbose", class_verbosity); + add_class_parameter("thermodynamics_verbose", class_verbosity); + add_class_parameter("perturbations_verbose", class_verbosity); + add_class_parameter("transfer_verbose", class_verbosity); + add_class_parameter("primordial_verbose", class_verbosity); + add_class_parameter("spectra_verbose", class_verbosity); + add_class_parameter("nonlinear_verbose", class_verbosity); + add_class_parameter("lensing_verbose", class_verbosity); + add_class_parameter("output_verbose", class_verbosity); + + // output parameters, only needed for the control CLASS .ini file that we output + std::stringstream zlist; + if (ztarget_ == zstart_) + zlist << ztarget_ << ((ztarget_!=0.0)? ", 0.0" : ""); + else + zlist << std::max(ztarget_, zstart_) << ", " << std::min(ztarget_, zstart_) << ", 0.0"; + add_class_parameter("z_pk", zlist.str()); + + music::ilog << "Computing transfer function via ClassEngine..." << std::endl; + double wtime = get_wtime(); + + the_ClassEngine_ = std::make_unique(pars_, false); + + wtime = get_wtime() - wtime; + music::ilog << "CLASS took " << wtime << " s." << std::endl; + } + + //! run ClassEngine with parameters set up + void run_ClassEngine(double z, std::vector &k, std::vector &dc, std::vector &tc, std::vector &db, std::vector &tb, + std::vector &dn, std::vector &tn, std::vector &dm, std::vector &tm) + { + k.clear(); + dc.clear(); db.clear(); dn.clear(); dm.clear(); + tc.clear(); tb.clear(); tn.clear(); tm.clear(); + + the_ClassEngine_->getTk(z, k, dc, db, dn, dm, tc, tb, tn, tm); + + const double h = cosmo_params_.get("h"); + + for (size_t i = 0; i < k.size(); ++i) + { + // convert to 'CAMB' format, since we interpolate loglog and + // don't want negative numbers... + auto ik2 = 1.0 / (k[i] * k[i]) * h * h; + dc[i] = -dc[i] * ik2; + db[i] = -db[i] * ik2; + dn[i] = -dn[i] * ik2; + dm[i] = -dm[i] * ik2; + tc[i] = -tc[i] * ik2; + tb[i] = -tb[i] * ik2; + tn[i] = -tn[i] * ik2; + tm[i] = -tm[i] * ik2; + } + } + +public: + explicit transfer_CLASS_plugin(config_file &cf, const cosmology::parameters& cosmo_params) + : transfer_function_plugin(cf,cosmo_params) + { + this->tf_isnormalised_ = true; + + ofs_class_input_.open("input_class_parameters.ini", std::ios::trunc); + + // all cosmological parameters need to be passed through the_cosmo_calc + + ztarget_ = pcf_->get_value_safe("cosmology", "ztarget", 0.0); + atarget_ = 1.0 / (1.0 + ztarget_); + zstart_ = pcf_->get_value("setup", "zstart"); + astart_ = 1.0 / (1.0 + zstart_); + + h_ = cosmo_params_["h"]; + + if (cosmo_params_["A_s"] > 0.0) { + music::ilog << "CLASS: Using A_s=" << cosmo_params_["A_s"] << " to normalise the transfer function." << std::endl; + }else{ + double sigma8 = cosmo_params_["sigma_8"]; + if( sigma8 < 0 ){ + throw std::runtime_error("Need to specify either A_s or sigma_8 for CLASS plugin..."); + } + music::ilog << "CLASS: Using sigma8_ =" << sigma8<< " to normalise the transfer function." << std::endl; + } + + // determine highest k we will need for the resolution selected + double lbox = pcf_->get_value("setup", "BoxLength"); + int nres = pcf_->get_value("setup", "GridRes"); + kmax_ = std::max(20.0, 2.0 * M_PI / lbox * nres / 2 * sqrt(3) * 2.0); // 120% of spatial diagonal, or k=10h Mpc-1 + + // initialise CLASS and get the normalisation + this->init_ClassEngine(); + double A_s_ = the_ClassEngine_->get_A_s(); // this either the input one, or the one computed from sigma8 + + // compute the normalisation to interface with MUSIC + double k_p = cosmo_params["k_p"] / cosmo_params["h"]; + tnorm_ = std::sqrt(2.0 * M_PI * M_PI * A_s_ * std::pow(1.0 / k_p, cosmo_params["n_s"] - 1) / std::pow(2.0 * M_PI, 3.0)); + + // compute the transfer function at z=0 using CLASS engine + std::vector k, dc, tc, db, tb, dn, tn, dm, tm; + this->run_ClassEngine(0.0, k, dc, tc, db, tb, dn, tn, dm, tm); + + delta_c0_.set_data(k, dc); + theta_c0_.set_data(k, tc); + delta_b0_.set_data(k, db); + theta_b0_.set_data(k, tb); + delta_n0_.set_data(k, dn); + theta_n0_.set_data(k, tn); + delta_m0_.set_data(k, dm); + theta_m0_.set_data(k, tm); + + // compute the transfer function at z=z_target using CLASS engine + this->run_ClassEngine(ztarget_, k, dc, tc, db, tb, dn, tn, dm, tm); + delta_c_.set_data(k, dc); + theta_c_.set_data(k, tc); + delta_b_.set_data(k, db); + theta_b_.set_data(k, tb); + delta_n_.set_data(k, dn); + theta_n_.set_data(k, tn); + delta_m_.set_data(k, dm); + theta_m_.set_data(k, tm); + + kmin_ = k[0]; + kmax_ = k.back(); + + music::ilog << "CLASS table contains k = " << this->get_kmin() << " to " << this->get_kmax() << " h Mpc-1." << std::endl; + + tf_distinct_ = true; + tf_withvel_ = true; + tf_withtotal0_ = true; + } + + ~transfer_CLASS_plugin() + { + } + + inline double compute(double k, tf_type type) const + { + k *= h_; + + if (k < kmin_ || k > kmax_) + { + return 0.0; + } + + real_t val(0.0); + switch (type) + { + // values at ztarget: + case delta_matter: + val = delta_m_(k); break; + case delta_cdm: + val = delta_c_(k); break; + case delta_baryon: + val = delta_b_(k); break; + case theta_matter: + val = theta_m_(k); break; + case theta_cdm: + val = theta_c_(k); break; + case theta_baryon: + val = theta_b_(k); break; + case delta_bc: + val = delta_b_(k)-delta_c_(k); break; + case theta_bc: + val = theta_b_(k)-theta_c_(k); break; + + // values at zstart: + case delta_matter0: + val = delta_m0_(k); break; + case delta_cdm0: + val = delta_c0_(k); break; + case delta_baryon0: + val = delta_b0_(k); break; + case theta_matter0: + val = theta_m0_(k); break; + case theta_cdm0: + val = theta_c0_(k); break; + case theta_baryon0: + val = theta_b0_(k); break; + default: + throw std::runtime_error("Invalid type requested in transfer function evaluation"); + } + return val * tnorm_; + } + + inline double get_kmin(void) const { return kmin_ / h_; } + inline double get_kmax(void) const { return kmax_ / h_; } +}; + +namespace +{ +transfer_function_plugin_creator_concrete creator("CLASS"); +} + +#endif // USE_CLASS diff --git a/src/plugins/transfer_bbks.cc b/src/plugins/transfer_bbks.cc index e26ef80..5bf0014 100644 --- a/src/plugins/transfer_bbks.cc +++ b/src/plugins/transfer_bbks.cc @@ -57,7 +57,7 @@ public: } //! computes the value of the BBKS transfer function for mode k (in h/Mpc) - inline double compute( double k, tf_type type ){ + inline double compute( double k, tf_type type ) const{ double q, f1, f2; if(k < 1e-7 ) @@ -71,11 +71,11 @@ public: } - inline double get_kmin( void ){ + inline double get_kmin( void ) const{ return 1e-4; } - inline double get_kmax( void ){ + inline double get_kmax( void ) const{ return 1.e4; } }; diff --git a/src/plugins/transfer_camb.cc b/src/plugins/transfer_camb.cc index abf4a93..793edd5 100644 --- a/src/plugins/transfer_camb.cc +++ b/src/plugins/transfer_camb.cc @@ -221,7 +221,7 @@ public: } // linear interpolation in log-log - inline double extrap_right(double k, const tf_type &type) { + inline double extrap_right(double k, const tf_type &type) const { int n = m_tab_k.size() - 1, n1 = n - 1; double v1(1.0), v2(1.0); @@ -267,7 +267,7 @@ public: return 0.0; } - inline double compute(double k, tf_type type) { + inline double compute(double k, tf_type type) const{ // use constant interpolation on the left side of the tabulated values if (k < m_kmin) { switch (type) { @@ -320,9 +320,9 @@ public: } } - inline double get_kmin(void) { return pow(10.0, m_tab_k[1]); } + inline double get_kmin(void) const { return pow(10.0, m_tab_k[1]); } - inline double get_kmax(void) { + inline double get_kmax(void) const { return pow(10.0, m_tab_k[m_tab_k.size() - 2]); } }; diff --git a/src/plugins/transfer_eisenstein.cc b/src/plugins/transfer_eisenstein.cc index 9832805..bec6cf1 100644 --- a/src/plugins/transfer_eisenstein.cc +++ b/src/plugins/transfer_eisenstein.cc @@ -102,7 +102,7 @@ struct eisenstein_transfer } //! private member function: computes transfer function for mode k (k in Mpc) - inline double TFfit_onek(double k, double *tf_baryon, double *tf_cdm) + inline double TFfit_onek(double k, double *tf_baryon, double *tf_cdm) const /* Input: k -- Wavenumber at which to calculate transfer function, in Mpc^-1. *tf_baryon, *tf_cdm -- Input value not used; replaced on output if the input was not NULL. */ @@ -173,7 +173,7 @@ struct eisenstein_transfer fc_ = (cp["Omega_m"] - cp["Omega_b"]) / cp["Omega_m"] ; } - inline double at_k(double k) + inline double at_k(double k) const { double tfb, tfcdm; TFfit_onek(k * m_h0, &tfb, &tfcdm); @@ -211,17 +211,17 @@ public: } //! Computes the transfer function for k in Mpc/h by calling TFfit_onek - inline double compute(double k, tf_type type) + inline double compute(double k, tf_type type) const { return etf_.at_k(k); } - inline double get_kmin(void) + inline double get_kmin(void) const { return 1e-5; } - inline double get_kmax(void) + inline double get_kmax(void) const { return 1.e5; } @@ -303,17 +303,17 @@ public: std::cerr << "WDM alpha = " << m_WDMalpha << std::endl; } - inline double compute(double k, tf_type type) + inline double compute(double k, tf_type type) const { return etf_.at_k(k) * pow(1.0 + pow(m_WDMalpha * k, 2.0 * wdmnu_), -5.0 / wdmnu_); } - inline double get_kmin(void) + inline double get_kmin(void) const { return 1e-4; } - inline double get_kmax(void) + inline double get_kmax(void) const { return 1.e4; } @@ -348,7 +348,7 @@ public: music::ilog.Print(" bino CDM: k_fs = %g, k_d = %g", kfs_, kd_); } - inline double compute(double k, tf_type type) + inline double compute(double k, tf_type type) const { double kkfs = k / kfs_; double kkfs2 = kkfs * kkfs; @@ -362,12 +362,12 @@ public: return 0.0; } - inline double get_kmin(void) + inline double get_kmin(void) const { return 1e-4; } - inline double get_kmax(void) + inline double get_kmax(void) const { return 1.e8; } diff --git a/src/plugins/transfer_linger++.cc b/src/plugins/transfer_linger++.cc index ce1b196..60d009d 100644 --- a/src/plugins/transfer_linger++.cc +++ b/src/plugins/transfer_linger++.cc @@ -197,7 +197,7 @@ public: gsl_interp_accel_free(acc_dtot0); } - inline double extrap_left(double k, const tf_type &type) + inline double extrap_left(double k, const tf_type &type) const { if (k < 1e-8) return 1.0; @@ -246,7 +246,7 @@ public: return pow(10.0, (v2 - v1) / dk * (delk) + v1); } - inline double extrap_right(double k, const tf_type &type) + inline double extrap_right(double k, const tf_type &type) const { double v1(1.0), v2(1.0); @@ -294,7 +294,7 @@ public: return pow(10.0, (v2 - v1) / dk * (delk) + v2); } - inline double compute(double k, tf_type type) + inline double compute(double k, tf_type type) const { double lk = log10(k); @@ -335,12 +335,12 @@ public: return 1.0; } - inline double get_kmin(void) + inline double get_kmin(void) const { return pow(10.0, m_tab_k[0]); } - inline double get_kmax(void) + inline double get_kmax(void) const { return pow(10.0, m_tab_k[m_tab_k.size() - 1]); } diff --git a/src/plugins/transfer_music.cc b/src/plugins/transfer_music.cc index 687c6e2..7164764 100644 --- a/src/plugins/transfer_music.cc +++ b/src/plugins/transfer_music.cc @@ -161,7 +161,7 @@ public: gsl_interp_accel_free (acc_theta_baryon); } - inline double extrap_left( double k, const tf_type& type ) + inline double extrap_left( double k, const tf_type& type ) const { if( k<1e-8 ) return 1.0; @@ -202,7 +202,7 @@ public: return pow(10.0,(v2-v1)/dk*(delk)+v1); } - inline double extrap_right( double k, const tf_type& type ) + inline double extrap_right( double k, const tf_type& type ) const { double v1(1.0), v2(1.0); @@ -242,7 +242,7 @@ public: return pow(10.0,(v2-v1)/dk*(delk)+v2); } - inline double compute( double k, tf_type type ){ + inline double compute( double k, tf_type type ) const{ double lk = log10(k); @@ -279,11 +279,11 @@ public: return 1.0; } - inline double get_kmin( void ){ + inline double get_kmin( void ) const{ return pow(10.0,m_tab_k[0]); } - inline double get_kmax( void ){ + inline double get_kmax( void ) const{ return pow(10.0,m_tab_k[m_tab_k.size()-1]); } diff --git a/src/transfer_function.hh b/src/transfer_function.hh index 3940b6e..0513557 100644 --- a/src/transfer_function.hh +++ b/src/transfer_function.hh @@ -74,13 +74,13 @@ public: virtual void intialise(void) {} //! compute value of transfer function at waven umber - virtual double compute(double k, tf_type type) = 0; + virtual double compute(double k, tf_type type) const = 0; //! return maximum wave number allowed - virtual double get_kmax(void) = 0; + virtual double get_kmax(void) const = 0; //! return minimum wave number allowed - virtual double get_kmin(void) = 0; + virtual double get_kmin(void) const = 0; //! return if density transfer function is distinct for baryons and DM bool tf_is_distinct(void) const