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