Source Code for Module ivs.sed.distance

 1  # -*- coding: utf-8 -*- 
 2  """ 
 3  Calculate the distance to a star with Lutz-Kelker bias 
 4  """ 
 5  import numpy as np 
 6  from numpy import exp,sqrt,sin,pi 
 7  import logging 
 8   
 9  logger = logging.getLogger("IVS.DIST") 
10   
11   
12 -def rho(z,z_sun=20.0,hd=31.8,hh=490.,sigma=1.91e-3,f=0.039): 
13      """ 
14      Stellar galactic density function. 
15       
16      Galactic coordinate z: self-gravitating isothermal disk plus a Gaussian halo 
17       
18      See, e.g., Maiz-Apellaniz, Alfaro and Sota 2007/2008 (Poster) 
19       
20      Other values we found in the literature: 
21       
22      z_sun,hd,sigma,f = 24.7,34.2,1.62e-3,0.058 
23       
24      @param z: galactic coordinate z (parsec) 
25      @type z: array or float 
26      @param z_sun: Sun's distance above the Galactic plane (20.0 +/- 2.9 pc) 
27      @type z_sun: float 
28      @param hd: disk scale height (31.8+/-1.6 pc) 
29      @type hd: float 
30      @param hh: halo half width (490+/-170 pc) 
31      @type hh: float 
32      @param sigma: number of stars per square parsec (total surface number density) 
33      (1.91e-3+/-0.11e-3) 
34      @type sigma: float 
35      @param f: fraction of stars in the halo population (0.039+/-0.015) 
36      @type f: float 
37      @return: halo+disk,halo and disk density function 
38      @rtype: 3-tuple 
39      """ 
40      disk = sigma* (1-f)/(4*hd) * np.cosh(0.5*(z+z_sun)/hd)**-2 
41      halo = sigma*f/(2*sqrt(2)*hh) * exp( -(0.5*(z+z_sun)/hh)**2) 
42      return halo+disk,halo,disk 
43   
44 -def probability_cd(r,plx,e_plx): 
45      """ 
46      Compute the probability for an object to be at distance r (pc), given 
47      its parallax (mas) and error on the parallax (mas) and a constant density 
48      function. 
49       
50      Unnormalised! 
51       
52      To obtain the probabilty, multiply with a stellar galactic density function. 
53       
54      @param r: distance (pc) 
55      @type r: float/array 
56      @param plx: parallax (mas) 
57      @type plx: float/array 
58      @param e_plx: error on parallax (mas) 
59      @type e_plx: float/array 
60      @return: probability function 
61      """ 
62      plx /= 1000. 
63      e_plx /= 1000. 
64      A = 1. 
65      prob1 = A*r**2*exp(-0.5*((1-r*plx)/(r*e_plx))**2) 
66      return prob1 
67   
68 -def distprob(r,theta,plx,**kwargs): 
69      """ 
70      Compute the probability for an object to be located at a distance r (pc), 
71      given its parallax and galactic lattitude. 
72       
73      theta in degrees 
74      plx is a tuple! 
75       
76      returns (unnormalised) probability density function. 
77      """ 
78      z = r*sin(theta/180.*pi) 
79      pcd = probability_cd(r,plx[0],plx[1]) 
80      galaxy,halo,disk = rho(z,**kwargs) 
81      dpr = pcd*galaxy 
82      logger.info("Calculate Lutz-Kelker bias for THETA=%.3fdeg"%(theta)) 
83      return dpr 
84