1 """module for performing common calculations
2
3 (c) 2007-2011 Matt Hilton
4
5 (c) 2013 Matt Hilton & Steven Boada
6
7 U{http://astlib.sourceforge.net}
8
9 The focus in this module is at present on calculations of distances in a given
10 cosmology. The parameters for the cosmological model are set using the
11 variables OMEGA_M0, OMEGA_L0, OMEGA_R0, H0 in the module namespace (see below for
12 details).
13
14 @var OMEGA_M0: The matter density parameter at z=0.
15 @type OMEGA_M0: float
16
17 @var OMEGA_L0: The dark energy density (in the form of a cosmological
18 constant) at z=0.
19 @type OMEGA_L0: float
20
21 @var OMEGA_R0: The radiation density at z=0 (note this is only used currently
22 in calculation of L{Ez}).
23 @type OMEGA_R0: float
24
25 @var H0: The Hubble parameter (in km/s/Mpc) at z=0.
26 @type H0: float
27
28 @var C_LIGHT: The speed of light in km/s.
29 @type C_LIGHT: float
30
31 """
32
33 OMEGA_M0 = 0.3
34 OMEGA_L0 = 0.7
35 OMEGA_R0 = 8.24E-5
36 H0 = 70.0
37
38 C_LIGHT = 3.0e5
39
40 import math
41 import sys
42
43
45 """Calculates the luminosity distance in Mpc at redshift z.
46
47 @type z: float
48 @param z: redshift
49 @rtype: float
50 @return: luminosity distance in Mpc
51
52 """
53
54 DM = dm(z)
55 DL = (1.0+z)*DM
56
57 return DL
58
59
61 """Calculates the angular diameter distance in Mpc at redshift z.
62
63 @type z: float
64 @param z: redshift
65 @rtype: float
66 @return: angular diameter distance in Mpc
67
68 """
69 DM = dm(z)
70 DA = DM/(1.0+z)
71
72 return DA
73
74
76 """Calculates the transverse comoving distance (proper motion distance) in
77 Mpc at redshift z.
78
79 @type z: float
80 @param z: redshift
81 @rtype: float
82 @return: transverse comoving distance (proper motion distance) in Mpc
83
84 """
85
86 OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0
87
88
89 N = 1000
90
91
92 xMax = 1.0
93 xMin = 1.0 / (1.0 + z)
94 xStep = float((xMax - xMin) / N)
95
96
97 yTotal = 0
98
99
100 for n in range(N):
101 x = xMin + xStep*n
102 yn = (1.0/math.sqrt(OMEGA_M0*x + OMEGA_L0*math.pow(x, 4) +
103 OMEGA_K*math.pow(x, 2)))
104
105
106 if n == 0 or n == N:
107 yTotal= yTotal + yn
108 else:
109 oddness = float(n/2.0);
110 fractPart = math.modf(oddness)[0]
111 if fractPart == 0.5:
112 yTotal = yTotal + (4*yn)
113 else:
114 yTotal = yTotal + (2*yn)
115
116 integralValue = (xMax-xMin)/(3*N) * yTotal
117
118 if OMEGA_K > 0.0:
119 DM = (C_LIGHT/H0 * math.pow(abs(OMEGA_K), -0.5) *
120 math.sinh(math.sqrt(abs(OMEGA_K)) * integralValue))
121 elif OMEGA_K == 0.0:
122 DM = C_LIGHT/H0 * integralValue
123 elif OMEGA_K < 0.0:
124 DM = (C_LIGHT/H0 * math.pow(abs(OMEGA_K), -0.5) *
125 math.sin(math.sqrt(abs(OMEGA_K)) * integralValue))
126
127 return DM
128
129
131 """Calculates the line of sight comoving distance in Mpc at redshift z.
132
133 @type z: float
134 @param z: redshift
135 @rtype: float
136 @return: transverse comoving distance (proper motion distance) in Mpc
137
138 """
139
140 OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0
141
142
143 N = 1000
144
145
146 xMax = 1.0
147 xMin = 1.0/(1.0+z)
148 xStep = float((xMax-xMin)/ N)
149
150
151 yTotal = 0
152
153
154 for n in range(N):
155 x = xMin + xStep*n
156 yn = (1.0 / math.sqrt(OMEGA_M0*x + OMEGA_L0*math.pow(x, 4) +
157 OMEGA_K * math.pow(x,2)))
158
159
160 if n == 0 or n == N:
161 yTotal = yTotal + yn
162 else:
163 oddness = float(n/2.0)
164 fractPart = math.modf(oddness)[0]
165 if fractPart == 0.5:
166 yTotal = yTotal + (4*yn)
167 else:
168 yTotal = yTotal + (2*yn)
169
170 integralValue = (xMax-xMin)/(3*N)*yTotal
171
172 DC= C_LIGHT/H0*integralValue
173
174 return DC
175
176
178 """Calculates the line of sight comoving volume element per steradian dV/dz
179 at redshift z.
180
181 @type z: float
182 @param z: redshift
183 @rtype: float
184 @return: comoving volume element per steradian
185
186 """
187
188 dH = C_LIGHT/H0
189 dVcdz=(dH*(da(z)**2)*((1+z)**2))/Ez(z)
190
191 return dVcdz
192
193
194 -def dl2z(distanceMpc):
195 """Calculates the redshift z corresponding to the luminosity distance given
196 in Mpc.
197
198 @type distanceMpc: float
199 @param distanceMpc: distance in Mpc
200 @rtype: float
201 @return: redshift
202
203 """
204
205 dTarget = distanceMpc
206
207 toleranceMpc = 0.1
208
209 zMin = 0.0
210 zMax = 10.0
211
212 diff = dl(zMax) - dTarget
213 while diff < 0:
214 zMax = zMax + 5.0
215 diff = dl(zMax) - dTarget
216
217 zTrial = zMin + (zMax-zMin)/2.0
218
219 dTrial = dl(zTrial)
220 diff = dTrial - dTarget
221 while abs(diff) > toleranceMpc:
222
223 if diff > 0:
224 zMax = zMax - (zMax-zMin)/2.0
225 else:
226 zMin = zMin + (zMax-zMin)/2.0
227
228 zTrial = zMin + (zMax-zMin)/2.0
229 dTrial = dl(zTrial)
230 diff = dTrial - dTarget
231
232 return zTrial
233
234
235 -def dc2z(distanceMpc):
236 """Calculates the redshift z corresponding to the comoving distance given
237 in Mpc.
238
239 @type distanceMpc: float
240 @param distanceMpc: distance in Mpc
241 @rtype: float
242 @return: redshift
243
244 """
245
246 dTarget = distanceMpc
247
248 toleranceMpc = 0.1
249
250 zMin = 0.0
251 zMax = 10.0
252
253 diff = dc(zMax) - dTarget
254 while diff < 0:
255 zMax = zMax + 5.0
256 diff = dc(zMax) - dTarget
257
258 zTrial = zMin + (zMax-zMin)/2.0
259
260 dTrial = dc(zTrial)
261 diff = dTrial - dTarget
262 while abs(diff) > toleranceMpc:
263
264 if diff > 0:
265 zMax = zMax - (zMax-zMin)/2.0
266 else:
267 zMin = zMin + (zMax-zMin)/2.0
268
269 zTrial = zMin + (zMax-zMin)/2.0
270 dTrial = dc(zTrial)
271 diff = dTrial - dTarget
272
273 return zTrial
274
275
277 """Calculates the age of the universe in Gyr at z=0 for the current set of
278 cosmological parameters.
279
280 @rtype: float
281 @return: age of the universe in Gyr at z=0
282
283 """
284
285 OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0
286
287
288 N = 1000
289
290
291 xMax = 1.0
292 xMin = 0.0
293 xStep = float((xMax-xMin)/N)
294
295
296 yTotal = 0
297
298
299 for n in range(1, N):
300 x = xMin + xStep*n
301 yn = (x/math.sqrt(OMEGA_M0*x + OMEGA_L0 * math.pow(x, 4) + OMEGA_K *
302 math.pow(x, 2)))
303
304
305 if n == 0 or n == N:
306 yTotal = yTotal + yn
307 else:
308 oddness = float(n/2.0);
309 fractPart = math.modf(oddness)[0]
310 if fractPart == 0.5:
311 yTotal = yTotal + (4*yn)
312 else:
313 yTotal = yTotal + (2*yn)
314
315 integralValue = (xMax-xMin)/(3*N) *yTotal
316
317 T0 = (1.0/H0*integralValue*3.08e19)/3.16e7/1e9
318
319 return T0
320
321
323 """ Calculates the lookback time in Gyr to redshift z for the current set
324 of cosmological parameters.
325
326 @type z: float
327 @param z: redshift
328 @rtype: float
329 @return: lookback time in Gyr to redshift z
330
331 """
332 OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0
333
334
335 N = 1000
336
337
338 xMax = 1.0
339 xMin = 1.0/(1.0+z)
340 xStep = float((xMax-xMin)/N)
341
342
343 yTotal = 0
344
345
346 for n in range(1, N):
347 x = xMin + xStep*n
348 yn = (x/math.sqrt(OMEGA_M0*x + OMEGA_L0 * math.pow(x, 4) + OMEGA_K *
349 math.pow(x, 2)))
350
351
352 if n == 0 or n == N:
353 yTotal = yTotal + yn
354 else:
355 oddness = float(n/2.0);
356 fractPart = math.modf(oddness)[0]
357 if fractPart == 0.5:
358 yTotal = yTotal +(4*yn)
359 else:
360 yTotal = yTotal + (2*yn)
361
362 integralValue = (xMax-xMin)/(3*N) *yTotal
363
364 T0 = (1.0/H0*integralValue*3.08e19)/3.16e7/1e9
365
366 return T0
367
368
370 """Calculates the age of the universe at redshift z for the current set of
371 cosmological parameters.
372
373 @type z: float
374 @param z: redshift
375 @rtype: float
376 @return: age of the universe in Gyr at redshift z
377
378 """
379
380 TZ = t0() - tl(z)
381
382 return TZ
383
384
386 """Calculates the redshift z corresponding to lookback time tlGyr given in
387 Gyr.
388
389 @type tlGyr: float
390 @param tlGyr: lookback time in Gyr
391 @rtype: float
392 @return: redshift
393
394 """
395
396 tTarget = tlGyr
397
398 toleranceGyr = 0.001
399
400 zMin = 0.0
401 zMax = 10.0
402
403 diff = tl(zMax) - tTarget
404 while diff < 0:
405 zMax = zMax + 5.0
406 diff = tl(zMax) - tTarget
407
408 zTrial = zMin + (zMax-zMin)/2.0
409
410 tTrial = tl(zTrial)
411 diff = tTrial - tTarget
412 while abs(diff) > toleranceGyr:
413
414 if diff > 0:
415 zMax = zMax - (zMax-zMin)/2.0
416 else:
417 zMin = zMin + (zMax-zMin)/2.0
418
419 zTrial = zMin + (zMax-zMin)/2.0
420 tTrial = tl(zTrial)
421 diff = tTrial - tTarget
422
423 return zTrial
424
425
427 """Calculates the redshift z corresponding to age of the universe tzGyr
428 given in Gyr.
429
430 @type tzGyr: float
431 @param tzGyr: age of the universe in Gyr
432 @rtype: float
433 @return: redshift
434
435 """
436 tl = t0() - tzGyr
437 z = tl2z(tl)
438
439 return z
440
441
443 """Calculates the absolute magnitude of an object at given luminosity
444 distance in Mpc.
445
446 @type appMag: float
447 @param appMag: apparent magnitude of object
448 @type distMpc: float
449 @param distMpc: distance to object in Mpc
450 @rtype: float
451 @return: absolute magnitude of object
452
453 """
454 absMag = appMag - (5.0*math.log10(distMpc*1.0e5))
455
456 return absMag
457
458
460 """Calculates the value of E(z), which describes evolution of the Hubble
461 parameter with redshift, at redshift z for the current set of cosmological
462 parameters. See, e.g., Bryan & Norman 1998 (ApJ, 495, 80).
463
464 @type z: float
465 @param z: redshift
466 @rtype: float
467 @return: value of E(z) at redshift z
468
469 """
470
471 Ez = math.sqrt(Ez2(z))
472
473 return Ez
474
475
477 """Calculates the value of E(z)^2, which describes evolution of the Hubble
478 parameter with redshift, at redshift z for the current set of cosmological
479 parameters. See, e.g., Bryan & Norman 1998 (ApJ, 495, 80).
480
481 @type z: float
482 @param z: redshift
483 @rtype: float
484 @return: value of E(z)^2 at redshift z
485
486 """
487
488
489
490
491 Ez2 = (OMEGA_R0 * math.pow(1.0+z, 4) +
492 OMEGA_M0* math.pow(1.0+z, 3) +
493 (1.0- OMEGA_M0- OMEGA_L0) *
494 math.pow(1.0+z, 2) + OMEGA_L0)
495
496 return Ez2
497
498
500 """Calculates the matter density of the universe at redshift z. See, e.g.,
501 Bryan & Norman 1998 (ApJ, 495, 80).
502
503 @type z: float
504 @param z: redshift
505 @rtype: float
506 @return: matter density of universe at redshift z
507
508 """
509 ez2 = Ez2(z)
510
511 Omega_Mz = (OMEGA_M0*math.pow(1.0+z, 3))/ez2
512
513 return Omega_Mz
514
515
517 """ Calculates the dark energy density of the universe at redshift z.
518
519 @type z: float
520 @param z: redshift
521 @rtype: float
522 @return: dark energy density of universe at redshift z
523
524 """
525 ez2 = Ez2(z)
526
527 return OMEGA_L0/ez2
528
529
531 """ Calculates the radiation density of the universe at redshift z.
532
533 @type z: float
534 @param z: redshift
535 @rtype: float
536 @return: radiation density of universe at redshift z
537
538 """
539 ez2 = Ez2(z)
540
541 return OMEGA_R0*math.pow(1+z, 4)/ez2
542
543
545 """Calculates the density contrast of a virialised region S{Delta}V(z),
546 assuming a S{Lambda}CDM-type flat cosmology. See, e.g., Bryan & Norman
547 1998 (ApJ, 495, 80).
548
549 @type z: float
550 @param z: redshift
551 @rtype: float
552 @return: density contrast of a virialised region at redshift z
553
554 @note: If OMEGA_M0+OMEGA_L0 is not equal to 1, this routine exits and
555 prints an error
556 message to the console.
557
558 """
559
560 OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0
561
562 if OMEGA_K == 0.0:
563 Omega_Mz = OmegaMz(z)
564 deltaVz = (18.0*math.pow(math.pi, 2)+82.0*(Omega_Mz-1.0)-39.0 *
565 math.pow(Omega_Mz-1, 2))
566 return deltaVz
567 else:
568 raise Exception("cosmology is NOT flat.")
569
570
572 """Calculates the virial radius (in Mpc) of a galaxy cluster at redshift z
573 with X-ray temperature kT, assuming self-similar evolution and a flat
574 cosmology. See Arnaud et al. 2002 (A&A, 389, 1) and Bryan & Norman 1998
575 (ApJ, 495, 80). A flat S{Lambda}CDM-type flat cosmology is assumed.
576
577 @type kT: float
578 @param kT: cluster X-ray temperature in keV
579 @type z: float
580 @param z: redshift
581 @type betaT: float
582 @param betaT: the normalisation of the virial relation, for which Evrard et
583 al. 1996 (ApJ,469, 494) find a value of 1.05
584 @rtype: float
585 @return: virial radius of cluster in Mpc
586
587 @note: If OMEGA_M0+OMEGA_L0 is not equal to 1, this routine exits and
588 prints an error message to the console.
589
590 """
591
592 OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0
593
594 if OMEGA_K == 0.0:
595 Omega_Mz = OmegaMz(z)
596 deltaVz = (18.0 * math.pow(math.pi, 2) + 82.0 * (Omega_Mz-1.0)- 39.0 *
597 math.pow(Omega_Mz-1, 2))
598 deltaz = (deltaVz*OMEGA_M0)/(18.0*math.pow(math.pi, 2)*Omega_Mz)
599
600
601
602 h50 = H0/50.0
603 Rv = (3.80*math.sqrt(betaT)*math.pow(deltaz, -0.5) *
604 math.pow(1.0+z, (-3.0/2.0)) * math.sqrt(kT/10.0)*(1.0/h50))
605
606 return Rv
607
608 else:
609 raise Exception("cosmology is NOT flat.")
610
611
612