1 """Cosmological distance measures.
2
3 Mostly follows David Hogg's pedagogical paper arXiv:astro-ph/9905116v4 .
4
5 Distance units are Mpc, time units are seconds.
6
7 """
8
9 import math
10
11 import numpy
12 import scipy
13 import scipy.integrate as si
14 import scipy.interpolate
15 import scipy.optimize
16
17 import constants as cc
18
20 """'Spatial curvature density' omega_k_0 for a cosmology (if needed).
21
22 If omega_k_0 is specified, return it. Otherwise return:
23
24 1.0 - omega_M_0 - omega_lambda_0
25
26 """
27
28 if 'omega_k_0' in cosmo:
29 omega_k_0 = cosmo['omega_k_0']
30 else:
31 omega_k_0 = 1. - cosmo['omega_M_0'] - cosmo['omega_lambda_0']
32 return omega_k_0
33
35 """Returns the cosmo dictionary with omega_k_0 set.
36 See get_omega_k_0.
37
38 Note that cosmo is not passed as \*\*cosmo for once. This function
39 modifies the dictionary in place and returns the result.
40
41 """
42 if 'omega_k_0' in cosmo:
43 return cosmo
44 else:
45 cosmo['omega_k_0'] = get_omega_k_0(**cosmo)
46 return cosmo
47
48
49
50
52 """The unitless Hubble expansion rate at redshift z.
53
54 In David Hogg's (arXiv:astro-ph/9905116v4) formalism, this is
55 equivalent to E(z), defined in his eq. 14.
56
57 Modified (JBJ, 29-Feb-2012) to include scalar w parameter
58
59 """
60
61 if 'w' in cosmo:
62
63 return (cosmo['omega_M_0'] * (1+z)**3. +
64 cosmo['omega_k_0'] * (1+z)**2. +
65 cosmo['omega_lambda_0'] * (1+z)**(1+cosmo['w']) )**0.5
66 else:
67 return (cosmo['omega_M_0'] * (1+z)**3. +
68 cosmo['omega_k_0'] * (1+z)**2. +
69 cosmo['omega_lambda_0'])**0.5
70
72 """The value of the Hubble constant at redshift z.
73
74 Units are s^-1
75
76 In David Hogg's (arXiv:astro-ph/9905116v4) formalism, this is
77 equivalent to H_0 * E(z) (see his eq. 14).
78
79 """
80 H_0 = cosmo['h'] * cc.H100_s
81
82 return H_0 * e_z(z, **cosmo)
83
85 """The value of the Hubble distance at redshift z.
86
87 Units are Mpc.
88
89 In David Hogg's (arXiv:astro-ph/9905116v4) formalism, this is
90 equivalent to D_H / E(z) = c / (H_0 E(z)) [see his eq. 14], which
91 appears in the definitions of many other distance measures.
92
93 """
94 H_0 = cosmo['h'] * cc.H100_s
95
96 return cc.c_light_Mpc_s / (H_0 * e_z(z, **cosmo))
97
99
100 e_z = (omega_M_0 * (1+z)**3. +
101 omega_k_0 * (1+z)**2. +
102 omega_lambda_0 * (1+z)**(1.+w))**0.5
103
104 H_0 = h * cc.H100_s
105
106 H_z = H_0 * e_z
107
108 return cc.c_light_Mpc_s / (H_z)
109
111 """The derivative of the comoving distance with redshift: dd_c/dz.
112
113 See equation 15 of David Hogg's arXiv:astro-ph/9905116v4
114
115 Units are Mpc.
116
117 """
118 if 'w' in cosmo:
119 w = cosmo['w']
120 else:
121 w = -1.
122
123 return _comoving_integrand(z,
124 cosmo['omega_M_0'],
125 cosmo['omega_lambda_0'],
126 cosmo['omega_k_0'],
127 cosmo['h'], w)
128
130 """The line-of-sight comoving distance (in Mpc) to redshift z.
131
132 See equation 15 of David Hogg's arXiv:astro-ph/9905116v4
133
134 Units are Mpc.
135
136 Optionally calculate the integral from z0 to z.
137
138 Returns
139 -------
140
141 d_co: ndarray
142 Comoving distance in Mpc.
143
144 Examples
145 --------
146
147 >>> import cosmolopy.distance as cd
148 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
149 >>> cosmo = cd.set_omega_k_0(cosmo)
150 >>> d_co = cd.comoving_distance(6., **cosmo)
151 >>> print "Comoving distance to z=6 is %.1f Mpc" % (d_co)
152 Comoving distance to z=6 is 8017.8 Mpc
153
154 """
155
156
157
158 if 'w' in cosmo:
159 w = cosmo['w']
160 else:
161 w = -1.
162
163 dc_func = \
164 numpy.vectorize(lambda z, z0, omega_M_0, omega_lambda_0, omega_k_0, h, w:
165 si.quad(_comoving_integrand, z0, z, limit=1000,
166 args=(omega_M_0, omega_lambda_0, omega_k_0, h, w)))
167 d_co, err = dc_func(z, z0,
168 cosmo['omega_M_0'],
169 cosmo['omega_lambda_0'],
170 cosmo['omega_k_0'],
171 cosmo['h'],
172 w
173 )
174 return d_co
175
177 """Returns comoving_distance_transverse.
178
179 Units are Mpc.
180
181 Examples
182 --------
183
184 >>> import cosmolopy.distance as cd
185 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
186 >>> cosmo = cd.set_omega_k_0(cosmo)
187 >>> d_M = cd.proper_motion_distance(6., **cosmo)
188 >>> print "Transverse comoving distance to z=6 is %.1f Mpc" % (d_M)
189 Transverse comoving distance to z=6 is 8017.8 Mpc
190
191 """
192 return comoving_distance_transverse(z, **cosmo)
193
195 """The transverse comoving distance (in Mpc) to redshift z.
196
197 This is also called the proper motion distance, D_M.
198
199 See equation 16 of David Hogg's arXiv:astro-ph/9905116v4
200
201 Units are Mpc.
202
203 This is the distance d_m, such that the comoving distance between
204 two events at the same redshift, but separated on the sky by some
205 angle delta_theta is d_m * delta_theta.
206
207 Examples
208 --------
209
210 >>> import cosmolopy.distance as cd
211 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
212 >>> cosmo = cd.set_omega_k_0(cosmo)
213 >>> d_M = cd.comoving_distance_transverse(6., **cosmo)
214 >>> print "Transverse comoving distance to z=6 is %.1f Mpc" % (d_M)
215 Transverse comoving distance to z=6 is 8017.8 Mpc
216
217 """
218
219 d_c = comoving_distance(z, 0.0, **cosmo)
220
221 omega_k_0 = get_omega_k_0(**cosmo)
222
223 if numpy.all(omega_k_0 == 0.0):
224 return d_c
225
226 d_h_0 = hubble_distance_z(0.0, **cosmo)
227 sqrt_ok0 = numpy.sqrt(numpy.abs(omega_k_0))
228 if not numpy.isscalar(omega_k_0):
229 sqrt_ok0[omega_k_0 == 0.0] = 1.0
230 argument = sqrt_ok0 * d_c / d_h_0
231 factor = d_h_0 * (1./sqrt_ok0)
232 d_m = ((omega_k_0 > 0.0) * (factor * numpy.sinh(argument)) +
233 (omega_k_0 == 0.0) * d_c +
234 (omega_k_0 < 0.0) * (factor * numpy.sin(argument)))
235
236 return d_m
237
239 """The angular-diameter distance (Mpc) to redshift z.
240
241 Optionally find the angular diameter distance between objects at
242 z0 and z (only implemented for omega_k_0 >= 0).
243
244 See equations 18-19 of David Hogg's arXiv:astro-ph/9905116v4
245
246 Units are Mpc.
247
248 Examples
249 --------
250
251 >>> import cosmolopy.distance as cd
252 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
253 >>> cosmo = cd.set_omega_k_0(cosmo)
254 >>> d_a = cd.angular_diameter_distance(6., **cosmo)
255 >>> print "Angular diameter distance = %.1f Mpc" % (d_a)
256 Angular diameter distance = 1145.4 Mpc
257
258 """
259
260 omega_k = numpy.atleast_1d(get_omega_k_0(**cosmo))
261 if (numpy.any(omega_k < 0) and not(z0 == 0)):
262 raise ValueError("Not implemented for Omega_k < 0 and z0 > 0")
263
264 dm2 = comoving_distance_transverse(z, **cosmo)
265 if z0 == 0:
266 return dm2 / (1. + z)
267
268 dm1 = comoving_distance_transverse(z0, **cosmo)
269
270 d_h_0 = hubble_distance_z(0.0, **cosmo)
271
272 term1 = dm1 * numpy.sqrt(1. + omega_k * (dm2/d_h_0)**2.)
273 term2 = dm2 * numpy.sqrt(1. + omega_k * (dm1/d_h_0)**2.)
274
275 da12 = (term2 - term1)/(1+z)
276
277 return da12
278
280 """The luminosity distance to redshift z.
281
282 Units are Mpc.
283
284 See, for example, David Hogg's arXiv:astro-ph/9905116v4
285
286 """
287 da = angular_diameter_distance(z, **cosmo)
288 return da * (1+z)**2.
289
291 """The differential comoving volume element dV_c/dz/dSolidAngle.
292
293 Dimensions are volume per unit redshift per unit solid angle.
294
295 Units are Mpc**3 Steradians^-1.
296
297 See David Hogg's arXiv:astro-ph/9905116v4, equation 28.
298
299 Examples
300 --------
301
302 >>> import cosmolopy.distance as cd
303 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
304 >>> cosmo = cd.set_omega_k_0(cosmo)
305 >>> dVc = cd.diff_comoving_volume(6.0, **cosmo)
306 >>> print "dV/dz/dSolidAngle at z=6 is %.3g Mpc**3" % (dVc)
307 dV/dz/dSolidAngle at z=6 is 2.63e+10 Mpc**3
308 """
309
310 d_h_0 = hubble_distance_z(0.0, **cosmo)
311 d_m = comoving_distance_transverse(z, **cosmo)
312 ez = e_z(z, **cosmo)
313 return d_h_0 * d_m**2. / ez
314
316 """The comoving volume out to redshift z.
317
318 See David Hogg's arXiv:astro-ph/9905116v4, equation 29.
319
320 Examples
321 --------
322
323 >>> import cosmolopy.distance as cd
324 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
325 >>> cosmo = cd.set_omega_k_0(cosmo)
326 >>> Vc = cd.comoving_volume(6.0, **cosmo)
327 >>> print "Vc = %.3g Mpc**3" % (Vc)
328 Vc = 2.16e+12 Mpc**3
329
330
331 >>> import cosmolopy.distance as cd
332 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.0, 'h' : 0.72}
333 >>> cosmo = cd.set_omega_k_0(cosmo)
334 >>> Vc = cd.comoving_volume(6.0, **cosmo)
335 >>> print "Vc = %.3g Mpc**3" % (Vc)
336 Vc = 1.68e+12 Mpc**3
337
338 """
339
340 dm = comoving_distance_transverse(z, **cosmo)
341
342 omega_k_0 = get_omega_k_0(**cosmo)
343
344 flat_volume = 4. * numpy.pi * dm**3. / 3.
345
346 if numpy.all(omega_k_0 == 0.0):
347 return flat_volume
348
349 d_h_0 = hubble_distance_z(0.0, **cosmo)
350
351 sqrt_ok0 = numpy.sqrt(numpy.abs(omega_k_0))
352 dmdh = dm/d_h_0
353 argument = sqrt_ok0 * dmdh
354 f1 = 4. * numpy.pi * d_h_0**3. / (2. * omega_k_0)
355 f2 = dmdh * numpy.sqrt(1. + omega_k_0 * (dmdh)**2.)
356 f3 = 1./sqrt_ok0
357
358 if numpy.isscalar(omega_k_0):
359 if omega_k_0 > 0.0:
360 return f1 * (f2 - f3 * numpy.arcsinh(argument))
361 elif omega_k_0 == 0.0:
362 return flat_volume
363 elif omega_k_0 < 0.0:
364 return f1 * (f2 - f3 * numpy.arcsin(argument))
365 else:
366 b = numpy.broadcast(omega_k_0,z,dm)
367 Vc = numpy.zeros(b.shape)
368 m1 = numpy.ones(b.shape, dtype=bool) * (omega_k_0 > 0.0)
369 Vc[m1] = (f1 * (f2 - f3 * numpy.arcsinh(argument)))[m1]
370
371 m1 = numpy.ones(b.shape, dtype=bool) * (omega_k_0 == 0.0)
372 Vc[m1] = flat_volume[m1]
373
374 m1 = numpy.ones(b.shape, dtype=bool) * (omega_k_0 < 0.0)
375 Vc[m1] = (f1 * (f2 - f3 * numpy.arcsin(argument)))[m1]
376 return Vc
377
379
380 e_z = (omega_M_0 * (1+z)**3. +
381 omega_k_0 * (1+z)**2. +
382 omega_lambda_0)**0.5
383
384 H_0 = h * cc.H100_s
385
386 H_z = H_0 * e_z
387
388 return 1./((1. + z) * H_z)
389
391 """The derivative of the lookback time with redshift: dt_L/dz.
392
393 See equation 30 of David Hogg's arXiv:astro-ph/9905116v4
394
395 Units are seconds.
396
397 """
398 return _lookback_integrand(z,
399 cosmo['omega_M_0'],
400 cosmo['omega_lambda_0'],
401 cosmo['omega_k_0'],
402 cosmo['h'])
403
405 """The lookback time (in s) to redshift z.
406
407 See equation 30 of David Hogg's arXiv:astro-ph/9905116v4
408
409 Units are s.
410
411 Optionally calculate the integral from z0 to z.
412
413 Returns
414 -------
415
416 t_look: ndarray
417 Lookback time in seconds.
418
419 """
420
421
422
423 lt_func = \
424 numpy.vectorize(lambda z, z0, omega_M_0, omega_lambda_0, omega_k_0, h:
425 si.quad(_lookback_integrand, z0, z, limit=1000,
426 args=(omega_M_0, omega_lambda_0, omega_k_0, h)))
427 t_look, err = lt_func(z, z0,
428 cosmo['omega_M_0'],
429 cosmo['omega_lambda_0'],
430 cosmo['omega_k_0'],
431 cosmo['h']
432 )
433 return t_look
434
435 -def age(z, use_flat=True, **cosmo):
436 """The age of the universe as seen at redshift z.
437
438 Age at z is lookback time at z'->Infinity minus lookback time at z.
439
440 See also: lookback_time.
441
442 Units are s.
443
444 Examples
445 --------
446
447 >>> import cosmolopy.distance as cd
448 >>> import cosmolopy.constants as cc
449 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
450 >>> cosmo = cd.set_omega_k_0(cosmo)
451 >>> t = cd.age(6.0, **cosmo)
452 >>> print "age at z=6.0 = %.3g Gyr" % (t/cc.Gyr_s)
453 age at z=6.0 = 0.892 Gyr
454
455 """
456 if use_flat and numpy.all(get_omega_k_0(**cosmo) == 0):
457 return age_flat(z, **cosmo)
458 fullage = lookback_time(numpy.Inf, **cosmo)
459 tl = lookback_time(z, **cosmo)
460 age = fullage - tl
461 return age
462
464 """The age of the universe assuming a flat cosmology.
465
466 Units are s.
467
468 Analytical formula from Peebles, p. 317, eq. 13.2.
469
470 Examples
471 --------
472
473 >>> import cosmolopy.distance as cd
474 >>> import cosmolopy.constants as cc
475 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
476 >>> cosmo = cd.set_omega_k_0(cosmo)
477 >>> t = cd.age_flat(6.0, **cosmo)
478 >>> print "age at z=6.0 is %.3g Gyr" % (t/cc.Gyr_s)
479 age at z=6.0 is 0.892 Gyr
480
481 """
482
483 omega_k = get_omega_k_0(**cosmo)
484 if (numpy.any(omega_k != 0)):
485
486 print "Warning: using lambda = 1 - omega_M for non-flat cosmology!"
487
488 om = cosmo['omega_M_0']
489 lam = 1. - cosmo['omega_M_0']
490 t_z = (2. *
491 numpy.arcsinh(numpy.sqrt(lam/om) * (1. + z)**(-3./2.)) /
492 (cc.H100_s * cosmo['h'] * 3. * numpy.sqrt(lam)))
493
494 return t_z
495
496 -def quick_distance_function(function, zmax = 20., zmin = 0., zstep = 0.001,
497 return_inverse=False, k=3,
498 **cosmo):
499 """Return an interpolation function that will give distance as a
500 funtion of z
501
502 If return_inverse is True, will also return a function giving z as
503 a function of distance.
504
505 Inputs
506 ------
507
508 function -- the distance function to interpolate (can be any
509 callable that takes a redshift argument plus cosmology keywords).
510
511 k -- spline order (`scipy.interpolate.InterpolatedUnivariateSpline`)
512
513 Returns
514 -------
515
516 distfunc
517
518 or
519
520 distfunc, zfunc
521
522 Examples
523 --------
524
525 >>> import cosmolopy.distance as cd
526 >>> import cosmolopy.constants as cc
527 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
528 >>> cosmo = cd.set_omega_k_0(cosmo)
529 >>> distfunc, redfunc = cd.quick_distance_function(cd.luminosity_distance, return_inverse=True, **cosmo)
530 >>> d = distfunc(6.3333)
531 >>> z = redfunc(d)
532 >>> "%.1g" % (distfunc(6.3333)/cd.luminosity_distance(6.3333, **cosmo) - 1.0)
533 '-2e-16'
534 >>> "%.1g" % (z/6.3333 - 1.0)
535 '0'
536
537 """
538 z = numpy.linspace(zmin, zmax, math.ceil((zmax-zmin)/zstep))
539 dists = function(z, **cosmo)
540 distfunc = scipy.interpolate.InterpolatedUnivariateSpline(z, dists, k=k)
541 if return_inverse:
542 redfunc = scipy.interpolate.InterpolatedUnivariateSpline(dists, z, k=k)
543 return distfunc, redfunc
544 else:
545 return distfunc
546
547 -def quick_age_function(zmax = 20., zmin = 0., zstep = 0.001,
548 return_inverse=False,
549 **cosmo):
550 """Return an interpolation function that will give age as a funtion of z
551
552 Units are s.
553
554 If return_inverse is True, will also return a function giving z as
555 a function of age.
556
557 Returns
558 -------
559
560 agefunc
561
562 or
563
564 agefunc, redfunc
565
566 Examples
567 --------
568
569 >>> import cosmolopy.distance as cd
570 >>> import cosmolopy.constants as cc
571 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
572 >>> cosmo = cd.set_omega_k_0(cosmo)
573 >>> agefunc = cd.quick_age_function(**cosmo)
574 >>> t = agefunc(6.0)
575 >>> print "age at z=6.0 is %.3g Gyr" % (t/cc.Gyr_s)
576 age at z=6.0 is 0.892 Gyr
577
578
579 """
580 z = numpy.arange(zmin, zmax, zstep)
581 ages = age(z, **cosmo)
582 agefunc = scipy.interpolate.interp1d(z, ages)
583 if return_inverse:
584 redfunc = scipy.interpolate.interp1d(ages[::-1], z[::-1])
585 return agefunc, redfunc
586 else:
587 return agefunc
588
590 """Return an interpolation function giving z as a funtion of age
591 of the universe.
592
593 Units of time are s.
594
595 Returns
596 -------
597
598 redfunc
599
600 Examples
601 --------
602
603 >>> import cosmolopy.distance as cd
604 >>> import cosmolopy.constants as cc
605 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
606 >>> cosmo = cd.set_omega_k_0(cosmo)
607 >>> redfunc = cd.quick_redshift_age_function(**cosmo)
608 >>> z = redfunc(1.0 * cc.Gyr_s)
609 >>> print "When age=1.0Gyr z=%.2f" % (z)
610 When age=1.0Gyr z=5.49
611
612 """
613 z = numpy.arange(zmin, zmax, zstep)
614 z = z[::-1]
615 ages = age(z, **cosmo)
616 return scipy.interpolate.interp1d(ages, z)
617
619 """The light travel distance to redshift z.
620
621 Units are Mpc.
622
623 Examples
624 --------
625
626 >>> import cosmolopy.distance as cd
627 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
628 >>> cosmo = cd.set_omega_k_0(cosmo)
629 >>> dlookback = cd.light_travel_distance(3.0, 2.0, **cosmo)
630 >>> print "Lookback distance from z=2 to 3 is %.2g Mpc" % (dlookback)
631 Lookback distance from z=2 to 3 is 3.3e+02 Mpc
632
633 """
634 t_look = lookback_time(z, z0, **cosmo)
635 return cc.c_light_Mpc_s * t_look
636
638 """The redshift corresponding to a given light travel distance.
639
640 Units are the same as light_travel_distance (Mpc).
641
642 Examples
643 --------
644
645 >>> import cosmolopy.distance as cd
646 >>> import cosmolopy.constants as cc
647 >>> cosmo = {'omega_M_0' : 0.3, 'omega_lambda_0' : 0.7, 'h' : 0.72}
648 >>> cosmo = cd.set_omega_k_0(cosmo)
649 >>> z = cd.redshift_d_light(10. * cc.c_light_Mpc_Gyr, **cosmo)
650 Optimization terminated successfully.
651 Current function value: 0.000112
652 Iterations: 26
653 Function evaluations: 52
654 >>> print "Redshift for a lookback time of 10Gyr is z=%.3f" % (z)
655 Redshift for a lookback time of 10Gyr is z=2.025
656
657 """
658
659 dl_diff = lambda z: abs(dl - light_travel_distance(z, **cosmo)[0])
660 z = scipy.optimize.fmin(dl_diff, z_guess, **fmin_args)
661 return z
662
663 if __name__ == "__main__":
664 import doctest
665 doctest.testmod()
666