Module distance

source code

Calculate the distance to a star with Lutz-Kelker bias

Section 1. Bias introduced by the shape of the Galaxy

In this section, we compare the distance calculation via the inverse of the parallax with a method that takes into account the distribution of stars in the Milky Way. The latter can be important when the error on the parallax is larger than 25%: then, it becomes more likely to have observed a star located further in the galaxy but which appears erronously closer, just because there are many more stars in the volume Delta_d (as the volume is much larger). This depends on the the location of the star in the Milky Way (in the Halo, the number of stars is low anyway) and thus also on the shape of the Milky Way. Here, the Galaxy is approximated as the sum of a self-gravitating isothermal disk and a Gaussian Halo (rho). For more details, see, e.g., Maiz-Apellaniz, Alfaro and Sota 2007/2008 (Poster).

To correctly compute the distance to a star given it's parallax, we need to know it's position wrt the galactic plane. So, we need information on the coordinates of the star:

>>> from ivs.catalogs import sesame

We take three examples: Vega, which has a well known parallax, HIP14, which has a mediocre value, and HIP15, which has a badly determined value. When retrieving the information from SIMBAD, make sure we have Van Leeuwen's (2007) parallax, and that we have the galactic coordinates of the target.

>>> vega = sesame.search('vega',fix=True)
>>> hip14= sesame.search('hip14',fix=True)
>>> hip15= sesame.search('hip15',fix=True)

We have the following parameters:

   vega : 130.23 +/- 0.36 (0.28%)
   HIP14: 4.86 +/- 0.67 (13.79%)
   HIP15: 1.91 +/- 1.14 (59.69%)

We can compute the distance probability functions via distprob, but first we have to define arrays of distances to compute the probability on (in parsec):

>>> r1 = np.linspace(7,8,1000)     # vega is closeby
>>> r2 = np.linspace(0,500,1000)   # this one is around 200 pc
>>> r3 = np.linspace(0,1.5e4,1000) # this one seems to have a heavy tail

We only need the galactic latitude and the parallax:

>>> prob1 = distprob(r1,vega['galpos'][1],(vega['plx']['v'],vega['plx']['e']))
>>> prob2 = distprob(r2,hip14['galpos'][1],(hip14['plx']['v'],hip14['plx']['e']))
>>> prob3 = distprob(r3,hip15['galpos'][1],(hip15['plx']['v'],hip15['plx']['e']))

It is useful to compare with the Gaussian approximation. First we need some extra modules to handle computations with uncertainties, and to evaluate the Gaussian function.

>>> from ivs.sigproc import evaluate
>>> from ivs.units.uncertainties import ufloat

Then we can approximate the distance to the star (in pc) as the inverse of the parallax (in arcsec).

>>> vega_ = 1./ufloat((vega['plx']['v']/1000.,vega['plx']['e']/1000.))
>>> hip14_ = 1./ufloat((hip14['plx']['v']/1000.,hip14['plx']['e']/1000.))
>>> hip15_ = 1./ufloat((hip15['plx']['v']/1000.,hip15['plx']['e']/1000.))
>>> prob1_ = evaluate.gauss(r1,[1.,vega_.nominal_value,vega_.std_dev()])
>>> prob2_ = evaluate.gauss(r2,[1.,hip14_.nominal_value,hip14_.std_dev()])
>>> prob3_ = evaluate.gauss(r3,[1.,hip15_.nominal_value,hip15_.std_dev()])

We summarize everything in a plot: up to a relative error of ~25%, the Gaussian approximation works pretty well. From then on, even the peak of the distribution shows a bias (the Lutz-Kelker bias), but also heavy tails can appear.

>>> p = pl.figure()
>>> p = pl.subplot(221)
>>> p = pl.plot(r1,prob1/prob1.max(),'k-',lw=2,label='Vega (LK)')
>>> p = pl.plot(r1,prob1_,'r--',lw=2,label='Vega (GA)')
>>> p = pl.legend()
>>> p = pl.xlabel('Distance (pc)');p = pl.ylabel('Unnormalised probability')
>>> p = pl.xlim(7.55,7.8)

>>> p = pl.subplot(222)
>>> p = pl.plot(r2,prob2/prob2.max(),'k-',lw=2,label='HIP14 (LK)')
>>> p = pl.plot(r2,prob2_,'r--',lw=2,label='HIP14 (GA)')
>>> p = pl.xlabel('Distance (pc)');p = pl.ylabel('Unnormalised probability')
>>> p = pl.legend()

>>> p = pl.subplot(212)
>>> p = pl.plot(r3,prob3/prob3.max(),'k-',lw=2,label='HIP15 (LK)')
>>> p = pl.plot(r3,prob3_,'r--',lw=2,label='HIP15 (GA)')
>>> p = pl.xlabel('Distance (pc)');p = pl.ylabel('Unnormalised probability')
>>> p = pl.legend()

]include figure]]ivs_observations_distance_LK.png]