CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
PkXiNonLinear.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2010 by Federico Marulli *
3  * federico.marulli3@unibo.it *
4  * *
5  * This program is free software; you can redistribute it and/or *
6  * modify it under the terms of the GNU General Public License as *
7  * published by the Free Software Foundation; either version 2 of *
8  * the License, or (at your option) any later version. *
9  * *
10  * This program is distributed in the hope that it will be useful, *
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13  * GNU General Public License for more details. *
14  * *
15  * You should have received a copy of the GNU General Public *
16  * License along with this program; if not, write to the Free *
17  * Software Foundation, Inc., *
18  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19  ********************************************************************/
20 
36 #include "Cosmology.h"
37 
38 using namespace std;
39 
40 using namespace cbl;
41 using namespace cosmology;
42 
43 
44 // =====================================================================================
45 
46 
47 double cbl::cosmology::Cosmology::f_k (const double kk, const shared_ptr<cbl::glob::FuncGrid> Pk, const double qmin, const double qmax, const double prec)
48 {
49  auto integrand = [&] (const double qq) {
50  double _Pk = Pk->operator()(qq);
51  double fact = 6.*pow(kk, 7)*qq-79.*pow(kk, 5)*pow(qq, 3)+50.*pow(kk,3)*pow(qq, 5)-21.*kk*pow(qq, 7)+3./4*pow(kk*kk-qq*qq, 3)*(2.*kk*kk+7.*qq*qq)*log(pow(fabs(kk-qq), 2)/pow(fabs(kk+qq), 2));
52  return pow(504.*pow(kk, 3)*pow(qq, 5), -1)*fact*_Pk*qq*qq;
53  };
54 
55  return 4.*par::pi*wrapper::gsl::GSL_integrate_cquad(integrand, qmin, qmax, prec);
56 }
57 
58 
59 // =====================================================================================
60 
61 
62 double cbl::cosmology::Cosmology::g_k (const double kk, const shared_ptr<cbl::glob::FuncGrid> Pk, const double qmin, const double qmax, const double prec)
63 {
64  auto integrand = [&] (const double qq) {
65  double _Pk = Pk->operator()(qq);
66  double fact = 6.*pow(kk, 7)*qq-41.*pow(kk, 5)*pow(qq, 3)+2.*pow(kk,3)*pow(qq, 5)-3.*kk*pow(qq, 7)+3./4*pow(kk*kk-qq*qq, 3)*(2.*kk*kk+qq*qq)*log(pow(fabs(kk-qq), 2)/pow(fabs(kk+qq), 2));
67  return pow(168.*pow(kk, 3)*pow(qq, 5), -1)*fact*_Pk*qq*qq;
68  };
69 
70  return 4.*par::pi*wrapper::gsl::GSL_integrate_cquad(integrand, qmin, qmax, prec);
71 }
72 
73 
74 // =====================================================================================
75 
76 
77 double cbl::cosmology::Cosmology::F2 (const double k, const double q, const double kq)
78 {
79  return 5./7. + kq/2. *(k/q+q/k) + 2./7.*kq*kq;
80 }
81 
82 
83 // =====================================================================================
84 
85 
86 double cbl::cosmology::Cosmology::G2 (const double k, const double q, const double kq)
87 {
88  return 3./7. + kq/2. *(k/q+q/k) + 4./7.*kq*kq;
89 }
90 
91 
92 // ============================================================================================
93 
94 
95 double cbl::cosmology::Cosmology::Pk_1loop (const double kk, const shared_ptr<cbl::glob::FuncGrid> Pk, const int corrtype, const double qmin, const double qmax, const double prec)
96 {
97  function<double(double, double, double)> func1, func2;
98 
99  if (corrtype==0) {
100  func1 = [&] (const double k, const double q, const double kq) {return F2(k, q, kq);};
101  func2 = [&] (const double k, const double q, const double kq) {return F2(k, q, kq);};
102  }
103  else if (corrtype==1) {
104  func1 = [&] (const double k, const double q, const double kq) {return F2(k, q, kq);};
105  func2 = [&] (const double k, const double q, const double kq) {return G2(k, q, kq);};
106  }
107  else if (corrtype==2) {
108  func1 = [&] (const double k, const double q, const double kq) {return G2(k, q, kq);};
109  func2 = [&] (const double k, const double q, const double kq) {return G2(k, q, kq);};
110  }
111  else
112  ErrorCBL("the input value of corrtype is not allowed!", "Pk_1loop", "PkXiNonLinear.cpp");
113 
114  auto integrand = [&] (const double qq)
115  {
116  auto integrand_intermediate = [&] (const double xx)
117  {
118  double kq = sqrt(kk*kk+qq*qq-2.*kk*qq*xx);
119  double akq = (kk*xx-qq)/kq;
120  return func1(kq, qq, akq)*func2(kq, qq, akq)*Pk->operator()(kq);
121  };
122  double res = wrapper::gsl::GSL_integrate_cquad(integrand_intermediate, -1, 1, prec);
123  return Pk->operator()(qq)*qq*qq*res;
124  };
125 
126  return 4.*par::pi*wrapper::gsl::GSL_integrate_cquad(integrand, qmin, qmax, prec);
127 }
128 
129 
130 // ============================================================================================
131 
132 
133 double cbl::cosmology::Cosmology::Pk_DeltaDelta (const double kk, const std::shared_ptr<cbl::glob::FuncGrid> Pk, const double qmin, const double qmax, const double prec)
134 {
135  double GG = pow(exp(f_k(kk, Pk, qmin, qmax, prec)), 2);
136  return pow(2*par::pi, 3)*GG*(Pk->operator()(kk)+Pk_1loop(kk, Pk, 0, qmin, qmax, prec));
137 }
138 
139 
140 // ============================================================================================
141 
142 
143 double cbl::cosmology::Cosmology::Pk_DeltaTheta (const double kk, const std::shared_ptr<cbl::glob::FuncGrid> Pk, const double qmin, const double qmax, const double prec)
144 {
145  double GG = pow(exp(f_k(kk, Pk, qmin, qmax, prec)), 2); // to be checked: possible bug!!!
146  return pow(2*par::pi, 3)*GG*(Pk->operator()(kk)+Pk_1loop(kk, Pk, 1, qmin, qmax, prec));
147 }
148 
149 
150 // ============================================================================================
151 
152 
153 double cbl::cosmology::Cosmology::Pk_ThetaTheta (const double kk, const std::shared_ptr<cbl::glob::FuncGrid> Pk, const double qmin, const double qmax, const double prec)
154 {
155  double GG = exp(g_k(kk, Pk, qmin, qmax, prec)); // to be checked: possible bug!!!
156  return pow(2*par::pi, 3)*GG*(Pk->operator()(kk)+Pk_1loop(kk, Pk, 2, qmin, qmax, prec));
157 }
158 
159 // ============================================================================================
160 
161 
162 std::vector<double> cbl::cosmology::Cosmology::Pk_DeltaDelta (const std::vector<double> kk, const double redshift, const std::string method_Pk, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec, const std::string file_par, const bool unit1)
163 {
164  vector<double> pkLin = Pk_matter(kk, method_Pk, false, redshift, store_output, output_root, norm, k_min, k_max, prec, file_par, unit1);
165 
166  for (size_t i=0; i<kk.size(); i++)
167  pkLin[i] /= pow(2*par::pi, 3);
168 
169  auto pkLinInterp = make_shared<glob::FuncGrid>(glob::FuncGrid(kk, pkLin, "Spline"));
170 
171  vector<double> pkDeltaDelta(kk.size(), 0);
172 
173  for (size_t i=0; i<kk.size(); i++)
174  pkDeltaDelta[i] = Pk_DeltaDelta(kk[i], pkLinInterp, k_min, k_max, prec);
175 
176  return pkDeltaDelta;
177 }
178 
179 
180 // ============================================================================================
181 
182 
183 std::vector<double> cbl::cosmology::Cosmology::Pk_DeltaTheta (const std::vector<double> kk, const double redshift, const std::string method_Pk, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec, const std::string file_par, const bool unit1)
184 {
185  vector<double> pkLin = Pk_matter(kk, method_Pk, false, redshift, store_output, output_root, norm, k_min, k_max, prec, file_par, unit1);
186 
187  for (size_t i=0; i<kk.size(); i++)
188  pkLin[i] /= pow(2*par::pi, 3);
189 
190  auto pkLinInterp = make_shared<glob::FuncGrid>(glob::FuncGrid(kk, pkLin, "Spline"));
191 
192  vector<double> pkDeltaTheta(kk.size(), 0);
193 
194  for (size_t i=0; i<kk.size(); i++)
195  pkDeltaTheta[i] = Pk_DeltaTheta(kk[i], pkLinInterp, k_min, k_max, prec);
196 
197  return pkDeltaTheta;
198 }
199 
200 // ============================================================================================
201 
202 
203 std::vector<double> cbl::cosmology::Cosmology::Pk_ThetaTheta (const std::vector<double> kk, const double redshift, const std::string method_Pk, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec, const std::string file_par, const bool unit1)
204 {
205  vector<double> pkLin = Pk_matter(kk, method_Pk, false, redshift, store_output, output_root, norm, k_min, k_max, prec, file_par, unit1);
206 
207  for (size_t i=0; i<kk.size(); i++)
208  pkLin[i] /= pow(2*par::pi, 3);
209 
210  auto pkLinInterp = make_shared<glob::FuncGrid>(glob::FuncGrid(kk, pkLin, "Spline"));
211 
212  vector<double> pkThetaTheta(kk.size(), 0);
213 
214  for (size_t i=0; i<kk.size(); i++)
215  pkThetaTheta[i] = Pk_ThetaTheta(kk[i], pkLinInterp, k_min, k_max, prec);
216 
217  return pkThetaTheta;
218 }
219 
220 
221 // =====================================================================================
222 
223 
224 std::vector<std::vector<double>> cbl::cosmology::Cosmology::Pk_TNS_AB_multipoles (std::vector<double> kk, const std::string method, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec)
225 {
226  cbl::Path path;
227  string dir = path.DirCosmo()+"/External/CPT_Library/";
228  string output_tmpCPT = dir+"tmpCPT/";
229  string MKout = "mkdir -p " + output_tmpCPT; if (system(MKout.c_str())) {}
230  double sigma8_z0 = sigma8_Pk(method, 0., store_output, output_root);
231 
232  // Pklin_z0
233  const vector<double> Pklin = Pk_matter(kk, method, false, 0., store_output, output_root, norm, k_min, k_max, prec);
234  string file = "Pklin.dat";
235  ofstream File_Pklin(output_tmpCPT + file);
236  for (size_t nn=0; nn<kk.size(); ++nn)
237  File_Pklin << kk[nn] << "\t" << Pklin[nn] << endl;
238  File_Pklin.close();
239 
240  // setting parameters for PkA and PkB
241  string File_par = "params_AB_termsTNS.ini";
242  ofstream fsAB(output_tmpCPT + File_par);
243  fsAB << redshift << "\n"
244  << "1 100. \n"
245  << "Pklin.dat \n"
246  << "1 \n" << conv(m_hh, par::fDP6) <<"\n"
247  << "2 \n" << par::TCMB << "\n"
248  << "3 \n" << n_spec() << "\n"
249  << "4 \n" << sigma8_z0 << "\n"
250  << "5 \n" << conv(m_Omega_matter, par::fDP6) << "\n"
251  << "6 \n" << conv(m_Omega_baryon, par::fDP6) << "\n"
252  << "7 \n" << conv(m_w0, par::fDP6) << "\n"
253  << "0 \n" << "1 \n" << "0. \n";
254  fsAB.close();
255  string calc_pk_correction = "cd " + output_tmpCPT + " && " + dir + "calc_pk_correction < " + File_par; if (system (calc_pk_correction.c_str())) {}
256  string calc_pk_correction2 = "cd " + output_tmpCPT + " && " + dir + "calc_pk_correction2 < " + File_par; if (system (calc_pk_correction2.c_str())) {}
257 
258  double KK, PK, PK0EH, PK0corr, PK2EH, PK2corr, PK4EH, PK4corr;
259  vector<double> kA, pkA0, pkA2, pkA4, kB, pkB0, pkB2, pkB4;
260 
261  const string filenameB = output_tmpCPT + "corr_pkred.dat"; // PkB terms
262  ifstream finB(filenameB.c_str());
263  while (finB >> KK >> PK >> PK0EH >> PK0corr >> PK2EH >> PK2corr >> PK4EH >> PK4corr){
264  kB.emplace_back(KK);
265  pkB0.emplace_back(PK0corr);
266  pkB2.emplace_back(PK2corr);
267  pkB4.emplace_back(PK4corr);
268  }
269  finB.clear();
270 
271  const string filenameA = output_tmpCPT + "corr_pkred2.dat"; // PkA terms
272  ifstream finA(filenameA.c_str());
273  while (finA >> KK >> PK >> PK0EH >> PK0corr >> PK2EH >> PK2corr >> PK4EH >> PK4corr){
274  kA.emplace_back(KK);
275  pkA0.emplace_back(PK0corr);
276  pkA2.emplace_back(PK2corr);
277  pkA4.emplace_back(PK4corr);
278  }
279  finA.clear();
280 
281  vector<double> pkA0_new(kk.size()), pkA2_new(kk.size()), pkA4_new(kk.size()), pkB0_new(kk.size()), pkB2_new(kk.size()), pkB4_new(kk.size());
282  glob::FuncGrid interp_PkA0(kA, pkA0, "Spline");
283  glob::FuncGrid interp_PkA2(kA, pkA2, "Spline");
284  glob::FuncGrid interp_PkA4(kA, pkA4, "Spline");
285  glob::FuncGrid interp_PkB0(kB, pkB0, "Spline");
286  glob::FuncGrid interp_PkB2(kB, pkB2, "Spline");
287  glob::FuncGrid interp_PkB4(kB, pkB4, "Spline");
288 
289  pkA0_new = interp_PkA0.eval_func(kk);
290  pkA2_new = interp_PkA2.eval_func(kk);
291  pkA4_new = interp_PkA4.eval_func(kk);
292  pkB0_new = interp_PkB0.eval_func(kk);
293  pkB2_new = interp_PkB2.eval_func(kk);
294  pkB4_new = interp_PkB4.eval_func(kk);
295 
296  if (system (("rm -rf "+output_tmpCPT).c_str())) {}
297 
298  return {pkA0_new, pkA2_new, pkA4_new, pkB0_new, pkB2_new, pkB4_new};
299 }
300 
301 
302 // =====================================================================================
303 
304 
305 std::vector<std::vector<double>> cbl::cosmology::Cosmology::Pk_TNS_AB_1loop (std::vector<double> kk, const double mu, const std::string method, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec)
306 {
307  vector<vector<double>> Pk_AB_multipoles = Pk_TNS_AB_multipoles(kk, method, redshift, store_output, output_root, norm, k_min, k_max, prec);
308  vector<double> Pk_A, Pk_B;
309 
310  for (size_t ii=0; ii < Pk_AB_multipoles[0].size(); ++ii) {
311  Pk_A.emplace_back(Pk_AB_multipoles[0][ii] + Pk_AB_multipoles[1][ii]*cbl::legendre_polynomial(mu, 2) + Pk_AB_multipoles[2][ii]*cbl::legendre_polynomial(mu, 4));
312  Pk_B.emplace_back(Pk_AB_multipoles[3][ii] + Pk_AB_multipoles[4][ii]*cbl::legendre_polynomial(mu, 2) + Pk_AB_multipoles[5][ii]*cbl::legendre_polynomial(mu, 4));
313  }
314 
315  return {Pk_A, Pk_B};
316 }
317 
318 
319 // =====================================================================================
320 
321 
322 std::vector<std::vector<double>> cbl::cosmology::Cosmology::Pk_TNS_AB_terms_1loop (std::vector<double> kk, const std::string method, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec)
323 {
324  cbl::Path path;
325  string dir = path.DirCosmo()+"/External/CPT_Library/";
326  string output_tmpCPT = dir+"tmpCPT/";
327  string MKout = "mkdir -p " + output_tmpCPT; if (system(MKout.c_str())) {}
328  double sigma8_z0 = sigma8_Pk(method, 0., store_output, output_root);
329 
330  // input Pklin_z0
331  const vector<double> Pklin = Pk_matter(kk, method, false, 0., store_output, output_root, norm, k_min, k_max, prec);
332  string file = "Pklin.dat";
333  ofstream File_Pklin(output_tmpCPT + file);
334  for (size_t nn=0; nn<kk.size(); ++nn)
335  File_Pklin << kk[nn] << "\t" << Pklin[nn] << endl;
336  File_Pklin.close();
337 
338  // setting parameters for PkA and PkB
339  string File_par = "params_AB_termsTNS.ini";
340  ofstream fsAB(output_tmpCPT + File_par);
341  fsAB << redshift << "\n"
342  << "1 100. \n"
343  << "Pklin.dat \n"
344  << "1 \n" << conv(m_hh, par::fDP6) <<"\n"
345  << "2 \n" << par::TCMB << "\n"
346  << "3 \n" << n_spec() << "\n"
347  << "4 \n" << sigma8_z0 << "\n"
348  << "5 \n" << conv(m_Omega_matter, par::fDP6) << "\n"
349  << "6 \n" << conv(m_Omega_baryon, par::fDP6) << "\n"
350  << "7 \n" << conv(m_w0, par::fDP6) << "\n"
351  << "0 \n" << "1 \n" << "0. \n";
352  fsAB.close();
353  string calc_pk_correction = "cd " + output_tmpCPT + " && " + dir + "calc_pk_correction < " + File_par; if (system (calc_pk_correction.c_str())) {}
354  string calc_pk_correction2 = "cd " + output_tmpCPT + " && " + dir + "calc_pk_correction2 < " + File_par; if (system (calc_pk_correction2.c_str())) {}
355 
356  double KK, A11, A12, A22, A23, A33, B12, B13, B14, B22, B23, B24, B33, B34, B44;
357  vector<double> kA, kB, pk_A11, pk_A12, pk_A22, pk_A23, pk_A33, pk_B12, pk_B13, pk_B14, pk_B22, pk_B23, pk_B24, pk_B33, pk_B34, pk_B44;
358 
359  const string filenameB = output_tmpCPT + "pkstd_corr_tree.dat"; // pk_B12, pk_B13, pk_B14, pk_B22, pk_B23, pk_B24, pk_B33, pk_B34, pk_B44
360 
361  ifstream finB(filenameB.c_str());
362  while (finB >> KK >> B12 >> B13 >> B14 >> B22 >> B23 >> B24 >> B33 >> B34 >> B44){
363  kB.emplace_back(KK);
364  pk_B12.emplace_back(B12);
365  pk_B13.emplace_back(B13);
366  pk_B14.emplace_back(B14);
367  pk_B22.emplace_back(B22);
368  pk_B23.emplace_back(B23);
369  pk_B24.emplace_back(B24);
370  pk_B33.emplace_back(B33);
371  pk_B34.emplace_back(B34);
372  pk_B44.emplace_back(B44);
373  }
374  finB.clear();
375 
376  const string filenameA = output_tmpCPT + "pkstd_corr2_tree.dat"; // pk_A11, pk_A12, pk_A22, pk_A23, pk_A33
377  ifstream finA(filenameA.c_str());
378  while (finA >> KK >> A11 >> A12 >> A22 >> A23 >> A33){
379  kA.emplace_back(KK);
380  pk_A11.emplace_back(A11);
381  pk_A12.emplace_back(A12);
382  pk_A22.emplace_back(A22);
383  pk_A23.emplace_back(A23);
384  pk_A33.emplace_back(A33);
385  }
386  finA.clear();
387 
388  vector<double> pk_A11_new(kk.size()), pk_A12_new(kk.size()), pk_A22_new(kk.size()), pk_A23_new(kk.size()), pk_A33_new(kk.size()), pk_B12_new(kk.size()), pk_B13_new(kk.size()), pk_B14_new(kk.size()), pk_B22_new(kk.size()), pk_B23_new(kk.size()), pk_B24_new(kk.size()), pk_B33_new(kk.size()), pk_B34_new(kk.size()), pk_B44_new(kk.size());
389 
390  glob::FuncGrid interp_A11(kA, pk_A11, "Spline");
391  glob::FuncGrid interp_A12(kA, pk_A12, "Spline");
392  glob::FuncGrid interp_A22(kA, pk_A22, "Spline");
393  glob::FuncGrid interp_A23(kA, pk_A23, "Spline");
394  glob::FuncGrid interp_A33(kA, pk_A33, "Spline");
395 
396  glob::FuncGrid interp_B12(kB, pk_B12, "Spline");
397  glob::FuncGrid interp_B13(kB, pk_B13, "Spline");
398  glob::FuncGrid interp_B14(kB, pk_B14, "Spline");
399  glob::FuncGrid interp_B22(kB, pk_B22, "Spline");
400  glob::FuncGrid interp_B23(kB, pk_B23, "Spline");
401  glob::FuncGrid interp_B24(kB, pk_B24, "Spline");
402  glob::FuncGrid interp_B33(kB, pk_B33, "Spline");
403  glob::FuncGrid interp_B34(kB, pk_B34, "Spline");
404  glob::FuncGrid interp_B44(kB, pk_B44, "Spline");
405 
406  pk_A11_new = interp_A11.eval_func(kk);
407  pk_A12_new = interp_A12.eval_func(kk);
408  pk_A22_new = interp_A22.eval_func(kk);
409  pk_A23_new = interp_A23.eval_func(kk);
410  pk_A33_new = interp_A33.eval_func(kk);
411  pk_B12_new = interp_B12.eval_func(kk);
412  pk_B13_new = interp_B13.eval_func(kk);
413  pk_B14_new = interp_B14.eval_func(kk);
414  pk_B22_new = interp_B22.eval_func(kk);
415  pk_B23_new = interp_B23.eval_func(kk);
416  pk_B24_new = interp_B24.eval_func(kk);
417  pk_B33_new = interp_B33.eval_func(kk);
418  pk_B34_new = interp_B34.eval_func(kk);
419  pk_B44_new = interp_B44.eval_func(kk);
420 
421  // define the normalization
422  int Norm = norm;
423  if (Norm==-1) Norm = (m_sigma8>0) ? 1 : 0;
424  if (Norm==1) {
425  double sigma8;
426  const double RR = 8.;
427  glob::FuncGrid interpPk(kk, Pklin, "Spline");
428  auto func_sigma = [&] (double _k) { return pow(TopHat_WF(_k*RR)*_k, 2)*interpPk(_k); };
429  sigma8 = sqrt(1./(2.*pow(par::pi, 2))*wrapper::gsl::GSL_integrate_qag (func_sigma, k_min, k_max, 1.e-5))/DN(redshift, 0.);
430  m_Pk0_CAMB = pow(m_sigma8/sigma8,2);
431  }
432  else { m_Pk0_CAMB = 1.;}
433 
434  for (size_t i=0; i<kk.size(); i++) {
435  pk_A11_new[i] *= m_Pk0_CAMB;
436  pk_A12_new[i] *= m_Pk0_CAMB;
437  pk_A22_new[i] *= m_Pk0_CAMB;
438  pk_A23_new[i] *= m_Pk0_CAMB;
439  pk_A33_new[i] *= m_Pk0_CAMB;
440  pk_B12_new[i] *= m_Pk0_CAMB;
441  pk_B13_new[i] *= m_Pk0_CAMB;
442  pk_B14_new[i] *= m_Pk0_CAMB;
443  pk_B22_new[i] *= m_Pk0_CAMB;
444  pk_B23_new[i] *= m_Pk0_CAMB;
445  pk_B24_new[i] *= m_Pk0_CAMB;
446  pk_B33_new[i] *= m_Pk0_CAMB;
447  pk_B34_new[i] *= m_Pk0_CAMB;
448  pk_B44_new[i] *= m_Pk0_CAMB;
449  }
450 
451  if (system (("rm -rf " + output_tmpCPT).c_str())) {}
452 
453  return {pk_A11_new, pk_A12_new, pk_A22_new, pk_A23_new, pk_A33_new, pk_B12_new, pk_B13_new, pk_B14_new, pk_B22_new, pk_B23_new, pk_B24_new, pk_B33_new, pk_B34_new, pk_B44_new};
454 }
455 
456 
457 // =====================================================================================
458 
459 
460 std::vector<std::vector<double>> cbl::cosmology::Cosmology::Pk_TNS_AB_1loop (std::vector<double> kk, const double mu, const double linear_growth_rate, const double bias, const std::string method, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec)
461 {
462  double beta = linear_growth_rate/bias;
463 
464  vector<vector<double>> Pk_AB = Pk_TNS_AB_terms_1loop(kk, method, redshift, store_output, output_root, norm, k_min, k_max, prec);
465  vector<double> Pk_A, Pk_B;
466 
467  for(size_t ii=0; ii < kk.size(); ++ii) {
468  Pk_A.emplace_back(beta*pow(mu,2)*Pk_AB[0][ii] + pow(beta,2)*(pow(mu,2)*Pk_AB[1][ii] + pow(mu,4)*Pk_AB[2][ii]) + pow(beta,3)*(pow(mu,4)*Pk_AB[3][ii] + pow(mu,6)*Pk_AB[4][ii]));
469  Pk_B.emplace_back(pow(mu,2)*(pow(beta,2)*Pk_AB[5][ii] + pow(beta,3)*Pk_AB[6][ii] + pow(beta,4)*Pk_AB[7][ii]) + pow(mu,4)*(pow(beta,2)*Pk_AB[8][ii] + pow(beta,3)*Pk_AB[9][ii] + pow(beta,4)*Pk_AB[10][ii]) + pow(mu,6)*(pow(beta,3)*Pk_AB[11][ii] + pow(beta,4)*Pk_AB[12][ii]) + pow(mu,8)*pow(beta,4)*Pk_AB[13][ii]);
470  }
471 
472  return {Pk_A, Pk_B};
473 }
474 
475 
476 // =====================================================================================
477 
478 
479 std::vector<std::vector<double>> cbl::cosmology::Cosmology::Pk_TNS_dd_dt_tt (std::vector<double> kk, const std::string method, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec)
480 {
481  cbl::Path path;
482  string dir = path.DirCosmo()+"/External/CPT_Library/";
483  string output_tmpCPT = dir+"tmpCPT/";
484  string MKout = "mkdir -p " + output_tmpCPT; if (system(MKout.c_str())) {}
485  double sigma8_z0 = sigma8_Pk(method, 0., store_output, output_root);
486 
487  // input Pklin_z0
488  const vector<double> Pklin = Pk_matter(kk, method, false, 0., store_output, output_root, norm, k_min, k_max, prec);
489  string file = "Pklin.dat";
490  ofstream File_Pklin(output_tmpCPT + file);
491  for (size_t nn=0; nn<kk.size(); ++nn)
492  File_Pklin << kk[nn] << "\t" << Pklin[nn] << endl;
493  File_Pklin.close();
494 
495  // setting input parameters
496  string File_par = "params_stdPT2.ini";
497  ofstream File_stdPT2(output_tmpCPT + File_par);
498  File_stdPT2 << sigma8_z0 << "\n Pklin.dat";
499  File_stdPT2.close();
500  string stdPT2 = "cd " + output_tmpCPT + " && " + dir + "stdPT2 < " + File_par; if (system (stdPT2.c_str())) {}
501 
502  File_par = "params_read_pk2.ini";
503  ofstream File_read(output_tmpCPT + File_par);
504  File_read << "0 \n"
505  << "1 \n" << conv(m_Omega_matter, par::fDP6) <<"\n"
506  << "2 \n" << conv(m_Omega_baryon, par::fDP6) << "\n"
507  << "3 \n" << conv(m_hh, par::fDP6) << "\n"
508  << "4 \n" << conv(m_n_spec, par::fDP6) << "\n"
509  << "5 \n" << sigma8_z0 << "\n"
510  << "6 \n" << conv(m_w0, par::fDP6) << "\n"
511  << "0 \n" << "0 " << redshift << "\n"
512  << "1 \n" << "0 \n" << "0";
513  File_read.close();
514 
515  string read_pk2 = "cd " + output_tmpCPT + " && " + dir + "read_pk2 < " + File_par; if (system (read_pk2.c_str())) {}
516 
517  double Kspt, PKlinNWspt, PKlinspt, PDDspt, PDTspt, PTTspt;
518  vector<double> k_spt, PklinNW_spt, Pklin_spt, Pdd_spt, Pdt_spt, Ptt_spt;
519 
520  const string filenamePK = output_tmpCPT + "pkstd2.dat"; // Pkdd, Pkdt, Pktt
521  ifstream finPK(filenamePK.c_str());
522 
523  while (finPK >> Kspt >> PKlinNWspt >> PKlinspt >> PDDspt >> PDTspt >> PTTspt)
524  if (Kspt>0 && PKlinNWspt>0 && PKlinspt>0 && PDDspt>0 && PDTspt>0 && PTTspt>0) {
525  k_spt.emplace_back(Kspt);
526  PklinNW_spt.emplace_back(PKlinNWspt);
527  Pklin_spt.emplace_back(PKlinspt);
528  Pdd_spt.emplace_back(PDDspt);
529  Pdt_spt.emplace_back(PDTspt);
530  Ptt_spt.emplace_back(PTTspt);
531  }
532  finPK.clear();
533 
534  vector<double> Pkdd_new(kk.size()), Pkdt_new(kk.size()), Pktt_new(kk.size()), Pklin_new(kk.size());
535  glob::FuncGrid interp_Pklin(k_spt, Pklin_spt, "Spline");
536  glob::FuncGrid interp_Pkdd(k_spt, Pdd_spt, "Spline");
537  glob::FuncGrid interp_Pkdt(k_spt, Pdt_spt, "Spline");
538  glob::FuncGrid interp_Pktt(k_spt, Ptt_spt, "Spline");
539 
540  Pklin_new = interp_Pklin.eval_func(kk);
541  Pkdd_new = interp_Pkdd.eval_func(kk);
542  Pkdt_new = interp_Pkdt.eval_func(kk);
543  Pktt_new = interp_Pktt.eval_func(kk);
544 
545  // define the normalization
546  int Norm = norm;
547  if (Norm==-1) Norm = (m_sigma8>0) ? 1 : 0;
548  if (Norm==1) {
549  double sigma8;
550  const double RR = 8.;
551  glob::FuncGrid interpPk(kk, Pklin_new, "Spline");
552  auto func_sigma = [&] (double _k) { return pow(TopHat_WF(_k*RR)*_k, 2)*interpPk(_k); };
553  sigma8 = sqrt(1./(2.*pow(par::pi, 2))*wrapper::gsl::GSL_integrate_qag (func_sigma, k_min, k_max, 1.e-5))/DN(redshift, 0.);
554  m_Pk0_CAMB = pow(m_sigma8/sigma8, 2);
555  }
556  else { m_Pk0_CAMB = 1.; }
557 
558  for (size_t i=0; i<kk.size(); i++) {
559  Pkdd_new[i] *= m_Pk0_CAMB;
560  Pkdt_new[i] *= m_Pk0_CAMB;
561  Pktt_new[i] *= m_Pk0_CAMB;
562  }
563 
564  if (system (("rm -rf "+output_tmpCPT).c_str())) {}
565 
566  return {Pkdd_new, Pkdt_new, Pktt_new};
567 }
568 
569 
570 // =====================================================================================
571 
572 
573 std::vector<std::vector<double>> cbl::cosmology::Cosmology::Pk_eTNS_terms_1loop (std::vector<double> kk, const std::string method, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec)
574 {
575  cbl::Path path;
576  string dir = path.DirCosmo()+"/External/CAMB_SPT_private/";
577  string output_tmpCPT = dir+"tmpCPT_eTNS/";
578  string MKout = "mkdir -p " + output_tmpCPT; if (system(MKout.c_str())) {}
579  double HH0 = m_hh*100.;
580 
581  // input Pklin_z0
582 
583  const vector<double> Pklin = Pk_matter(kk, method, false, 0., store_output, output_root, norm, k_min, k_max, prec);
584  int Norm = norm;
585  if (Norm==-1) Norm = (m_sigma8>0) ? 1 : 0;
586  if (Norm==1) Pk_0(method, redshift, store_output, output_root, k_min, k_max, prec);
587 
588  double PP0 = 1.;
589  if (method=="EisensteinHu") PP0 = m_Pk0_EH;
590  else if (method=="CAMB" || method=="MGCAMB") PP0 = m_Pk0_CAMB;
591  else if (method=="MPTbreeze-v1") PP0 = m_Pk0_MPTbreeze;
592  else if (method=="CLASS") PP0 = m_Pk0_CLASS;
593 
594  string file = "Pklin.dat";
595  ofstream File_Pklin(output_tmpCPT + file);
596  for (size_t nn=0; nn<kk.size(); ++nn)
597  File_Pklin << kk[nn] << "\t" << Pklin[nn] << endl;
598  File_Pklin.close();
599 
600  // setting parameters for Pk_xy
601  string File_par = "SPT_NLB_params.ini";
602  ofstream fsPkt(output_tmpCPT + File_par);
603  fsPkt << "output_root = SPT_NLB \n" // OutputRoot
604  << "get_scalar_cls = F \n"
605  << "get_vector_cls = F \n"
606  << "get_tensor_cls = F \n"
607  << "COBE_normalize = F \n"
608  << "CMB_outputscale = 7.4311e12\n"
609  << "get_transfer = T \n"
610  << "do_nonlinear = 3 \n"
611  << "DoNonLocalBias = T \n"
612  << "RSDmodel = 1 \n"
613  << "w = " << conv(m_w0, par::fDP6) <<"\n"
614  << "cs2_lam = 1 \n"
615  << "hubble = " << conv(HH0, par::fDP6) <<"\n"
616  << "use_physical = T \n"
617  << "ombh2 = " << conv(m_Omega_baryon*m_hh*m_hh, par::fDP6) <<"\n"
618  << "omch2 = " << conv(m_Omega_CDM*m_hh*m_hh, par::fDP6) <<"\n"
619  << "omnuh2 = " << conv(m_Omega_neutrinos*m_hh*m_hh, par::fDP6) <<"\n"
620  << "omk = " << conv(m_Omega_k, par::fDP6) <<"\n"
621  << "temp_cmb = " << cbl::par::TCMB <<"\n"
622  << "helium_fraction = 0.24 \n"
623  << "massless_neutrinos = " << conv(m_massless_neutrinos, par::fDP6) <<"\n"
624  << "massive_neutrinos = " << conv(m_massive_neutrinos, par::fINT) <<"\n"
625  << "nu_mass_eigenstates = 1 \n"
626  << "nu_mass_degeneracies = 0 \n"
627  << "nu_mass_fractions = 1 \n"
628  << "transfer_high_precision = T \n"
629  << "transfer_kmax = " << k_max <<"\n"
630  << "transfer_k_per_logint = 20 \n"
631  << "transfer_num_redshifts = 1 \n"
632  << "transfer_interp_matterpower = T \n"
633  << "transfer_power_var = 7 \n"
634  << "transfer_redshift(1) = " << cbl::conv(redshift, cbl::par::fDP1) <<"\n"
635  << "transfer_filename(1) = Tk_z.dat \n"
636  << "transfer_matterpower(1) = Pk_z.dat \n"
637  << "reionization = T \n"
638  << "re_use_optical_depth = T \n"
639  << "re_optical_depth = " << conv(m_tau, par::fDP6) <<"\n"
640  << "re_delta_redshift = 1.5 \n"
641  << "re_ionization_frac = -1 \n"
642  << "pivot_scalar = " << conv(m_scalar_pivot, par::fDP6) <<"\n"
643  << "pivot_tensor = 0.002 \n"
644  << "initial_power_num = 1 \n"
645  << "scalar_spectral_index(1) = " << conv(m_n_spec, par::fDP6) <<"\n"
646  << "scalar_nrun(1) = 0 \n"
647  << "scalar_amp(1) = " << conv(m_scalar_amp, par::ee3) <<"\n"
648  << "RECFAST_fudge_He = 0.86 \n"
649  << "RECFAST_Heswitch = 6 \n"
650  << "RECFAST_Hswitch = T \n"
651  << "RECFAST_fudge = 1.14 \n"
652  << "do_lensing_bispectrum = F \n"
653  << "do_primordial_bispectrum = F \n"
654  << "initial_condition = 1 \n"
655  << "accurate_polarization = T \n"
656  << "accurate_reionization = T \n"
657  << "accurate_BB = F \n"
658  << "do_late_rad_truncation = T \n"
659  << "do_tensor_neutrinos = T \n"
660  << "feedback_level = 1 \n"
661  << "massive_nu_approx = 1 \n"
662  << "number_of_threads = 0 \n"
663  << "accuracy_boost = 2 \n"
664  << "l_accuracy_boost = 1 \n"
665  << "high_accuracy_default = F \n"
666  << "l_sample_boost = 1 ";
667  fsPkt.close();
668  string calc_pk_eTNScorrection = "cd " + output_tmpCPT + " && " + dir + "camb " + File_par; if (system (calc_pk_eTNScorrection.c_str())) {}
669 
670  double Kspt, PKlin, PDD, PDV, PVV, PB2D, PB2V, PB22, PBS2D, PBS2V, PB2S2, PBS22, sigma32PKlin, BB1, BB2, BBS2;
671  vector<double> k_spt, Pdd, Pdv, Pvv, Pb2d, Pb2v, Pb22, Pbs2d, Pbs2v, Pb2s2, Pbs22, sigma32Pklin, Bb1, Bb2, Bbs2;
672 
673  const string filenamePK = output_tmpCPT + "SPT_NLB_Pk_z.dat"; // k, Pklin, Pdd, Pdv, Pvv, Pb2d, Pb2v, Pb22, Pbs2d, Pbs2v, Pb2s2, Pbs22, sigma3^2Pklin, Bb1, Bb2, Bbs2
674 
675  ifstream finPK(filenamePK.c_str());
676  while (finPK >> Kspt >> PKlin >> PDD >> PDV >> PVV >> PB2D >> PB2V >> PB22 >> PBS2D >> PBS2V >> PB2S2 >> PBS22 >> sigma32PKlin >> BB1 >> BB2 >> BBS2)
677  if (Kspt>0 && PKlin>0 && PDD>0 && PDV>0 && PVV>0) {
678  k_spt.emplace_back(Kspt);
679  Pdd.emplace_back(PDD*PP0);
680  Pdv.emplace_back(PDV*PP0);
681  Pvv.emplace_back(PVV*PP0);
682  Pb2d.emplace_back(PB2D*PP0);
683  Pb2v.emplace_back(PB2V*PP0);
684  Pb22.emplace_back(PB22*PP0);
685  Pbs2d.emplace_back(PBS2D*PP0);
686  Pbs2v.emplace_back(PBS2V*PP0);
687  Pb2s2.emplace_back(PB2S2*PP0);
688  Pbs22.emplace_back(PBS22*PP0);
689  sigma32Pklin.emplace_back(sigma32PKlin*PP0);
690  Bb1.emplace_back(BB1*PP0);
691  Bb2.emplace_back(BB2*PP0);
692  Bbs2.emplace_back(BBS2*PP0);
693  }
694  finPK.clear();
695 
696  vector<double> Pdd_new(kk.size()), Pdv_new(kk.size()), Pvv_new(kk.size()), Pb2d_new(kk.size()), Pb2v_new(kk.size()), Pb22_new(kk.size()), Pbs2d_new(kk.size()), Pbs2v_new(kk.size()), Pb2s2_new(kk.size()), Pbs22_new(kk.size()), sigma32Pklin_new(kk.size()), Bb1_new(kk.size()), Bb2_new(kk.size()), Bbs2_new(kk.size());
697 
698  glob::FuncGrid interp_Pdd(k_spt, Pdd, "Spline");
699  glob::FuncGrid interp_Pdv(k_spt, Pdv, "Spline");
700  glob::FuncGrid interp_Pvv(k_spt, Pvv, "Spline");
701  glob::FuncGrid interp_Pb2d(k_spt, Pb2d, "Spline");
702  glob::FuncGrid interp_Pb2v(k_spt, Pb2v, "Spline");
703  glob::FuncGrid interp_Pb22(k_spt, Pb22, "Spline");
704  glob::FuncGrid interp_Pbs2d(k_spt, Pbs2d, "Spline");
705  glob::FuncGrid interp_Pbs2v(k_spt, Pbs2v, "Spline");
706  glob::FuncGrid interp_Pb2s2(k_spt, Pb2s2, "Spline");
707  glob::FuncGrid interp_Pbs22(k_spt, Pbs22, "Spline");
708  glob::FuncGrid interp_sigma32Pklin(k_spt, sigma32Pklin, "Spline");
709  glob::FuncGrid interp_Bb1(k_spt, Bb1, "Spline");
710  glob::FuncGrid interp_Bb2(k_spt, Bb2, "Spline");
711  glob::FuncGrid interp_Bbs2(k_spt, Bbs2, "Spline");
712 
713  Pdd_new = interp_Pdd.eval_func(kk);
714  Pdv_new = interp_Pdv.eval_func(kk);
715  Pvv_new = interp_Pvv.eval_func(kk);
716  Pb2d_new = interp_Pb2d.eval_func(kk);
717  Pb2v_new = interp_Pb2v.eval_func(kk);
718  Pb22_new = interp_Pb22.eval_func(kk);
719  Pbs2d_new = interp_Pbs2d.eval_func(kk);
720  Pbs2v_new = interp_Pbs2v.eval_func(kk);
721  Pb2s2_new = interp_Pb2s2.eval_func(kk);
722  Pbs22_new = interp_Pbs22.eval_func(kk);
723  sigma32Pklin_new = interp_sigma32Pklin.eval_func(kk);
724  Bb1_new = interp_Bb1.eval_func(kk);
725  Bb2_new = interp_Bb2.eval_func(kk);
726  Bbs2_new = interp_Bbs2.eval_func(kk);
727 
728  if (system (("rm -rf "+output_tmpCPT).c_str())) {}
729 
730  return {Pdd_new, Pdv_new, Pvv_new, Pb2d_new, Pb2v_new, Pb22_new, Pbs2d_new, Pbs2v_new, Pb2s2_new, Pbs22_new, sigma32Pklin_new, Bb1_new, Bb2_new, Bbs2_new};
731 }
Cosmology.h
The class Cosmology.
cbl::Path
The class Path.
Definition: Path.h:145
cbl::Path::DirCosmo
std::string DirCosmo()
get the directory where the CosmoBolognaLbi are stored
Definition: Path.h:178
cbl::cosmology::Cosmology::Pk_DeltaDelta
double Pk_DeltaDelta(const double kk, const std::shared_ptr< cbl::glob::FuncGrid > Pk, const double qmin, const double qmax, const double prec=1.e-3)
the real-space matter non-linear power spectrum , computed at 1-loop
Definition: PkXiNonLinear.cpp:133
cbl::cosmology::Cosmology::Pk_1loop
double Pk_1loop(const double kk, const std::shared_ptr< cbl::glob::FuncGrid > PkLin, const int corrtype, const double qmin, const double qmax, const double prec=1.e-3)
the one-loop power spectrum
Definition: PkXiNonLinear.cpp:95
cbl::cosmology::Cosmology::F2
double F2(const double k, const double q, const double kq)
function used to estimate the non-linear power spectrum
Definition: PkXiNonLinear.cpp:77
cbl::cosmology::Cosmology::G2
double G2(const double k, const double q, const double kq)
function used to estimate the non-linear power spectrum
Definition: PkXiNonLinear.cpp:86
cbl::cosmology::Cosmology::Pk_TNS_AB_1loop
std::vector< std::vector< double > > Pk_TNS_AB_1loop(std::vector< double > kk, const double mu, const std::string method, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec)
the A and B correction terms for the TNS model at 1-loop from the multipole expansion
Definition: PkXiNonLinear.cpp:305
cbl::cosmology::Cosmology::Pk_TNS_dd_dt_tt
std::vector< std::vector< double > > Pk_TNS_dd_dt_tt(std::vector< double > kk, const std::string method, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min=0.001, const double k_max=100., const double prec=1.e-2)
the non-linear , , matter power spectra
Definition: PkXiNonLinear.cpp:479
cbl::cosmology::Cosmology::Pk_ThetaTheta
double Pk_ThetaTheta(const double kk, const std::shared_ptr< cbl::glob::FuncGrid > Pk, const double qmin, const double qmax, const double prec=1.e-3)
the real-space matter non-linear power spectrum , computed at 1-loop
Definition: PkXiNonLinear.cpp:153
cbl::cosmology::Cosmology::Pk_eTNS_terms_1loop
std::vector< std::vector< double > > Pk_eTNS_terms_1loop(std::vector< double > kk, const std::string method, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min=0.001, const double k_max=100., const double prec=1.e-2)
The expanded correction terms for the extended TNS model (eTNS)
Definition: PkXiNonLinear.cpp:573
cbl::cosmology::Cosmology::Pk_TNS_AB_terms_1loop
std::vector< std::vector< double > > Pk_TNS_AB_terms_1loop(std::vector< double > kk, const std::string method, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min=0.001, const double k_max=100., const double prec=1.e-2)
the expanded A and B correction terms for the TNS model
Definition: PkXiNonLinear.cpp:322
cbl::cosmology::Cosmology::f_k
double f_k(const double k, const std::shared_ptr< cbl::glob::FuncGrid > PkLin, const double qmin, const double qmax, const double prec=1.e-3)
function used to estimate the non-linear power spectrum
Definition: PkXiNonLinear.cpp:47
cbl::cosmology::Cosmology::Pk_TNS_AB_multipoles
std::vector< std::vector< double > > Pk_TNS_AB_multipoles(std::vector< double > kk, const std::string method, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec)
the multipoles of the A and B correction terms for the TNS model
Definition: PkXiNonLinear.cpp:224
cbl::cosmology::Cosmology::Pk_DeltaTheta
double Pk_DeltaTheta(const double kk, const std::shared_ptr< cbl::glob::FuncGrid > Pk, const double qmin, const double qmax, const double prec=1.e-3)
the real-space matter non-linear power spectrum , computed at 1-loop
Definition: PkXiNonLinear.cpp:143
cbl::cosmology::Cosmology::g_k
double g_k(const double k, const std::shared_ptr< cbl::glob::FuncGrid > PkLin, const double qmin, const double qmax, const double prec=1.e-3)
function used to estimate the non-linear power spectrum
Definition: PkXiNonLinear.cpp:62
cbl::glob::FuncGrid
The class FuncGrid.
Definition: FuncGrid.h:55
cbl::glob::FuncGrid::eval_func
std::vector< double > eval_func(const std::vector< double > xx) const
evaluate the function at the xx points
Definition: FuncGrid.cpp:153
cbl::par::fDP6
static const char fDP6[]
conversion symbol for: double -> std::string
Definition: Constants.h:145
cbl::par::fDP1
static const char fDP1[]
conversion symbol for: double -> std::string
Definition: Constants.h:130
cbl::par::ee3
static const char ee3[]
conversion symbol for: double -> std::string
Definition: Constants.h:160
cbl::par::fINT
static const char fINT[]
conversion symbol for: int -> std::string
Definition: Constants.h:121
cbl::par::pi
static const double pi
: the ratio of a circle's circumference to its diameter
Definition: Constants.h:199
cbl::par::kB
static const double kB
: the Boltzmann constant, the conversion factor relating temperature and energy [eV K-1]
Definition: Constants.h:224
cbl::par::TCMB
static const double TCMB
: the present day CMB temperature [K]
Definition: Constants.h:275
cbl::modelling::twopt::bias
double bias(const double Mmin, const double sigmalgM, const double M0, const double M1, const double alpha, const std::shared_ptr< void > inputs)
the mean galaxy bias
Definition: ModelFunction_TwoPointCorrelation1D_monopole.cpp:1026
cbl::wrapper::gsl::GSL_integrate_cquad
double GSL_integrate_cquad(gsl_function func, const double a, const double b, const double rel_err=1.e-3, const double abs_err=0, const int nevals=100)
integral, using the gsl cquad method
Definition: GSLwrapper.cpp:159
cbl::wrapper::gsl::GSL_integrate_qag
double GSL_integrate_qag(gsl_function Func, const double a, const double b, const double rel_err=1.e-3, const double abs_err=0, const int limit_size=1000, const int rule=6)
integral, computed using the GSL qag method
Definition: GSLwrapper.cpp:180
cbl
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
cbl::conv
std::string conv(const T val, const char *fact)
convert a number to a std::string
Definition: Kernel.h:903
cbl::TopHat_WF
T TopHat_WF(const T kR)
the top-hat window function
Definition: Func.h:1195
cbl::ErrorCBL
int ErrorCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL, const cbl::glob::ExitCode exitCode=cbl::glob::ExitCode::_error_)
throw an exception: it is used for handling exceptions inside the CosmoBolognaLib
Definition: Kernel.h:780
cbl::legendre_polynomial
double legendre_polynomial(const double mu, const int l)
the order l Legendre polynomial
Definition: Func.cpp:1786