1
2 """
3 Calculate the distance to a star with Lutz-Kelker bias
4
5 Section 1. Bias introduced by the shape of the Galaxy
6 =====================================================
7
8 In this section, we compare the distance calculation via the inverse of the
9 parallax with a method that takes into account the distribution of stars in the
10 Milky Way. The latter can be important when the error on the parallax is larger
11 than 25%: then, it becomes more likely to have observed a star located further
12 in the galaxy but which appears erronously closer, just because there are many
13 more stars in the volume Delta_d (as the volume is much larger). This depends
14 on the the location of the star in the Milky Way (in the Halo, the number of
15 stars is low anyway) and thus also on the shape of the Milky Way. Here, the
16 Galaxy is approximated as the sum of a self-gravitating isothermal disk and a
17 Gaussian Halo (L{rho}). For more details, see, e.g., Maiz-Apellaniz, Alfaro and
18 Sota 2007/2008 (Poster).
19
20 To correctly compute the distance to a star given it's parallax, we need to
21 know it's position wrt the galactic plane. So, we need information on the
22 coordinates of the star:
23
24 >>> from ivs.catalogs import sesame
25
26 We take three examples: Vega, which has a well known parallax, HIP14, which has a
27 mediocre value, and HIP15, which has a badly determined value. When retrieving
28 the information from SIMBAD, make sure we have Van Leeuwen's (2007) parallax,
29 and that we have the galactic coordinates of the target.
30
31 >>> vega = sesame.search('vega',fix=True)
32 >>> hip14= sesame.search('hip14',fix=True)
33 >>> hip15= sesame.search('hip15',fix=True)
34
35 We have the following parameters::
36
37 vega : 130.23 +/- 0.36 (0.28%)
38 HIP14: 4.86 +/- 0.67 (13.79%)
39 HIP15: 1.91 +/- 1.14 (59.69%)
40
41 We can compute the distance probability functions via L{distprob}, but first
42 we have to define arrays of distances to compute the probability on (in parsec):
43
44 >>> r1 = np.linspace(7,8,1000) # vega is closeby
45 >>> r2 = np.linspace(0,500,1000) # this one is around 200 pc
46 >>> r3 = np.linspace(0,1.5e4,1000) # this one seems to have a heavy tail
47
48 We only need the galactic latitude and the parallax:
49
50 >>> prob1 = distprob(r1,vega['galpos'][1],(vega['plx']['v'],vega['plx']['e']))
51 >>> prob2 = distprob(r2,hip14['galpos'][1],(hip14['plx']['v'],hip14['plx']['e']))
52 >>> prob3 = distprob(r3,hip15['galpos'][1],(hip15['plx']['v'],hip15['plx']['e']))
53
54 It is useful to compare with the Gaussian approximation. First we need some
55 extra modules to handle computations with uncertainties, and to evaluate the
56 Gaussian function.
57
58 >>> from ivs.sigproc import evaluate
59 >>> from ivs.units.uncertainties import ufloat
60
61 Then we can approximate the distance to the star (in pc) as the inverse of the
62 parallax (in arcsec).
63
64 >>> vega_ = 1./ufloat((vega['plx']['v']/1000.,vega['plx']['e']/1000.))
65 >>> hip14_ = 1./ufloat((hip14['plx']['v']/1000.,hip14['plx']['e']/1000.))
66 >>> hip15_ = 1./ufloat((hip15['plx']['v']/1000.,hip15['plx']['e']/1000.))
67 >>> prob1_ = evaluate.gauss(r1,[1.,vega_.nominal_value,vega_.std_dev()])
68 >>> prob2_ = evaluate.gauss(r2,[1.,hip14_.nominal_value,hip14_.std_dev()])
69 >>> prob3_ = evaluate.gauss(r3,[1.,hip15_.nominal_value,hip15_.std_dev()])
70
71 We summarize everything in a plot: up to a relative error of ~25%, the Gaussian
72 approximation works pretty well. From then on, even the peak of the distribution
73 shows a bias (the Lutz-Kelker bias), but also heavy tails can appear.
74
75 >>> p = pl.figure()
76 >>> p = pl.subplot(221)
77 >>> p = pl.plot(r1,prob1/prob1.max(),'k-',lw=2,label='Vega (LK)')
78 >>> p = pl.plot(r1,prob1_,'r--',lw=2,label='Vega (GA)')
79 >>> p = pl.legend()
80 >>> p = pl.xlabel('Distance (pc)');p = pl.ylabel('Unnormalised probability')
81 >>> p = pl.xlim(7.55,7.8)
82
83
84 >>> p = pl.subplot(222)
85 >>> p = pl.plot(r2,prob2/prob2.max(),'k-',lw=2,label='HIP14 (LK)')
86 >>> p = pl.plot(r2,prob2_,'r--',lw=2,label='HIP14 (GA)')
87 >>> p = pl.xlabel('Distance (pc)');p = pl.ylabel('Unnormalised probability')
88 >>> p = pl.legend()
89
90
91 >>> p = pl.subplot(212)
92 >>> p = pl.plot(r3,prob3/prob3.max(),'k-',lw=2,label='HIP15 (LK)')
93 >>> p = pl.plot(r3,prob3_,'r--',lw=2,label='HIP15 (GA)')
94 >>> p = pl.xlabel('Distance (pc)');p = pl.ylabel('Unnormalised probability')
95 >>> p = pl.legend()
96
97 ]include figure]]ivs_observations_distance_LK.png]
98
99 """
100 import numpy as np
101 from numpy import exp,sqrt,sin,pi
102 import logging
103
104 logger = logging.getLogger("IVS.DIST")
105
106
107 -def rho(z,z_sun=20.0,hd=31.8,hh=490.,sigma=1.91e-3,f=0.039):
108 """
109 Stellar galactic density function.
110
111 Galactic coordinate z: self-gravitating isothermal disk plus a Gaussian halo
112
113 See, e.g., Maiz-Apellaniz, Alfaro and Sota 2007/2008 (Poster)
114
115 Other values we found in the literature:
116
117 z_sun,hd,sigma,f = 24.7,34.2,1.62e-3,0.058
118
119 @param z: galactic coordinate z (parsec)
120 @type z: array or float
121 @param z_sun: Sun's distance above the Galactic plane (20.0 +/- 2.9 pc)
122 @type z_sun: float
123 @param hd: disk scale height (31.8+/-1.6 pc)
124 @type hd: float
125 @param hh: halo half width (490+/-170 pc)
126 @type hh: float
127 @param sigma: number of stars per square parsec (total surface number density)
128 (1.91e-3+/-0.11e-3)
129 @type sigma: float
130 @param f: fraction of stars in the halo population (0.039+/-0.015)
131 @type f: float
132 @return: halo+disk,halo and disk density function
133 @rtype: 3-tuple
134 """
135 disk = sigma* (1-f)/(4*hd) * np.cosh(0.5*(z+z_sun)/hd)**-2
136 halo = sigma*f/(2*sqrt(2)*hh) * exp( -(0.5*(z+z_sun)/hh)**2)
137 return halo+disk,halo,disk
138
140 """
141 Compute the probability for an object to be at distance r (pc), given
142 its parallax (mas) and error on the parallax (mas) and a constant density
143 function.
144
145 Unnormalised!
146
147 To obtain the probabilty, multiply with a stellar galactic density function.
148
149 @param r: distance (pc)
150 @type r: float/array
151 @param plx: parallax (mas) and error
152 @type plx: tuple of float/array
153 @return: probability function
154 @rtype: float/array
155 """
156 plx,e_plx = plx[0]/1000,plx[1]/1000.
157 A = 1.
158 prob1 = A*r**2*exp(-0.5*((1-r*plx)/(r*e_plx))**2)
159 return prob1
160
162 """
163 Compute the probability for an object to be located at a distance r (pc),
164 given its parallax and galactic lattitude.
165
166 theta in degrees
167 plx is a tuple!
168
169 returns (unnormalised) probability density function.
170
171 @param r: distance (pc)
172 @type r: float/array
173 @param theta: galactic latitude (deg)
174 @type theta: float
175 @param plx: parallax (mas) and error
176 @type plx: tuple of float/array
177 @return: probability function
178 @rtype: float/array
179 """
180 z = r*sin(theta/180.*pi)
181 pcd = probability_cd(r,plx)
182 galaxy,halo,disk = rho(z,**kwargs)
183 dpr = pcd*galaxy
184 logger.info("Calculate Lutz-Kelker bias for THETA=%.3fdeg"%(theta))
185 return dpr
186
187
188 if __name__ == "__main__":
189 import doctest
190 import pylab as pl
191 doctest.testmod()
192 pl.show()
193