Python Tutorials I

In this section we will review some of the concepts seen in the introduction, such as units, coordinate systems, etc. and work some tutorials in Python

This is not intended to be a complete course on programing in python .The idea is mor to put python in the context of astro-article physics, and address some of the questions seen in the following tutorials.

Working with units

In this first case or tutorial we are going to show a n example of how to work with units using the module units that exist for example in astropy

Code

import astropy.units as u
from astropy import constants as const

M_Earth = 5.97E24 * u.kg
M_Sun = 1.99E30 * u.kg
M_MW  = 1E12 * M_Sun
#By adding the quantity u.kg you can print directly the mass in Kg
print ("Mass Earth is: %s" %M_Earth)

Output

There are many other modules in python that have the units utility

Note that the variables defined above already have their units attached to them, so when you make a print statement it will provide as well the units as a string. That's the reason of the %s format. More examples:

R_Earth = 6.371E6 * u.m # meters
R_Sun = 6.955E8 * u.m # meters
AU = 1.496E11 * u.m # meters

With the radius we can calculate the mean density of Earth and the Sun. We will show that units are preserved along calculations:

Code

from numpy import pi
vol_sphere = lambda r: 4*pi/3*r**3
rho_Sun = M_Sun / vol_sphere(R_Sun)
rho_Earth = M_Earth / vol_sphere(R_Earth)

#A unit can be changed calling the .to(u.unit) method
print ("Density of Earth: %s" %(rho_Earth.to(u.g/u.cm**3)))
print ("Density of Sun: %s" %rho_Sun)

Output

We can use this module to make different transformations of units, for example from light years to meters:

Code

ly = 1 * u.lyr

print ("Number of seconds for light to travel from Sun to Earth:", 1./const.c.to(u.AU/u.s))
print ("Meters in a light year: %s", ly.to(u.m))

Output

Assuming that the Galaxy is roughly a disk 50 kpc in diameter and 500 pc thick we can now calculate its density:

Code

V_Gal =  pi * (25000*u.pc)**2 * 500*u.pc

print "Volume of the Milky Way is approximately: %s " % V_Gal.to(u.m**3)
M_Gal = 1E12 * M_Sun
rho_Gal = M_Gal / V_Gal
print "Average density of Milky Way is %s" % rho_Gal.to(u.g/u.cm**3)

Output

Exercise: How long does the light travel from the Sun to the Earth in minutes? How long does the light travel from the Galactic center (assume a distance of 8 kpc) in years?

Coordinates

The coordinates package of astropy provides classes for representing a variety of celestial/spatial coordinates and their velocity components, as well as tools for converting between common coordinate systems in a uniform way. The simplest way to use coordinates is to use the SkyCoord class. SkyCoord objects are instantiated by passing in positions (and optional velocities) with specified units and a coordinate frame.

from astropy import units as u
from astropy.coordinates import SkyCoord
c = SkyCoord(ra=10.625*u.degree, dec=41.2*u.degree, frame='icrs')
print(c)
print("ra = {:.2f}, {:.2f} rad".format(c.ra,c.ra.radian))
print("dec = {:.2f}, {:.2f} rad".format(c.dec,c.dec.radian))

With SkyCoord you can handle more than one set of coordinates:

c = SkyCoord([1, 2, 3], [-30, 45, 8], frame="icrs", unit="deg")  # 3 coords
print(c)

You can use from_name method of SkyCoord to retrieve the coordinates of known sources. It uses Sesame (http://cds.u-strasbg.fr/cgi-bin/Sesame) to retrieve coordinates for a particular named object:

from astropy.coordinates import SkyCoord
SkyCoord.from_name("PSR J1012+5307")  

Exercise: Define the sky coordinate for your favorite object and find the distance to the crab nebula and Galactic center.

Working with Tables

Astropy.table (http://docs.astropy.org/en/stable/table/) provides functionality for storing and manipulating heterogeneous tables of data.

Let's create a simple table with three columns of data named a, b, and c. These columns have integer, float, and string values respectively:

from astropy.table import Table
a = [6, 4, 5]
b = [2.0, 5.0, 8.2]
c = ['x', 'y', 'z']
t = Table([a, b, c], names=('a', 'b', 'c'), meta={'name': 'first table'})
print(t)

To check which columns are available you can use Table.colnames:

t.colnames

And access individual columns just by their name:

t['a']

And also a subset of columns:

t[['a', 'b']]

Rows can be accessed using numpy indexing:

t[0:1]

Or by using a boolean numpy array for indexing:

selection = t['a'] < 6
t[selection]

Astropy tables can be serialized into many formats. Some examples below:

t.write('example.fits', overwrite=True, format='fits')
t.write('example.ecsv', overwrite=True, format='ascii.ecsv')
Table.read('example.fits')

To sort the table by column 'a':

t.sort('b')

🖋️ Exercise: The (Lon, Lat) positions of Crab, Sagittarius A* and Cassiopeia A in the Galactic Coordinate system are:

• (184.55754381,-5.78427369) Crab

• (359.94422947,-0.04615714) Sag A*

• (111.74169477,-2.13544151) Cas A

Put in a Table the positions (with the first two columns for Lon, Lat) of the three objects (rows), and add two columns with the RA and DEC coordinates of the objects

The Astroquery package

Another package that allows one to retrieve information from astronomy databases for further processing is astroquery. It can query in different databases (SIMBAD , SDSS). You can import this package to explore the SDSS database by inserting the following line at the top of your script:

from astroquery.sdss import SDSS 
from astroquery.simbad import Simbad

🖋️ Exercise: Using the astroquery.sdss imported above, call the function query_region to get the following output [’ra’, ’dec’, ’u’, ’g’, ’r’, ’i’, ’z’, ’type’] in 1 arc minutes around ’pos’. Call this output ’xid’. Hint: To specify the output, use the option photoobj_fields If you check the type of this output, you see that it is of type Table.

• Select only sources that are identified as galaxies and put them in a table called ’galaxy_table’

• Select stars and put them in ’a table called ’stars_table’

• How many objects of each type do we have? The objects types in SDSS can be found in this table: http://www.sdss3.org/dr8/algorithms/classify.php#photo_sb

Working with coordinates

Now let's do a little tutorial. We are going to work. Let's do some "representation" of the galactic plane. We generate some random points using numpy following a 2-dimensional gaussian in the $\ell$: $-\pi, +\pi$ and $b$:$-\pi/2, +\pi/2$ space. Now we are going to use matplotlib to make plots. If you are using jupyter notebook, don't forget to use the magic command %matplotlib inline if you want to make the plots appear inline inside the notebook:

Code

%matplotlib inline

#We call random.multivariate_normal 
#to generate randon normal points at 0 
import numpy as np

#import matplotlib for plots
import matplotlib.pylab as plt
#Lets use the inline figure format as svg
%config InlineBackend.figure_format = 'svg'

plt.style.use("seaborn-poster")

disk = np.random.multivariate_normal(mean=[0,0], 
cov=np.diag([1,0.001]), size=5000)
#disk is a list of pairs [l, b] in radians

f = plt.figure(figsize=(9,7))
ax = plt.subplot(111, projection='aitoff')
#There are several projections: Aitoff, Hammer, Lambert, 
#MollWeide
ax.set_title("Galactic\n Coordinates")
ax.grid(True)
ll = disk[:,0]
bb = disk[:,1]
ax.plot(ll, bb, 'o', markersize=2, alpha=0.3)
plt.show()

Output

Now let's plot it in equatorial coordinates (right ascension, declination). Because matplotlib needs the coordinates in radians and between −$\pi$ and $\pi$, not between 0 and $2\pi$, we have to convert them.

Code

from astropy import units as u
from astropy.coordinates import SkyCoord

c_gal = SkyCoord(l=ll*u.radian, b=bb*u.radian, frame='galactic')
c_gal_icrs = c_gal.icrs

ra_rad = c_gal_icrs.ra.wrap_at(180 * u.deg).radian
dec_rad = c_gal_icrs.dec.radian 

plt.figure(figsize=(9,7))
ax = plt.subplot(111, projection="mollweide")
ax.set_title("Equatorial coordinates")
plt.grid(True)
ax.plot(ra_rad, dec_rad, 'o', markersize=2, alpha=0.3)
plt.show()

Output

Usually right ascension is plotted from right (0°) to left (360°), but it can also be expressed in hours (1 hr = 15°).

Calculating the Age of the Universe

We are going to use python to numerically solve the look-back time for any given set of cosmological parameters. For that we are going to rewrite the loopback formula in terms of the scale factor $a$ since we have the redshift scale relation:

$$(1+z) = \frac{1}{a}$$

we can prove thatL

$$\frac{{\rm d}z}{1 +z} = -\frac{{\rm d}a}{a}$$

$${\rm d}t = \frac{{\rm d}a}{H(a)a}$$

with

$$H^2(a) = H^2_0 [\Omega^0_M a^{-3} + \Omega^0_r a^{-4} + \Omega^0_k a^{-2} +\Omega^0_\Lambda]$$

Code

First we are going to use the astropy.cosmology module to retrieve the values of the cosmological constants using Planck data from 2013.

from astropy.cosmology import Planck13
H0 = Planck13.H(0).value

#We define the friedman equation ignoring the radiation component omega_r <<
def friedman(a, omega_M, omega_lambda):
    omega_k = 1 - omega_M - omega_lambda
    H2 = H0**2 * (omega_M * a**-3 + omega_k * a**-2 + omega_lambda)
    return H2

import scipy.integrate as integrate

def calculate_times(omega_m, omega_lambda):
    #Lets check if there is a maximal, ie, if H^2(a) = 0
    astep = 0.001
    scales = np.arange(0.1, 3, astep)
    amax = 1e9
    for a in scales:
        if( friedman(a, omega_M, omega_lambda) < 0):
            amax = a
            scales = np.arange(0.1, 2*amax, 0.01) #we go from a = 0 to a= 2*amax
            print ("Found a-max in %.2f" %amax)
            break

    #time at a = a_max
    tmax = integrate.quad(lambda x: 1/x * 1/np.sqrt(friedman(x,omega_m, omega_lambda)), 0, amax - astep)[0]
    #time today a = 1
    t0 = integrate.quad(lambda x: 1/x * 1/np.sqrt(friedman(x,omega_m, omega_lambda)), 0, 1)[0]

    #Empty x,y for the plots
    times = []
    scale = []
    for a in scales:
        if a < amax: #If a < amax we do the typical integral 0 -> a
            time = integrate.quad(lambda x: 1/x * 1/np.sqrt(friedman(x,omega_m, omega_lambda)), 0, a)[0]
            times.append(H0*(time - t0)) #We calibrate at -t0 to place all curves together
            scfrom IPython.core.display import HTML
        if a >= amax: #If a > amax we are out the domain we integrate backwards
            time = 2*tmax - integrate.quad(lambda x: 1/x * 1/np.sqrt(friedman(x,omega_m, omega_lambda)), 0, 2*amax -a)[0]
            times.append(H0*(time - t0))
            scale.append(2*amax - a)
    return np.array(times), np.array(scale)

fig, ax = plt.subplots(1, 1, figsize=(8,5))

omega_M = 0.3
omega_lambda = 0.7
x, y = calculate_times(omega_M, omega_lambda)
ax.plot(x, y, label="$\Omega_M = %.1f, \Omega_\Lambda = %.1f$" %(omega_M, omega_lambda))

omega_M = 1.0
omega_lambda = 0
x, y = calculate_times(omega_M, omega_lambda)
ax.plot(x, y, label="$\Omega_M = %.1f, \Omega_\Lambda = %.1f$" %(omega_M, omega_lambda))

omega_M = 0.3
omega_lambda = 0
x, y = calculate_times(omega_M, omega_lambda)
ax.plot(x, y, label="$\Omega_M = %.1f, \Omega_\Lambda = %.1f$" %(omega_M, omega_lambda))

omega_M = 4
omega_lambda = 0
x, y = calculate_times(omega_M, omega_lambda)
ax.plot(x, y, label="$\Omega_M = %.1f, \Omega_\Lambda = %.1f$" %(omega_M, omega_lambda))


plt.legend(loc="upper left")
ax.grid()
ax.set_xlim(-1,1.5)
ax.set_ylim(0.2,2)
ax.set_ylabel("$a(t)$")
ax.set_xlabel("$H_0 (t - t_0)$")
plt.show()

Output