LAL  7.5.0.1-8083555
LALCosmologyCalculator.c
Go to the documentation of this file.
1 /* Copyright (C) 2012 Walter Del Pozzo, Tjonnie Li
2  *
3  * This program is free software; you can redistribute it and/or modify
4  * it under the terms of the GNU General Public License as published by
5  * the Free Software Foundation; either version 2 of the License, or
6  * (at your option) any later version.
7  *
8  * This program is distributed in the hope that it will be useful,
9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  * GNU General Public License for more details.
12  *
13  * You should have received a copy of the GNU General Public License
14  * along with with program; see the file COPYING. If not, write to the
15  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
16  * MA 02110-1301 USA
17  */
18 #include <lal/LALMalloc.h>
19 #include <lal/LALCosmologyCalculator.h>
20 
21 /**
22  * The set of functions in this module implements the standard cosmological distance measures
23  * defined a Friedmann-Robertson-Walker-Lemaitre metric.
24  * For a detailed reference to the various formulae, see Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 ).
25  */
26 
27 /**
28  * Computes the luminosity distance as the angular distance divided by (1+z)^2
29  * Eq. 21 in Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 )
30  */
31 
34  double z)
35 {
36  double x=1.0+z;
37  return XLALComovingTransverseDistance(omega, z)*x;
38 }
39 
40 /**
41  * Computes the angular size distance by dividing the comoving transverse distance by 1+z
42  * Eq. 18 in Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 )
43  */
44 
47  double z)
48 {
49  double x=1.0+z;
50  return XLALComovingTransverseDistance(omega, z)/x;
51 }
52 
53 /**
54  * Computes the comoving line-of-sight distance (usually called the proper distance) by integrating the Hubble parameter.
55  * Eq. 15 in Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 ).
56  */
57 
60  double z)
61 {
62  double dH = XLALHubbleDistance(omega);
63  return dH*XLALIntegrateHubbleParameter(omega,z);
64 }
65 /**
66  * Computes the comoving transverse distance from the comoving line-of-sight distance (proper distance) and
67  * Eq. 16 in Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 )
68  */
71  double z)
72 {
73  double ok = omega->ok;
74  double dH = XLALHubbleDistance(omega);
75  double dC = XLALComovingLOSDistance(omega,z);
76 
77  if (fabs(ok)<1e-7)
78  {
79 
80  return dC;
81  }
82  else if (ok>1e-7)
83  {
84  return dH*sinh(sqrt(ok)*dC/dH)/sqrt(ok);
85  }
86  else if (ok<1e-7)
87  {
88  return dH*sin(sqrt(fabs(ok))*dC/dH)/sqrt(fabs(ok));
89  }
90  else
91  {
92  fprintf(stderr,"Something funny happened. Aborting.\n");
93  exit(-1);
94  }
95 
96 }
97 /**
98  * Computes the Hubble distance, independent of the redshift. This is the distance with sets the units for all others.
99  * It returns the distance in Mpc.
100  * Eq. 4 in Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 )
101  */
104  )
105 {
106  double x = GSL_CONST_MKSA_SPEED_OF_LIGHT/(100.0*omega->h);
107  return 1e-3*x;
108 }
109 
110 /**
111  * Computes the inverse of the Hubble parameter at redshift z
112  * Eq. 14 in Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 )
113  */
114 
115 double XLALHubbleParameter(double z,
116  void *omega
117  )
118 {
120 
121  double E=0.0;
122  double x = 1.0+z;
123  double om=p->om;
124  double ol=p->ol;
125  double ok=p->ok;
126  double w0=p->w0;
127  double w1=p->w1;
128  double w2=p->w2;
129  E = sqrt(om*x*x*x+ok*x*x+ol*pow(x,3.*(1.0+w0+w1+w2))*exp(-3.0*((w1+w2)*z/x + w2*z*z/(2.0*x*x))));
130  return 1.0/E;
131 }
132 /**
133  * Computes the integral of inverse of the Hubble parameter at redshift z.
134  * The integral is computed using the built-in function gsl_integration_qag of the gsl library.
135  * The integral is performed using a Gauss-Kronrod adaptive method (see http://www.gnu.org/software/gsl/manual/html_node/QAG-adaptive-integration.html )
136  */
138 {
139  double result = 0.0;
140  double error;
141  double epsabs = 5e-5;
142  double epsrel = 1e-5;
143  size_t limit = 1024;
144  int key = 1;
145 
146  gsl_function F;
147  F.function = &XLALHubbleParameter;
148  F.params = omega;
149 
150  gsl_integration_workspace * w
151  = gsl_integration_workspace_alloc (1024);
152 
153  gsl_integration_qag (&F, 0.0, z, epsabs, epsrel,
154  limit, key, w, &result, &error);
155 
156  gsl_integration_workspace_free (w);
157 
158  return result;
159 }
160 /**
161  * This function computes the value of a uniform probability distribution over the comoving volume.
162  * Details of the derivation of these formulae can be found in Coward, Burman 2005 ( http://arxiv.org/abs/astro-ph/0505181 )
163  */
166  double z,
167  double zmax)
168 {
169  double norm = XLALIntegrateComovingVolumeDensity(omega, zmax);
170  double unnorm_density = XLALUniformComovingVolumeDensity(z,(void *)omega);
171  return unnorm_density/norm;
172 }
173 
174 /**
175  * This function computes the value of a uniform probability density distribution over the comoving volume at redshift z.
176  */
177 
179  double z,
180  void *omega)
181 {
182 
184 
185  double x = 1.0+z;
186  double unnorm_density = XLALComovingVolumeElement(z,p)/x;
187  return unnorm_density;
188 }
189 /**
190  * This function computes the comoving volume between 0 and z
191  */
194  double z)
195 {
196  double V = XLALIntegrateComovingVolume(omega, z);
197  return V;
198 }
199 
200 /**
201  * This function computes the comoving volume element (a 4&pi shell) between redshift z and z+dz.
202  */
203 
205  double z,
206  void *omega)
207 {
208 
210 
211  double dm = XLALComovingTransverseDistance(omega,z);
212  double E = XLALHubbleParameter(z,omega);
213  double dVdz = 4.0*M_PI*dm*dm*E*XLALHubbleDistance(p);
214  return dVdz;
215 }
216 
217 /**
218  * Function that integrates the comoving volume element.
219  * The integration is performed using gsl_integration_qagiu from the gsl library (see http://www.gnu.org/software/gsl/manual/html_node/QAGI-adaptive-integration-on-infinite-intervals.html )
220  */
222 {
223  double result = 0.0;
224  double error;
225  double epsabs = 5e-5;
226  double epsrel = 1e-5;
227  size_t limit = 1024;
228  int key = 1;
229 
230  gsl_function F;
231  F.function = &XLALComovingVolumeElement;
232  F.params = omega;
233 
234  gsl_integration_workspace * w
235  = gsl_integration_workspace_alloc (1024);
236 
237  if (z<0.0) gsl_integration_qagiu (&F, 0.0, epsabs, epsrel,
238  limit, w, &result, &error);
239 
240  else gsl_integration_qag (&F, 0.0, z, epsabs, epsrel,
241  limit, key, w, &result, &error);
242 
243  gsl_integration_workspace_free (w);
244 
245  return result;
246 }
247 
248 
249 /**
250  * Function that integrates the uniform in comoving volume density to compute the normalisation factor.
251  * Consistently with Coward, Burman 2005 ( http://arxiv.org/abs/astro-ph/0505181 ) the integral should be always performed up to infinity.
252  * However, if the user specifies a value for zmax, the probability density will be normalised by integrating up to that value.
253  * To let the upper limit of integration go to infinity, specify zmax < 0.
254  * The integration is performed using gsl_integration_qagiu from the gsl library (see http://www.gnu.org/software/gsl/manual/html_node/QAGI-adaptive-integration-on-infinite-intervals.html )
255  */
257 {
258  double result = 0.0;
259  double error;
260  double epsabs = 5e-5;
261  double epsrel = 1e-5;
262  size_t limit = 1024;
263  int key = 1;
264 
265  gsl_function F;
266  F.function = &XLALUniformComovingVolumeDensity;
267  F.params = omega;
268 
269  gsl_integration_workspace * w
270  = gsl_integration_workspace_alloc (1024);
271 
272  if (z<0.0) gsl_integration_qagiu (&F, 0.0, epsabs, epsrel,
273  limit, w, &result, &error);
274 
275  else gsl_integration_qag (&F, 0.0, z, epsabs, epsrel,
276  limit, key, w, &result, &error);
277 
278  gsl_integration_workspace_free (w);
279 
280  return result;
281 }
282 /**
283  * Creates a LALCosmologicalParameters structure from the values of the cosmological parameters.
284  * Note that the boundary condition \f$\Omega_m + \Omega_k + \Omega_\Lambda = 1\f$ is imposed here.
285  */
286 LALCosmologicalParameters *XLALCreateCosmologicalParameters(double h, double om, double ol, double w0, double w1, double w2)
287 {
289  p->h = h;
290  p->om=om;
291  p->ol=ol;
292  p->ok=1.0-om-ol;
293  p->w0=w0;
294  p->w1=w1;
295  p->w2=w2;
296  return p;
297 }
298 /**
299  * Destroys a LALCosmologicalParameters structure.
300  */
301 
303 {
304  LALFree(omega);
305 }
306 /**
307  * The next set of functions return specific parameters from the LALCosmologicalParameters structure
308  */
309 
311 {
312  return omega->h;
313 }
315 {
316  return omega->om;
317 }
319 {
320  return omega->ol;
321 }
323 {
324  return omega->ok;
325 }
327 {
328  return omega->w0;
329 }
331 {
332  return omega->w1;
333 }
335 {
336  return omega->w2;
337 }
338 /**
339  * Function to set a LALCosmologicalParameters structure to the default LambdaCDM value (see http://arxiv.org/abs/1303.5076 )
340  */
342 {
343  omega->h = 0.671;
344  omega->om=0.3175;
345  omega->ok=1.0-0.3175-0.68251;
346  omega->ol=0.68251;
347  omega->w0=-1.0;
348  omega->w1=0.0;
349  omega->w2=0.0;
350 }
351 /**
352  * Function to create a LALCosmologicalRateParameters structure
353  */
355 {
357  p->r0=r0;
358  p->Q=Q;
359  p->W=W;
360  p->R=R;
361  return p;
362 }
363 /**
364  * Destroys a LALCosmologicalParameters structure.
365  */
366 
368 {
369  LALFree(rate);
370 }
371 /**
372  * Returns the local rate
373  */
375 {
376  return rate->r0;
377 }
378 /**
379  * Implements the fit to the SFR in Eq.7 of Coward, Burman 2005 ( http://arxiv.org/abs/astro-ph/0505181 )
380  * See also Porciani & Madau ( http://arxiv.org/pdf/astro-ph/0008294v2.pdf ) and references therein
381  */
382 double XLALStarFormationDensity(double z, void *params)
383 {
385  double hz = XLALHubbleParameter(z,p->omega);
386  double x = 1.0/sqrt(1.+z);
387  double hz0 = x*x*x;
388  return hz0/hz*p->rate->r0*(1.0+p->rate->W)*exp(p->rate->Q*z)/(exp(p->rate->R*z)+p->rate->W);
389 }
390 /**
391  * Returns the Rate weighted uniform comoving volume density
392  */
393 double XLALRateWeightedUniformComovingVolumeDensity(double z, void *params)
394 {
396  double dvdz = XLALUniformComovingVolumeDensity(z,p->omega);
397  double ez = XLALStarFormationDensity(z,p);
398  return ez*dvdz;
399 }
400 
401 /**
402  * Function that integrates the uniform in comoving volume density to compute the normalisation factor.
403  * Consistently with Coward, Burman 2005 ( http://arxiv.org/abs/astro-ph/0505181 ) the integral should be always performed up to infinity.
404  * However, if the user specifies a value for zmax, the probability density will be normalised by integrating up to that value.
405  * To let the upper limit of integration go to infinity, specify zmax < 0.
406  * The integration is performed using gsl_integration_qagiu from the gsl library (see http://www.gnu.org/software/gsl/manual/html_node/QAGI-adaptive-integration-on-infinite-intervals.html )
407  */
408 
410 {
411  double result = 0.0;
412  double error;
413  double epsabs = 5e-5;
414  double epsrel = 1e-5;
415  size_t limit = 1024;
416  int key = 1;
417 
418  gsl_function F;
420  F.params = p;
421 
422  gsl_integration_workspace * w
423  = gsl_integration_workspace_alloc (1024);
424 
425  if (z<0.0) gsl_integration_qagiu (&F, 0.0, epsabs, epsrel,
426  limit, w, &result, &error);
427 
428  else gsl_integration_qag (&F, 0.0, z, epsabs, epsrel,
429  limit, key, w, &result, &error);
430 
431  gsl_integration_workspace_free (w);
432 
433  return result;
434 }
435 /**
436  * Returns the source redshifts distribution function obtained by normalizing the differential rate dR/dz integrated throughout the cosmos, as seen in our frame.
437  * It gives the probability density distribution function of an event occurring in the redshift range 0 to z, see Eq.12 in Coward & Burman ( http://arxiv.org/pdf/astro-ph/0505181v1.pdf )
438  */
440 {
442  double unnorm_density = XLALRateWeightedUniformComovingVolumeDensity(z,(void *)p);
443  return unnorm_density/norm;
444 }
445 /**
446  * Creates a cosmological parameters and rate structure
447  */
449 {
451  p->omega=XLALCreateCosmologicalParameters(0.0,0.0,0.0,0.0,0.0,0.0);
452  p->rate=XLALCreateCosmologicalRateParameters(0.0,0.0,0.0,0.0);
453  return p;
454 }
455 /**
456  * Destroys a cosmological parameters and rate structure
457  */
459 {
462  LALFree(p);
463 }
464 /**
465  * Function to set a LALCosmologicalRateParameters structure to the default independent of z value
466  */
468 {
469  params->r0=5E-12;
470  params->W= 0.0;
471  params->Q= 0.0;
472  params->R= 0.0;
473 }
double XLALGetOmegaLambda(LALCosmologicalParameters *omega)
double XLALComovingVolume(LALCosmologicalParameters *omega, double z)
This function computes the comoving volume between 0 and z.
double XLALGetW2(LALCosmologicalParameters *omega)
double XLALLuminosityDistance(LALCosmologicalParameters *omega, double z)
The set of functions in this module implements the standard cosmological distance measures defined a ...
double XLALGetLocalRate(LALCosmologicalRateParameters *rate)
Returns the local rate.
double XLALGetHubbleConstant(LALCosmologicalParameters *omega)
The next set of functions return specific parameters from the LALCosmologicalParameters structure.
double XLALIntegrateComovingVolume(LALCosmologicalParameters *omega, double z)
Function that integrates the comoving volume element.
double XLALGetOmegaK(LALCosmologicalParameters *omega)
double XLALRateWeightedUniformComovingVolumeDensity(double z, void *params)
Returns the Rate weighted uniform comoving volume density.
double XLALIntegrateComovingVolumeDensity(LALCosmologicalParameters *omega, double z)
Function that integrates the uniform in comoving volume density to compute the normalisation factor.
double XLALComovingTransverseDistance(LALCosmologicalParameters *omega, double z)
Computes the comoving transverse distance from the comoving line-of-sight distance (proper distance) ...
double XLALHubbleParameter(double z, void *omega)
Computes the inverse of the Hubble parameter at redshift z Eq.
LALCosmologicalRateParameters * XLALCreateCosmologicalRateParameters(double r0, double W, double Q, double R)
Function to create a LALCosmologicalRateParameters structure.
double XLALUniformComovingVolumeDistribution(LALCosmologicalParameters *omega, double z, double zmax)
This function computes the value of a uniform probability distribution over the comoving volume.
double XLALGetOmegaMatter(LALCosmologicalParameters *omega)
double XLALComovingVolumeElement(double z, void *omega)
This function computes the comoving volume element (a 4&pi shell) between redshift z and z+dz.
double XLALIntegrateHubbleParameter(LALCosmologicalParameters *omega, double z)
Computes the integral of inverse of the Hubble parameter at redshift z.
LALCosmologicalParametersAndRate * XLALCreateCosmologicalParametersAndRate(void)
Creates a cosmological parameters and rate structure.
double XLALRateWeightedComovingVolumeDistribution(LALCosmologicalParametersAndRate *p, double z, double zmax)
Returns the source redshifts distribution function obtained by normalizing the differential rate dR/d...
double XLALGetW0(LALCosmologicalParameters *omega)
void XLALDestroyCosmologicalParametersAndRate(LALCosmologicalParametersAndRate *p)
Destroys a cosmological parameters and rate structure.
double XLALComovingLOSDistance(LALCosmologicalParameters *omega, double z)
Computes the comoving line-of-sight distance (usually called the proper distance) by integrating the ...
double XLALAngularDistance(LALCosmologicalParameters *omega, double z)
Computes the angular size distance by dividing the comoving transverse distance by 1+z Eq.
void XLALDestroyCosmologicalParameters(LALCosmologicalParameters *omega)
Destroys a LALCosmologicalParameters structure.
double XLALUniformComovingVolumeDensity(double z, void *omega)
This function computes the value of a uniform probability density distribution over the comoving volu...
double XLALStarFormationDensity(double z, void *params)
Implements the fit to the SFR in Eq.7 of Coward, Burman 2005 ( http://arxiv.org/abs/astro-ph/0505181 ...
double XLALGetW1(LALCosmologicalParameters *omega)
double XLALHubbleDistance(LALCosmologicalParameters *omega)
Computes the Hubble distance, independent of the redshift.
void XLALSetCosmologicalParametersDefaultValue(LALCosmologicalParameters *omega)
Function to set a LALCosmologicalParameters structure to the default LambdaCDM value (see http://arxi...
void XLALSetCosmologicalRateParametersDefaultValue(LALCosmologicalRateParameters *params)
Function to set a LALCosmologicalRateParameters structure to the default independent of z value.
void XLALDestroyCosmologicalRateParameters(LALCosmologicalRateParameters *rate)
Destroys a LALCosmologicalParameters structure.
LALCosmologicalParameters * XLALCreateCosmologicalParameters(double h, double om, double ol, double w0, double w1, double w2)
Creates a LALCosmologicalParameters structure from the values of the cosmological parameters.
double XLALIntegrateRateWeightedComovingVolumeDensity(LALCosmologicalParametersAndRate *p, double z)
Function that integrates the uniform in comoving volume density to compute the normalisation factor.
#define LALFree(p)
Definition: LALMalloc.h:96
#define fprintf
static REAL8TimeSeries * error(const REAL8TimeSeries *s1, const REAL8TimeSeries *s0)
static double Q(double mu, double x)
Definition: XLALMarcumQ.c:99