mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
changed CAMB plugin to use new interpolation interface rather than directly GSL
This commit is contained in:
Oliver Hahn
parent
bf6dd7e0ee
commit
e72547f81b
3 changed files with 105 additions and 293 deletions
|
|
@ -29,6 +29,7 @@ public:
|
|||
interpolated_function_1d(const std::vector<double> &data_x, const std::vector<double> &data_y)
|
||||
: isinit_(false)
|
||||
{
|
||||
static_assert(!(logx & periodic),"Class \'interpolated_function_1d\' cannot both be periodic and logarithmic in x!");
|
||||
this->set_data( data_x, data_y );
|
||||
}
|
||||
|
||||
|
|
@ -44,7 +45,6 @@ public:
|
|||
|
||||
assert(data_x_.size() == data_y_.size());
|
||||
assert(data_x_.size() > 5);
|
||||
assert(!(logx & periodic));
|
||||
|
||||
if (logx) for (auto &d : data_x_) d = std::log(d);
|
||||
if (logy) for (auto &d : data_y_) d = std::log(d);
|
||||
|
|
@ -61,8 +61,8 @@ public:
|
|||
double operator()(double x) const noexcept
|
||||
{
|
||||
assert( isinit_ && !(logx&&x<=0.0) );
|
||||
double xa = logx ? std::log(x) : x;
|
||||
double y(gsl_spline_eval(gsl_sp_, xa, gsl_ia_));
|
||||
const double xa = logx ? std::log(x) : x;
|
||||
const double y(gsl_spline_eval(gsl_sp_, xa, gsl_ia_));
|
||||
return logy ? std::exp(y) : y;
|
||||
}
|
||||
};
|
||||
|
|
@ -15,14 +15,9 @@
|
|||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
#include <gsl/gsl_errno.h>
|
||||
#include <gsl/gsl_spline.h>
|
||||
|
||||
#include <vector>
|
||||
|
||||
#include <transfer_function_plugin.hh>
|
||||
|
||||
const double tiny = 1e-30;
|
||||
#include <math/interpolate.hh>
|
||||
|
||||
class transfer_CAMB_file_plugin : public TransferFunction_plugin
|
||||
{
|
||||
|
|
@ -31,49 +26,31 @@ private:
|
|||
|
||||
using TransferFunction_plugin::cosmo_params_;
|
||||
|
||||
std::string m_filename_Pk, m_filename_Tk;
|
||||
std::vector<double> m_tab_k, m_tab_Tk_tot, m_tab_Tk_cdm, m_tab_Tk_baryon;
|
||||
std::vector<double> m_tab_Tvk_tot, m_tab_Tvk_cdm, m_tab_Tvk_baryon;
|
||||
gsl_interp_accel *acc_tot, *acc_cdm, *acc_baryon;
|
||||
gsl_interp_accel *acc_vtot, *acc_vcdm, *acc_vbaryon;
|
||||
gsl_spline *spline_tot, *spline_cdm, *spline_baryon;
|
||||
gsl_spline *spline_vtot, *spline_vcdm, *spline_vbaryon;
|
||||
interpolated_function_1d<true, true, false> delta_c_, delta_b_, delta_n_, delta_m_, theta_c_, theta_b_, theta_n_, theta_m_;
|
||||
|
||||
double m_kmin, m_kmax;
|
||||
unsigned m_nlines;
|
||||
|
||||
bool m_linbaryoninterp;
|
||||
// bool m_linbaryoninterp;
|
||||
|
||||
void read_table(void)
|
||||
void read_table( const std::string& filename )
|
||||
{
|
||||
|
||||
m_nlines = 0;
|
||||
m_linbaryoninterp = false;
|
||||
size_t nlines{0};
|
||||
// m_linbaryoninterp = false;
|
||||
|
||||
#ifdef WITH_MPI
|
||||
if (MPI::COMM_WORLD.Get_rank() == 0)
|
||||
std::vector<double> k, dc, tc, db, tb, dn, tn, dm, tm;
|
||||
|
||||
if( CONFIG::MPI_task_rank == 0 )
|
||||
{
|
||||
#endif
|
||||
music::ilog.Print("Reading tabulated transfer function data from file \n \'%s\'", m_filename_Tk.c_str());
|
||||
|
||||
music::ilog << "Reading tabulated transfer function data from file:" << std::endl
|
||||
<< " \'" << filename << "\'" << std::endl;
|
||||
|
||||
std::string line;
|
||||
std::ifstream ifs(m_filename_Tk.c_str());
|
||||
std::ifstream ifs(filename.c_str());
|
||||
|
||||
if (!ifs.good())
|
||||
throw std::runtime_error("Could not find transfer function file \'" + m_filename_Tk + "\'");
|
||||
|
||||
m_tab_k.clear();
|
||||
m_tab_Tk_tot.clear();
|
||||
m_tab_Tk_cdm.clear();
|
||||
m_tab_Tk_baryon.clear();
|
||||
m_tab_Tvk_tot.clear();
|
||||
m_tab_Tvk_cdm.clear(); //>[150609SH: add]
|
||||
m_tab_Tvk_baryon.clear(); //>[150609SH: add]
|
||||
|
||||
m_kmin = 1e30;
|
||||
m_kmax = -1e30;
|
||||
std::ofstream ofs("dump_transfer.txt");
|
||||
|
||||
throw std::runtime_error("Could not find transfer function file \'" + filename + "\'");
|
||||
|
||||
while (!ifs.eof())
|
||||
{
|
||||
getline(ifs, line);
|
||||
|
|
@ -86,100 +63,95 @@ private:
|
|||
|
||||
std::stringstream ss(line);
|
||||
|
||||
double k, Tkc, Tkb, Tktot, Tkvtot, Tkvc, Tkvb, dummy;
|
||||
double Tk, Tdc, Tdb, Tdn, Tdm, Tvb, Tvc, dummy;
|
||||
|
||||
ss >> k;
|
||||
ss >> Tkc; // cdm
|
||||
ss >> Tkb; // baryon
|
||||
ss >> Tk; // k
|
||||
ss >> Tdc; // cdm
|
||||
ss >> Tdb; // baryon
|
||||
ss >> dummy; // photon
|
||||
ss >> dummy; // nu
|
||||
ss >> dummy; // mass_nu
|
||||
ss >> Tktot; // total
|
||||
ss >> Tdn; // mass_nu
|
||||
ss >> Tdm; // total
|
||||
ss >> dummy; // no_nu
|
||||
ss >> dummy; // total_de
|
||||
ss >> dummy; // Weyl
|
||||
ss >> Tkvc; // v_cdm
|
||||
ss >> Tkvb; // v_b
|
||||
ss >> Tvc; // v_cdm
|
||||
ss >> Tvb; // v_b
|
||||
ss >> dummy; // v_b-v_cdm
|
||||
|
||||
if (ss.bad() || ss.fail())
|
||||
{
|
||||
music::elog.Print("error reading the transfer function file (corrupt or not in expected format)!");
|
||||
throw std::runtime_error("error reading transfer function file \'" + m_filename_Tk + "\'");
|
||||
throw std::runtime_error("error reading transfer function file \'" + filename + "\'");
|
||||
}
|
||||
|
||||
if (cosmo_params_["Omega_b"] < 1e-6)
|
||||
Tkvtot = Tktot;
|
||||
else
|
||||
Tkvtot = cosmo_params_["f_c"] * Tkvc + cosmo_params_["f_b"]* Tkvb;
|
||||
// if (cosmo_params_["Omega_b"] < 1e-6)
|
||||
// Tkvtot = Tktot;
|
||||
// else
|
||||
// Tkvtot = cosmo_params_["f_c"] * Tkvc + cosmo_params_["f_b"]* Tkvb;
|
||||
|
||||
m_linbaryoninterp |= Tkb < 0.0 || Tkvb < 0.0;
|
||||
// m_linbaryoninterp |= Tkb < 0.0 || Tkvb < 0.0;
|
||||
|
||||
m_tab_k.push_back(log10(k));
|
||||
|
||||
m_tab_Tk_tot.push_back(Tktot);
|
||||
m_tab_Tk_baryon.push_back(Tkb);
|
||||
m_tab_Tk_cdm.push_back(Tkc);
|
||||
m_tab_Tvk_tot.push_back(Tkvtot);
|
||||
m_tab_Tvk_baryon.push_back(Tkvb);
|
||||
m_tab_Tvk_cdm.push_back(Tkvc);
|
||||
|
||||
++m_nlines;
|
||||
|
||||
if (k < m_kmin)
|
||||
m_kmin = k;
|
||||
if (k > m_kmax)
|
||||
m_kmax = k;
|
||||
}
|
||||
|
||||
for (size_t i = 0; i < m_tab_k.size(); ++i)
|
||||
{
|
||||
m_tab_Tk_tot[i] = log10(m_tab_Tk_tot[i]);
|
||||
m_tab_Tk_cdm[i] = log10(m_tab_Tk_cdm[i]);
|
||||
m_tab_Tvk_cdm[i] = log10(m_tab_Tvk_cdm[i]);
|
||||
m_tab_Tvk_tot[i] = log10(m_tab_Tvk_tot[i]);
|
||||
|
||||
if (!m_linbaryoninterp)
|
||||
{
|
||||
m_tab_Tk_baryon[i] = log10(m_tab_Tk_baryon[i]);
|
||||
m_tab_Tvk_baryon[i] = log10(m_tab_Tvk_baryon[i]);
|
||||
}
|
||||
k.push_back(Tk);
|
||||
dc.push_back(Tdc);
|
||||
db.push_back(Tdb);
|
||||
dn.push_back(Tdn);
|
||||
dm.push_back(Tdm);
|
||||
tc.push_back(Tvc);
|
||||
tb.push_back(Tvb);
|
||||
tn.push_back(Tdm);
|
||||
tm.push_back(Tdm);
|
||||
++nlines;
|
||||
}
|
||||
|
||||
ifs.close();
|
||||
music::ilog.Print("Read CAMB transfer function table with %d rows", nlines);
|
||||
|
||||
music::ilog.Print("Read CAMB transfer function table with %d rows", m_nlines);
|
||||
// if (m_linbaryoninterp)
|
||||
// music::ilog.Print("Using log-lin interpolation for baryons\n (TF is not positive definite)");
|
||||
|
||||
if (m_linbaryoninterp)
|
||||
music::ilog.Print("Using log-lin interpolation for baryons\n (TF is not "
|
||||
"positive definite)");
|
||||
|
||||
#ifdef WITH_MPI
|
||||
}
|
||||
|
||||
unsigned n = m_tab_k.size();
|
||||
MPI::COMM_WORLD.Bcast(&n, 1, MPI_UNSIGNED, 0);
|
||||
#if defined(USE_MPI)
|
||||
unsigned n = k.size();
|
||||
MPI_Bcast(&n, 1, MPI_UNSIGNED, 0, MPI_COMM_WORLD);
|
||||
|
||||
if (MPI::COMM_WORLD.Get_rank() > 0)
|
||||
if (CONFIG::MPI_task_rank > 0)
|
||||
{
|
||||
m_tab_k.assign(n, 0);
|
||||
m_tab_Tk_tot.assign(n, 0);
|
||||
m_tab_Tk_cdm.assign(n, 0);
|
||||
m_tab_Tk_baryon.assign(n, 0);
|
||||
m_tab_Tvk_tot.assign(n, 0);
|
||||
m_tab_Tvk_cdm.assign(n, 0);
|
||||
m_tab_Tvk_baryon.assign(n, 0);
|
||||
k.assign(n, 0);
|
||||
dc.assign(n, 0);
|
||||
tc.assign(n, 0);
|
||||
db.assign(n, 0);
|
||||
tb.assign(n, 0);
|
||||
dn.assign(n, 0);
|
||||
tn.assign(n, 0);
|
||||
dm.assign(n, 0);
|
||||
tm.assign(n, 0);
|
||||
}
|
||||
|
||||
MPI::COMM_WORLD.Bcast(&m_tab_k[0], n, MPI_DOUBLE, 0);
|
||||
MPI::COMM_WORLD.Bcast(&m_tab_Tk_tot[0], n, MPI_DOUBLE, 0);
|
||||
MPI::COMM_WORLD.Bcast(&m_tab_Tk_cdm[0], n, MPI_DOUBLE, 0);
|
||||
MPI::COMM_WORLD.Bcast(&m_tab_Tk_baryon[0], n, MPI_DOUBLE, 0);
|
||||
MPI::COMM_WORLD.Bcast(&m_tab_Tvk_tot[0], n, MPI_DOUBLE, 0);
|
||||
MPI::COMM_WORLD.Bcast(&m_tab_Tvk_cdm[0], n, MPI_DOUBLE, 0);
|
||||
MPI::COMM_WORLD.Bcast(&m_tab_Tvk_baryon[0], n, MPI_DOUBLE, 0);
|
||||
|
||||
MPI_Bcast(&k[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
|
||||
MPI_Bcast(&dc[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
|
||||
MPI_Bcast(&tc[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
|
||||
MPI_Bcast(&db[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
|
||||
MPI_Bcast(&tb[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
|
||||
MPI_Bcast(&dn[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
|
||||
MPI_Bcast(&tn[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
|
||||
MPI_Bcast(&dm[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
|
||||
MPI_Bcast(&tm[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
|
||||
#endif
|
||||
|
||||
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);
|
||||
|
||||
// do not use first and last value since interpolation becomes lower order
|
||||
m_kmin = k[1];
|
||||
m_kmax = k[k.size()-2];
|
||||
}
|
||||
|
||||
public:
|
||||
|
|
@ -188,37 +160,11 @@ public:
|
|||
{
|
||||
music::wlog << "The CAMB file plugin is not well tested! Proceed with checks of correctness of output before running a simulation!" << std::endl;
|
||||
|
||||
m_filename_Tk = pcf_->get_value<std::string>("cosmology", "transfer_file");
|
||||
std::string filename = pcf_->get_value<std::string>("cosmology", "transfer_file");
|
||||
|
||||
read_table();
|
||||
|
||||
acc_tot = gsl_interp_accel_alloc();
|
||||
acc_cdm = gsl_interp_accel_alloc();
|
||||
acc_baryon = gsl_interp_accel_alloc();
|
||||
acc_vtot = gsl_interp_accel_alloc();
|
||||
acc_vcdm = gsl_interp_accel_alloc();
|
||||
acc_vbaryon = gsl_interp_accel_alloc();
|
||||
|
||||
spline_tot = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
|
||||
spline_cdm = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
|
||||
spline_baryon = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
|
||||
spline_vtot = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
|
||||
spline_vcdm =
|
||||
gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
|
||||
spline_vbaryon =
|
||||
gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
|
||||
|
||||
gsl_spline_init(spline_tot, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size());
|
||||
gsl_spline_init(spline_cdm, &m_tab_k[0], &m_tab_Tk_cdm[0], m_tab_k.size());
|
||||
gsl_spline_init(spline_baryon, &m_tab_k[0], &m_tab_Tk_baryon[0],
|
||||
m_tab_k.size());
|
||||
gsl_spline_init(spline_vtot, &m_tab_k[0], &m_tab_Tvk_tot[0],
|
||||
m_tab_k.size());
|
||||
gsl_spline_init(spline_vcdm, &m_tab_k[0], &m_tab_Tvk_cdm[0],
|
||||
m_tab_k.size());
|
||||
gsl_spline_init(spline_vbaryon, &m_tab_k[0], &m_tab_Tvk_baryon[0],
|
||||
m_tab_k.size());
|
||||
this->read_table( filename );
|
||||
|
||||
// set properties of this transfer function plugin:
|
||||
tf_distinct_ = true; // different density between CDM v.s. Baryon
|
||||
tf_withvel_ = true; // using velocity transfer function
|
||||
tf_withtotal0_ = false; // only have 1 file for the moment
|
||||
|
|
@ -226,182 +172,54 @@ public:
|
|||
|
||||
~transfer_CAMB_file_plugin()
|
||||
{
|
||||
gsl_spline_free(spline_tot);
|
||||
gsl_spline_free(spline_cdm);
|
||||
gsl_spline_free(spline_baryon);
|
||||
gsl_spline_free(spline_vtot);
|
||||
gsl_spline_free(spline_vcdm);
|
||||
gsl_spline_free(spline_vbaryon);
|
||||
|
||||
gsl_interp_accel_free(acc_tot);
|
||||
gsl_interp_accel_free(acc_cdm);
|
||||
gsl_interp_accel_free(acc_baryon);
|
||||
gsl_interp_accel_free(acc_vtot);
|
||||
gsl_interp_accel_free(acc_vcdm);
|
||||
gsl_interp_accel_free(acc_vbaryon);
|
||||
}
|
||||
|
||||
// linear interpolation in log-log
|
||||
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);
|
||||
|
||||
double lk = log10(k);
|
||||
double dk = m_tab_k[n] - m_tab_k[n1];
|
||||
double delk = lk - m_tab_k[n];
|
||||
double dc{0.0}, db{0.0};
|
||||
|
||||
switch (type)
|
||||
{
|
||||
case total0:
|
||||
case total:
|
||||
v1 = m_tab_Tk_tot[n1];
|
||||
v2 = m_tab_Tk_tot[n];
|
||||
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
||||
|
||||
case cdm0:
|
||||
case cdm:
|
||||
v1 = m_tab_Tk_cdm[n1];
|
||||
v2 = m_tab_Tk_cdm[n];
|
||||
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
||||
|
||||
case baryon0:
|
||||
case baryon:
|
||||
v1 = m_tab_Tk_baryon[n1];
|
||||
v2 = m_tab_Tk_baryon[n];
|
||||
if (m_linbaryoninterp)
|
||||
return std::max((v2 - v1) / dk * (delk) + v2, tiny);
|
||||
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
||||
|
||||
case deltabc:
|
||||
v1 = m_tab_Tk_cdm[n1];
|
||||
v2 = m_tab_Tk_cdm[n];
|
||||
dc = pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
||||
v1 = m_tab_Tk_baryon[n1];
|
||||
v2 = m_tab_Tk_baryon[n];
|
||||
db = pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
||||
return db - dc;
|
||||
|
||||
case vtotal0:
|
||||
case vtotal:
|
||||
v1 = m_tab_Tvk_tot[n1];
|
||||
v2 = m_tab_Tvk_tot[n];
|
||||
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
||||
|
||||
case vcdm0:
|
||||
case vcdm:
|
||||
v1 = m_tab_Tvk_cdm[n1];
|
||||
v2 = m_tab_Tvk_cdm[n];
|
||||
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
||||
|
||||
case vbaryon0:
|
||||
case vbaryon:
|
||||
v1 = m_tab_Tvk_baryon[n1];
|
||||
v2 = m_tab_Tvk_baryon[n];
|
||||
if (m_linbaryoninterp)
|
||||
return std::max((v2 - v1) / dk * (delk) + v2, tiny);
|
||||
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
||||
|
||||
default:
|
||||
throw std::runtime_error(
|
||||
"Invalid type requested in transfer function evaluation");
|
||||
}
|
||||
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
//!< return log-log interpolated values for transfer funtion '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)
|
||||
{
|
||||
case total0:
|
||||
case total:
|
||||
return pow(10.0, m_tab_Tk_tot[0]);
|
||||
// use constant interpolation on the left side of the tabulated values, i.e.
|
||||
// set values k<k_min to value at k_min! (since transfer functions asymptote to constant)
|
||||
k = std::max(k,m_kmin);
|
||||
|
||||
case cdm0:
|
||||
case cdm:
|
||||
return pow(10.0, m_tab_Tk_cdm[0]);
|
||||
|
||||
case baryon0:
|
||||
case baryon:
|
||||
if (m_linbaryoninterp)
|
||||
return m_tab_Tk_baryon[0];
|
||||
return pow(10.0, m_tab_Tk_baryon[0]);
|
||||
|
||||
case deltabc:
|
||||
return pow(10.0, m_tab_Tk_baryon[0]) - pow(10.0, m_tab_Tk_cdm[0]);
|
||||
|
||||
case vtotal0:
|
||||
case vtotal:
|
||||
return pow(10.0, m_tab_Tvk_tot[0]);
|
||||
|
||||
case vcdm0:
|
||||
case vcdm:
|
||||
return pow(10.0, m_tab_Tvk_cdm[0]);
|
||||
|
||||
case vbaryon0:
|
||||
case vbaryon:
|
||||
if (m_linbaryoninterp)
|
||||
return m_tab_Tvk_baryon[0];
|
||||
return pow(10.0, m_tab_Tvk_baryon[0]);
|
||||
|
||||
default:
|
||||
throw std::runtime_error(
|
||||
"Invalid type requested in transfer function evaluation");
|
||||
}
|
||||
}
|
||||
// use linear interpolation on the right side of the tabulated values
|
||||
else if (k > m_kmax)
|
||||
return extrap_right(k, type);
|
||||
|
||||
double lk = log10(k);
|
||||
switch (type)
|
||||
{
|
||||
case total0:
|
||||
case total:
|
||||
return pow(10.0, gsl_spline_eval(spline_tot, lk, acc_tot));
|
||||
case total:
|
||||
return delta_m_(k);
|
||||
|
||||
case cdm0:
|
||||
case cdm:
|
||||
return pow(10.0, gsl_spline_eval(spline_cdm, lk, acc_cdm));
|
||||
return delta_c_(k);
|
||||
|
||||
case baryon0:
|
||||
case baryon:
|
||||
if (m_linbaryoninterp)
|
||||
return gsl_spline_eval(spline_baryon, lk, acc_baryon);
|
||||
return pow(10.0, gsl_spline_eval(spline_baryon, lk, acc_baryon));
|
||||
|
||||
case deltabc:
|
||||
return pow(10.0, gsl_spline_eval(spline_baryon, lk, acc_baryon)) - pow(10.0, gsl_spline_eval(spline_cdm, lk, acc_cdm));
|
||||
return delta_b_(k);
|
||||
|
||||
case vtotal0:
|
||||
case vtotal:
|
||||
return pow(10.0, gsl_spline_eval(spline_vtot, lk, acc_vtot));
|
||||
return theta_m_(k);
|
||||
|
||||
case vcdm0:
|
||||
case vcdm:
|
||||
return pow(10.0, gsl_spline_eval(spline_vcdm, lk, acc_vcdm));
|
||||
return theta_c_(k);
|
||||
|
||||
case vbaryon0:
|
||||
case vbaryon:
|
||||
if (m_linbaryoninterp)
|
||||
return gsl_spline_eval(spline_vbaryon, lk, acc_vbaryon);
|
||||
return pow(10.0, gsl_spline_eval(spline_vbaryon, lk, acc_vbaryon));
|
||||
return theta_b_(k);
|
||||
|
||||
case deltabc:
|
||||
return delta_b_(k)-delta_c_(k);
|
||||
|
||||
default:
|
||||
throw std::runtime_error(
|
||||
"Invalid type requested in transfer function evaluation");
|
||||
throw std::runtime_error("Invalid type requested in transfer function evaluation");
|
||||
}
|
||||
}
|
||||
|
||||
inline double get_kmin(void) const { return pow(10.0, m_tab_k[1]); }
|
||||
//!< Return minimum k for which we can interpolate
|
||||
inline double get_kmin(void) const { return m_kmin; }
|
||||
|
||||
inline double get_kmax(void) const { return pow(10.0, m_tab_k[m_tab_k.size() - 2]); }
|
||||
//!< Return maximum k for which we can interpolate
|
||||
inline double get_kmax(void) const { return m_kmax; }
|
||||
};
|
||||
|
||||
namespace
|
||||
|
|
|
|||
|
|
@ -41,12 +41,6 @@ private:
|
|||
interpolated_function_1d<true, true, false> delta_c_, delta_b_, delta_n_, delta_m_, theta_c_, theta_b_, theta_n_, theta_m_;
|
||||
interpolated_function_1d<true, true, false> delta_c0_, delta_b0_, delta_n0_, delta_m0_, theta_c0_, theta_b0_, theta_n0_, theta_m0_;
|
||||
|
||||
// single fluid growing/decaying mode decomposition
|
||||
// gsl_interp_accel *gsl_ia_Cplus_, *gsl_ia_Cminus_;
|
||||
// gsl_spline *gsl_sp_Cplus_, *gsl_sp_Cminus_;
|
||||
// std::vector<double> tab_Cplus_, tab_Cminus_;
|
||||
|
||||
//double Omega_m_, Omega_b_, N_ur_, zstart_, ztarget_, kmax_, kmin_, h_, astart_, atarget_, A_s_, n_s_, sigma8_, Tcmb_, tnorm_;
|
||||
double zstart_, ztarget_, astart_, atarget_, kmax_, kmin_, h_, tnorm_;
|
||||
|
||||
ClassParams pars_;
|
||||
|
|
@ -175,9 +169,9 @@ private:
|
|||
|
||||
the_ClassEngine_->getTk(z, k, dc, db, dn, dm, tc, tb, tn, tm);
|
||||
|
||||
double h = cosmo_params_.get("h");
|
||||
double fc = cosmo_params_.get("f_c");
|
||||
double fb = cosmo_params_.get("f_b");
|
||||
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)
|
||||
{
|
||||
|
|
|
|||
Loading…
Reference in a new issue