CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
SizeFunction.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2010 by Federico Marulli and Tommaso Ronconi *
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 
35 #include "Cosmology.h"
36 
37 using namespace std;
38 using namespace cbl;
39 
40 
41 // =====================================================================================
42 
43 
44 double cbl::cosmology::Cosmology::deltav_L (const double deltav_NL, const double b_eff, double slope, double offset) const
45 {
46  if (b_eff==1.)
47  {slope = 1.; offset = 0.;}
48 
49  return 1.594*(1.-pow(1+deltav_NL/(slope*b_eff+offset), (-1./1.594)));
50 }
51 
52 
53 // =====================================================================================
54 
55 
56 double cbl::cosmology::Cosmology::deltav_NL (const double deltav_L) const
57 {
58  return pow((1.-deltav_L/1.594), -1.594) - 1.;
59 }
60 
61 
62 // =====================================================================================
63 
64 
65 double cbl::cosmology::Cosmology::r_rL (const double deltav) const
66 {
67  return pow(pow((1.-deltav/1.594), -1.594), -1./3.);
68 }
69 
70 
71 // =====================================================================================
72 
73 
74 double cbl::cosmology::Cosmology::f_nu (const double SS, const double del_v, const double del_c) const
75 {
76  double radnu = fabs(del_v)/SS;
77  double nu = pow(radnu, 2.);
78  double DDD = fabs(del_v)/(del_c+fabs(del_v));
79  double xx = DDD/radnu;
80 
81  if (xx<=0.276)
82  return sqrt(2./par::pi)*radnu*exp(-0.5*nu);
83 
84  else {
85  double ff = 0.;
86  for (int j = 1; j <= 4; j++)
87  ff += 2.*exp(-pow((j*par::pi*xx), 2.)/2.)*j*par::pi*pow(xx, 2.)*sin(j*par::pi*DDD);
88  return ff;
89  }
90 }
91 
92 
93 // =====================================================================================
94 
95 
96 vector<double> cbl::cosmology::Cosmology::AP_corr(const cbl::cosmology::Cosmology cosm_true, const std::vector<double> redshift)
97 {
98  std::vector<double> APcorr_vect(redshift.size(),0.);
99 
100  for (size_t ii=0; ii<redshift.size(); ii++) {
101  double q_par = HH(redshift[ii])/cosm_true.HH(redshift[ii]);
102  double q_perp = cosm_true.D_M(redshift[ii])/D_M(redshift[ii]);
103  APcorr_vect[ii] = pow(q_par*q_perp*q_perp,(-1./3.));
104  }
105  return APcorr_vect;
106 }
107 
108 
109 // =====================================================================================
110 
111 
112 double cbl::cosmology::Cosmology::size_function (const double RV, const double redshift, const std::string model, const double b_eff, double slope, double offset, const double deltav_NL, const double del_c, const std::string method_Pk, const bool store_output, const std::string output_root, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file) const
113 {
114 
115  double del_v = deltav_L(deltav_NL, b_eff, slope, offset);
116  double RL;
117 
118  if ((model == "Vdn") || (model == "SvdW"))
119  RL = RV/r_rL(del_v);
120 
121  else if (model == "linear")
122  RL = RV;
123 
124  else
125  { ErrorCBL("the model name is not allowed; the allowed names are: SvdW (Sheth and van de Weygaert, 2004), linear/Vdn (Jennings, Li and Hu, 2013)", "size_function", "SizeFunction.cpp"); return 0; }
126 
127  double fact = DN(redshift, 0.);
128 
129  double sigmaR = sqrt(sigma2R(RL, method_Pk, 0., true, output_root, interpType, k_max, input_file, is_parameter_file));
130  double sigmaRz = sigmaR*fact;
131  double SSSR = sigmaRz*sigmaRz;
132 
133  double Dln_SigmaR = dnsigma2R(1, RL, method_Pk, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file)*(RL/(2.*SSSR))*pow(fact, 2.);
134 
135  if (model == "Vdn") return f_nu(sigmaRz, del_v, del_c)/volume_sphere(RV)*fabs(Dln_SigmaR);
136  else return f_nu(sigmaRz, del_v, del_c)/volume_sphere(RL)*fabs(Dln_SigmaR);
137 
138 }
139 
140 
141 // =====================================================================================
142 
143 
144 std::vector<double> cbl::cosmology::Cosmology::size_function (const std::vector<double> RV, const double redshift, const std::string model, const double b_eff, double slope, double offset, const double deltav_NL, const double del_c, const std::string method_Pk, const bool store_output, const std::string output_root, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file) const
145 {
146 
147  double del_v = deltav_L(deltav_NL, b_eff, slope, offset);
148 
149  std::vector<double> RL (RV.size());
150 
151  if ((model == "Vdn") || (model == "SvdW"))
152  for (size_t i=0; i<RV.size(); i++) RL[i] = RV[i]/r_rL(del_v);
153 
154  else if (model == "linear")
155  for (size_t i=0; i<RV.size(); i++) RL[i] = RV[i];
156 
157  else
158  ErrorCBL("the model name is not allowed; the allowed names are: SvdW (Sheth and van de Weygaert, 2004), linear/Vdn (Jennings, Li and Hu, 2013)", "size_function", "SizeFunction.cpp");
159 
160  double fact = DN(redshift);
161 
162  std::vector<double> sigmaR(RV.size());
163  std::vector<double> sigmaRz(RV.size());
164  std::vector<double> SSSR(RV.size());
165  std::vector<double> Dln_SigmaR(RV.size());
166 
167  for (size_t i=0; i<RV.size(); i++) {
168 
169  sigmaR[i] = sqrt(sigma2R(RL[i], method_Pk, 0., true, output_root, interpType, k_max, input_file, is_parameter_file));
170  sigmaRz[i] = sigmaR[i]*fact;
171  SSSR[i] = sigmaRz[i]*sigmaRz[i];
172  Dln_SigmaR[i] = dnsigma2R(1, RL[i], method_Pk, 0., true, output_root, interpType, k_max, input_file, is_parameter_file)*(RL[i]/(2.*SSSR[i]))*pow(fact, 2.);
173 
174  }
175 
176  if (!store_output) {
177  remove_output_Pk_tables(method_Pk, false, 0., output_root);
178  remove_output_Pk_tables(method_Pk, false, redshift, output_root);
179  }
180 
181  std::vector<double> result(RV.size());
182 
183  if (model == "Vdn") for (size_t i=0; i<RV.size(); i++) result[i] = f_nu(sigmaRz[i], del_v, del_c)/volume_sphere(RV[i])*fabs(Dln_SigmaR[i]);
184  else for (size_t i=0; i<RV.size(); i++) result[i] = f_nu(sigmaRz[i], del_v, del_c)/volume_sphere(RL[i])*fabs(Dln_SigmaR[i]);
185 
186  return result;
187 
188 }
189 
190 
191 // =====================================================================================
192 
193 
194 std::vector<std::vector<double>> cbl::cosmology::Cosmology::Nvoids (const double min_r, const double max_r, const int num_bins, const double mean_z, const double Volume, const std::string model, const double b_eff, double slope, double offset, const double deltav_NL, const double del_c, const std::string method_Pk, const bool store_output, const std::string output_root, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file) const
195 {
196  double del_v = deltav_L(deltav_NL, b_eff, slope, offset);
197  double NLcorr;
198 
199  double fact = DN(mean_z);
200 
201  std::vector<double> sigmaR(num_bins);
202  std::vector<double> sigmaRz(num_bins);
203  std::vector<double> SSSR(num_bins);
204  std::vector<double> Dln_SigmaR(num_bins);
205 
206  std::vector<double> r_bins = cbl::logarithmic_bin_vector(num_bins+1, min_r, max_r);
207 
208  std::vector<double> RV(num_bins);
209  std::vector<double> RL (num_bins);
210 
211  if ((model == "Vdn") || (model == "SvdW")) NLcorr = r_rL(del_v);
212  else if (model == "linear") NLcorr = 1.;
213  else ErrorCBL("the model name is not allowed; the allowed names are: SvdW (Sheth and van de Weygaert, 2004), linear/Vdn (Jennings, Li and Hu, 2013)", "size_function", "SizeFunction.cpp");
214 
215  for (int i=0; i<num_bins; i++) {
216 
217  RV[i] = pow(10., log10(r_bins[i])+0.5*(log10(r_bins[i+1]/r_bins[i])));
218  RL[i] = RV[i]/NLcorr;
219  sigmaR[i] = sqrt(sigma2R(RL[i], method_Pk, 0., true, output_root, interpType, k_max, input_file, is_parameter_file));
220  sigmaRz[i] = sigmaR[i]*fact;
221  SSSR[i] = sigmaRz[i]*sigmaRz[i];
222  Dln_SigmaR[i] = dnsigma2R(1, RL[i], method_Pk, 0., true, output_root, interpType, k_max, input_file, is_parameter_file)*(RL[i]/(2.*SSSR[i]))*pow(fact, 2.);
223 
224  }
225 
226  if (!store_output) {
227  remove_output_Pk_tables(method_Pk, false, 0., output_root);
228  remove_output_Pk_tables(method_Pk, false, mean_z, output_root);
229  }
230 
231  std::vector<std::vector<double>> result(2,std::vector<double>(num_bins));
232 
233  if (model == "Vdn") for (int i=0; i<num_bins; i++) {
234  result[0][i] = RV[i];
235  result[1][i] = f_nu(sigmaRz[i], del_v, del_c)/volume_sphere(RV[i])*fabs(Dln_SigmaR[i])*(r_bins[i+1]-r_bins[i])/RV[i]*Volume;
236  }
237 
238  else for (int i=0; i<num_bins; i++) {
239  result[0][i] = RV[i];
240  result[1][i] = f_nu(sigmaRz[i], del_v, del_c)/volume_sphere(RL[i])*fabs(Dln_SigmaR[i])*(r_bins[i+1]-r_bins[i])/RV[i]*Volume;
241  }
242  return result;
243 
244 }
245 
246 // =====================================================================================
247 
248 
249 std::vector<std::vector<double>> cbl::cosmology::Cosmology::Nvoids (const double min_r, const double max_r, const int num_bins, const double min_z, const double max_z, const double mean_z, const double Area, const std::string model, const double b_eff, double slope, double offset, const double deltav_NL, const double del_c, const std::string method_Pk, const bool store_output, const std::string output_root, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file) const
250 {
251  double del_v = deltav_L(deltav_NL, b_eff, slope, offset);
252  double NLcorr;
253 
254  double fact = DN(mean_z);
255 
256  std::vector<double> sigmaR(num_bins);
257  std::vector<double> sigmaRz(num_bins);
258  std::vector<double> SSSR(num_bins);
259  std::vector<double> Dln_SigmaR(num_bins);
260 
261  std::vector<double> r_bins = cbl::logarithmic_bin_vector(num_bins+1, min_r, max_r);
262 
263  std::vector<double> RV(num_bins);
264  std::vector<double> RL (num_bins);
265 
266  if ((model == "Vdn") || (model == "SvdW")) NLcorr = r_rL(del_v);
267  else if (model == "linear") NLcorr = 1.;
268  else ErrorCBL("the model name is not allowed; the allowed names are: SvdW (Sheth and van de Weygaert, 2004), linear/Vdn (Jennings, Li and Hu, 2013)", "size_function", "SizeFunction.cpp");
269 
270  for (int i=0; i<num_bins; i++) {
271 
272  RV[i] = pow(10., log10(r_bins[i])+0.5*(log10(r_bins[i+1]/r_bins[i])));
273  RL[i] = RV[i]/NLcorr;
274  sigmaR[i] = sqrt(sigma2R(RL[i], method_Pk, 0., true, output_root, interpType, k_max, input_file, is_parameter_file));
275  sigmaRz[i] = sigmaR[i]*fact;
276  SSSR[i] = sigmaRz[i]*sigmaRz[i];
277  Dln_SigmaR[i] = dnsigma2R(1, RL[i], method_Pk, 0., true, output_root, interpType, k_max, input_file, is_parameter_file)*(RL[i]/(2.*SSSR[i]))*pow(fact, 2.);
278 
279  }
280 
281  if (!store_output) {
282  remove_output_Pk_tables(method_Pk, false, 0., output_root);
283  remove_output_Pk_tables(method_Pk, false, mean_z, output_root);
284  }
285 
286  std::vector<std::vector<double>> result(2,std::vector<double>(num_bins));
287 
288  double volume = Volume(min_z,max_z,Area);
289  if (model == "Vdn") for (int i=0; i<num_bins; i++) {
290  result[0][i] = RV[i];
291  result[1][i] = f_nu(sigmaRz[i], del_v, del_c)/volume_sphere(RV[i])*fabs(Dln_SigmaR[i])*(r_bins[i+1]-r_bins[i])/RV[i]*volume;
292  }
293 
294  else for (int i=0; i<num_bins; i++) {
295  result[0][i] = RV[i];
296  result[1][i] = f_nu(sigmaRz[i], del_v, del_c)/volume_sphere(RL[i])*fabs(Dln_SigmaR[i])*(r_bins[i+1]-r_bins[i])/RV[i]*volume;
297  }
298  return result;
299 
300 }
301 
302 
303 // =====================================================================================
304 
305 
306 double cbl::cosmology::Cosmology::size_function (const double RV, const double redshift, const std::string model_mf, const double del_v, const std::string model_sf, const std::string method_Pk, const bool store_output, const std::string output_root, const double Delta, const std::string interpType, const int norm, const double k_min, const double k_max, const double prec, const std::string input_file, const bool is_parameter_file)
307 {
308  double RL;
309 
310  if ((model_sf == "Vdn") || (model_sf == "SvdW"))
311  RL = RV/r_rL(del_v);
312 
313  else if (model_sf == "linear")
314  RL = RV;
315 
316  else
317  { ErrorCBL("the model name is not allowed; the allowed names are: SvdW (Sheth and van de Weygaert, 2004), linear/Vdn (Jennings, Li and Hu, 2013)", "size_function", "SizeFunction.cpp"); return 0; }
318 
319  double RHO = rho_m(redshift, true);
320  double MM = Mass(RL, RHO);
321 
322  if (model_sf == "Vdn")
323  return 3.*MM*cosmology::Cosmology::mass_function(MM, redshift, model_mf, method_Pk, store_output, output_root, Delta, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file, false, del_v)*deltav_NL(del_v);
324 
325  else if ((model_sf == "SvdW") || (model_sf == "linear"))
326  return 3.*MM*cosmology::Cosmology::mass_function(MM, redshift, model_mf, method_Pk, store_output, output_root, Delta, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file, false, del_v);
327 
328  else
329  { ErrorCBL("the model name is not allowed; the allowed names are: SvdW (Sheth and van de Weygaert, 2004), linear/Vdn (Jennings, Li and Hu, 2013)", "size_function", "SizeFunction.cpp"); return 0; }
330 
331 }
The class Cosmology.
The class Cosmology.
Definition: Cosmology.h:277
double deltav_NL(const double deltav=-2.71) const
Non-Linear (under)density contrast.
double r_rL(const double deltav=-2.71) const
expansion factor
double size_function(const double RV, const double redshift, const std::string model, const double b_eff, double slope=0.854, double offset=0.420, const double deltav_NL=-0.795, const double del_c=1.69, const std::string method_Pk="EisensteinHu", const bool store_output=true, const std::string output_root="test", const std::string interpType="Linear", const double k_max=100., const std::string input_file=par::defaultString, const bool is_parameter_file=true) const
the void size function
double deltav_L(const double deltav_NL, const double b_eff, double slope=0.854, double offset=0.420) const
Linear (under)density contrast.
double HH(const double redshift=0.) const
the Hubble function
Definition: Cosmology.cpp:570
std::vector< double > AP_corr(const cbl::cosmology::Cosmology cosm_true, const std::vector< double > redshift)
Supplementary function to compute a correction factor to apply to the void size function,...
double f_nu(const double SS, const double del_v, const double del_c) const
(approximation)
double D_M(const double redshift) const
the comoving transverse distance at a given redshift
Definition: Cosmology.cpp:832
std::vector< std::vector< double > > Nvoids(const double min_r, const double max_r, const int num_bins, const double mean_z, const double Volume, const std::string model, const double b_eff, double slope=0.854, double offset=0.420, const double deltav_NL=-0.795, const double del_c=1.69, const std::string method_Pk="EisensteinHu", const bool store_output=true, const std::string output_root="test", const std::string interpType="Linear", const double k_max=100., const std::string input_file=par::defaultString, const bool is_parameter_file=true) const
number of voids computed from the void size function model for bins of radii spaced in log scale and ...
static const double pi
: the ratio of a circle's circumference to its diameter
Definition: Constants.h:199
double mass_function(const double mass, cosmology::Cosmology cosmology, const double redshift, const std::string model_MF, const bool store_output, const double Delta, const bool isDelta_critical, const cbl::glob::FuncGrid interp_Pk, const double kmax)
compute the mass function
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
double sigmaR(const double RR, const int corrType, const std::vector< double > rr, const std::vector< double > corr)
the rms mass fluctuation within a given radius
Definition: FuncXi.cpp:227
T Mass(const T RR, const T Rho)
the mass of a sphere of a given radius and density
Definition: Func.h:1243
T volume_sphere(const T RR)
the volume of a sphere of a given radius
Definition: Func.h:1231
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
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
double volume(const double boxSize, const int frac, const double Bord, const double mean_redshift, cosmology::Cosmology &real_cosm)
get the volume of a simulation box
Definition: Func.cpp:78