# Python Tutorials I

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`](http://astropy.readthedocs.org/en/latest/units/)

{% tabs %}
{% tab title="Code" %}

```python
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)
```

{% endtab %}

{% tab title="Output" %}

```
Mass Earth is: 5.97e+24 kg
```

{% endtab %}
{% endtabs %}

{% hint style="info" %}
There are many other modules in python that have the units utility
{% endhint %}

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:<br>

```python
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:

{% tabs %}
{% tab title="Code" %}

```python
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)
```

{% endtab %}

{% tab title="Output" %}
`Density of Earth: 5.51141236929 g / cm3 Density of Sun: 1412.12474755 kg / m3`
{% endtab %}
{% endtabs %}

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

{% tabs %}
{% tab title="Code" %}

```python
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))
```

{% endtab %}

{% tab title="Output" %}

```
Number of seconds for light to travel from Sun to Earth: 499.004783836 s / AU
Meters in a light year: 9.46073047258e+15 m
```

{% endtab %}
{% endtabs %}

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

{% tabs %}
{% tab title="Code" %}

```python
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)
```

{% endtab %}

{% tab title="Output" %}

```python
Volume of the Milky Way is approximately: 2.88437372034e+61 m3 
```

```
Average density of Milky Way is 6.89924466433e-23 g / cm3
```

{% endtab %}
{% endtabs %}

{% hint style="warning" %}
**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?
{% endhint %}

## Coordinates &#x20;

The [coordinates package ](http://docs.astropy.org/en/stable/coordinates/)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.

```python
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:

```python
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:

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

{% hint style="warning" %}
**Exercise:** Define the sky coordinate for your favorite object and find the distance to the crab nebula and Galactic center.
{% endhint %}

## 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:

```python
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`:

```python
t.colnames
```

And access individual columns just by their name:

```python
t['a']
```

And also a subset of columns:

```python
t[['a', 'b']]
```

Rows can be accessed using `numpy` indexing:

```python
t[0:1]
```

Or by using a boolean `numpy` array for indexing:

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

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

```python
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':

```python
t.sort('b')
```

{% hint style="info" %}
:pen\_fountain: **Exercise:** The (Lon, Lat) positions of Crab, Sagittarius A\* and Cassiopeia A in the Galactic Coordinate system are:

* (184.55754381,-5.78427369) Crab&#x20;
* (359.94422947,-0.04615714) Sag A\*&#x20;
* (111.74169477,-2.13544151) Cas A&#x20;

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
{% endhint %}

## 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](http://simbad.u-strasbg.fr/simbad/) , [SDSS](http://skyserver.sdss.org/dr14/en/tools/chart/navi.aspx)). You can import this package to explore the SDSS database by inserting the following line at the top of your script:

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

{% hint style="info" %}
:pen\_fountain: **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.

* &#x20;Select only sources that are identified as galaxies and put them in a table called ’galaxy\_table’&#x20;
* Select stars and put them in ’a table called ’stars\_table’&#x20;
* 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>
  {% endhint %}

## Working with coordinates&#x20;

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:

{% tabs %}
{% tab title="Code" %}

```python
%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()
```

{% endtab %}

{% tab title="Output" %}
![Toy representation of the galactic plane in galactic coordinates](/files/-LVxTj_8RmgfniCLDREH)
{% endtab %}
{% endtabs %}

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.

{% tabs %}
{% tab title="Code" %}

```python
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()
```

{% endtab %}

{% tab title="Output" %}
![Toy representation of the galactic plane in galactic coordinates ](/files/-LVxURcpgrk2JHxhaKGb)
{% endtab %}
{% endtabs %}

{% hint style="warning" %}
Usually right ascension is plotted from right (0°) to left (360°), but it can also be expressed in hours (1 hr = 15°).
{% endhint %}

## 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]
$$

{% tabs %}
{% tab title="Code" %}
First we are going to use the `astropy.cosmology` module to retrieve the values of the cosmological constants using Planck data from 2013.&#x20;

```python

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


```

```python
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)
```

```python
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()
```

{% endtab %}

{% tab title="Output" %}
![](/files/-M23TnelWF7l6m00ovrA)
{% endtab %}
{% endtabs %}


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://astroparticle.gitbook.io/docs/basic-concepts/tutorial-1.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.