diff --git a/example.conf b/example.conf index 6b7d706..6b57a7c 100644 --- a/example.conf +++ b/example.conf @@ -14,6 +14,7 @@ zstart = 24.0 # starting redshift LPTorder = 3 # order of the LPT to be used (1,2 or 3) DoBaryons = no # also do baryon ICs? +DoBaryonVrel = no # if doing baryons, incl. also relative velocity to linear order? DoFixing = yes # do mode fixing à la Angulo&Pontzen (https://arxiv.org/abs/1603.05253) DoInversion = no # invert phases (for paired simulations) diff --git a/include/cosmology_calculator.hh b/include/cosmology_calculator.hh index bf0ef1b..21c9376 100644 --- a/include/cosmology_calculator.hh +++ b/include/cosmology_calculator.hh @@ -38,7 +38,8 @@ namespace cosmology * @brief provides functions to compute cosmological quantities * * This class provides member functions to compute cosmological quantities - * related to the Friedmann equations and linear perturbation theory + * related to the Friedmann equations and linear perturbation theory, it also + * provides the functionality to work with back-scaled cosmological fields */ class calculator { @@ -54,34 +55,42 @@ private: interpolated_function_1d D_of_a_, f_of_a_, a_of_D_; double Dnow_, Dplus_start_, Dplus_target_, astart_, atarget_; + //! wrapper for GSL adaptive integration routine, do not use if many integrations need to be done as it allocates and deallocates memory + //! set to 61-point Gauss-Kronrod and large workspace, used for sigma_8 normalisation real_t integrate(double (*func)(double x, void *params), double a, double b, void *params) const { + constexpr size_t wspace_size{100000}; + + double result{0.0}; + double error{0.0}; + gsl_function F; F.function = func; F.params = params; - - double result; - double error; - - gsl_set_error_handler_off(); - gsl_integration_workspace *w = gsl_integration_workspace_alloc(100000); - gsl_integration_qag(&F, a, b, 0, REL_PRECISION, 100000, 6, w, &result, &error); - - gsl_integration_workspace_free(w); - - gsl_set_error_handler(NULL); + + auto errh = gsl_set_error_handler_off(); + gsl_integration_workspace *wspace = gsl_integration_workspace_alloc(wspace_size); + gsl_integration_qag(&F, a, b, 0, REL_PRECISION, wspace_size, GSL_INTEG_GAUSS61, wspace, &result, &error); + gsl_integration_workspace_free(wspace); + gsl_set_error_handler(errh); if (error / result > REL_PRECISION) music::wlog << "no convergence in function 'integrate', rel. error=" << error / result << std::endl; - return (real_t)result; + return static_cast(result); } + //! compute the linear theory growth factor D+ by solving the single fluid ODE, returns tables D(a), f(a) + /*! + * @param tab_a reference to STL vector for values of a at which table entries exist + * @param tab_D reference to STL vector for values D(a) with a from tab_a + * @param tab_f reference to STL vector for values f(a)=dlog D / dlog a with a from tab_a + */ void compute_growth( std::vector& tab_a, std::vector& tab_D, std::vector& tab_f ) { using v_t = vec_t<3, double>; - // set ICs + // set ICs, very deep in radiation domination const double a0 = 1e-10; const double D0 = a0; const double Dprime0 = 2.0 * D0 * H_of_a(a0) / std::pow(phys_const::c_SI, 2); @@ -125,7 +134,7 @@ private: tab_a.push_back(yy[0]); tab_D.push_back(yy[1]); - tab_f.push_back(yy[2]); + tab_f.push_back(yy[2]); // temporarily store D' in table dt = dtnext; } @@ -162,7 +171,7 @@ public: a_of_D_.set_data(tab_D,tab_a); Dnow_ = D_of_a_(1.0); - Dplus_start_ = D_of_a_( astart_ ) / Dnow_; + Dplus_start_ = D_of_a_( astart_ ) / Dnow_; Dplus_target_ = D_of_a_( atarget_ ) / Dnow_; music::ilog << "Linear growth factors: D+_target = " << Dplus_target_ << ", D+_start = " << Dplus_start_ << std::endl; @@ -185,9 +194,7 @@ public: << " : " << transfer_function_->get_kmax() << " h/Mpc" << std::endl; } - ~calculator() - { - } + ~calculator() { } //! Write out a correctly scaled power spectrum at time a void write_powerspectrum(real_t a, std::string fname) const @@ -213,22 +220,23 @@ public: << std::setw(20) << ("P_dbar(k,a=1)") << std::setw(20) << ("P_tcdm(k,a=1)") << std::setw(20) << ("P_tbar(k,a=1)") - << std::setw(20) << ("P_dtot(K,a=1)") << std::endl; for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01) { ofs << std::setw(20) << std::setprecision(10) << k - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, total)*Dplus_start_, 2.0) - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, cdm)*Dplus_start_, 2.0) - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, baryon)*Dplus_start_, 2.0) - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, vcdm)*Dplus_start_, 2.0) - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, vbaryon)*Dplus_start_, 2.0) - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, total0), 2.0) - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, cdm0), 2.0) - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, baryon0), 2.0) - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, vcdm0), 2.0) - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, vbaryon0), 2.0) - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, vtotal), 2.0) + << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter)*Dplus_start_, 2.0) + << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm)*Dplus_start_, 2.0) + << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon)*Dplus_start_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter)*Dplus_start_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm)*Dplus_start_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon)*Dplus_start_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_cdm)*Dplus_start_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_baryon)*Dplus_start_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter0)* Dplus_start_ / Dplus_target_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm0)* Dplus_start_ / Dplus_target_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon0)* Dplus_start_ / Dplus_target_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_cdm0)* Dplus_start_ / Dplus_target_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_baryon0)* Dplus_start_ / Dplus_target_, 2.0) << std::endl; } } @@ -257,15 +265,24 @@ public: double fb = cosmo_param_["f_b"], fc = cosmo_param_["f_c"]; for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01) { - double dm = this->get_amplitude(k, total) * Dplus_start_ / Dplus_target_; - double dbc = this->get_amplitude(k, baryon) - this->get_amplitude(k, cdm); - double db = dm + fc * dbc; - double dc = dm - fb * dbc; + const double dm = this->get_amplitude(k, delta_matter) * Dplus_start_ / Dplus_target_; + const double dbc = this->get_amplitude(k, delta_bc); + const double db = dm + fc * dbc; + const double dc = dm - fb * dbc; + const double tm = this->get_amplitude(k, delta_matter) * Dplus_start_ / Dplus_target_; + const double tbc = this->get_amplitude(k, theta_bc); + const double tb = dm + fc * dbc; + const double tc = dm - fb * dbc; + ofs << std::setw(20) << std::setprecision(10) << k << std::setw(20) << std::setprecision(10) << dc << std::setw(20) << std::setprecision(10) << db << std::setw(20) << std::setprecision(10) << dm - << std::setw(20) << std::setprecision(10) << dbc + << std::setw(20) << std::setprecision(10) << dbc + 2 * tbc * (std::sqrt( Dplus_target_ / Dplus_start_ ) - 1.0) + << std::setw(20) << std::setprecision(10) << tc / std::pow( Dplus_start_ / Dplus_target_, 0.5 ) + << std::setw(20) << std::setprecision(10) << tb / std::pow( Dplus_start_ / Dplus_target_, 0.5 ) + << std::setw(20) << std::setprecision(10) << tm / std::pow( Dplus_start_ / Dplus_target_, 0.5 ) + << std::setw(20) << std::setprecision(10) << tbc / std::pow( Dplus_start_ / Dplus_target_, 0.5 ) << std::endl; } } @@ -327,9 +344,9 @@ public: const double x = k * 8.0; const double w = (x < 0.001)? 1.0-0.1*x*x : 3.0 * (std::sin(x) - x * std::cos(x)) / (x * x * x); - + static double nspect = (double)pcc->cosmo_param_["n_s"]; - double tf = pcc->transfer_function_->compute(k, total); + double tf = pcc->transfer_function_->compute(k, delta_matter); //... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1 return k * k * w * w * pow((double)k, (double)nspect) * tf * tf; @@ -346,47 +363,42 @@ public: const double w = (x < 0.001)? 1.0-0.1*x*x : 3.0 * (std::sin(x) - x * std::cos(x)) / (x * x * x); static double nspect = static_cast(pcc->cosmo_param_["n_s"]); - double tf = pcc->transfer_function_->compute(k, total0); + double tf = pcc->transfer_function_->compute(k, delta_matter0); //... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1 return k * k * w * w * std::pow(k, nspect) * tf * tf; } - // static double dSigma_bc(double k, void *pParams) - // { - // if (k <= 0.0) - // return 0.0f; - - // cosmology::calculator *pcc = reinterpret_cast(pParams); - - // static double nspect = static_cast(pcc->cosmo_param_.nspect); - // double tf = pcc->transfer_function_->compute(k, deltabc); - - // //... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1 - // return k * k * std::pow(k, nspect) * tf * tf; // * cosmo_param_.sqrtpnorm * Dplus_target_; - // } - //! Computes the amplitude of a mode from the power spectrum - /*! Function evaluates the supplied transfer function ptransfer_fun_ - * and returns the amplitude of fluctuations at wave number k at z=0 + /*! Function evaluates the supplied transfer function transfer_function_ + * and returns the amplitude of fluctuations at wave number k (in h/Mpc) back-scaled to z=z_start * @param k wave number at which to evaluate + * @param type one of the species: {delta,theta}_{matter,cdm,baryon,neutrino} */ inline real_t get_amplitude( const real_t k, const tf_type type) const { - return std::pow(k, 0.5 * cosmo_param_["n_s"]) * transfer_function_->compute(k, type) * cosmo_param_["sqrtpnorm"];// * ((type!=deltabc)? 1.0 : 1.0/Dplus_target_); + return std::pow(k, 0.5 * cosmo_param_["n_s"]) * transfer_function_->compute(k, type) * cosmo_param_["sqrtpnorm"]; } - inline real_t get_amplitude_phibc( const real_t k ) const + //! Compute amplitude of the back-scaled delta_bc mode, with decaying velocity v_bc included or not (in which case delta_bc=const) + inline real_t get_amplitude_delta_bc( const real_t k, bool withvbc ) const { + const real_t Dratio = Dplus_target_ / Dplus_start_; + const real_t dbc = transfer_function_->compute(k, delta_bc) + (withvbc? 2 * transfer_function_->compute(k, theta_bc) * (std::sqrt(Dratio) - 1.0) : 0.0); // need to multiply with Dplus_target since sqrtpnorm rescales like that - return -std::pow(k, 0.5 * cosmo_param_["n_s"]-2.0) * transfer_function_->compute(k, deltabc) * cosmo_param_["sqrtpnorm"] * Dplus_target_; + return std::pow(k, 0.5 * cosmo_param_["n_s"]) * dbc * (cosmo_param_["sqrtpnorm"] * Dplus_target_); } - inline real_t get_amplitude_rhobc( const real_t k ) const + + //! Compute amplitude of the back-scaled relative velocity theta_bc mode if withvbc==true, otherwise return zero + inline real_t get_amplitude_theta_bc( const real_t k, bool withvbc ) const { + const real_t Dratio = Dplus_target_ / Dplus_start_; + const real_t tbc = transfer_function_->compute(k, theta_bc) * std::sqrt(Dratio); // need to multiply with Dplus_target since sqrtpnorm rescales like that - return std::pow(k, 0.5 * cosmo_param_["n_s"]) * transfer_function_->compute(k, deltabc) * cosmo_param_["sqrtpnorm"] * Dplus_target_; + return withvbc ? std::pow(k, 0.5 * cosmo_param_["n_s"]) * tbc * (cosmo_param_["sqrtpnorm"] * Dplus_target_) : 0.0; } + //! Computes the normalization for the power spectrum /*! * integrates the power spectrum to fix the normalization to that given @@ -407,19 +419,6 @@ public: return std::sqrt(sigma0); } - // real_t compute_sigma_bc( void ) - // { - // real_t sigma0, kmin, kmax; - // kmax = 100.0; //transfer_function_->get_kmax(); - // kmin = transfer_function_->get_kmin(); - - // sigma0 = 4.0 * M_PI * integrate(&dSigma_bc, static_cast(kmin), static_cast(kmax), this); - // sigma0 = std::sqrt(sigma0); - // sigma0 *= cosmo_param_.sqrtpnorm * Dplus_target_; - // std::cout << "kmin = " << kmin << ", kmax = " << kmax << ", sigma_bc = " << sigma0 << std::endl; - // return sigma0; - // } - //! Computes the normalization for the power spectrum /*! * integrates the power spectrum to fix the normalization to that given diff --git a/include/cosmology_parameters.hh b/include/cosmology_parameters.hh index 35523db..0505893 100644 --- a/include/cosmology_parameters.hh +++ b/include/cosmology_parameters.hh @@ -31,7 +31,7 @@ namespace cosmology std::map pmap_; //!< All parameters are stored here as key-value pairs public: - //!< get routine for cosmological parameter key-value pairs + //! get routine for cosmological parameter key-value pairs double get(const std::string &key) const { auto it = pmap_.find(key); @@ -44,7 +44,7 @@ namespace cosmology return it->second; } - //!< set routine for cosmological parameter key-value pairs + //! set routine for cosmological parameter key-value pairs void set(const std::string &key, const double value) { auto it = pmap_.find(key); @@ -60,16 +60,16 @@ namespace cosmology } } - //!< shortcut get routine for cosmological parameter key-value pairs through bracket operator + //! shortcut get routine for cosmological parameter key-value pairs through bracket operator inline double operator[](const std::string &key) const { return this->get(key); } - //!< no default constructor + //! no default constructor parameters() = delete; - //!< default copy constructor + //! default copy constructor parameters(const parameters &) = default; - //!< main constructor for explicit construction from input config file + //! main constructor for explicit construction from input config file explicit parameters( config_file &cf ) { // CMB @@ -154,7 +154,7 @@ namespace cosmology music::ilog << "A_s = " << std::setw(16) << this->get("A_s"); else music::ilog << "sigma_8 = " << std::setw(16) << this->get("sigma_8"); - music::ilog << "nspec = " << std::setw(16) << this->get("n_s") << std::endl; + music::ilog << "n_s = " << std::setw(16) << this->get("n_s") << std::endl; music::ilog << " Omega_c = " << std::setw(16) << this->get("Omega_c") << "Omega_b = " << std::setw(16) << this->get("Omega_b") << "Omega_m = " << std::setw(16) << this->get("Omega_m") << std::endl; music::ilog << " Omega_r = " << std::setw(16) << this->get("Omega_r") << "Omega_nu = " << std::setw(16) << this->get("Omega_nu_massive") << "∑m_nu = " << std::setw(11) << sum_m_nu << "eV" << std::endl; music::ilog << " Omega_DE = " << std::setw(16) << this->get("Omega_DE") << "w_0 = " << std::setw(16) << this->get("w_0") << "w_a = " << std::setw(16) << this->get("w_a") << std::endl; diff --git a/include/transfer_function_plugin.hh b/include/transfer_function_plugin.hh index bfd362d..fc9550f 100644 --- a/include/transfer_function_plugin.hh +++ b/include/transfer_function_plugin.hh @@ -25,19 +25,20 @@ enum tf_type { - total, - cdm, - baryon, - vtotal, - vcdm, - vbaryon, - deltabc, - total0, - cdm0, - baryon0, - vtotal0, - vcdm0, - vbaryon0, + delta_matter, + delta_cdm, + delta_baryon, + theta_matter, + theta_cdm, + theta_baryon, + delta_bc, + theta_bc, + delta_matter0, + delta_cdm0, + delta_baryon0, + theta_matter0, + theta_cdm0, + theta_baryon0, }; class TransferFunction_plugin diff --git a/src/ic_generator.cc b/src/ic_generator.cc index 7e0dc31..4a63648 100644 --- a/src/ic_generator.cc +++ b/src/ic_generator.cc @@ -45,7 +45,7 @@ int initialise( config_file& the_config ) the_random_number_generator = std::move(select_RNG_plugin(the_config)); the_cosmo_calc = std::make_unique(the_config); the_output_plugin = std::move(select_output_plugin(the_config, the_cosmo_calc)); - + return 0; } @@ -96,6 +96,9 @@ int run( config_file& the_config ) //-------------------------------------------------------------------------------------------------------- //! do baryon ICs? const bool bDoBaryons = the_config.get_value_safe("setup", "DoBaryons", false ); + //! enable also back-scaled decaying relative velocity mode? only first order! + const bool bDoLinearBCcorr = the_config.get_value_safe("cosmology", "DoBaryonVrel", false); + // compute mass fractions std::map< cosmo_species, double > Omega; if( bDoBaryons ){ double Om = the_config.get_value("cosmology", "Omega_m"); @@ -164,8 +167,6 @@ int run( config_file& the_config ) const double g1 = -Dplus0; const double g2 = ((LPTorder>1)? -3.0/7.0*Dplus0*Dplus0 : 0.0); const double g3 = ((LPTorder>2)? 1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0); - // const double g3a = ((LPTorder>2)? 1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0); - // const double g3b = ((LPTorder>2)? -10.0/21.*Dplus0*Dplus0*Dplus0 : 0.0); const double g3c = ((LPTorder>2)? 1.0/7.0*Dplus0*Dplus0*Dplus0 : 0.0); // vfac = d log D+ / dt @@ -328,7 +329,7 @@ int run( config_file& the_config ) phi.FourierTransformForward(false); phi.assign_function_of_grids_kdep([&](auto k, auto wn) { real_t kmod = k.norm(); - ccomplex_t delta = wn * the_cosmo_calc->get_amplitude(kmod, total); + ccomplex_t delta = wn * the_cosmo_calc->get_amplitude(kmod, delta_matter); return -delta / (kmod * kmod); }, wnoise); @@ -525,7 +526,7 @@ int run( config_file& the_config ) wnoise.FourierTransformForward(); rho.FourierTransformForward(false); rho.assign_function_of_grids_kdep( [&]( auto k, auto wn ){ - return wn * the_cosmo_calc->get_amplitude_rhobc(k.norm());; + return wn * the_cosmo_calc->get_amplitude_delta_bc(k.norm(),bDoLinearBCcorr); }, wnoise ); rho.zero_DC_mode(); rho.FourierTransformBackward(); @@ -556,7 +557,7 @@ int run( config_file& the_config ) wnoise.FourierTransformForward(); rho.FourierTransformForward(false); rho.assign_function_of_grids_kdep( [&]( auto k, auto wn ){ - return wn * the_cosmo_calc->get_amplitude_rhobc(k.norm());; + return wn * the_cosmo_calc->get_amplitude_delta_bc(k.norm(), false); }, wnoise ); rho.zero_DC_mode(); rho.FourierTransformBackward(); @@ -682,6 +683,7 @@ int run( config_file& the_config ) A3[1]->FourierTransformForward(); A3[2]->FourierTransformForward(); } + wnoise.FourierTransformForward(); // write out positions for( int idim=0; idim<3; ++idim ){ @@ -769,6 +771,12 @@ int run( config_file& the_config ) tmp.kelem(idx) += vfac3 * (lg.gradient(idimp,tmp.get_k3(i,j,k)) * A3[idimpp]->kelem(idx) - lg.gradient(idimpp,tmp.get_k3(i,j,k)) * A3[idimp]->kelem(idx)); } + // if multi-species, then add vbc component backwards + if( bDoBaryons & bDoLinearBCcorr ){ + real_t knorm = wnoise.get_k(i,j,k).norm(); + tmp.kelem(idx) -= vfac1 * C_species * the_cosmo_calc->get_amplitude_theta_bc(knorm, bDoLinearBCcorr) * wnoise.kelem(i,j,k) * lg.gradient(idim,tmp.get_k3(i,j,k)) / (knorm*knorm); + } + // correct with interpolation kernel if we used interpolation to read out the positions (for glasses) if( the_output_plugin->write_species_as( this_species ) == output_type::particles && lattice_type == particle::lattice_glass){ tmp.kelem(idx) *= interp.compensation_kernel( tmp.get_k(i,j,k) ); diff --git a/src/plugins/transfer_CAMB_file.cc b/src/plugins/transfer_CAMB_file.cc index 7869e16..53b2316 100644 --- a/src/plugins/transfer_CAMB_file.cc +++ b/src/plugins/transfer_CAMB_file.cc @@ -44,7 +44,7 @@ private: { music::ilog << "Reading tabulated transfer function data from file:" << std::endl << " \'" << filename << "\'" << std::endl; - + std::string line; std::ifstream ifs(filename.c_str()); @@ -183,32 +183,35 @@ public: switch (type) { - case total0: - case total: + case delta_matter0: + case delta_matter: return delta_m_(k); - case cdm0: - case cdm: + case delta_cdm0: + case delta_cdm: return delta_c_(k); - case baryon0: - case baryon: + case delta_baryon0: + case delta_baryon: return delta_b_(k); - case vtotal0: - case vtotal: + case theta_matter0: + case theta_matter: return theta_m_(k); - case vcdm0: - case vcdm: + case theta_cdm0: + case theta_cdm: return theta_c_(k); - case vbaryon0: - case vbaryon: + case theta_baryon0: + case theta_baryon: return theta_b_(k); - case deltabc: + case delta_bc: return delta_b_(k)-delta_c_(k); + + case theta_bc: + return theta_b_(k)-theta_c_(k); default: throw std::runtime_error("Invalid type requested in transfer function evaluation"); diff --git a/src/plugins/transfer_CLASS.cc b/src/plugins/transfer_CLASS.cc index b558000..71f320a 100644 --- a/src/plugins/transfer_CLASS.cc +++ b/src/plugins/transfer_CLASS.cc @@ -178,8 +178,6 @@ private: the_ClassEngine_->getTk(z, k, dc, db, dn, dm, tc, tb, tn, tm); const double h = cosmo_params_.get("h"); - const double fc = cosmo_params_.get("f_c"); - const double fb = cosmo_params_.get("f_b"); for (size_t i = 0; i < k.size(); ++i) { @@ -288,33 +286,35 @@ public: switch (type) { // values at ztarget: - case total: + case delta_matter: val = delta_m_(k); break; - case cdm: + case delta_cdm: val = delta_c_(k); break; - case baryon: + case delta_baryon: val = delta_b_(k); break; - case vtotal: + case theta_matter: val = theta_m_(k); break; - case vcdm: + case theta_cdm: val = theta_c_(k); break; - case vbaryon: + case theta_baryon: val = theta_b_(k); break; - case deltabc: + 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 total0: + case delta_matter0: val = delta_m0_(k); break; - case cdm0: + case delta_cdm0: val = delta_c0_(k); break; - case baryon0: + case delta_baryon0: val = delta_b0_(k); break; - case vtotal0: + case theta_matter0: val = theta_m0_(k); break; - case vcdm0: + case theta_cdm0: val = theta_c0_(k); break; - case vbaryon0: + case theta_baryon0: val = theta_b0_(k); break; default: throw std::runtime_error("Invalid type requested in transfer function evaluation"); diff --git a/src/plugins/transfer_eisenstein.cc b/src/plugins/transfer_eisenstein.cc index baf52b8..1c4fc5a 100644 --- a/src/plugins/transfer_eisenstein.cc +++ b/src/plugins/transfer_eisenstein.cc @@ -237,7 +237,8 @@ public: //! Computes the transfer function for k in Mpc/h by calling TFfit_onek inline double compute(double k, tf_type type) const { - return (type!=deltabc)? etf_.at_k(k) : 0.0; + if( type == theta_bc || type == delta_bc ) return 0.0; + return etf_.at_k(k); } inline double get_kmin(void) const @@ -332,6 +333,7 @@ public: inline double compute(double k, tf_type type) const { + if( type == theta_bc || type == delta_bc ) return 0.0; return etf_.at_k(k) * pow(1.0 + pow(m_WDMalpha * k, 2.0 * wdmnu_), -5.0 / wdmnu_); } @@ -384,6 +386,8 @@ public: double kkfs2 = kkfs * kkfs; double kkd2 = (k / kd_) * (k / kd_); + if( type == theta_bc || type == delta_bc ) return 0.0; + // in principle the Green et al. (2004) works only up to k/k_fs < 1 // the fit crosses zero at (k/k_fs)**2 = 3/2, we just zero it there... if (kkfs2 < 1.5)