Loading [MathJax]/extensions/TeX/AMSsymbols.js
LAL 7.7.0.1-3a66518
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
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/LALConstants.h>
20#include <lal/LALCosmologyCalculator.h>
21
22#include <stdio.h>
23#include <stdlib.h>
24#include <string.h>
25#include <math.h>
26#include <gsl/gsl_const_mksa.h>
27#include <gsl/gsl_integration.h>
28
29/**
30 * The set of functions in this module implements the standard cosmological distance measures
31 * defined a Friedmann-Robertson-Walker-Lemaitre metric.
32 * For a detailed reference to the various formulae, see Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 ).
33 */
34
35/**
36 * Computes the luminosity distance as the angular distance divided by (1+z)^2
37 * Eq. 21 in Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 )
38 */
39
42 double z)
43{
44 double x=1.0+z;
45 return XLALComovingTransverseDistance(omega, z)*x;
46}
47
48/**
49 * Computes the angular size distance by dividing the comoving transverse distance by 1+z
50 * Eq. 18 in Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 )
51 */
52
55 double z)
56{
57 double x=1.0+z;
58 return XLALComovingTransverseDistance(omega, z)/x;
59}
60
61/**
62 * Computes the comoving line-of-sight distance (usually called the proper distance) by integrating the Hubble parameter.
63 * Eq. 15 in Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 ).
64 */
65
68 double z)
69{
70 double dH = XLALHubbleDistance(omega);
71 return dH*XLALIntegrateHubbleParameter(omega,z);
72}
73/**
74 * Computes the comoving transverse distance from the comoving line-of-sight distance (proper distance) and
75 * Eq. 16 in Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 )
76 */
79 double z)
80{
81 double ok = omega->ok;
82 double dH = XLALHubbleDistance(omega);
83 double dC = XLALComovingLOSDistance(omega,z);
84
85 if (fabs(ok)<1e-7)
86 {
87
88 return dC;
89 }
90 else if (ok>1e-7)
91 {
92 return dH*sinh(sqrt(ok)*dC/dH)/sqrt(ok);
93 }
94 else if (ok<1e-7)
95 {
96 return dH*sin(sqrt(fabs(ok))*dC/dH)/sqrt(fabs(ok));
97 }
98 else
99 {
100 fprintf(stderr,"Something funny happened. Aborting.\n");
101 exit(-1);
102 }
103
104}
105/**
106 * Computes the Hubble distance, independent of the redshift. This is the distance with sets the units for all others.
107 * It returns the distance in Mpc.
108 * Eq. 4 in Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 )
109 */
112 )
113{
114 double x = GSL_CONST_MKSA_SPEED_OF_LIGHT/(100.0*omega->h);
115 return 1e-3*x;
116}
117
118/**
119 * Computes the inverse of the Hubble parameter at redshift z
120 * Eq. 14 in Hogg 1999 ( http://arxiv.org/abs/astro-ph/9905116 )
121 */
122
123double XLALHubbleParameter(double z,
124 void *omega
125 )
126{
128
129 double E=0.0;
130 double x = 1.0+z;
131 double om=p->om;
132 double ol=p->ol;
133 double ok=p->ok;
134 double w0=p->w0;
135 double w1=p->w1;
136 double w2=p->w2;
137 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))));
138 return 1.0/E;
139}
140/**
141 * Computes the integral of inverse of the Hubble parameter at redshift z.
142 * The integral is computed using the built-in function gsl_integration_qag of the gsl library.
143 * The integral is performed using a Gauss-Kronrod adaptive method (see http://www.gnu.org/software/gsl/manual/html_node/QAG-adaptive-integration.html )
144 */
146{
147 double result = 0.0;
148 double error;
149 double epsabs = 5e-5;
150 double epsrel = 1e-5;
151 size_t limit = 1024;
152 int key = 1;
153
154 gsl_function F;
155 F.function = &XLALHubbleParameter;
156 F.params = omega;
157
158 gsl_integration_workspace * w
159 = gsl_integration_workspace_alloc (1024);
160
161 gsl_integration_qag (&F, 0.0, z, epsabs, epsrel,
162 limit, key, w, &result, &error);
163
164 gsl_integration_workspace_free (w);
165
166 return result;
167}
168/**
169 * This function computes the value of a uniform probability distribution over the comoving volume.
170 * Details of the derivation of these formulae can be found in Coward, Burman 2005 ( http://arxiv.org/abs/astro-ph/0505181 )
171 */
174 double z,
175 double zmax)
176{
177 double norm = XLALIntegrateComovingVolumeDensity(omega, zmax);
178 double unnorm_density = XLALUniformComovingVolumeDensity(z,(void *)omega);
179 return unnorm_density/norm;
180}
181
182/**
183 * This function computes the value of a uniform probability density distribution over the comoving volume at redshift z.
184 */
185
187 double z,
188 void *omega)
189{
190
192
193 double x = 1.0+z;
194 double unnorm_density = XLALComovingVolumeElement(z,p)/x;
195 return unnorm_density;
196}
197/**
198 * This function computes the comoving volume between 0 and z
199 */
202 double z)
203{
204 double V = XLALIntegrateComovingVolume(omega, z);
205 return V;
206}
207
208/**
209 * This function computes the comoving volume element (a 4&pi shell) between redshift z and z+dz.
210 */
211
213 double z,
214 void *omega)
215{
216
218
219 double dm = XLALComovingTransverseDistance(omega,z);
220 double E = XLALHubbleParameter(z,omega);
221 double dVdz = 4.0*M_PI*dm*dm*E*XLALHubbleDistance(p);
222 return dVdz;
223}
224
225/**
226 * Function that integrates the comoving volume element.
227 * 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 )
228 */
230{
231 double result = 0.0;
232 double error;
233 double epsabs = 5e-5;
234 double epsrel = 1e-5;
235 size_t limit = 1024;
236 int key = 1;
237
238 gsl_function F;
239 F.function = &XLALComovingVolumeElement;
240 F.params = omega;
241
242 gsl_integration_workspace * w
243 = gsl_integration_workspace_alloc (1024);
244
245 if (z<0.0) gsl_integration_qagiu (&F, 0.0, epsabs, epsrel,
246 limit, w, &result, &error);
247
248 else gsl_integration_qag (&F, 0.0, z, epsabs, epsrel,
249 limit, key, w, &result, &error);
250
251 gsl_integration_workspace_free (w);
252
253 return result;
254}
255
256
257/**
258 * Function that integrates the uniform in comoving volume density to compute the normalisation factor.
259 * Consistently with Coward, Burman 2005 ( http://arxiv.org/abs/astro-ph/0505181 ) the integral should be always performed up to infinity.
260 * However, if the user specifies a value for zmax, the probability density will be normalised by integrating up to that value.
261 * To let the upper limit of integration go to infinity, specify zmax < 0.
262 * 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 )
263 */
265{
266 double result = 0.0;
267 double error;
268 double epsabs = 5e-5;
269 double epsrel = 1e-5;
270 size_t limit = 1024;
271 int key = 1;
272
273 gsl_function F;
275 F.params = omega;
276
277 gsl_integration_workspace * w
278 = gsl_integration_workspace_alloc (1024);
279
280 if (z<0.0) gsl_integration_qagiu (&F, 0.0, epsabs, epsrel,
281 limit, w, &result, &error);
282
283 else gsl_integration_qag (&F, 0.0, z, epsabs, epsrel,
284 limit, key, w, &result, &error);
285
286 gsl_integration_workspace_free (w);
287
288 return result;
289}
290
291/**
292 * Creates a LALCosmologicalParameters structure from the values of the cosmological parameters.
293 * Note that the boundary condition \f$\Omega_m + \Omega_k + \Omega_\Lambda = 1\f$ is imposed here.
294 */
295LALCosmologicalParameters *XLALCreateCosmologicalParameters(double h, double om, double ol, double w0, double w1, double w2)
296{
298 p->h = h;
299 p->om=om;
300 p->ol=ol;
301 p->ok=1.0-om-ol;
302 p->w0=w0;
303 p->w1=w1;
304 p->w2=w2;
305 return p;
306}
307
309{
312 return p;
313}
314
315/**
316 * Destroys a LALCosmologicalParameters structure.
317 */
318
320{
321 LALFree(omega);
322}
323/**
324 * The next set of functions return specific parameters from the LALCosmologicalParameters structure
325 */
326
328{
329 return omega->h;
330}
332{
333 return omega->om;
334}
336{
337 return omega->ol;
338}
340{
341 return omega->ok;
342}
344{
345 return omega->w0;
346}
348{
349 return omega->w1;
350}
352{
353 return omega->w2;
354}
355/**
356 * Function to set a LALCosmologicalParameters structure to the default FlatLambdaCDM
357 */
359{
360 /* Flat LambdaCDM */
361 omega->h = LAL_H0_DIMENSIONLESS;
362 omega->om = LAL_OMEGA_M;
363 omega->ok = 0.0; /* flat */
364 omega->ol = 1.0 - omega->om - omega->ok; /* closed */
365 omega->w0 = -1.0; /* cosmological constant */
366 omega->w1 = 0.0;
367 omega->w2 = 0.0;
368}
369/**
370 * Function to create a LALCosmologicalRateParameters structure
371 */
373{
375 p->r0=r0;
376 p->Q=Q;
377 p->W=W;
378 p->R=R;
379 return p;
380}
381/**
382 * Destroys a LALCosmologicalParameters structure.
383 */
384
386{
387 LALFree(rate);
388}
389/**
390 * Returns the local rate
391 */
393{
394 return rate->r0;
395}
396/**
397 * Implements the fit to the SFR in Eq.7 of Coward, Burman 2005 ( http://arxiv.org/abs/astro-ph/0505181 )
398 * See also Porciani & Madau ( http://arxiv.org/pdf/astro-ph/0008294v2.pdf ) and references therein
399 */
400double XLALStarFormationDensity(double z, void *params)
401{
403 double hz = XLALHubbleParameter(z,p->omega);
404 double x = 1.0/sqrt(1.+z);
405 double hz0 = x*x*x;
406 return hz0/hz*p->rate->r0*(1.0+p->rate->W)*exp(p->rate->Q*z)/(exp(p->rate->R*z)+p->rate->W);
407}
408/**
409 * Returns the Rate weighted uniform comoving volume density
410 */
412{
414 double dvdz = XLALUniformComovingVolumeDensity(z,p->omega);
415 double ez = XLALStarFormationDensity(z,p);
416 return ez*dvdz;
417}
418
419/**
420 * Function that integrates the uniform in comoving volume density to compute the normalisation factor.
421 * Consistently with Coward, Burman 2005 ( http://arxiv.org/abs/astro-ph/0505181 ) the integral should be always performed up to infinity.
422 * However, if the user specifies a value for zmax, the probability density will be normalised by integrating up to that value.
423 * To let the upper limit of integration go to infinity, specify zmax < 0.
424 * 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 )
425 */
426
428{
429 double result = 0.0;
430 double error;
431 double epsabs = 5e-5;
432 double epsrel = 1e-5;
433 size_t limit = 1024;
434 int key = 1;
435
436 gsl_function F;
438 F.params = p;
439
440 gsl_integration_workspace * w
441 = gsl_integration_workspace_alloc (1024);
442
443 if (z<0.0) gsl_integration_qagiu (&F, 0.0, epsabs, epsrel,
444 limit, w, &result, &error);
445
446 else gsl_integration_qag (&F, 0.0, z, epsabs, epsrel,
447 limit, key, w, &result, &error);
448
449 gsl_integration_workspace_free (w);
450
451 return result;
452}
453/**
454 * Returns the source redshifts distribution function obtained by normalizing the differential rate dR/dz integrated throughout the cosmos, as seen in our frame.
455 * 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 )
456 */
458{
460 double unnorm_density = XLALRateWeightedUniformComovingVolumeDensity(z,(void *)p);
461 return unnorm_density/norm;
462}
463/**
464 * Creates a cosmological parameters and rate structure
465 */
467{
469 p->omega=XLALCreateCosmologicalParameters(0.0,0.0,0.0,0.0,0.0,0.0);
470 p->rate=XLALCreateCosmologicalRateParameters(0.0,0.0,0.0,0.0);
471 return p;
472}
473/**
474 * Destroys a cosmological parameters and rate structure
475 */
477{
480 LALFree(p);
481}
482/**
483 * Function to set a LALCosmologicalRateParameters structure to the default independent of z value
484 */
486{
487 params->r0=5E-12;
488 params->W= 0.0;
489 params->Q= 0.0;
490 params->R= 0.0;
491}
double XLALGetOmegaLambda(LALCosmologicalParameters *omega)
double XLALComovingVolume(LALCosmologicalParameters *omega, double z)
This function computes the comoving volume between 0 and z.
double XLALGetW2(LALCosmologicalParameters *omega)
LALCosmologicalRateParameters * XLALCreateCosmologicalRateParameters(double r0, double W, double Q, double R)
Function to create a LALCosmologicalRateParameters structure.
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.
double XLALUniformComovingVolumeDistribution(LALCosmologicalParameters *omega, double z, double zmax)
This function computes the value of a uniform probability distribution over the comoving volume.
LALCosmologicalParameters * XLALCreateDefaultCosmologicalParameters(void)
LALCosmologicalParametersAndRate * XLALCreateCosmologicalParametersAndRate(void)
Creates a cosmological parameters and rate structure.
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.
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 FlatLambdaCDM.
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.
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.
double XLALIntegrateRateWeightedComovingVolumeDensity(LALCosmologicalParametersAndRate *p, double z)
Function that integrates the uniform in comoving volume density to compute the normalisation factor.
#define LALMalloc(n)
Definition: LALMalloc.h:93
#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
#define LAL_OMEGA_M
Default dimensionless Hubble constant, dimensionless.
Definition: LALConstants.h:650
#define LAL_H0_DIMENSIONLESS
Default dimensionless Hubble constant, dimensionless.
Definition: LALConstants.h:649