CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
3PCF.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::denominator_Q (const double r1, const double r2, const double theta, const vector<double> rr, const vector<double> xi_matter) const
48 {
49  const double r3 = sqrt(r1*r1+r2*r2-2.*r1*r2*cos(theta));
50  const glob::FuncGrid interp_xi_matter(rr, xi_matter, "Spline");
51 
52  const double xi1 = interp_xi_matter(r1);
53  const double xi2 = interp_xi_matter(r2);
54  const double xi3 = interp_xi_matter(r3);
55 
56  return xi1*xi2+xi2*xi3+xi3*xi1;
57 }
58 
59 
60 // =====================================================================================
61 
62 
63 void cbl::cosmology::Cosmology::integrals_Q_nonLocal (vector<double> &xi_matter, vector<double> &Phi, const vector<double> rr, const vector<double> kk, const vector<double> Pk_matter, const double prec) const
64 {
65  xi_matter = wrapper::fftlog::transform_FFTlog(rr, 1, kk, Pk_matter, 0, 0, 1, 1);
66 
67  const int nk = kk.size();
68  glob::FuncGrid interp_Pk = glob::FuncGrid(kk, Pk_matter, "Spline", BinType::_logarithmic_);
69 
70  Phi.erase(Phi.begin(), Phi.end());
71  Phi.resize(nk, 0);
72 
73  for (size_t i=0; i<rr.size(); i++) {
74 
75  auto integrand = [&] (const double _k)
76  {
77  const double kr = _k*rr[i];
78  return pow(TopHat_WF(kr), 2)*interp_Pk(_k)*sin(kr)/kr;
79  };
80 
81  Phi[i] = 1./(2.*pow(par::pi, 2))*wrapper::gsl::GSL_integrate_qag(integrand, 0, 100, prec);
82  }
83 
84 }
85 
86 
87 // =====================================================================================
88 
89 
90 double cbl::cosmology::Cosmology::Gamma_3PCF (const double r1, const double r2, const double theta, const vector<double> xi, const vector<double> dPhi) const
91 {
92  return ((xi[0]+3*dPhi[0]/r1)*(xi[1]+3*dPhi[1]/r2))*legendre_polynomial(cos(theta), 2);
93 }
94 
95 
96 // =====================================================================================
97 
98 
99 double cbl::cosmology::Cosmology::Q_nonLocal (const double r1, const double r2, const double theta, std::vector<double> &rr, std::vector<double> &xi_matter, std::vector<double> &Phi, const std::vector<double> kk, const std::vector<double> Pk_matter) const
100 {
101  if (rr.size()==0) {
102  rr = linear_bin_vector(200, 1., 300.);
103  integrals_Q_nonLocal(xi_matter, Phi, rr, kk, Pk_matter, 1.e-3);
104  }
105 
106  glob::FuncGrid interp_xiDM(rr, xi_matter, "Spline");
107  glob::FuncGrid interp_Phi(rr, Phi, "Spline");
108 
109  const double mu12 = cos(theta);
110 
111  const double r3 = sqrt(r1*r1+r2*r2-2.*r1*r2*mu12);
112  double a12 = theta;
113 
114  double a23 = 0.;
115  double a13 = 0.;
116  if (r1<=r2) {
117  a23 = asin(sin(a12)*r1/r3);
118  a13 = par::pi-a12-a23;
119  }
120  else {
121  a13 = asin(sin(a12)*r2/r3);
122  a23 = par::pi-a12-a13;
123  }
124 
125  const double xi1 = interp_xiDM(r1);
126  const double xi2 = interp_xiDM(r2);
127  const double xi3 = interp_xiDM(r3);
128 
129  const double dPhi1 = interp_Phi.D1v(r1);
130  const double dPhi2 = interp_Phi.D1v(r2);
131  const double dPhi3 = interp_Phi.D1v(r3);
132 
133  double gamma = 0.;
134  gamma += Gamma_3PCF(r1, r2, a12, {xi1, xi2}, {dPhi1, dPhi2});
135  gamma += Gamma_3PCF(r2, r3, a23, {xi2, xi3}, {dPhi2, dPhi3});
136  gamma += Gamma_3PCF(r3, r1, a13, {xi3, xi1}, {dPhi3, dPhi1});
137 
138  return 2./3*(gamma/denominator_Q(r1, r2, theta, rr, xi_matter)-1);
139 }
140 
141 
142 // =====================================================================================
143 
144 
145 std::vector<double> cbl::cosmology::Cosmology::Q_nonLocal (const double r1, const double r2, const std::vector<double> theta, const std::vector<double> kk, const std::vector<double> Pk_matter) const
146 {
147  const int ntheta = theta.size();
148  vector<double> rr, xi_matter, Phi;
149  vector<double> qNL(ntheta, 0);
150 
151  for (int i=0; i<ntheta; i++)
152  qNL[i] = Q_nonLocal(r1, r2, theta[i], rr, xi_matter, Phi, kk, Pk_matter);
153 
154  return qNL;
155 }
156 
157 
158 // =====================================================================================
159 
160 
161 void cbl::cosmology::Cosmology::integrals_zeta_Slepian (std::vector<double> &xi_matter, std::vector<double> &xi_matter_m1, std::vector<double> &xi_matter_p1, std::vector<double> &xi_matter_2, const std::vector<double> rr, const std::vector<double> kk, const std::vector<double> Pk_matter) const
162 {
163  vector<double> Pk_matter_m1 = Pk_matter, Pk_matter_p1 = Pk_matter;
164  const int nk = kk.size();
165 
166  for (int i=0; i< nk; i++) {
167  Pk_matter_m1[i] *= pow(kk[i], -1);
168  Pk_matter_p1[i] *= kk[i];
169  }
170 
171  xi_matter = wrapper::fftlog::transform_FFTlog(rr, 1, kk, Pk_matter, 0, 0, 1, 1);
172  xi_matter_m1 = wrapper::fftlog::transform_FFTlog(rr, 1, kk, Pk_matter_m1, 1, 0, 1, 1);
173  xi_matter_p1 = wrapper::fftlog::transform_FFTlog(rr, 1, kk, Pk_matter_p1, 1, 0, 1, 1);
174  xi_matter_2 = wrapper::fftlog::transform_FFTlog(rr, 1, kk, Pk_matter, 2, 0, 1, 1);
175 }
176 
177 
178 // =====================================================================================
179 
180 
181 double cbl::cosmology::Cosmology::zeta_precyclic_Slepian (const double r1, const double r2, const double mu, const double b1, const double b2, const glob::FuncGrid interp_xi_matter, const glob::FuncGrid interp_xi_matter_m1, const glob::FuncGrid interp_xi_matter_p1, const glob::FuncGrid interp_xi_matter_2) const
182 {
183  const double mu12 = mu;
184 
185  if (abs(mu12)>1)
186  return 0;
187 
188  const double r3 = sqrt(r1*r1+r2*r2-2.*r1*r2*mu12);
189  const double a12 = acos(mu12);
190 
191  double a23 = 0.;
192  double a13 = 0.;
193  if (r1<=r2) {
194  a23 = asin(sin(a12)*r1/r3);
195  a13 = par::pi-a12-a23;
196  }
197  else {
198  a13 = asin(sin(a12)*r2/r3);
199  a23 = par::pi-a12-a13;
200  }
201 
202  const double mu13 = cos(a13);
203  const double mu23 = cos(a23);
204 
205  const double b1_3 = b1*b1*b1;
206  const double zeta_pc_0_12 = (2.*b1*b1*b2+(34./21)*b1_3)*interp_xi_matter(r1)*interp_xi_matter(r2);
207  const double zeta_pc_1_12 = -b1_3*(interp_xi_matter_m1(r1)*interp_xi_matter_p1(r2)+interp_xi_matter_m1(r2)*interp_xi_matter_p1(r1))*legendre_polynomial(mu12, 1);
208  const double zeta_pc_2_12 = (8./21)*b1_3*interp_xi_matter_2(r1)*interp_xi_matter_2(r2)*legendre_polynomial(mu12, 2);
209 
210  const double zeta_pc_0_13 = (2.*b1*b1*b2+(34./21)*b1_3)*interp_xi_matter(r1)*interp_xi_matter(r3);
211  const double zeta_pc_1_13 = -b1_3*(interp_xi_matter_m1(r1)*interp_xi_matter_p1(r3)+interp_xi_matter_m1(r3)*interp_xi_matter_p1(r1))*legendre_polynomial(mu13, 1);
212  const double zeta_pc_2_13 = 8./21*b1_3*interp_xi_matter_2(r1)*interp_xi_matter_2(r3)*legendre_polynomial(mu13, 2);
213 
214  const double zeta_pc_0_23 = (2.*b1*b1*b2+(34./21)*b1_3)*interp_xi_matter(r2)*interp_xi_matter(r3);
215  const double zeta_pc_1_23 = -b1_3*(interp_xi_matter_m1(r2)*interp_xi_matter_p1(r3)+interp_xi_matter_m1(r3)*interp_xi_matter_p1(r2))*legendre_polynomial(mu23, 1);
216  const double zeta_pc_2_23 = (8./21)*b1_3*interp_xi_matter_2(r2)*interp_xi_matter_2(r3)*legendre_polynomial(mu23, 2);
217 
218  return zeta_pc_0_12+zeta_pc_1_12+zeta_pc_2_12+zeta_pc_0_13+zeta_pc_1_13+zeta_pc_2_13+zeta_pc_0_23+zeta_pc_1_23+zeta_pc_2_23;
219 }
220 
221 
222 // =====================================================================================
223 
224 
225 double cbl::cosmology::Cosmology::zeta_precyclic_Slepian (const double r1, const double r2, const double r3, const double deltaR, const double b1, const double b2, const glob::FuncGrid interp_xi_matter, const glob::FuncGrid interp_xi_matter_m1, const glob::FuncGrid interp_xi_matter_p1, const glob::FuncGrid interp_xi_matter_2) const
226 {
227  double r1Min = r1-0.5*deltaR;
228  double r1Max = r1+0.5*deltaR;
229  double r2Min = r2-0.5*deltaR;
230  double r2Max = r2+0.5*deltaR;
231  double r3Min = r3-0.5*deltaR;
232  double r3Max = r3+0.5*deltaR;
233 
234  // Numerator
235  auto integrandNum_r1 = [&] (const double _r1) {
236 
237  auto integrandNum_r2 = [&] (const double _r2) {
238 
239  double a = max(r3Min, _r2-_r1);
240  double b = min(r3Max, _r2+_r1);
241 
242  if (a>b)
243  return 0.;
244  else {
245  auto integrandNum_r3 = [&] (const double _r3) {
246  double mu = (_r1*_r1+_r2*_r2-_r3*_r3)/(2*_r1*_r2);
247  return zeta_precyclic_Slepian (_r1, _r2, mu, b1, b2, interp_xi_matter, interp_xi_matter_m1, interp_xi_matter_p1, interp_xi_matter_2);
248  };
249 
250  return wrapper::gsl::GSL_integrate_cquad(integrandNum_r3, a, b, 1.e-4);
251  }
252  };
253 
254  return wrapper::gsl::GSL_integrate_cquad(integrandNum_r2, r2Min, r2Max, 1.e-4);
255  };
256 
257  double num = wrapper::gsl::GSL_integrate_cquad(integrandNum_r1, r1Min, r1Max, 1.e-4);
258 
259  // Denominator
260  auto integrandDen_r1 = [&] (const double _r1) {
261 
262  auto integrandDen_r2 = [&] (const double _r2) {
263 
264  double a = max(r3Min, _r2-_r1);
265  double b = min(r3Max, _r2+_r1);
266 
267  if (a>b)
268  return 0.;
269  else {
270  auto integrandDen_r3 = [&] (const double _r3) {
271  (void)_r3;
272  return 1.;
273  };
274 
275  return wrapper::gsl::GSL_integrate_cquad(integrandDen_r3, a, b, 1.e-4);
276  }
277  };
278 
279  return wrapper::gsl::GSL_integrate_cquad(integrandDen_r2, r2Min, r2Max, 1.e-4);
280  };
281 
282  double den = wrapper::gsl::GSL_integrate_cquad(integrandDen_r1, r1Min, r1Max, 1.e-4);
283 
284  return num/den;
285 }
286 
287 
288 // =====================================================================================
289 
290 
291 std::vector<double> cbl::cosmology::Cosmology::zeta_expansion_Slepian (const double r1, const double r2, const double b1, const double b2, std::vector<double> &rr, std::vector<double> &xi_matter, std::vector<double> &xi_matter_m1, std::vector<double> &xi_matter_p1, std::vector<double> &xi_matter_2, const int norders, const double prec) const
292 {
293  glob::FuncGrid interp_xi_matter(rr, xi_matter, "Spline");
294  glob::FuncGrid interp_xi_matter_m1(rr, xi_matter_m1, "Spline");
295  glob::FuncGrid interp_xi_matter_p1(rr, xi_matter_p1, "Spline");
296  glob::FuncGrid interp_xi_matter_2(rr, xi_matter_2, "Spline");
297 
298  vector<double> zeta_r1_r2(norders, 0);
299 
300  for (int i=0; i<norders; i++) {
301  auto integrand = [&] (const double mu12) {
302  return Cosmology::zeta_precyclic_Slepian(r1, r2, mu12, b1, b2, interp_xi_matter, interp_xi_matter_m1, interp_xi_matter_p1, interp_xi_matter_2)*legendre_polynomial (mu12, i);
303  };
304  zeta_r1_r2[i] = 0.5*(2*i+1)*wrapper::gsl::GSL_integrate_qag(integrand, -1, 1, prec);
305  }
306 
307  return zeta_r1_r2;
308 }
309 
310 
311 // =====================================================================================
312 
313 
314 double cbl::cosmology::Cosmology::zeta_DM_Slepian (const double r1, const double r2, const double theta, std::vector<double> &rr, std::vector<double> &xi_matter, std::vector<double> &xi_matter_m1, std::vector<double> &xi_matter_p1, std::vector<double> &xi_matter_2, const std::vector<double> kk, const std::vector<double> Pk_matter, const int norders, const double prec) const
315 {
316  if (rr.size()==0) {
317  rr = linear_bin_vector(200, 1., 300.);
318  integrals_zeta_Slepian(xi_matter, xi_matter_m1, xi_matter_p1, xi_matter_2, rr, kk, Pk_matter);
319  }
320 
321  const double mu = cos(theta);
322  vector<double> z_r1_r2 = Cosmology::zeta_expansion_Slepian(r1, r2, 1, 0, rr, xi_matter, xi_matter_m1, xi_matter_p1, xi_matter_2, norders, prec);
323 
324  double zeta_r1_r2_theta = 0.;
325  for (size_t i=0; i<z_r1_r2.size(); i++)
326  zeta_r1_r2_theta += z_r1_r2[i]*legendre_polynomial(mu, i);
327 
328  return zeta_r1_r2_theta;
329 }
330 
331 
332 // =====================================================================================
333 
334 
335 double cbl::cosmology::Cosmology::Q_DM_Slepian (const double r1, const double r2, const double theta, std::vector<double> &rr, std::vector<double> &xi_matter, std::vector<double> &xi_matter_m1, std::vector<double> &xi_matter_p1, std::vector<double> &xi_matter_2, const std::vector<double> kk, const std::vector<double> Pk_matter, const int norders, const double prec) const
336 {
337  const double zeta_DM = zeta_DM_Slepian(r1, r2, theta, rr, xi_matter, xi_matter_m1, xi_matter_p1, xi_matter_2, kk, Pk_matter, norders, prec);
338 
339  return zeta_DM/denominator_Q(r1, r2, theta, rr, xi_matter);
340 }
341 
342 
343 // =====================================================================================
344 
345 
346 void cbl::cosmology::Cosmology::integrals_zeta_BarrigaGatzanaga (std::vector<double> &xi_matter, std::vector<double> &Phi, const std::vector<double> rr, const std::vector<double> kk, const std::vector<double> Pk_matter) const
347 {
348  const int nk = kk.size();
349 
350  vector<double> Pk_matter_m4 = Pk_matter;
351  for (int i=0; i< nk; i++)
352  Pk_matter_m4[i] *= pow(kk[i], -2);
353 
354  xi_matter = wrapper::fftlog::transform_FFTlog(rr, 1, kk, Pk_matter, 0, 0, 1, 1);
355  Phi = wrapper::fftlog::transform_FFTlog(rr, 1, kk, Pk_matter_m4, 0, 0, 1, 1);
356 }
357 
358 
359 // =====================================================================================
360 
361 
362 double cbl::cosmology::Cosmology::zeta_single_BarrigaGatzanaga (const double r1, const double r2, const double theta, const std::vector<double> xi, const std::vector<double> dxi, const std::vector<double> dPhi) const
363 {
364  const double mu = cos(theta);
365 
366  const double xi1 = xi[0];
367  const double xi2 = xi[1];
368 
369  const double dxi1 = dxi[0];
370  const double dxi2 = dxi[1];
371 
372  const double dPhi1 = dPhi[0];
373  const double dPhi2 = dPhi[1];
374 
375  const double t1 = (10./7)*xi1*xi2;
376  const double t2 = -3*dPhi1*dPhi2/(r1*r2)-xi1*dPhi2/r2-xi2*dPhi1/r1;
377  const double t3 = mu*mu*( (xi1+3*dPhi1/r1)*(xi2+3*dPhi2/r2));
378  const double t4 = -mu*(dxi1*dPhi2+dxi2*dPhi1);
379 
380  return t1+(4./7)*(t2+t3)+t4;
381 }
382 
383 
384 // =====================================================================================
385 
386 
387 double cbl::cosmology::Cosmology::zeta_DM_BarrigaGatzanaga (const double r1, const double r2, const double theta, std::vector<double> &rr, std::vector<double> &xi_matter, std::vector<double> &Phi, const std::vector<double> kk, const std::vector<double> Pk_matter) const
388 {
389  if (rr.size()==0) {
390  rr = linear_bin_vector(200, 1., 300.);
391  integrals_zeta_BarrigaGatzanaga (xi_matter, Phi, rr, kk, Pk_matter);
392  }
393 
394  glob::FuncGrid interp_xiDM(rr, xi_matter, "Spline");
395  glob::FuncGrid interp_Phi(rr, Phi, "Spline");
396 
397  const double mu12 = cos(theta);
398 
399  const double r3 = sqrt(r1*r1+r2*r2-2.*r1*r2*mu12);
400  const double a12 = theta;
401 
402  double a23 = 0.;
403  double a13 = 0.;
404  if (r1<=r2) {
405  a23 = asin(sin(a12)*r1/r3);
406  a13 = par::pi-a12-a23;
407  }
408  else {
409  a13 = asin(sin(a12)*r2/r3);
410  a23 = par::pi-a12-a13;
411  }
412 
413  const double xi1 = interp_xiDM(r1);
414  const double xi2 = interp_xiDM(r2);
415  const double xi3 = interp_xiDM(r3);
416 
417  const double dxi1 = interp_xiDM.D1v(r1);
418  const double dxi2 = interp_xiDM.D1v(r2);
419  const double dxi3 = interp_xiDM.D1v(r3);
420 
421  const double dPhi1 = interp_Phi.D1v(r1);
422  const double dPhi2 = interp_Phi.D1v(r2);
423  const double dPhi3 = interp_Phi.D1v(r3);
424 
425  const double t1 = zeta_single_BarrigaGatzanaga(r1, r2, a12, {xi1, xi2}, {dxi1, dxi2}, {dPhi1, dPhi2});
426  const double t2 = zeta_single_BarrigaGatzanaga(r2, r3, a23, {xi2, xi3}, {dxi2, dxi3}, {dPhi2, dPhi3});
427  const double t3 = zeta_single_BarrigaGatzanaga(r3, r1, a13, {xi3, xi1}, {dxi3, dxi1}, {dPhi3, dPhi1});
428 
429  const double zeta = t1+t2+t3;
430  return zeta;
431 }
432 
433 
434 // =====================================================================================
435 
436 
437 double cbl::cosmology::Cosmology::Q_DM_BarrigaGatzanaga (const double r1, const double r2, const double theta, std::vector<double> &rr, std::vector<double> &xi_matter, std::vector<double> &Phi, const std::vector<double> kk, const std::vector<double> Pk_matter) const
438 {
439  return zeta_DM_BarrigaGatzanaga(r1, r2, theta, rr, xi_matter, Phi, kk, Pk_matter)/denominator_Q(r1, r2, theta, rr, xi_matter);
440 }
441 
442 
443 // =====================================================================================
444 
445 
446 std::vector<double> cbl::cosmology::Cosmology::zeta_DM (const double r1, const double r2, const std::vector<double> theta, const string model, const std::vector<double> kk, const std::vector<double> Pk_matter) const
447 {
448  const int ntheta = theta.size();
449  vector<double> rr, xi_matter;
450  vector<double> zDM(ntheta, 0.);
451 
452  if (model=="Slepian") {
453  vector<double> xi_matter_m1, xi_matter_p1, xi_matter_2;
454  for (int i=0; i<ntheta; i++)
455  zDM[i] = zeta_DM_Slepian(r1, r2, theta[i], rr, xi_matter, xi_matter_m1, xi_matter_p1, xi_matter_2, kk, Pk_matter, 9, 1.e-3);
456  }
457  else if (model=="BarrigaGatzanaga") {
458  vector<double> Phi;
459  for (int i=0; i<ntheta; i++)
460  zDM[i] = zeta_DM_BarrigaGatzanaga(r1, r2, theta[i], rr, xi_matter, Phi, kk, Pk_matter);
461  }
462  else
463  ErrorCBL("the chosen model is not implemented!", "zeta_DM", "3PCF.cpp");
464 
465  return zDM;
466 }
467 
468 // =====================================================================================
469 
470 
471 std::vector<double> cbl::cosmology::Cosmology::Q_DM (const double r1, const double r2, const std::vector<double> theta, const string model, const std::vector<double> kk, const std::vector<double> Pk_matter) const
472 {
473  const int ntheta = theta.size();
474  vector<double> rr, xi_matter;
475  vector<double> qDM(ntheta, 0);
476 
477  if (model=="Slepian") {
478  vector<double> xi_matter_m1, xi_matter_p1, xi_matter_2;
479  for (int i=0; i<ntheta; i++)
480  qDM[i] = Q_DM_Slepian(r1, r2, theta[i], rr, xi_matter, xi_matter_m1, xi_matter_p1, xi_matter_2, kk, Pk_matter, 9, 1.e-3);
481  }
482  else if (model=="BarrigaGatzanaga")
483  {
484  vector<double> Phi;
485  for (int i=0; i<ntheta; i++)
486  qDM[i] = Q_DM_BarrigaGatzanaga(r1, r2, theta[i], rr, xi_matter, Phi, kk, Pk_matter);
487  }
488  else
489  ErrorCBL("the chosen model is not implemented!", "Q_DM", "3PCF.cpp");
490 
491  return qDM;
492 
493 }
494 
495 
496 // =====================================================================================
497 
498 
499 std::vector<double> cbl::cosmology::Cosmology::zeta_halo (const double r1, const double r2, const std::vector<double> theta, const double b1, const double b2, const string model, const std::vector<double> kk, const std::vector<double> Pk_matter) const
500 {
501  const int ntheta = theta.size();
502  vector<double> rr, xi_matter;
503  vector<double> zH(ntheta, 0);
504 
505  if (model=="Slepian") {
506  vector<double> xi_matter_m1, xi_matter_p1, xi_matter_2;
507  for (int i=0; i<ntheta; i++)
508  zH[i] = b1*b1*b1*zeta_DM_Slepian(r1, r2, theta[i], rr, xi_matter, xi_matter_m1, xi_matter_p1, xi_matter_2, kk, Pk_matter, 9, 1.e-3)+b1*b1*b2*denominator_Q(r1, r2, theta[i], rr, xi_matter);
509  }
510  else if (model=="BarrigaGatzanaga")
511  {
512  vector<double> Phi;
513  for (int i=0; i<ntheta; i++)
514  zH[i] = b1*b1*b1*zeta_DM_BarrigaGatzanaga(r1, r2, theta[i], rr, xi_matter, Phi, kk, Pk_matter)+b1*b1*b2*denominator_Q(r1, r2, theta[i], rr, xi_matter);
515  }
516  else
517  ErrorCBL("the chosen model is not implemented!", "z_halo", "3PCF.cpp");
518 
519  return zH;
520 }
521 
522 
523 // =====================================================================================
524 
525 
526 std::vector<double> cbl::cosmology::Cosmology::Q_halo (const double r1, const double r2, const std::vector<double> theta, const double b1, const double b2, const std::string model, const std::vector<double> kk, const std::vector<double> Pk_matter) const
527 {
528  const int ntheta = theta.size();
529  vector<double> qDM = Cosmology::Q_DM(r1, r2, theta, model, kk, Pk_matter);
530  vector<double> qH(ntheta, 0);
531 
532  for (int i=0; i<ntheta; i++)
533  qH[i] = qDM[i]/b1+b2/(b1*b1);
534 
535  return qH;
536 }
537 
538 
539 // =====================================================================================
540 
541 
542 std::vector<double> cbl::cosmology::Cosmology::Q_halo (const double r1, const double r2, const std::vector<double> theta, const double b1, const double b2, const double g2, const std::string model, const std::vector<double> kk, const std::vector<double> Pk_matter) const
543 {
544  const int ntheta = theta.size();
545  vector<double> qH = Cosmology::Q_halo(r1, r2, theta, b1, b2, model, kk, Pk_matter);
546  vector<double> qNL = Cosmology::Q_nonLocal(r1, r2, theta, kk, Pk_matter);
547 
548  for (int i=0; i<ntheta; i++)
549  qH[i] += g2/b1*qNL[i];
550 
551  return qH;
552 }
553 
554 
555 // =====================================================================================
556 
557 
558 std::vector<double> cbl::cosmology::Cosmology::zeta_DM_eq (const std::vector<double> rr, const std::string model, const std::vector<double> kk, const std::vector<double> Pk_matter) const
559 {
560  const int nr = rr.size();
561  const double theta = par::pi/3;
562  vector<double> _rr, xi_matter;
563  vector<double> zDM(nr, 0);
564 
565  if (model=="Slepian") {
566  vector<double> xi_matter_m1, xi_matter_p1, xi_matter_2;
567  for (int i=0; i<nr; i++)
568  zDM[i] = zeta_DM_Slepian(rr[i], rr[i], theta, _rr, xi_matter, xi_matter_m1, xi_matter_p1, xi_matter_2, kk, Pk_matter, 9, 1.e-3);
569  }
570  else if (model=="BarrigaGatzanaga")
571  {
572  vector<double> Phi;
573  for (int i=0; i<nr; i++)
574  zDM[i] = zeta_DM_BarrigaGatzanaga(rr[i], rr[i], theta, _rr, xi_matter, Phi, kk, Pk_matter);
575  }
576  else
577  ErrorCBL("the chosen model is not implemented!", "zeta_DM_eq", "3PCF.cpp");
578 
579  return zDM;
580 }
581 
582 
583 // =====================================================================================
584 
585 
586 std::vector<double> cbl::cosmology::Cosmology::Q_DM_eq (const std::vector<double> rr, const std::string model, const std::vector<double> kk, const std::vector<double> Pk_matter) const
587 {
588  const int nr = rr.size();
589  const double theta = par::pi/3;
590  vector<double> _rr, xi_matter;
591  vector<double> Q_DM(nr, 0);
592 
593  if (model=="Slepian") {
594  vector<double> xi_matter_m1, xi_matter_p1, xi_matter_2;
595  for (int i=0; i<nr; i++)
596  Q_DM[i] = Q_DM_Slepian(rr[i], rr[i], theta, _rr, xi_matter, xi_matter_m1, xi_matter_p1, xi_matter_2, kk, Pk_matter, 9, 1.e-3);
597  }
598  else if (model=="BarrigaGatzanaga")
599  {
600  vector<double> Phi;
601  for (int i=0; i<nr; i++)
602  Q_DM[i] = Q_DM_BarrigaGatzanaga(rr[i], rr[i], theta, _rr, xi_matter, Phi, kk, Pk_matter);
603  }
604  else
605  ErrorCBL("the chosen model is not implemented!", "Q_DM_eq", "3PCF.cpp");
606 
607  return Q_DM;
608 }
609 
610 
611 // =====================================================================================
612 
613 
614 double cbl::cosmology::Cosmology::zeta_multipoles_covariance (const double Volume, const double nObjects, const int l, const int l_prime, const double r1, const double r2, const double r1_prime, const double r2_prime, const double deltaR, const std::vector<double> kk, const std::vector<double> Pk, const std::vector<double> rr, const std::vector<double> Xi, const double prec)
615 {
616  (void)prec;
617 
618  const double kr0 = 1.; //check
619 
620  size_t nk = static_cast<int>(kk.size());
621  size_t nr = static_cast<int>(rr.size());
622  double dlogr = log(rr[1])-log(rr[0]);
623 
624  const double density = nObjects/Volume;
625 
626  // r binning
627 
628  function<double(double, int, double)> func1;
629  function<double(double, int, int, double, double)> func2;
630  function<double(double, double)> func1_noise;
631 
632  if (deltaR<0) {
633  func1 = [&] (const double kk, const int l, const double rr) {
634  return jl(kk*rr, l); };
635  func2 = [&] (const double kk, const int l, const int lp, const double rr, const double rrp) {
636  return jl(kk*rr, l)*jl(kk*rrp, lp); };
637  } else {
638  func1 = [&] (const double kk, const int l, const double rr) {
639  return jl_distance_average(kk, l, rr-deltaR*0.5, rr+deltaR*0.5); };
640  func2 = [&] (const double kk, const int l, const int lp, const double rr, const double rrp) {
641  return jl_distance_average(kk, l, rr-deltaR*0.5, rr+deltaR*0.5)*jl_distance_average(kk, lp, rrp-deltaR*0.5, rrp+deltaR*0.5); };
642  }
643 
644  func1_noise = [&] (const double rr, const double bin) {
645  double rmin = bin-deltaR/2;
646  double rmax = bin+deltaR/2;
647  if (rr>=rmin && rr<=rmax && deltaR>0)
648  return 1./(4*par::pi/3*(pow(rmax, 3)-pow(rmin, 3)));
649  else
650  return 0.;
651  };
652 
653  // integrand
654 
655  // Terms of type I_ll = int [(Pk+1/n) jl] jl
656 
657  vector<double> pk_r1 = Pk, pk_r2 = Pk;
658  vector<double> pk_r1p = Pk, pk_r2p = Pk;
659 
660  for (size_t i=0; i<nk; i++) {
661  pk_r1[i] *= func1(kk[i], l, r1);
662  pk_r2[i] *= func1(kk[i], l, r2);
663  pk_r1p[i] *= func1(kk[i], l_prime, r1_prime);
664  pk_r2p[i] *= func1(kk[i], l_prime, r2_prime);
665  }
666 
667  vector<double> I1_r_r1 = wrapper::fftlog::transform_FFTlog(rr, 1, kk, pk_r1, l, 0, kr0, 1);
668  vector<double> I1_r_r2 = wrapper::fftlog::transform_FFTlog(rr, 1, kk, pk_r2, l, 0, kr0, 1);
669  vector<double> I1_r_r1p = wrapper::fftlog::transform_FFTlog(rr, 1, kk, pk_r1p, l_prime, 0, kr0, 1);
670  vector<double> I1_r_r2p = wrapper::fftlog::transform_FFTlog(rr, 1, kk, pk_r2p, l_prime, 0, kr0, 1);
671 
672  for (size_t i=0; i<nr; i++) {
673  I1_r_r1[i] += func1_noise(rr[i], r1)/density;
674  I1_r_r2[i] += func1_noise(rr[i], r2)/density;
675  I1_r_r1p[i] += func1_noise(rr[i], r1_prime)/density;
676  I1_r_r2p[i] += func1_noise(rr[i], r2_prime)/density;
677  }
678 
679  // l2
680  vector<int> l2;
681  for (int ll=abs(l-l_prime); ll<l+l_prime+1; ll++)
682  l2.push_back(ll);
683 
684  const size_t n_l2 = l2.size();
685 
686  // Terms of type I_l1_l2_l = int [(Pk+1/n) jl1 jl2] jl
687 
688  vector<double> pk_r1_r1p = Pk, pk_r2_r2p = Pk, pk_r2_r1p = Pk, pk_r1_r2p = Pk;
689  for (size_t i=0; i<nk; i++) {
690  pk_r1_r1p[i] = (pk_r1_r1p[i]+1./density)*func2(kk[i], l, l_prime, r1, r1_prime);
691  pk_r2_r2p[i] = (pk_r2_r2p[i]+1./density)*func2(kk[i], l, l_prime, r2, r2_prime);
692  pk_r2_r1p[i] = (pk_r2_r1p[i]+1./density)*func2(kk[i], l, l_prime, r2, r1_prime);
693  pk_r1_r2p[i] = (pk_r1_r2p[i]+1./density)*func2(kk[i], l, l_prime, r1, r2_prime);
694  }
695 
696  vector<vector<double>> I2_r1_r1p(n_l2, vector<double>(nr, 0.));
697  vector<vector<double>> I2_r2_r2p(n_l2, vector<double>(nr, 0.));
698  vector<vector<double>> I2_r2_r1p(n_l2, vector<double>(nr, 0.));
699  vector<vector<double>> I2_r1_r2p(n_l2, vector<double>(nr, 0.));
700 
701  for (size_t ll=0; ll<n_l2; ll++) {
702  double wig = gsl_sf_coupling_3j(2*l, 2*l_prime, 2*l2[ll], 0, 0, 0);
703  if (wig!=0) {
704  I2_r1_r1p[ll] = wrapper::fftlog::transform_FFTlog(rr, 1, kk, pk_r1_r1p, l2[ll], 0, kr0, 1);
705  I2_r2_r2p[ll] = wrapper::fftlog::transform_FFTlog(rr, 1, kk, pk_r2_r2p, l2[ll], 0, kr0, 1);
706  I2_r2_r1p[ll] = wrapper::fftlog::transform_FFTlog(rr, 1, kk, pk_r2_r1p, l2[ll], 0, kr0, 1);
707  I2_r1_r2p[ll] = wrapper::fftlog::transform_FFTlog(rr, 1, kk, pk_r1_r2p, l2[ll], 0, kr0, 1);
708  }
709  }
710 
711  // integrand
712 
713  double Int = 0;
714 
715  for (size_t i=0; i<nr; i++) {
716  double sum = 0;
717 
718  double Xi_r = Xi[i];
719 
720  double f_r_r1 = I1_r_r1[i];
721  double f_r_r2 = I1_r_r2[i];
722  double f_r_r1p = I1_r_r1p[i];
723  double f_r_r2p = I1_r_r2p[i];
724 
725  for (size_t ll=0; ll<n_l2; ll++) {
726  int ell2 = l2[ll];
727  double wig = gsl_sf_coupling_3j(2*l, 2*l_prime, 2*ell2, 0, 0, 0);
728  if (wig!=0) {
729  double t1 = (2*ell2+1)*pow(wig,2);
730  double f_r_r1_r1p = I2_r1_r1p[ll][i];
731  double f_r_r2_r2p = I2_r2_r2p[ll][i];
732  double f_r_r2_r1p = I2_r2_r1p[ll][i];
733  double f_r_r1_r2p = I2_r1_r2p[ll][i];
734 
735  double t2 = pow(-1, l2[ll])*Xi_r*(f_r_r1_r1p*f_r_r2_r2p+f_r_r2_r1p*f_r_r1_r2p);
736  double t3 = pow(-1, 0.5*(l+l_prime+ell2))*(f_r_r1*f_r_r1p*f_r_r2_r2p+f_r_r1*f_r_r2p*f_r_r2_r1p+f_r_r2*f_r_r1p*f_r_r1_r2p+f_r_r2*f_r_r2p*f_r_r1_r1p);
737  sum += t1*(t2+t3);
738  }
739  }
740 
741  Int += rr[i]*rr[i]*sum*rr[i];
742  }
743 
744  double fact = 4.*par::pi/Volume*(2*l+1)*(2*l_prime+1)*pow(-1, l+l_prime);
745 
746  return fact*Int*dlogr;
747 }
748 
749 
750 // =====================================================================================
751 
752 
753 std::vector<std::vector<double>> cbl::cosmology::Cosmology::zeta_covariance (const double Volume, const double nObjects, const std::vector<double> theta, const double r1, const double r2, const double deltaR, const std::vector<double> kk, const std::vector<double> Pk, const int norders, const double prec, const bool method, const int nExtractions, std::vector<double> mean, const int seed)
754 {
755  (void)method;
756  (void)nExtractions;
757  (void)mean;
758  (void)seed;
759 
760  vector<double> rr, Xi;
761  wrapper::fftlog::transform_FFTlog(rr, Xi, 1, kk, Pk, 0);
762 
763  vector<vector<double>> zeta_l1l2_covariance(norders, vector<double>(norders, 0.));
764 
765  for (int i=0; i<norders; i++)
766  for (int j=i; j<norders; j++) {
767  zeta_l1l2_covariance[i][j] = zeta_multipoles_covariance(Volume, nObjects, i, j, r1, r2, r1, r2, deltaR, kk, Pk, rr, Xi, prec);
768  zeta_l1l2_covariance[j][i] = zeta_l1l2_covariance[i][j];
769  }
770 
771  const int ntheta = int(theta.size());
772 
773  vector<vector<double>> Pl_theta(ntheta, vector<double>(norders, 0));
774  for (int i=0; i<ntheta; i++) {
775  for (int j=0; j<norders; j++)
776  Pl_theta[i][j] = legendre_polynomial (cos(theta[i]), j);
777  }
778 
779  vector<vector<double>> zeta_covariance(ntheta, vector<double>(ntheta, 0));
780  for (int i=0; i<ntheta; i++)
781  for (int j=0; j<ntheta; j++)
782  for (int l1=0; l1<norders; l1++)
783  for (int l2=0; l2<norders; l2++)
784  zeta_covariance[i][j] += zeta_l1l2_covariance[l1][l2]*Pl_theta[i][l1]*Pl_theta[j][l2];
785 
786  return zeta_covariance;
787 }
788 
789 
790 // =====================================================================================a
791 
792 
793 void cbl::cosmology::Cosmology::xi_r_n (std::vector<double> &xi_n, const std::vector<double> rr, const int nn, const std::vector<double> kk, const std::vector<double> Pk)
794 {
795  xi_n = wrapper::fftlog::transform_FFTlog (rr, 1, kk, Pk, nn, 0, par::pi, 1);
796 }
797 
798 
799 // =====================================================================================
800 
801 
802 void cbl::cosmology::Cosmology::xi_r_n_pm (std::vector<double> &xi_n_p, std::vector<double> &xi_n_m, const std::vector<double> rr, const int nn, const std::vector<double> kk, const std::vector<double> Pk)
803 {
804  vector<double> pk_p(Pk.size(), 0), pk_m(Pk.size(), 0);
805 
806  for (size_t i=0; i<Pk.size(); i++){
807  pk_p[i] = kk[i]*Pk[i];
808  pk_m[i] = Pk[i]/kk[i];
809  }
810 
811  xi_n_p = wrapper::fftlog::transform_FFTlog (rr, 1, kk, pk_p, nn, 0, par::pi, 1);
812  xi_n_m = wrapper::fftlog::transform_FFTlog (rr, 1, kk, pk_m, nn, 0, par::pi, 1);
813 }
814 
815 
816 // =====================================================================================
817 
818 
819 void cbl::cosmology::Cosmology::eff_l_l1 (std::vector<std::vector<double>> &eff, const std::vector<double> rr, const int l, const int l1, const std::vector<double> kk, const std::vector<double> Pk)
820 {
821  double min_rr = Min(rr);
822  double max_rr = Max(rr);
823  vector<double> new_r = linear_bin_vector(rr.size(), min_rr, max_rr);
824  eff.resize(rr.size());
825 
826  for (size_t i=0; i<rr.size(); i++)
827  {
828  vector<double> _pk(Pk.size(), 0);
829 
830  for (size_t j=0; j<Pk.size(); j++)
831  _pk[j] = kk[j]*Pk[j]*jl(kk[j]*rr[i], l);
832 
833  eff[i] = wrapper::fftlog::transform_FFTlog (new_r, 1, kk, _pk, l1, 0, par::pi, 1);
834  }
835 
836 }
837 
838 
839 // =====================================================================================
840 
841 
842 void cbl::cosmology::Cosmology::I_ELL_ell (std::vector<std::vector<double>> &II, const std::vector<double> rr, const int ll, const int LL, const std::vector<double> kk, const std::vector<double> Pk)
843 {
844  II.resize(rr.size(), vector<double>(rr.size(), 0));
845  double min_rr = Min(rr);
846  double max_rr = Max(rr);
847  vector<double> new_r = linear_bin_vector(rr.size(), min_rr, max_rr);
848 
849  for (int l1 = 0; l1<=LL+ll; l1++)
850  {
851  if ( (LL>= fabs(l1-ll)) && (LL <= l1+ll)){
852  double fact = pow(-1., l1+ll)*(2.*l1+1)*(2.*ll+1)*pow(gsl_sf_coupling_3j(2*l1, 2*ll, 2*LL, 0, 0, 0),2);
853 
854  if(fact!=0){
855  vector<vector<double>> eff;
856  eff_l_l1 (eff, rr, ll, l1, kk, Pk);
857  for (size_t r1=0; r1<rr.size(); r1++){
858  for (size_t r2=r1; r2<rr.size(); r2++){
859 
860  glob::FuncGrid interp_r1_eff(new_r, eff[r1], "Spline");
861  glob::FuncGrid interp_r2_eff(new_r, eff[r2], "Spline");
862 
863  auto integrand = [&] ( const double _r) {
864  return interp_r1_eff(_r)*interp_r2_eff(_r)*_r;
865  };
866  II[r1][r2] += fact*wrapper::gsl::GSL_integrate_qag(integrand, min_rr, max_rr, 1.e-3); //Check the integral limits
867  if(r1!=r2)
868  II[r2][r1] += II[r1][r2];
869 
870  }
871  }
872  }
873  }
874  }
875 }
876 
877 
878 // =====================================================================================
879 
880 
881 void cbl::cosmology::Cosmology::k_ell (std::vector<std::vector<double>> &KK, const std::vector<double> rr, const int ll, const std::vector<double> kk, const std::vector<double> Pk)
882 {
883 
884  vector<vector<double>> I1l, I3l, I5l;
885 
886  I_ELL_ell (I1l, rr, ll, 1, kk, Pk);
887  I_ELL_ell (I3l, rr, ll, 3, kk, Pk);
888  I_ELL_ell (I5l, rr, ll, 5, kk, Pk);
889 
890  KK.resize(rr.size(), vector<double>(rr.size(), 0));
891 
892  const double fact = 64./77175;
893  for (size_t r1=0; r1<rr.size(); r1++)
894  for (size_t r2=r1; r2<rr.size(); r2++) {
895  KK[r1][r2] = fact*(9.*I1l[r1][r2]-14.*I3l[r1][r2]+5.*I5l[r1][r2]);
896 
897  if(r1!=r2)
898  KK[r2][r1] = KK[r1][r2];
899  }
900 
901 }
902 
903 
904 // =====================================================================================
905 
906 
907 double cbl::cosmology::Cosmology::zeta_ell_0_factor (const double b1, const double gamma, const double beta)
908 {
909  return pow(b1, 3)*(34./21*(1.+4.*beta/3+1154.*beta*beta/1275+936*pow(beta, 3)/2975+21*pow(beta, 4)/425)+gamma*(1+2.*beta/3+beta*beta/9));
910 }
911 
912 
913 // =====================================================================================
914 
915 
916 double cbl::cosmology::Cosmology::zeta_ell_0_factor_tidal (const double gamma_t, const double beta)
917 {
918  return 16.*beta*beta*gamma_t/675;
919 }
920 
921 
922 // =====================================================================================
923 
924 
925 double cbl::cosmology::Cosmology::zeta_ell_1_factor (const double b1, const double beta)
926 {
927  return -pow(b1, 3)*(1.+4.*beta/3+82*beta*beta/75+12.*pow(beta, 3)/25+3.*pow(beta, 4)/35);
928 }
929 
930 
931 // =====================================================================================
932 
933 
934 double cbl::cosmology::Cosmology::zeta_ell_2_factor (const double b1, const double gamma, const double beta)
935 {
936  return pow(b1, 3)*(8./21*(1.+4.*beta/3+52*beta*beta/21+81.*pow(beta, 3)/49+12.*pow(beta, 4)/35)+32*gamma/945*beta*beta);
937 }
938 
939 
940 // =====================================================================================
941 
942 
943 double cbl::cosmology::Cosmology::zeta_ell_2_factor_tidal (const double gamma_t, const double beta)
944 {
945  return 2.5*(8./15+16*beta/45+344*beta*beta/4725)*gamma_t;
946 }
947 
948 
949 // =====================================================================================
950 
951 
952 double cbl::cosmology::Cosmology::zeta_ell_3_factor (const double b1, const double beta)
953 {
954  return -pow(b1, 3)*(8*beta*beta/75+16.*pow(beta, 3)/175+8.*pow(beta, 4)/315);
955 }
956 
957 
958 // =====================================================================================
959 
960 
961 double cbl::cosmology::Cosmology::zeta_ell_4_factor (const double b1, const double beta)
962 {
963  return pow(b1, 3)*(-32.*beta*beta/3675+32.*pow(beta, 3)/8575+128.*pow(beta, 4)/11025);
964 }
965 
966 
967 // =====================================================================================
968 
969 
970 double cbl::cosmology::Cosmology::zeta_ell_4_factor_tidal (const double gamma_t, const double beta)
971 {
972  return 32.*beta*beta*gamma_t/525;
973 }
974 
975 
976 // =====================================================================================
977 
978 
979 double cbl::cosmology::Cosmology::zeta_ell_k_factor (const double b1, const double beta)
980 {
981  return pow(b1, 3)*(7.*beta*beta+3*pow(beta, 3));
982 }
983 
984 
985 // =====================================================================================
986 
987 
988 double cbl::cosmology::Cosmology::zeta_ell_precyclic (const double r1, const double r2, const int ell, const double b1, const double b2, const double bt, const double beta, std::vector<std::shared_ptr<cbl::glob::FuncGrid>> interp_xi_ell, const bool use_k, std::shared_ptr<cbl::glob::FuncGrid2D> interp_k_ell)
989 {
990  const double gamma = 2.*b2/b1;
991  const double gamma_t = 2.*bt/b1;
992 
993  double fact=0;
994 
995  if(ell==0){
996  fact = (zeta_ell_0_factor( b1, gamma, beta)+zeta_ell_0_factor_tidal(gamma_t, beta))*(interp_xi_ell[0]->operator()(r1)*interp_xi_ell[0]->operator()(r2));
997  }
998  else if (ell==1){
999  fact = zeta_ell_1_factor(b1, beta)*(interp_xi_ell[0]->operator()(r1)*interp_xi_ell[1]->operator()(r2)+interp_xi_ell[0]->operator()(r2)*interp_xi_ell[1]->operator()(r1));
1000  }
1001  else if (ell==2){
1002  fact = (zeta_ell_2_factor( b1, gamma, beta)+zeta_ell_2_factor_tidal(gamma_t, beta))*(interp_xi_ell[0]->operator()(r1)*interp_xi_ell[0]->operator()(r2));
1003  }
1004  else if (ell==3){
1005  fact = zeta_ell_3_factor(b1, beta)*(interp_xi_ell[0]->operator()(r1)*interp_xi_ell[1]->operator()(r2)+interp_xi_ell[0]->operator()(r2)*interp_xi_ell[1]->operator()(r1));
1006  }
1007  else if (ell==4){
1008  fact = (zeta_ell_4_factor(b1, beta) + zeta_ell_4_factor_tidal(gamma_t, beta))*(interp_xi_ell[0]->operator()(r1)*interp_xi_ell[0]->operator()(r2));
1009  }
1010 
1011  return ((use_k) ? fact+zeta_ell_k_factor (b1, beta)*interp_k_ell->operator()(r1, r2) : fact);
1012 }
1013 
1014 
1015 // =====================================================================================
1016 
1017 
1018 std::vector<double> cbl::cosmology::Cosmology::zeta_RSD (const double r1, const double r2, const int ntheta, const double b1, const double b2, const double bt, const double beta, const std::vector<double> rr, const std::vector<double> kk, const std::vector<double> Pk, const bool include_limits, const int max_ll, const bool use_k)
1019 {
1020  (void)max_ll; (void)use_k;
1021  vector<vector<double>> Kl (rr.size(), vector<double>(rr.size(), 0));
1022  auto interp_Kl= make_shared< glob::FuncGrid2D > (glob::FuncGrid2D(rr, rr, Kl, "Linear"));
1023 
1024  const int nbins = (include_limits) ? ntheta+1 : ntheta;
1025  vector<double> zeta(ntheta, 0);
1026  vector<double> r3(ntheta), cos_a12(ntheta), cos_a23(ntheta), cos_a31(ntheta);
1027 
1028  double theta_binSize = (include_limits) ? 0. : par::pi/ntheta*0.5;
1029  for(int i=0; i<nbins; i++) {
1030  double a12 = double(i)*par::pi/ntheta+theta_binSize;
1031  cos_a12[i] = cos(a12);
1032 
1033  r3[i] = sqrt(r1*r1+r2*r2-2.*r1*r2*cos_a12[i]);
1034 
1035  double a23, a31;
1036  if (r1<=r2) {
1037  a23 = asin(sin(a12)*r1/r3[i]);
1038  a31 = par::pi-a12-a23;
1039  }
1040  else {
1041  a31 = asin(sin(a12)*r2/r3[i]);
1042  a23 = par::pi-a12-a31;
1043  }
1044  cos_a23[i] = cos(a23);
1045  cos_a31[i] = cos(a31);
1046  }
1047 
1048  //Even multipoles (0, 2, 4)
1049 
1050  for (int i=0; i<3; i++) {
1051  double ll = 2*i;
1052  vector<double> xil; xi_r_n(xil, rr, ll, kk, Pk);
1053  auto interp_xil = make_shared<glob::FuncGrid> (glob::FuncGrid (rr, xil, "Spline") );
1054 
1055  for(int t=0; t<nbins; t++){
1056  zeta[t] += zeta_ell_precyclic (r1, r2, ll, b1, b2, bt, beta, {interp_xil}, use_k, interp_Kl)*legendre_polynomial(cos_a12[t], ll);
1057  zeta[t] += zeta_ell_precyclic (r2, r3[t], ll, b1, b2, bt, beta, {interp_xil}, use_k, interp_Kl)*legendre_polynomial(cos_a23[t], ll);
1058  zeta[t] += zeta_ell_precyclic (r3[t], r1, ll, b1, b2, bt, beta, {interp_xil}, use_k, interp_Kl)*legendre_polynomial(cos_a31[t], ll);
1059  }
1060  }
1061 
1062  // Odd Multipoles (1, 3)
1063  for (int i=0; i<2; i++) {
1064  double ll = 2*i+1;
1065  vector<double> xil_p, xil_m; xi_r_n_pm (xil_p, xil_m, rr, ll, kk, Pk);
1066  auto interp_xil_p = make_shared<glob::FuncGrid> (glob::FuncGrid (rr, xil_p, "Spline") );
1067  auto interp_xil_m = make_shared<glob::FuncGrid> (glob::FuncGrid (rr, xil_m, "Spline") );
1068 
1069  for(int t=0; t<nbins; t++){
1070  zeta[t] += zeta_ell_precyclic (r1, r2, ll, b1, b2, bt, beta, {interp_xil_p, interp_xil_m}, use_k, interp_Kl)*legendre_polynomial(cos_a12[t], ll);
1071  zeta[t] += zeta_ell_precyclic (r2, r3[t], ll, b1, b2, bt, beta, {interp_xil_p, interp_xil_m}, use_k, interp_Kl)*legendre_polynomial(cos_a23[t], ll);
1072  zeta[t] += zeta_ell_precyclic (r3[t], r1, ll, b1, b2, bt, beta, {interp_xil_p, interp_xil_m}, use_k, interp_Kl)*legendre_polynomial(cos_a31[t], ll);
1073  }
1074  }
1075 
1076  return zeta;
1077 }
1078 
1079 
1080 // =====================================================================================
1081 
1082 
1083 std::vector<double> cbl::cosmology::Cosmology::zeta_RSD (const double r1, const double r2, const int ntheta, const double b1, const double b2, const double bt, const double redshift, const std::string method_Pk, const int step_r, const int step_k, const bool store_output, const std::string output_root, const bool force_RealSpace, const bool include_limits, const int max_ll, const bool use_k)
1084 {
1085  double rmax = r1+r2;
1086 
1087  double beta = (force_RealSpace) ? 0 : linear_growth_rate(redshift)/b1;
1088  vector<double> rr = linear_bin_vector(step_r, 1., rmax);
1089  vector<double> kk = logarithmic_bin_vector(step_k, 1.e-4, 10.);
1090  vector<double> _Pk = Pk_matter(kk, method_Pk, false, redshift, store_output, output_root);
1091 
1092  return zeta_RSD (r1, r2, ntheta, b1, b2, bt, beta, rr, kk, _Pk, include_limits, max_ll, use_k);
1093 }
1094 
Cosmology.h
The class Cosmology.
cbl::cosmology::Cosmology::integrals_zeta_BarrigaGatzanaga
void integrals_zeta_BarrigaGatzanaga(std::vector< double > &xi_matter, std::vector< double > &Phi, const std::vector< double > rr, const std::vector< double > kk, const std::vector< double > Pk_matter) const
integrals used to compute the Barriga & Gatzanaga al. 2002 three-point correlation function model
Definition: 3PCF.cpp:346
cbl::cosmology::Cosmology::I_ELL_ell
void I_ELL_ell(std::vector< std::vector< double >> &II, const std::vector< double > rr, const int ll, const int LL, const std::vector< double > kk, const std::vector< double > Pk)
compute the quantity
Definition: 3PCF.cpp:842
cbl::cosmology::Cosmology::zeta_ell_1_factor
double zeta_ell_1_factor(const double b1, const double beta)
the multiplicative factor for , with local bias
Definition: 3PCF.cpp:925
cbl::cosmology::Cosmology::integrals_zeta_Slepian
void integrals_zeta_Slepian(std::vector< double > &xi_matter, std::vector< double > &xi_matter_m1, std::vector< double > &xi_matter_p1, std::vector< double > &xi_matter_2, const std::vector< double > rr, const std::vector< double > kk, const std::vector< double > Pk_matter) const
integrals used to compute the Slepian et al. 2015 three-point correlation function model
Definition: 3PCF.cpp:161
cbl::cosmology::Cosmology::zeta_RSD
std::vector< double > zeta_RSD(const double r1, const double r2, const int ntheta, const double b1, const double b2, const double bt, const double beta, const std::vector< double > rr, const std::vector< double > kk, const std::vector< double > Pk, const bool include_limits=false, const int max_ll=4, const bool use_k=false)
the
Definition: 3PCF.cpp:1018
cbl::cosmology::Cosmology::zeta_multipoles_covariance
double zeta_multipoles_covariance(const double Volume, const double nObjects, const int l, const int l_prime, const double r1, const double r2, const double r1_prime, const double r2_prime, const double deltaR, const std::vector< double > kk, const std::vector< double > Pk, const std::vector< double > rr, const std::vector< double > Xi, const double prec=1.e-3)
the dark matter three-point correlation function multipoles covariance model, by Slepian et al....
Definition: 3PCF.cpp:614
cbl::cosmology::Cosmology::zeta_ell_k_factor
double zeta_ell_k_factor(const double b1, const double beta)
the multiplicative factor for , with local bias
Definition: 3PCF.cpp:979
cbl::cosmology::Cosmology::zeta_ell_precyclic
double zeta_ell_precyclic(const double r1, const double r2, const int ell, const double b1, const double b2, const double bt, const double beta, std::vector< std::shared_ptr< glob::FuncGrid >> interp_xi_ell, const bool use_k, std::shared_ptr< glob::FuncGrid2D > interp_k_ell)
the pre-cyclic
Definition: 3PCF.cpp:988
cbl::cosmology::Cosmology::Q_DM_eq
std::vector< double > Q_DM_eq(const std::vector< double > rr, const std::string model, const std::vector< double > kk, const std::vector< double > Pk_matter) const
the dark matter equilateral reduced three-point correlation function
Definition: 3PCF.cpp:586
cbl::cosmology::Cosmology::k_ell
void k_ell(std::vector< std::vector< double >> &KK, const std::vector< double > rr, const int ll, const std::vector< double > kk, const std::vector< double > Pk)
compute the quantity
Definition: 3PCF.cpp:881
cbl::cosmology::Cosmology::zeta_halo
std::vector< double > zeta_halo(const double r1, const double r2, const std::vector< double > theta, const double b1, const double b2, const std::string model, const std::vector< double > kk, const std::vector< double > Pk_matter) const
the local-bias model of the three-point correlation function of dark matter haloes
Definition: 3PCF.cpp:499
cbl::cosmology::Cosmology::zeta_DM_BarrigaGatzanaga
double zeta_DM_BarrigaGatzanaga(const double r1, const double r2, const double theta, std::vector< double > &rr, std::vector< double > &xi_matter, std::vector< double > &Phi, const std::vector< double > kk, const std::vector< double > Pk_matter) const
the dark matter three-point correlation function model by Barriga & Gatzanaga et al....
Definition: 3PCF.cpp:387
cbl::cosmology::Cosmology::Q_halo
std::vector< double > Q_halo(const double r1, const double r2, const std::vector< double > theta, const double b1, const double b2, const std::string model, const std::vector< double > kk, const std::vector< double > Pk_matter) const
the local-bias model of the reduced three-point correlation function of dark matter haloes
Definition: 3PCF.cpp:526
cbl::cosmology::Cosmology::zeta_ell_3_factor
double zeta_ell_3_factor(const double b1, const double beta)
the multiplicative factor for , with local bias
Definition: 3PCF.cpp:952
cbl::cosmology::Cosmology::zeta_DM_Slepian
double zeta_DM_Slepian(const double r1, const double r2, const double theta, std::vector< double > &rr, std::vector< double > &xi_matter, std::vector< double > &xi_matter_m1, std::vector< double > &xi_matter_p1, std::vector< double > &xi_matter_2, const std::vector< double > kk, const std::vector< double > Pk_matter, const int norders=9, const double prec=1.e-3) const
the dark matter three-point correlation function model by Slepian et al. 2015
Definition: 3PCF.cpp:314
cbl::cosmology::Cosmology::zeta_expansion_Slepian
std::vector< double > zeta_expansion_Slepian(const double r1, const double r2, const double b1, const double b2, std::vector< double > &rr, std::vector< double > &xi_matter, std::vector< double > &xi_matter_m1, std::vector< double > &xi_matter_p1, std::vector< double > &xi_matter_2, const int norders=9, const double prec=1.e-3) const
the terms of the expansion
Definition: 3PCF.cpp:291
cbl::cosmology::Cosmology::Q_nonLocal
double Q_nonLocal(const double r1, const double r2, const double theta, std::vector< double > &rr, std::vector< double > &xi_matter, std::vector< double > &Phi, const std::vector< double > kk, const std::vector< double > Pk_matter) const
the non-local contribution to the reduced dark matter three-point correlation function
Definition: 3PCF.cpp:99
cbl::cosmology::Cosmology::zeta_ell_0_factor_tidal
double zeta_ell_0_factor_tidal(const double gamma_t, const double beta)
the multiplicative factor for , with non-local bias
Definition: 3PCF.cpp:916
cbl::cosmology::Cosmology::zeta_ell_2_factor
double zeta_ell_2_factor(const double b1, const double gamma, const double beta)
the multiplicative factor for , with local bias
Definition: 3PCF.cpp:934
cbl::cosmology::Cosmology::zeta_precyclic_Slepian
double zeta_precyclic_Slepian(const double r1, const double r2, const double mu, const double b1, const double b2, const glob::FuncGrid interp_xi_matter, const glob::FuncGrid interp_xi_matter_m1, const glob::FuncGrid interp_xi_matter_p1, const glob::FuncGrid interp_xi_matter_2) const
the pre-cyclic three-point correlation function as described in Slepian et al. 2015
Definition: 3PCF.cpp:181
cbl::cosmology::Cosmology::xi_r_n_pm
void xi_r_n_pm(std::vector< double > &xi_n_p, std::vector< double > &xi_n_m, const std::vector< double > rr, const int nn, const std::vector< double > kk, const std::vector< double > Pk)
compute the power spectrum integral transform
Definition: 3PCF.cpp:802
cbl::cosmology::Cosmology::zeta_ell_4_factor_tidal
double zeta_ell_4_factor_tidal(const double gamma_t, const double beta)
the multiplicative factor for , with non-local bias
Definition: 3PCF.cpp:970
cbl::cosmology::Cosmology::zeta_ell_2_factor_tidal
double zeta_ell_2_factor_tidal(const double gamma_t, const double beta)
the multiplicative factor for , with non-local bias
Definition: 3PCF.cpp:943
cbl::cosmology::Cosmology::zeta_single_BarrigaGatzanaga
double zeta_single_BarrigaGatzanaga(const double r1, const double r2, const double theta, const std::vector< double > xi, const std::vector< double > dxi, const std::vector< double > dPhi) const
the single term of the dark matter three-point correlation function model by Barriga & Gatzanaga et a...
Definition: 3PCF.cpp:362
cbl::cosmology::Cosmology::Gamma_3PCF
double Gamma_3PCF(const double r1, const double r2, const double theta, const std::vector< double > xi, const std::vector< double > dPhi) const
function to compute non-local contribution to three-point correlation function; specifically,...
Definition: 3PCF.cpp:90
cbl::cosmology::Cosmology::integrals_Q_nonLocal
void integrals_Q_nonLocal(std::vector< double > &xi_matter, std::vector< double > &Phi, const std::vector< double > rr, const std::vector< double > kk, const std::vector< double > Pk_matter, const double prec) const
integral functions for the three-point correlation model
Definition: 3PCF.cpp:63
cbl::cosmology::Cosmology::Q_DM_Slepian
double Q_DM_Slepian(const double r1, const double r2, const double theta, std::vector< double > &rr, std::vector< double > &xi_matter, std::vector< double > &xi_matter_m1, std::vector< double > &xi_matter_p1, std::vector< double > &xi_matter_2, const std::vector< double > kk, const std::vector< double > Pk_matter, const int norders=9, const double prec=1.e-3) const
the dark matter reduced three-point correlation function model by Slepian et al. 2015
Definition: 3PCF.cpp:335
cbl::cosmology::Cosmology::eff_l_l1
void eff_l_l1(std::vector< std::vector< double >> &eff, const std::vector< double > rr, const int l, const int l1, const std::vector< double > kk, const std::vector< double > Pk)
compute the power spectrum integral transform
Definition: 3PCF.cpp:819
cbl::cosmology::Cosmology::denominator_Q
double denominator_Q(const double r1, const double r2, const double theta, const std::vector< double > rr, const std::vector< double > xi_matter) const
the normalization factor for reduced three-point correlation function
Definition: 3PCF.cpp:47
cbl::cosmology::Cosmology::zeta_covariance
std::vector< std::vector< double > > zeta_covariance(const double Volume, const double nObjects, const std::vector< double > theta, const double r1, const double r2, const double deltaR, const std::vector< double > kk, const std::vector< double > Pk, const int norders=10, const double prec=1.e-3, const bool method=false, const int nExtractions=10000, const std::vector< double > mean={}, const int seed=543)
the dark matter three-point correlation function covariance model
Definition: 3PCF.cpp:753
cbl::cosmology::Cosmology::zeta_ell_0_factor
double zeta_ell_0_factor(const double b1, const double gamma, const double beta)
the multiplicative factor for , with local bias
Definition: 3PCF.cpp:907
cbl::cosmology::Cosmology::Q_DM_BarrigaGatzanaga
double Q_DM_BarrigaGatzanaga(const double r1, const double r2, const double theta, std::vector< double > &rr, std::vector< double > &xi_matter, std::vector< double > &Phi, const std::vector< double > kk, const std::vector< double > Pk_matter) const
the dark matter reduced three-point correlation function model by Barriga & Gatzanaga et al....
Definition: 3PCF.cpp:437
cbl::cosmology::Cosmology::Q_DM
std::vector< double > Q_DM(const double r1, const double r2, const std::vector< double > theta, const std::string model, const std::vector< double > kk, const std::vector< double > Pk_matter) const
the dark matter reduced three-point correlation function
Definition: 3PCF.cpp:471
cbl::cosmology::Cosmology::zeta_DM_eq
std::vector< double > zeta_DM_eq(const std::vector< double > rr, const std::string model, const std::vector< double > kk, const std::vector< double > Pk_matter) const
the dark matter equilateral three-point correlation function
Definition: 3PCF.cpp:558
cbl::cosmology::Cosmology::zeta_ell_4_factor
double zeta_ell_4_factor(const double b1, const double beta)
the multiplicative factor for , with local bias
Definition: 3PCF.cpp:961
cbl::cosmology::Cosmology::zeta_DM
std::vector< double > zeta_DM(const double r1, const double r2, const std::vector< double > theta, const std::string model, const std::vector< double > kk, const std::vector< double > Pk_matter) const
the dark matter three-point correlation function
Definition: 3PCF.cpp:446
cbl::cosmology::Cosmology::xi_r_n
void xi_r_n(std::vector< double > &xi_n, const std::vector< double > rr, const int nn, const std::vector< double > kk, const std::vector< double > Pk)
compute the power spectrum integral transform
Definition: 3PCF.cpp:793
cbl::glob::FuncGrid2D
The class FuncGrid2D.
Definition: FuncGrid.h:302
cbl::glob::FuncGrid
The class FuncGrid.
Definition: FuncGrid.h:55
cbl::glob::FuncGrid::D1v
double D1v(const double xx) const
compute the first derivative at xx
Definition: FuncGrid.cpp:167
cbl::par::pi
static const double pi
: the ratio of a circle's circumference to its diameter
Definition: Constants.h:199
cbl::modelling::threept::zeta_RSD
std::vector< double > zeta_RSD(const std::vector< double > theta, const std::shared_ptr< void > inputs, std::vector< double > &parameter)
model for the connected three-point correlation function
Definition: ModelFunction_ThreePointCorrelation_comoving_connected.cpp:49
cbl::wrapper::fftlog::transform_FFTlog
std::vector< double > transform_FFTlog(const std::vector< double > yy, const int dir, const std::vector< double > xx, const std::vector< double > fx, const double mu=0, const double q=0, const double kr=1, const int kropt=0)
wrapper of the FFTlog to compute the discrete Fourier sine or cosine transform of a logarithmically...
Definition: FFTlog.cpp:43
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::Min
T Min(const std::vector< T > vect)
minimum element of a std::vector
Definition: Kernel.h:1324
cbl::logarithmic_bin_vector
std::vector< T > logarithmic_bin_vector(const size_t nn, const T min, const T max)
fill a std::vector with logarithmically spaced values
Definition: Kernel.h:1621
cbl::linear_bin_vector
std::vector< T > linear_bin_vector(const size_t nn, const T min, const T max)
fill a std::vector with linearly spaced values
Definition: Kernel.h:1604
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::Max
T Max(const std::vector< T > vect)
maximum element of a std::vector
Definition: Kernel.h:1336
cbl::jl_distance_average
double jl_distance_average(const double kk, const int order, const double r_down, const double r_up)
the distance average for the order l-th spherical Bessel function
Definition: Func.cpp:2109
cbl::legendre_polynomial
double legendre_polynomial(const double mu, const int l)
the order l Legendre polynomial
Definition: Func.cpp:1786
cbl::jl
double jl(const double xx, const int order)
the order l spherical Bessel function
Definition: Func.cpp:2066