Package astLib :: Module astCalc
[hide private]
[frames] | no frames]

Source Code for Module astLib.astCalc

  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  #------------------------------------------------------------------------------ 
44 -def dl(z):
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 #------------------------------------------------------------------------------
60 -def da(z):
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 #------------------------------------------------------------------------------
75 -def dm(z):
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 # Assume 1000 intervals 89 N = 1000 90 91 # Step size through integration 92 xMax = 1.0 93 xMin = 1.0 / (1.0 + z) 94 xStep = float((xMax - xMin) / N) 95 96 # Running total of y as we integrate 97 yTotal = 0 98 99 # Simpson's rule integration 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 # Ends 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: # Odd 112 yTotal = yTotal + (4*yn) 113 else: # Even 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 #------------------------------------------------------------------------------
130 -def dc(z):
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 # Assume 1000 intervals 143 N = 1000 144 145 # Step size through integration 146 xMax = 1.0 147 xMin = 1.0/(1.0+z) 148 xStep = float((xMax-xMin)/ N) 149 150 # Running total of y as we integrate 151 yTotal = 0 152 153 # Simpson's rule integration 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 # Ends 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: # Odd 166 yTotal = yTotal + (4*yn) 167 else: # Even 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 #------------------------------------------------------------------------------
177 -def dVcdz(z):
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 #------------------------------------------------------------------------------
276 -def t0():
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 # Assume 1000 intervals 288 N = 1000 289 290 # Step size through integration 291 xMax = 1.0 292 xMin = 0.0 293 xStep = float((xMax-xMin)/N) 294 295 # Running total of y as we integrate 296 yTotal = 0 297 298 # Simpson's rule integration 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 # Ends 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: # Odd 311 yTotal = yTotal + (4*yn) 312 else: # Even 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 #------------------------------------------------------------------------------
322 -def tl(z):
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 # Assume 1000 intervals 335 N = 1000 336 337 # Step size through integration 338 xMax = 1.0 339 xMin = 1.0/(1.0+z) 340 xStep = float((xMax-xMin)/N) 341 342 # Running total of y as we integrate 343 yTotal = 0 344 345 # Simpson's rule integration 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 # Ends 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: # Odd 358 yTotal = yTotal +(4*yn) 359 else: # Even 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 #------------------------------------------------------------------------------
369 -def tz(z):
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 #------------------------------------------------------------------------------
385 -def tl2z(tlGyr):
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 #------------------------------------------------------------------------------
426 -def tz2z(tzGyr):
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 #------------------------------------------------------------------------------
442 -def absMag(appMag, distMpc):
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 #------------------------------------------------------------------------------
459 -def Ez(z):
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 #------------------------------------------------------------------------------
476 -def Ez2(z):
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 # This form of E(z) is more reliable at high redshift. It is basically the 488 # same for all redshifts below 10. But above that, the radiation term 489 # begins to d ominate. From Peebles 1993. 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 #------------------------------------------------------------------------------
499 -def OmegaMz(z):
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 #------------------------------------------------------------------------------
516 -def OmegaLz(z):
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 #------------------------------------------------------------------------------
530 -def OmegaRz(z):
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 #------------------------------------------------------------------------------
544 -def DeltaVz(z):
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 #------------------------------------------------------------------------------
571 -def RVirialXRayCluster(kT, z, betaT):
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 # The equation quoted in Arnaud, Aghanim & Neumann is for h50, so need 601 # to scale it 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