1 """module for performing common calculations
2
3 (c) 2007-2011 Matt Hilton
4
5 (c) 2013-2014 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
12 for 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 import math
34 try:
35 from scipy import integrate
36 except ImportError:
37 print("WARNING: astCalc failed to import scipy modules - ", )
38 print("some functions will not work")
39
40 OMEGA_M0 = 0.3
41 OMEGA_L0 = 0.7
42 OMEGA_R0 = 8.24E-5
43 H0 = 70.0
44
45 C_LIGHT = 3.0e5
46
47
49 """Calculates the luminosity distance in Mpc at redshift z.
50
51 @type z: float
52 @param z: redshift
53 @rtype: float
54 @return: luminosity distance in Mpc
55
56 """
57
58 DM = dm(z)
59 DL = (1.0 + z) * DM
60
61 return DL
62
63
64
66 """Calculates the angular diameter distance in Mpc at redshift z.
67
68 @type z: float
69 @param z: redshift
70 @rtype: float
71 @return: angular diameter distance in Mpc
72
73 """
74 DM = dm(z)
75 DA = DM / (1.0 + z)
76
77 return DA
78
79
80
82 """Calculates the transverse comoving distance (proper motion distance) in
83 Mpc at redshift z.
84
85 @type z: float
86 @param z: redshift
87 @rtype: float
88 @return: transverse comoving distance (proper motion distance) in Mpc
89
90 """
91
92 def _yn(x):
93
94 return (1.0 / math.sqrt(OMEGA_M0 * x + OMEGA_L0 * math.pow(x, 4) +
95 OMEGA_K * math.pow(x, 2)))
96
97 OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0
98
99
100 xMax = 1.0
101 xMin = 1.0 / (1.0 + z)
102
103 integralValue, integralError = integrate.quad(_yn, xMin, xMax)
104
105 if OMEGA_K > 0.0:
106 DM = (C_LIGHT / H0 * math.pow(
107 abs(OMEGA_K), -0.5) *
108 math.sinh(math.sqrt(abs(OMEGA_K)) * integralValue))
109 elif OMEGA_K == 0.0:
110 DM = C_LIGHT / H0 * integralValue
111 elif OMEGA_K < 0.0:
112 DM = (C_LIGHT / H0 * math.pow(
113 abs(OMEGA_K), -0.5) *
114 math.sin(math.sqrt(abs(OMEGA_K)) * integralValue))
115
116 return DM
117
118
119
121 """Calculates the line of sight comoving distance in Mpc at redshift z.
122
123 @type z: float
124 @param z: redshift
125 @rtype: float
126 @return: transverse comoving distance (proper motion distance) in Mpc
127
128 """
129 def _yn(x):
130
131 return (1.0 / math.sqrt(OMEGA_M0 * x + OMEGA_L0 * math.pow(x, 4) +
132 OMEGA_K * math.pow(x, 2)))
133
134 OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0
135
136
137 xMax = 1.0
138 xMin = 1.0 / (1.0 + z)
139
140 integralValue, integralError = integrate.quad(_yn, xMin, xMax)
141
142 DC = C_LIGHT / H0 * integralValue
143
144 return DC
145
146
147
149 """Calculates the line of sight comoving volume element per steradian dV/dz
150 at redshift z. The unit is Mpc^3/sr. To find the total volume of the
151 universe integrate this function between two redshifts and multiply by the
152 total area of the sky, 4pi steradians.
153
154 @type z: float
155 @param z: redshift
156 @rtype: float
157 @return: comoving volume element per steradian
158
159 """
160
161 dH = C_LIGHT / H0
162 dVcdz = (dH * (math.pow(da(z), 2)) * (math.pow(1 + z, 2)) / Ez(z))
163
164 return dVcdz
165
166
167
168 -def dl2z(distanceMpc):
169 """Calculates the redshift z corresponding to the luminosity distance given
170 in Mpc.
171
172 @type distanceMpc: float
173 @param distanceMpc: distance in Mpc
174 @rtype: float
175 @return: redshift
176
177 """
178
179 dTarget = distanceMpc
180
181 toleranceMpc = 0.1
182
183 zMin = 0.0
184 zMax = 10.0
185
186 diff = dl(zMax) - dTarget
187 while diff < 0:
188 zMax = zMax + 5.0
189 diff = dl(zMax) - dTarget
190
191 zTrial = zMin + (zMax - zMin) / 2.0
192
193 dTrial = dl(zTrial)
194 diff = dTrial - dTarget
195 while abs(diff) > toleranceMpc:
196
197 if diff > 0:
198 zMax = zMax - (zMax - zMin) / 2.0
199 else:
200 zMin = zMin + (zMax - zMin) / 2.0
201
202 zTrial = zMin + (zMax - zMin) / 2.0
203 dTrial = dl(zTrial)
204 diff = dTrial - dTarget
205
206 return zTrial
207
208
209
210 -def dc2z(distanceMpc):
211 """Calculates the redshift z corresponding to the comoving distance given
212 in Mpc.
213
214 @type distanceMpc: float
215 @param distanceMpc: distance in Mpc
216 @rtype: float
217 @return: redshift
218
219 """
220
221 dTarget = distanceMpc
222
223 toleranceMpc = 0.1
224
225 zMin = 0.0
226 zMax = 10.0
227
228 diff = dc(zMax) - dTarget
229 while diff < 0:
230 zMax = zMax + 5.0
231 diff = dc(zMax) - dTarget
232
233 zTrial = zMin + (zMax - zMin) / 2.0
234
235 dTrial = dc(zTrial)
236 diff = dTrial - dTarget
237 while abs(diff) > toleranceMpc:
238
239 if diff > 0:
240 zMax = zMax - (zMax - zMin) / 2.0
241 else:
242 zMin = zMin + (zMax - zMin) / 2.0
243
244 zTrial = zMin + (zMax - zMin) / 2.0
245 dTrial = dc(zTrial)
246 diff = dTrial - dTarget
247
248 return zTrial
249
250
251
253 """Calculates the age of the universe in Gyr at z=0 for the current set of
254 cosmological parameters.
255
256 @rtype: float
257 @return: age of the universe in Gyr at z=0
258
259 """
260
261 def _yn(x):
262
263 return (1.0 / math.sqrt(OMEGA_M0 * x + OMEGA_L0 * math.pow(x, 4) +
264 OMEGA_K * math.pow(x, 2)))
265
266 OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0
267
268
269 xMax = 1.0
270 xMin = 0
271
272 integralValue, integralError = integrate.quad(_yn, xMin, xMax)
273
274 T0 = (1.0 / H0 * integralValue * 3.08e19) / 3.16e7 / 1e9
275
276 return T0
277
278
279
281 """ Calculates the lookback time in Gyr to redshift z for the current set
282 of cosmological parameters.
283
284 @type z: float
285 @param z: redshift
286 @rtype: float
287 @return: lookback time in Gyr to redshift z
288
289 """
290
291 def _yn(x):
292
293 return (1.0 / math.sqrt(OMEGA_M0 * x + OMEGA_L0 * math.pow(x, 4) +
294 OMEGA_K * math.pow(x, 2)))
295
296 OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0
297
298
299 xMax = 1.0
300 xMin = 1. / (1. + z)
301
302 integralValue, integralError = integrate.quad(_yn, xMin, xMax)
303
304 T0 = (1.0 / H0 * integralValue * 3.08e19) / 3.16e7 / 1e9
305
306 return T0
307
308
309
311 """Calculates the age of the universe at redshift z for the current set of
312 cosmological parameters.
313
314 @type z: float
315 @param z: redshift
316 @rtype: float
317 @return: age of the universe in Gyr at redshift z
318
319 """
320
321 TZ = t0() - tl(z)
322
323 return TZ
324
325
326
328 """Calculates the redshift z corresponding to lookback time tlGyr given in
329 Gyr.
330
331 @type tlGyr: float
332 @param tlGyr: lookback time in Gyr
333 @rtype: float
334 @return: redshift
335
336 @note: Raises ValueError if tlGyr is not positive.
337
338 """
339 if tlGyr < 0.:
340 raise ValueError('Lookback time must be positive')
341
342 tTarget = tlGyr
343
344 toleranceGyr = 0.001
345
346 zMin = 0.0
347 zMax = 10.0
348
349 diff = tl(zMax) - tTarget
350 while diff < 0:
351 zMax = zMax + 5.0
352 diff = tl(zMax) - tTarget
353
354 zTrial = zMin + (zMax - zMin) / 2.0
355
356 tTrial = tl(zTrial)
357 diff = tTrial - tTarget
358 while abs(diff) > toleranceGyr:
359
360 if diff > 0:
361 zMax = zMax - (zMax - zMin) / 2.0
362 else:
363 zMin = zMin + (zMax - zMin) / 2.0
364
365 zTrial = zMin + (zMax - zMin) / 2.0
366 tTrial = tl(zTrial)
367 diff = tTrial - tTarget
368
369 return zTrial
370
371
372
374 """Calculates the redshift z corresponding to age of the universe tzGyr
375 given in Gyr.
376
377 @type tzGyr: float
378 @param tzGyr: age of the universe in Gyr
379 @rtype: float
380 @return: redshift
381
382 @note: Raises ValueError if Universe age not positive
383
384 """
385 if tzGyr <= 0:
386 raise ValueError('Universe age must be positive.')
387 tl = t0() - tzGyr
388 z = tl2z(tl)
389
390 return z
391
392
393
395 """Calculates the absolute magnitude of an object at given luminosity
396 distance in Mpc.
397
398 @type appMag: float
399 @param appMag: apparent magnitude of object
400 @type distMpc: float
401 @param distMpc: distance to object in Mpc
402 @rtype: float
403 @return: absolute magnitude of object
404
405 """
406 absMag = appMag - (5.0 * math.log10(distMpc * 1.0e5))
407
408 return absMag
409
410
411
413 """Calculates the value of E(z), which describes evolution of the Hubble
414 parameter with redshift, at redshift z for the current set of cosmological
415 parameters. See, e.g., Bryan & Norman 1998 (ApJ, 495, 80).
416
417 @type z: float
418 @param z: redshift
419 @rtype: float
420 @return: value of E(z) at redshift z
421
422 """
423
424 Ez = math.sqrt(Ez2(z))
425
426 return Ez
427
428
429
431 """Calculates the value of E(z)^2, which describes evolution of the Hubble
432 parameter with redshift, at redshift z for the current set of cosmological
433 parameters. See, e.g., Bryan & Norman 1998 (ApJ, 495, 80).
434
435 @type z: float
436 @param z: redshift
437 @rtype: float
438 @return: value of E(z)^2 at redshift z
439
440 """
441
442
443
444
445 Ez2 = (OMEGA_R0 * math.pow(1.0 + z, 4) + OMEGA_M0 * math.pow(1.0 + z, 3) +
446 (1.0 - OMEGA_M0 - OMEGA_L0) * math.pow(1.0 + z, 2) + OMEGA_L0)
447
448 return Ez2
449
450
451
453 """Calculates the matter density of the universe at redshift z. See, e.g.,
454 Bryan & Norman 1998 (ApJ, 495, 80).
455
456 @type z: float
457 @param z: redshift
458 @rtype: float
459 @return: matter density of universe at redshift z
460
461 """
462 ez2 = Ez2(z)
463
464 Omega_Mz = (OMEGA_M0 * math.pow(1.0 + z, 3)) / ez2
465
466 return Omega_Mz
467
468
469
471 """ Calculates the dark energy density of the universe at redshift z.
472
473 @type z: float
474 @param z: redshift
475 @rtype: float
476 @return: dark energy density of universe at redshift z
477
478 """
479 ez2 = Ez2(z)
480
481 return OMEGA_L0 / ez2
482
483
484
486 """ Calculates the radiation density of the universe at redshift z.
487
488 @type z: float
489 @param z: redshift
490 @rtype: float
491 @return: radiation density of universe at redshift z
492
493 """
494 ez2 = Ez2(z)
495
496 return OMEGA_R0 * math.pow(1 + z, 4) / ez2
497
498
499
501 """Calculates the density contrast of a virialised region S{Delta}V(z),
502 assuming a S{Lambda}CDM-type flat cosmology. See, e.g., Bryan & Norman
503 1998 (ApJ, 495, 80).
504
505 @type z: float
506 @param z: redshift
507 @rtype: float
508 @return: density contrast of a virialised region at redshift z
509
510 @note: If OMEGA_M0+OMEGA_L0 is not equal to 1, this routine exits and
511 prints an error
512 message to the console.
513
514 """
515
516 OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0
517
518 if OMEGA_K == 0.0:
519 Omega_Mz = OmegaMz(z)
520 deltaVz = (18.0 * math.pow(math.pi, 2) + 82.0 *
521 (Omega_Mz - 1.0) - 39.0 * math.pow(Omega_Mz - 1, 2))
522 return deltaVz
523 else:
524 raise Exception("cosmology is NOT flat.")
525
526
527
529 """Calculates the virial radius (in Mpc) of a galaxy cluster at redshift z
530 with X-ray temperature kT, assuming self-similar evolution and a flat
531 cosmology. See Arnaud et al. 2002 (A&A, 389, 1) and Bryan & Norman 1998
532 (ApJ, 495, 80). A flat S{Lambda}CDM-type flat cosmology is assumed.
533
534 @type kT: float
535 @param kT: cluster X-ray temperature in keV
536 @type z: float
537 @param z: redshift
538 @type betaT: float
539 @param betaT: the normalisation of the virial relation, for which Evrard et
540 al. 1996 (ApJ,469, 494) find a value of 1.05
541 @rtype: float
542 @return: virial radius of cluster in Mpc
543
544 @note: If OMEGA_M0+OMEGA_L0 is not equal to 1, this routine exits and
545 prints an error message to the console.
546
547 """
548
549 OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0
550
551 if OMEGA_K == 0.0:
552 Omega_Mz = OmegaMz(z)
553 deltaVz = (18.0 * math.pow(math.pi, 2) + 82.0 *
554 (Omega_Mz - 1.0) - 39.0 * math.pow(Omega_Mz - 1, 2))
555 deltaz = (deltaVz * OMEGA_M0) / (18.0 * math.pow(math.pi, 2) *
556 Omega_Mz)
557
558
559
560 h50 = H0 / 50.0
561 Rv = (3.80 * math.sqrt(betaT) * math.pow(deltaz, -0.5) *
562 math.pow(1.0 + z,
563 (-3.0 / 2.0)) * math.sqrt(kT / 10.0) * (1.0 / h50))
564
565 return Rv
566
567 else:
568 raise Exception("cosmology is NOT flat.")
569
570
571