#### Get Python

First make sure that you have Python set up on your system. If you do not yet have a favorite Python environment, we suggest Anaconda.

#### Install Python packages

We will need Healpy for reading the 3D localization file. For this example, we will also use Astropy/Astroquery for retrieving archival data, Matplotlib for plotting, and Scipy for statistics functions. Open a termianl and run the following commands:

$ conda install astropy scipy matplotlib
$ pip install astroquery healpy

#### Try it Out: Compute Probability in the Direction of PGC 23997

Fire up a Python interpreter and import packages:

$ python
>>> import healpy as hp
>>> import numpy as np
>>> from matplotlib import pyplot as plt
>>> from scipy.stats import norm
>>> from astroquery.ned import Ned
>>> from astropy.utils.data import download_file
>>> from astropy.cosmology import default_cosmology
>>> cosmo = default_cosmology.get()

Next, download an example 3D sky map:

>>> url = 'http://asd.gsfc.nasa.gov/Leo.Singer/going-the-distance/18951_bayestar.fits.gz'
>>> filename = download_file(url, cache=True)

Read the four HEALPix layers with Healpy:

>>> prob, distmu, distsigma, distnorm = hp.read_map(filename, field=[0, 1, 2, 3])

Look up PGC 23997 in NED:

>>> pgc = Ned.query_object('PGC 23997')

Convert its RA and Dec to physicists' spherical polar coordinates (𝜃, 𝜙) with 𝜃=0 at the North celestial pole, and its redshift to a luminosity distance in Mpc:

>>> theta = 0.5 * np.pi - np.deg2rad(pgc[0]['DEC(deg)'])
>>> phi = np.deg2rad(pgc[0]['RA(deg)'])
>>> dist = cosmo.luminosity_distance(pgc[0]['Redshift']).value

Look up the HEALPix pixel that contains the galaxy:

>>> npix = len(prob)
>>> nside = hp.npix2nside(npix)
>>> ipix = hp.ang2pix(nside, theta, phi)

Finally, plot the conditional distance distribution in the of the galaxy:

>>> r = np.linspace(0, 120)
>>> p = norm(distmu[ipix], distsigma[ipix]).pdf(r) * distnorm[ipix] * r**2
>>> plt.plot(r, p)
>>> plt.axvline(dist)
>>> plt.xlabel('Distance (Mpc)')
>>> plt.ylabel('Probability / Mpc')
>>> plt.savefig('pgc23997.png')