2020-08-21 17:03:33 +02:00
|
|
|
// 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 <http://www.gnu.org/licenses/>.
|
2019-05-07 01:05:16 +02:00
|
|
|
#pragma once
|
|
|
|
|
|
|
|
#include <array>
|
2020-03-29 14:49:17 +02:00
|
|
|
#include <vec.hh>
|
2019-05-07 01:05:16 +02:00
|
|
|
|
|
|
|
#include <cosmology_parameters.hh>
|
2020-03-29 14:49:17 +02:00
|
|
|
#include <physical_constants.hh>
|
2019-05-07 01:05:16 +02:00
|
|
|
#include <transfer_function_plugin.hh>
|
2020-04-04 23:59:13 +02:00
|
|
|
#include <math/ode_integrate.hh>
|
2019-05-07 01:05:16 +02:00
|
|
|
#include <logger.hh>
|
|
|
|
|
2020-04-04 23:59:13 +02:00
|
|
|
#include <math/interpolate.hh>
|
2020-04-02 12:48:52 +02:00
|
|
|
|
2019-05-07 01:05:16 +02:00
|
|
|
#include <gsl/gsl_integration.h>
|
|
|
|
#include <gsl/gsl_errno.h>
|
|
|
|
|
2020-03-29 14:49:17 +02:00
|
|
|
namespace cosmology
|
|
|
|
{
|
|
|
|
|
2019-05-07 01:05:16 +02:00
|
|
|
/*!
|
2020-03-29 14:49:17 +02:00
|
|
|
* @class cosmology::calculator
|
2019-05-07 01:05:16 +02:00
|
|
|
* @brief provides functions to compute cosmological quantities
|
|
|
|
*
|
|
|
|
* This class provides member functions to compute cosmological quantities
|
2020-09-10 05:02:43 +02:00
|
|
|
* related to the Friedmann equations and linear perturbation theory, it also
|
|
|
|
* provides the functionality to work with back-scaled cosmological fields
|
2019-05-07 01:05:16 +02:00
|
|
|
*/
|
2020-03-29 14:49:17 +02:00
|
|
|
class calculator
|
2019-05-07 01:05:16 +02:00
|
|
|
{
|
2020-03-29 14:49:17 +02:00
|
|
|
public:
|
|
|
|
//! data structure to store cosmological parameters
|
|
|
|
cosmology::parameters cosmo_param_;
|
|
|
|
|
|
|
|
//! pointer to an instance of a transfer function plugin
|
|
|
|
std::unique_ptr<TransferFunction_plugin> transfer_function_;
|
|
|
|
|
2019-05-10 04:48:55 +02:00
|
|
|
private:
|
2020-04-04 01:24:05 +02:00
|
|
|
static constexpr double REL_PRECISION = 1e-10;
|
2022-05-23 14:53:43 +02:00
|
|
|
interpolated_function_1d<true,true,false> D_of_a_, f_of_a_, a_of_D_, E_of_a_, g_of_a_, Fa_of_a_, Fb_of_a_, dotFa_of_a_, dotFb_of_a_, Fc_of_a_, dotFc_of_a_; // toma
|
2020-05-02 21:02:24 +02:00
|
|
|
double Dnow_, Dplus_start_, Dplus_target_, astart_, atarget_;
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2020-09-11 01:25:45 +02:00
|
|
|
double m_n_s_, m_sqrtpnorm_;
|
|
|
|
|
2020-09-10 05:02:43 +02:00
|
|
|
//! 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
|
2019-05-24 11:32:15 +02:00
|
|
|
real_t integrate(double (*func)(double x, void *params), double a, double b, void *params) const
|
2019-05-07 01:05:16 +02:00
|
|
|
{
|
2020-09-10 05:02:43 +02:00
|
|
|
constexpr size_t wspace_size{100000};
|
|
|
|
|
|
|
|
double result{0.0};
|
|
|
|
double error{0.0};
|
|
|
|
|
2019-05-07 01:05:16 +02:00
|
|
|
gsl_function F;
|
|
|
|
F.function = func;
|
|
|
|
F.params = params;
|
2020-09-10 05:02:43 +02:00
|
|
|
|
|
|
|
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);
|
2019-05-07 01:05:16 +02:00
|
|
|
|
|
|
|
if (error / result > REL_PRECISION)
|
2020-04-04 20:27:51 +02:00
|
|
|
music::wlog << "no convergence in function 'integrate', rel. error=" << error / result << std::endl;
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2020-09-10 05:02:43 +02:00
|
|
|
return static_cast<real_t>(result);
|
2019-05-07 01:05:16 +02:00
|
|
|
}
|
|
|
|
|
2020-09-10 05:02:43 +02:00
|
|
|
//! 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
|
|
|
|
*/
|
2022-05-20 10:32:52 +02:00
|
|
|
void compute_growth( std::vector<double>& tab_a, std::vector<double>& tab_D, std::vector<double>& tab_f,
|
|
|
|
std::vector<double>& tab_E, std::vector<double>& tab_g,
|
2022-05-23 14:53:43 +02:00
|
|
|
std::vector<double>& tab_Fa, std::vector<double>& tab_dotFa,
|
|
|
|
std::vector<double>& tab_Fb, std::vector<double>& tab_dotFb,
|
|
|
|
std::vector<double>& tab_Fc, std::vector<double>& tab_dotFc
|
2022-05-20 10:32:52 +02:00
|
|
|
)
|
2020-03-29 14:49:17 +02:00
|
|
|
{
|
2022-05-23 14:53:43 +02:00
|
|
|
using v_t = vec_t<10, double>;
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2020-09-10 05:02:43 +02:00
|
|
|
// set ICs, very deep in radiation domination
|
2020-03-29 14:49:17 +02:00
|
|
|
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);
|
|
|
|
const double t0 = 1.0 / (a0 * H_of_a(a0));
|
|
|
|
|
2022-05-17 16:36:59 +02:00
|
|
|
// toma
|
2022-05-20 10:32:52 +02:00
|
|
|
// second order growth
|
2022-05-17 16:36:59 +02:00
|
|
|
// initialisation in EdS
|
2022-05-20 10:32:52 +02:00
|
|
|
const double E0 = -3.0/7.0*a0*a0;
|
|
|
|
const double Eprime0 = -12.0/7.0*std::pow(a0,3.0/2.0) ;
|
2022-05-17 16:36:59 +02:00
|
|
|
|
2022-05-20 10:32:52 +02:00
|
|
|
// third order growth
|
|
|
|
const double Fa0 = -1.0/3.0*a0*a0*a0 ;
|
|
|
|
const double Faprime0 = -2.0*std::pow(a0,5.0/2.0);
|
|
|
|
|
|
|
|
const double Fb0 = 10./21*a0*a0*a0;
|
|
|
|
const double Fbprime0 = 60/21*std::pow(a0,5.0/2.0);
|
|
|
|
|
2022-05-23 14:53:43 +02:00
|
|
|
const double Fc0 = -1.0/7.0*a0*a0*a0;
|
|
|
|
|
|
|
|
v_t y0({a0, D0, Dprime0, E0, Eprime0, Fa0, Faprime0, Fb0, Fbprime0, Fc0});
|
2020-03-29 14:49:17 +02:00
|
|
|
|
|
|
|
// set up integration
|
|
|
|
double dt = 1e-9;
|
|
|
|
double dtdid, dtnext;
|
|
|
|
const double amax = 2.0;
|
|
|
|
|
|
|
|
v_t yy(y0);
|
|
|
|
double t = t0;
|
|
|
|
const double eps = 1e-10;
|
|
|
|
|
2020-09-01 15:29:40 +02:00
|
|
|
const double Omega_m( cosmo_param_["Omega_m"] ), H0( cosmo_param_["H0"] );
|
|
|
|
|
2020-03-29 14:49:17 +02:00
|
|
|
while (yy[0] < amax)
|
|
|
|
{
|
|
|
|
// RHS of ODEs
|
|
|
|
auto rhs = [&](double t, v_t y) -> v_t {
|
|
|
|
auto a = y[0];
|
|
|
|
auto D = y[1];
|
|
|
|
auto Dprime = y[2];
|
2022-05-20 10:32:52 +02:00
|
|
|
|
|
|
|
auto E = y[3];
|
|
|
|
auto Eprime = y[4];
|
|
|
|
|
|
|
|
auto Fa = y[5];
|
|
|
|
auto Faprime = y[6];
|
|
|
|
|
|
|
|
auto Fb = y[7];
|
|
|
|
auto Fbprime = y[8];
|
|
|
|
|
2020-03-29 14:49:17 +02:00
|
|
|
v_t dy;
|
|
|
|
// da/dtau = a^2 H(a)
|
|
|
|
dy[0] = a * a * H_of_a(a);
|
|
|
|
// d D/dtau
|
|
|
|
dy[1] = Dprime;
|
|
|
|
// d^2 D / dtau^2
|
2020-09-01 15:29:40 +02:00
|
|
|
dy[2] = -a * H_of_a(a) * Dprime + 3.0 / 2.0 * Omega_m * std::pow(H0, 2) * D / a;
|
2022-05-17 16:36:59 +02:00
|
|
|
|
|
|
|
// toma
|
2022-05-20 10:32:52 +02:00
|
|
|
dy[3] = Eprime;
|
|
|
|
dy[4] = -a * H_of_a(a) * Eprime + 3.0 / 2.0 * Omega_m * std::pow(H0, 2) * (E-D*D) / a;
|
|
|
|
|
|
|
|
dy[5] = Faprime;
|
|
|
|
dy[6] = -a * H_of_a(a) * Faprime + 3.0 / 2.0 * Omega_m * std::pow(H0, 2) * (Fa-2.0*D*D*D) / a;
|
|
|
|
|
|
|
|
dy[7] = Fbprime;
|
|
|
|
dy[8] = -a * H_of_a(a) * Fbprime + 3.0 / 2.0 * Omega_m * std::pow(H0, 2) * (Fb + 2.0*D*D*D - 2.0*E*D) / a;
|
2022-05-23 14:53:43 +02:00
|
|
|
|
|
|
|
dy[9] = D*Eprime-E*Dprime;
|
|
|
|
|
2020-03-29 14:49:17 +02:00
|
|
|
return dy;
|
|
|
|
};
|
|
|
|
|
|
|
|
// scale by predicted value to get approx. constant fractional errors
|
|
|
|
v_t yyscale = yy.abs() + dt * rhs(t, yy).abs();
|
2022-05-20 10:32:52 +02:00
|
|
|
|
2020-03-29 14:49:17 +02:00
|
|
|
|
|
|
|
// call integrator
|
|
|
|
ode_integrate::rk_step_qs(dt, t, yy, yyscale, rhs, eps, dtdid, dtnext);
|
|
|
|
|
2022-05-20 10:32:52 +02:00
|
|
|
|
2020-04-02 12:48:52 +02:00
|
|
|
tab_a.push_back(yy[0]);
|
|
|
|
tab_D.push_back(yy[1]);
|
2020-09-10 05:02:43 +02:00
|
|
|
tab_f.push_back(yy[2]); // temporarily store D' in table
|
2022-05-20 10:32:52 +02:00
|
|
|
tab_E.push_back(yy[3]); // toma
|
|
|
|
tab_g.push_back(yy[4]); // temporarily store E' in table
|
|
|
|
|
|
|
|
tab_Fa.push_back(yy[5]); //
|
2022-05-23 14:53:43 +02:00
|
|
|
tab_dotFa.push_back(yy[6]); // temporarily store Fa' in table
|
2022-05-20 10:32:52 +02:00
|
|
|
|
|
|
|
tab_Fb.push_back(yy[7]); //
|
2022-05-23 14:53:43 +02:00
|
|
|
tab_dotFb.push_back(yy[8]); // temporarily store Fb' in table
|
|
|
|
|
|
|
|
tab_Fc.push_back(yy[9]); //
|
|
|
|
tab_dotFc.push_back(yy[9]); // temporarily store un-important random value
|
2020-03-29 14:49:17 +02:00
|
|
|
|
|
|
|
dt = dtnext;
|
|
|
|
}
|
|
|
|
|
2022-05-17 16:36:59 +02:00
|
|
|
std::ofstream output_file;
|
|
|
|
output_file.open("GrowthFactors.txt");
|
2020-03-29 14:49:17 +02:00
|
|
|
// compute f, before we stored here D'
|
2022-05-23 14:53:43 +02:00
|
|
|
output_file << "#" << "a" <<" "<< "D" << " " << "f" << " " << "E" << " " << "g" << " " << "Fa" << " " << "dotFa" << " " << "Fb" << " " << "dotFb" << " " <<"Fc" << " " << "dotFc" <<"\n";
|
2020-04-02 12:48:52 +02:00
|
|
|
for (size_t i = 0; i < tab_a.size(); ++i)
|
2020-03-29 14:49:17 +02:00
|
|
|
{
|
2022-05-20 10:32:52 +02:00
|
|
|
|
2022-05-23 14:53:43 +02:00
|
|
|
//tab_D[i] = tab_D[i];
|
|
|
|
//tab_E[i] = tab_E[i]; // toma
|
|
|
|
//tab_Fa[i] = tab_Fa[i]; // toma
|
|
|
|
//tab_Fb[i] = tab_Fb[i]; // toma
|
|
|
|
//tab_Fc[i] = tab_Fc[i]; // toma
|
|
|
|
//tab_a[i] = tab_a[i];
|
|
|
|
tab_dotFc[i] = (tab_D[i]*tab_g[i] - tab_E[i]*tab_f[i]); // toma
|
|
|
|
tab_f[i] = tab_f[i] / (tab_a[i] * H_of_a(tab_a[i]) * tab_D[i]);
|
2022-05-20 10:32:52 +02:00
|
|
|
|
2022-05-17 16:36:59 +02:00
|
|
|
|
|
|
|
//toma
|
2022-05-23 14:53:43 +02:00
|
|
|
output_file << tab_a[i] <<" "<< tab_D[i] << " " << tab_f[i] << " " << tab_E[i] << " " << tab_g[i] << " " << tab_Fa[i] << " " << tab_dotFa[i] << " " << tab_Fb[i] << " " << tab_dotFb[i] << " " << tab_Fc[i] << " " << tab_dotFc[i] << "\n";
|
2020-03-29 14:49:17 +02:00
|
|
|
}
|
2022-05-17 16:36:59 +02:00
|
|
|
output_file.close();
|
2020-03-29 14:49:17 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
public:
|
2020-09-01 11:04:26 +02:00
|
|
|
|
2020-04-04 01:24:05 +02:00
|
|
|
calculator() = delete;
|
2020-09-01 11:04:26 +02:00
|
|
|
|
2020-04-04 01:24:05 +02:00
|
|
|
calculator(const calculator& c) = delete;
|
2020-09-01 11:04:26 +02:00
|
|
|
|
2019-05-07 01:05:16 +02:00
|
|
|
//! constructor for a cosmology calculator object
|
|
|
|
/*!
|
|
|
|
* @param acosmo a cosmological parameters structure
|
|
|
|
* @param pTransferFunction pointer to an instance of a transfer function object
|
|
|
|
*/
|
2020-04-04 20:55:24 +02:00
|
|
|
explicit calculator(config_file &cf)
|
2020-05-02 21:02:24 +02:00
|
|
|
: cosmo_param_(cf), astart_( 1.0/(1.0+cf.get_value<double>("setup","zstart")) ),
|
2020-09-07 22:55:39 +02:00
|
|
|
atarget_( 1.0/(1.0+cf.get_value_safe<double>("cosmology","ztarget",0.0)) )
|
2020-03-29 14:49:17 +02:00
|
|
|
{
|
2020-04-02 12:48:52 +02:00
|
|
|
// pre-compute growth factors and store for interpolation
|
2022-05-23 14:53:43 +02:00
|
|
|
std::vector<double> tab_a, tab_D, tab_f, tab_E, tab_g, tab_Fa, tab_dotFa, tab_Fb, tab_dotFb, tab_Fc, tab_dotFc; // toma
|
|
|
|
this->compute_growth(tab_a, tab_D, tab_f, tab_E, tab_g, tab_Fa, tab_dotFa, tab_Fb, tab_dotFb, tab_Fc, tab_dotFc); // toma
|
2020-04-02 12:48:52 +02:00
|
|
|
D_of_a_.set_data(tab_a,tab_D);
|
|
|
|
f_of_a_.set_data(tab_a,tab_f);
|
|
|
|
a_of_D_.set_data(tab_D,tab_a);
|
2022-05-17 16:36:59 +02:00
|
|
|
|
|
|
|
// toma
|
2022-05-20 10:32:52 +02:00
|
|
|
E_of_a_.set_data(tab_a,tab_E);
|
|
|
|
g_of_a_.set_data(tab_a,tab_g);
|
|
|
|
|
|
|
|
|
|
|
|
Fa_of_a_.set_data(tab_a,tab_Fa);
|
|
|
|
Fb_of_a_.set_data(tab_a,tab_Fb);
|
2022-05-23 14:53:43 +02:00
|
|
|
Fc_of_a_.set_data(tab_a,tab_Fc);
|
|
|
|
dotFa_of_a_.set_data(tab_a,tab_dotFa);
|
|
|
|
dotFb_of_a_.set_data(tab_a,tab_dotFb);
|
|
|
|
dotFc_of_a_.set_data(tab_a,tab_dotFc);
|
2022-05-17 16:36:59 +02:00
|
|
|
|
2020-04-02 12:48:52 +02:00
|
|
|
Dnow_ = D_of_a_(1.0);
|
|
|
|
|
2020-09-10 05:02:43 +02:00
|
|
|
Dplus_start_ = D_of_a_( astart_ ) / Dnow_;
|
2020-05-02 21:02:24 +02:00
|
|
|
Dplus_target_ = D_of_a_( atarget_ ) / Dnow_;
|
2020-04-02 19:25:54 +02:00
|
|
|
|
2020-08-16 20:53:32 +02:00
|
|
|
music::ilog << "Linear growth factors: D+_target = " << Dplus_target_ << ", D+_start = " << Dplus_start_ << std::endl;
|
2020-05-22 18:44:15 +02:00
|
|
|
|
2020-04-02 12:48:52 +02:00
|
|
|
// set up transfer functions and compute normalisation
|
2020-09-01 15:29:40 +02:00
|
|
|
transfer_function_ = std::move(select_TransferFunction_plugin(cf, cosmo_param_));
|
2019-05-24 08:13:58 +02:00
|
|
|
transfer_function_->intialise();
|
2020-09-07 22:55:39 +02:00
|
|
|
if( !transfer_function_->tf_isnormalised_ ){
|
2020-09-01 15:29:40 +02:00
|
|
|
cosmo_param_.set("pnorm", this->compute_pnorm_from_sigma8() );
|
2020-09-07 22:55:39 +02:00
|
|
|
}else{
|
2020-09-01 15:29:40 +02:00
|
|
|
cosmo_param_.set("pnorm", 1.0/Dplus_target_/Dplus_target_);
|
2020-04-04 01:24:05 +02:00
|
|
|
auto sigma8 = this->compute_sigma8();
|
2020-04-04 20:27:51 +02:00
|
|
|
music::ilog << "Measured sigma_8 for given PS normalisation is " << sigma8 << std::endl;
|
2020-04-02 12:48:52 +02:00
|
|
|
}
|
2020-09-01 15:29:40 +02:00
|
|
|
cosmo_param_.set("sqrtpnorm", std::sqrt(cosmo_param_["pnorm"]));
|
2020-04-02 12:48:52 +02:00
|
|
|
|
2020-04-04 20:27:51 +02:00
|
|
|
music::ilog << std::setw(32) << std::left << "TF supports distinct CDM+baryons"
|
2020-03-29 14:49:17 +02:00
|
|
|
<< " : " << (transfer_function_->tf_is_distinct() ? "yes" : "no") << std::endl;
|
2020-04-04 20:27:51 +02:00
|
|
|
music::ilog << std::setw(32) << std::left << "TF maximum wave number"
|
2020-03-29 14:49:17 +02:00
|
|
|
<< " : " << transfer_function_->get_kmax() << " h/Mpc" << std::endl;
|
2020-09-11 01:25:45 +02:00
|
|
|
|
|
|
|
m_n_s_ = cosmo_param_["n_s"];
|
|
|
|
m_sqrtpnorm_ = cosmo_param_["sqrtpnorm"];
|
2020-03-29 14:49:17 +02:00
|
|
|
}
|
|
|
|
|
2020-09-10 05:02:43 +02:00
|
|
|
~calculator() { }
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-05-24 11:32:15 +02:00
|
|
|
//! Write out a correctly scaled power spectrum at time a
|
2020-03-29 14:49:17 +02:00
|
|
|
void write_powerspectrum(real_t a, std::string fname) const
|
2019-05-24 11:32:15 +02:00
|
|
|
{
|
2020-04-02 12:48:52 +02:00
|
|
|
// const real_t Dplus0 = this->get_growth_factor(a);
|
2019-05-24 11:32:15 +02:00
|
|
|
|
2020-03-29 14:49:17 +02:00
|
|
|
if (CONFIG::MPI_task_rank == 0)
|
2019-05-24 11:32:15 +02:00
|
|
|
{
|
2020-03-29 14:49:17 +02:00
|
|
|
double kmin = std::max(1e-4, transfer_function_->get_kmin());
|
2019-05-29 05:16:56 +02:00
|
|
|
|
2019-05-24 11:32:15 +02:00
|
|
|
// write power spectrum to a file
|
|
|
|
std::ofstream ofs(fname.c_str());
|
2020-03-29 14:49:17 +02:00
|
|
|
std::stringstream ss;
|
2020-04-02 12:48:52 +02:00
|
|
|
ss << " ,ap=" << a << "";
|
2019-05-24 11:32:15 +02:00
|
|
|
ofs << "# " << std::setw(18) << "k [h/Mpc]"
|
2020-04-02 12:48:52 +02:00
|
|
|
<< std::setw(20) << ("P_dtot(k,a=ap)")
|
|
|
|
<< std::setw(20) << ("P_dcdm(k,a=ap)")
|
|
|
|
<< std::setw(20) << ("P_dbar(k,a=ap)")
|
|
|
|
<< std::setw(20) << ("P_tcdm(k,a=ap)")
|
|
|
|
<< std::setw(20) << ("P_tbar(k,a=ap)")
|
|
|
|
<< std::setw(20) << ("P_dtot(k,a=1)")
|
|
|
|
<< std::setw(20) << ("P_dcdm(k,a=1)")
|
|
|
|
<< std::setw(20) << ("P_dbar(k,a=1)")
|
|
|
|
<< std::setw(20) << ("P_tcdm(k,a=1)")
|
|
|
|
<< std::setw(20) << ("P_tbar(k,a=1)")
|
2020-03-29 14:49:17 +02:00
|
|
|
<< std::endl;
|
2020-09-07 23:55:56 +02:00
|
|
|
for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01)
|
2020-03-29 14:49:17 +02:00
|
|
|
{
|
2019-05-24 11:32:15 +02:00
|
|
|
ofs << std::setw(20) << std::setprecision(10) << k
|
2020-09-10 05:02:43 +02:00
|
|
|
<< 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)
|
2019-05-24 11:32:15 +02:00
|
|
|
<< std::endl;
|
|
|
|
}
|
|
|
|
}
|
2020-04-04 20:27:51 +02:00
|
|
|
music::ilog << "Wrote power spectrum at a=" << a << " to file \'" << fname << "\'" << std::endl;
|
2019-05-24 11:32:15 +02:00
|
|
|
}
|
|
|
|
|
2020-08-14 19:09:29 +02:00
|
|
|
//! Write out a correctly scaled power spectrum at starting time
|
|
|
|
void write_transfer( std::string fname ) const
|
|
|
|
{
|
|
|
|
// const real_t Dplus0 = this->get_growth_factor(a);
|
|
|
|
|
|
|
|
if (CONFIG::MPI_task_rank == 0)
|
|
|
|
{
|
|
|
|
double kmin = std::max(1e-4, transfer_function_->get_kmin());
|
|
|
|
|
|
|
|
// write power spectrum to a file
|
|
|
|
std::ofstream ofs(fname.c_str());
|
|
|
|
std::stringstream ss;
|
|
|
|
ss << " ,ap=" << astart_ << "";
|
|
|
|
ofs << "# " << std::setw(18) << "k [h/Mpc]"
|
|
|
|
<< std::setw(20) << ("delta_c(k,a=ap)")
|
|
|
|
<< std::setw(20) << ("delta_b(k,a=ap)")
|
|
|
|
<< std::setw(20) << ("delta_m(k,a=ap)")
|
|
|
|
<< std::setw(20) << ("delta_bc(k,a=ap)")
|
|
|
|
<< std::endl;
|
2020-09-01 15:29:40 +02:00
|
|
|
double fb = cosmo_param_["f_b"], fc = cosmo_param_["f_c"];
|
2020-09-07 23:55:56 +02:00
|
|
|
for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01)
|
2020-08-14 19:09:29 +02:00
|
|
|
{
|
2020-09-10 05:02:43 +02:00
|
|
|
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;
|
|
|
|
|
2020-08-14 19:09:29 +02:00
|
|
|
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
|
2020-09-10 05:02:43 +02:00
|
|
|
<< 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 )
|
2020-08-14 19:09:29 +02:00
|
|
|
<< std::endl;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
music::ilog << "Wrote input transfer functions at a=" << astart_ << " to file \'" << fname << "\'" << std::endl;
|
|
|
|
}
|
|
|
|
|
2020-04-02 12:48:52 +02:00
|
|
|
const cosmology::parameters &get_parameters(void) const noexcept
|
2019-05-10 04:48:55 +02:00
|
|
|
{
|
2019-05-07 01:05:16 +02:00
|
|
|
return cosmo_param_;
|
|
|
|
}
|
|
|
|
|
2020-03-29 14:49:17 +02:00
|
|
|
//! return the value of the Hubble function H(a) = dloga/dt
|
2020-04-02 12:48:52 +02:00
|
|
|
inline double H_of_a(double a) const noexcept
|
2019-05-07 01:05:16 +02:00
|
|
|
{
|
2020-03-29 14:49:17 +02:00
|
|
|
double HH2 = 0.0;
|
2022-05-23 14:53:43 +02:00
|
|
|
HH2 += cosmo_param_["Omega_r"] / (a * a * a * a);
|
2020-09-01 15:29:40 +02:00
|
|
|
HH2 += cosmo_param_["Omega_m"] / (a * a * a);
|
|
|
|
HH2 += cosmo_param_["Omega_k"] / (a * a);
|
|
|
|
HH2 += cosmo_param_["Omega_DE"] * std::pow(a, -3. * (1. + cosmo_param_["w_0"] + cosmo_param_["w_a"])) * exp(-3. * (1.0 - a) * cosmo_param_["w_a"]);
|
|
|
|
return cosmo_param_["H0"] * std::sqrt(HH2);
|
2019-05-07 01:05:16 +02:00
|
|
|
}
|
|
|
|
|
2020-03-29 14:49:17 +02:00
|
|
|
//! Computes the linear theory growth factor D+, normalised to D+(a=1)=1
|
2020-04-02 12:48:52 +02:00
|
|
|
real_t get_growth_factor(real_t a) const noexcept
|
2019-11-15 23:19:57 +01:00
|
|
|
{
|
2020-04-02 12:48:52 +02:00
|
|
|
return D_of_a_(a) / Dnow_;
|
|
|
|
}
|
|
|
|
|
|
|
|
//! Computes the inverse of get_growth_factor
|
|
|
|
real_t get_a( real_t Dplus ) const noexcept
|
|
|
|
{
|
|
|
|
return a_of_D_( Dplus * Dnow_ );
|
2019-11-15 23:19:57 +01:00
|
|
|
}
|
|
|
|
|
2020-03-29 14:49:17 +02:00
|
|
|
//! Computes the linear theory growth rate f
|
|
|
|
/*! Function computes (by interpolating on precalculated table)
|
|
|
|
* f = dlog D+ / dlog a
|
|
|
|
*/
|
2020-04-02 12:48:52 +02:00
|
|
|
real_t get_f(real_t a) const noexcept
|
2019-05-07 01:05:16 +02:00
|
|
|
{
|
2020-04-02 12:48:52 +02:00
|
|
|
return f_of_a_(a);
|
2019-05-07 01:05:16 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
//! Compute the factor relating particle displacement and velocity
|
|
|
|
/*! Function computes
|
2020-03-29 14:49:17 +02:00
|
|
|
* vfac = a * (H(a)/h) * dlogD+ / dlog a
|
|
|
|
*/
|
2020-04-02 12:48:52 +02:00
|
|
|
real_t get_vfact(real_t a) const noexcept
|
2019-05-07 01:05:16 +02:00
|
|
|
{
|
2020-09-01 15:29:40 +02:00
|
|
|
return f_of_a_(a) * a * H_of_a(a) / cosmo_param_["h"];
|
2019-05-07 01:05:16 +02:00
|
|
|
}
|
|
|
|
|
2022-05-17 16:36:59 +02:00
|
|
|
|
|
|
|
// toma
|
2022-05-20 10:32:52 +02:00
|
|
|
//! get E
|
2022-05-17 16:36:59 +02:00
|
|
|
real_t get_2growth_factor(real_t a) const noexcept
|
|
|
|
{
|
2022-05-20 10:32:52 +02:00
|
|
|
return E_of_a_(a);
|
2022-05-17 16:36:59 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
// toma
|
|
|
|
//! compute the factors w
|
2022-05-20 10:32:52 +02:00
|
|
|
real_t get_g(real_t a) const noexcept
|
2022-05-17 16:36:59 +02:00
|
|
|
{
|
2022-05-20 10:32:52 +02:00
|
|
|
return g_of_a_(a);
|
2022-05-17 16:36:59 +02:00
|
|
|
}
|
|
|
|
|
2022-05-20 10:32:52 +02:00
|
|
|
real_t get_3growthA_factor(real_t a) const noexcept
|
|
|
|
{
|
|
|
|
return Fa_of_a_(a);
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
real_t get_3growthB_factor(real_t a) const noexcept
|
|
|
|
{
|
|
|
|
return Fb_of_a_(a);
|
|
|
|
}
|
|
|
|
|
2022-05-23 14:53:43 +02:00
|
|
|
real_t get_3growthC_factor(real_t a) const noexcept
|
2022-05-20 10:32:52 +02:00
|
|
|
{
|
2022-05-23 14:53:43 +02:00
|
|
|
return Fc_of_a_(a);
|
2022-05-20 10:32:52 +02:00
|
|
|
}
|
|
|
|
|
2022-05-23 14:53:43 +02:00
|
|
|
|
|
|
|
real_t get_dotFa(real_t a) const noexcept
|
2022-05-20 10:32:52 +02:00
|
|
|
{
|
2022-05-23 14:53:43 +02:00
|
|
|
return dotFa_of_a_(a);
|
2022-05-20 10:32:52 +02:00
|
|
|
}
|
|
|
|
|
2022-05-23 14:53:43 +02:00
|
|
|
real_t get_dotFb(real_t a) const noexcept
|
|
|
|
{
|
|
|
|
return dotFb_of_a_(a);
|
|
|
|
}
|
|
|
|
|
|
|
|
real_t get_dotFc(real_t a) const noexcept
|
|
|
|
{
|
|
|
|
return dotFc_of_a_(a);
|
|
|
|
}
|
|
|
|
|
|
|
|
|
2022-05-20 10:32:52 +02:00
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
2019-05-07 01:05:16 +02:00
|
|
|
//! Integrand for the sigma_8 normalization of the power spectrum
|
|
|
|
/*! Returns the value of the primordial power spectrum multiplied with
|
|
|
|
the transfer function and the window function of 8 Mpc/h at wave number k */
|
|
|
|
static double dSigma8(double k, void *pParams)
|
|
|
|
{
|
2020-03-29 14:49:17 +02:00
|
|
|
cosmology::calculator *pcc = reinterpret_cast<cosmology::calculator *>(pParams);
|
|
|
|
|
2020-09-07 22:55:39 +02:00
|
|
|
const double x = k * 8.0;
|
2020-09-07 23:55:56 +02:00
|
|
|
const double w = (x < 0.001)? 1.0-0.1*x*x : 3.0 * (std::sin(x) - x * std::cos(x)) / (x * x * x);
|
2020-09-10 05:02:43 +02:00
|
|
|
|
2020-09-01 15:29:40 +02:00
|
|
|
static double nspect = (double)pcc->cosmo_param_["n_s"];
|
2020-09-10 05:02:43 +02:00
|
|
|
double tf = pcc->transfer_function_->compute(k, delta_matter);
|
2019-05-07 01:05:16 +02:00
|
|
|
|
|
|
|
//... 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;
|
|
|
|
}
|
|
|
|
|
|
|
|
//! Integrand for the sigma_8 normalization of the power spectrum
|
|
|
|
/*! Returns the value of the primordial power spectrum multiplied with
|
|
|
|
the transfer function and the window function of 8 Mpc/h at wave number k */
|
|
|
|
static double dSigma8_0(double k, void *pParams)
|
|
|
|
{
|
2020-03-29 14:49:17 +02:00
|
|
|
cosmology::calculator *pcc = reinterpret_cast<cosmology::calculator *>(pParams);
|
|
|
|
|
2020-09-07 22:55:39 +02:00
|
|
|
const double x = k * 8.0;
|
2020-09-07 23:55:56 +02:00
|
|
|
const double w = (x < 0.001)? 1.0-0.1*x*x : 3.0 * (std::sin(x) - x * std::cos(x)) / (x * x * x);
|
2020-09-07 22:55:39 +02:00
|
|
|
|
2020-09-01 15:29:40 +02:00
|
|
|
static double nspect = static_cast<double>(pcc->cosmo_param_["n_s"]);
|
2020-09-10 05:02:43 +02:00
|
|
|
double tf = pcc->transfer_function_->compute(k, delta_matter0);
|
2019-05-07 01:05:16 +02:00
|
|
|
|
|
|
|
//... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1
|
2020-05-22 18:44:15 +02:00
|
|
|
return k * k * w * w * std::pow(k, nspect) * tf * tf;
|
2019-05-07 01:05:16 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
//! Computes the amplitude of a mode from the power spectrum
|
2020-09-10 05:02:43 +02:00
|
|
|
/*! 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
|
2019-05-07 01:05:16 +02:00
|
|
|
* @param k wave number at which to evaluate
|
2020-09-10 05:02:43 +02:00
|
|
|
* @param type one of the species: {delta,theta}_{matter,cdm,baryon,neutrino}
|
2019-05-07 01:05:16 +02:00
|
|
|
*/
|
2020-05-22 18:44:15 +02:00
|
|
|
inline real_t get_amplitude( const real_t k, const tf_type type) const
|
|
|
|
{
|
2020-09-11 01:25:45 +02:00
|
|
|
return std::pow(k, 0.5 * m_n_s_) * transfer_function_->compute(k, type) * m_sqrtpnorm_;
|
2020-05-22 18:44:15 +02:00
|
|
|
}
|
|
|
|
|
2020-09-10 05:02:43 +02:00
|
|
|
//! 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
|
2019-05-07 01:05:16 +02:00
|
|
|
{
|
2020-09-10 05:02:43 +02:00
|
|
|
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);
|
2020-05-23 15:12:51 +02:00
|
|
|
// need to multiply with Dplus_target since sqrtpnorm rescales like that
|
2020-09-11 01:25:45 +02:00
|
|
|
return std::pow(k, 0.5 * m_n_s_) * dbc * (m_sqrtpnorm_ * Dplus_target_);
|
2020-05-23 16:14:37 +02:00
|
|
|
}
|
2020-09-10 05:02:43 +02:00
|
|
|
|
|
|
|
//! 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
|
2020-05-23 16:14:37 +02:00
|
|
|
{
|
2020-09-10 05:02:43 +02:00
|
|
|
const real_t Dratio = Dplus_target_ / Dplus_start_;
|
|
|
|
const real_t tbc = transfer_function_->compute(k, theta_bc) * std::sqrt(Dratio);
|
2020-05-23 16:14:37 +02:00
|
|
|
// need to multiply with Dplus_target since sqrtpnorm rescales like that
|
2020-09-11 01:25:45 +02:00
|
|
|
return withvbc ? std::pow(k, 0.5 * m_n_s_) * tbc * (m_sqrtpnorm_ * Dplus_target_) : 0.0;
|
2019-05-07 01:05:16 +02:00
|
|
|
}
|
|
|
|
|
2020-09-10 05:02:43 +02:00
|
|
|
|
2019-05-07 01:05:16 +02:00
|
|
|
//! Computes the normalization for the power spectrum
|
|
|
|
/*!
|
|
|
|
* integrates the power spectrum to fix the normalization to that given
|
|
|
|
* by the sigma_8 parameter
|
|
|
|
*/
|
2020-04-02 12:48:52 +02:00
|
|
|
real_t compute_sigma8(void)
|
2019-05-07 01:05:16 +02:00
|
|
|
{
|
|
|
|
real_t sigma0, kmin, kmax;
|
2019-05-24 08:13:58 +02:00
|
|
|
kmax = transfer_function_->get_kmax();
|
|
|
|
kmin = transfer_function_->get_kmin();
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2019-05-24 08:13:58 +02:00
|
|
|
if (!transfer_function_->tf_has_total0())
|
2020-05-22 18:44:15 +02:00
|
|
|
sigma0 = 4.0 * M_PI * integrate(&dSigma8, static_cast<double>(kmin), static_cast<double>(kmax), this);
|
2020-04-02 12:48:52 +02:00
|
|
|
else{
|
2020-05-22 18:44:15 +02:00
|
|
|
sigma0 = 4.0 * M_PI * integrate(&dSigma8_0, static_cast<double>(kmin), static_cast<double>(kmax), this);
|
2020-04-02 12:48:52 +02:00
|
|
|
}
|
2019-05-07 01:05:16 +02:00
|
|
|
|
2020-04-02 12:48:52 +02:00
|
|
|
return std::sqrt(sigma0);
|
|
|
|
}
|
|
|
|
|
|
|
|
//! Computes the normalization for the power spectrum
|
|
|
|
/*!
|
|
|
|
* integrates the power spectrum to fix the normalization to that given
|
|
|
|
* by the sigma_8 parameter
|
|
|
|
*/
|
|
|
|
real_t compute_pnorm_from_sigma8(void)
|
|
|
|
{
|
|
|
|
auto measured_sigma8 = this->compute_sigma8();
|
2020-09-01 15:29:40 +02:00
|
|
|
return cosmo_param_["sigma_8"] * cosmo_param_["sigma_8"] / (measured_sigma8 * measured_sigma8);
|
2019-05-07 01:05:16 +02:00
|
|
|
}
|
|
|
|
};
|
|
|
|
|
|
|
|
//! compute the jeans sound speed
|
|
|
|
/*! given a density in g/cm^-3 and a mass in g it gives back the sound
|
|
|
|
* speed in cm/s for which the input mass is equal to the jeans mass
|
|
|
|
* @param rho density
|
|
|
|
* @param mass mass scale
|
|
|
|
* @returns jeans sound speed
|
|
|
|
*/
|
2020-09-01 15:29:40 +02:00
|
|
|
// inline double jeans_sound_speed(double rho, double mass)
|
|
|
|
// {
|
|
|
|
// const double G = 6.67e-8;
|
|
|
|
// return pow(6.0 * mass / M_PI * std::sqrt(rho) * std::pow(G, 1.5), 1.0 / 3.0);
|
|
|
|
// }
|
2020-03-29 14:49:17 +02:00
|
|
|
|
2022-05-17 16:36:59 +02:00
|
|
|
} // namespace cosmology
|