Import standard modules:
In [4]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import HTML
HTML('../style/course.css') #apply general CSS
Out[4]:
In [7]:
from IPython.display import HTML
HTML('../style/code_toggle.html')
import healpy as hp
%pylab inline
pylab.rcParams['figure.figsize'] = (15, 10)
import matplotlib
import ephem
We can use a geographical coordinate system to uniquely identify a position on earth. We normally use the coordinates latitude $L_a$ (to measure north and south) and longitude $L_o$ (to measure east and west) to accomplish this. The equatorial coordinate system is depicted in Fig. 3.1.1 ⤵.
Figure 3.1.1: The geographical coordinates latitude $L_a$ and longitude $L_o$.
We also require a coordinate system to map the celestial objects. For all intents and purposes we may think of our universe as being projected onto a sphere of arbitrary radius. This sphere surrounds the Earth and is known as the celestial sphere. This is not a true representation of our universe, but it is a very useful approximate astronomical construct. The celestial equator is obtained by projecting the equator of the earth onto the celestial sphere. The stars themselves do not move on the celestial sphere and therefore have a unique location on it. The Sun is an exception, it changes position in a periodic fashion during the year (as the Earth orbits around the Sun). The path it traverses on the celestial sphere is known as the ecliptic.
The north celestial pole (NCP) is an important location on the celestial sphere and is obtained by projecting the north pole of the earth onto the celestial sphere. The star Polaris is very close to the NCP and serves as a reference when positioning a telescope.
The south celestial pole (SCP) is obtained in a similar way. The imaginary circle known as the celestial equator is in the same plane as the equator of the earth and is obtained by projecting the equator of the earth onto the celestial sphere. The southern hemisphere counterpart of Polaris is KT:GM: Do you want to add Sigma Octanis to the Glossary? Sigma Octanis.
We use a specific point on the celestial equator from which we measure the location of all other celestial objects. This point is known as the first point of Aries ($\gamma$) or the vernal equinox. The vernal equinox is the point where the ecliptic intersects the celestial equator (south to north). We discuss the vernal equinox in more detail in $\S$ 3.2.2 ➞ .
We use the equatorial coordinates to uniquely identify the location of celestial objects rotating with the celestial sphere around the SCP/NCP axis.
The Right Ascension $\alpha$ - We define the hour circle of an object as the circle on the celestial sphere that crosses the NCP and the object itself, while also perpendicularly intersecting with the celestial equator. The right ascension of an object is the angular distance between the vernal equinox and the hour circle of a celestial object measured along the celestial equator and is measured eastward. It is measured in Hours Minutes Seconds (e.g. $\alpha = 03^\text{h}13^\text{m}32.5^\text{s}$) and spans $360^\circ$ on the celestial sphere from $\alpha = 00^\text{h}00^\text{m}00^\text{s}$ (the coordinates of $\gamma$) to $\alpha = 23^\text{h}59^\text{m}59^\text{s}$.
The Declination $\delta$ - the declination of an object is the angular distance from the celestial equator measured along its hour circle (it is positive in the northern celestial hemisphere and negative in the southern celestial hemisphere). It is measured in Degrees Arcmin Arcsec (e.g. $\delta = -15^\circ23'44''$) which spans from $\delta = -90^\circ00'00''$ (SCP) to $+\delta = 90^\circ00'00''$ (NCP).
The equatorial coordinates are presented graphically in Fig. 3.1.2 ⤵ .
Figure 3.1.2: The equatorial coordinates $\alpha$ and $\delta$. The vernal equinox $\gamma$, the equatorial reference point is also depicted. The vernal equinox is the point where the ecliptic (the path the sun traverses over one year) intersects the celestial equator. KT:XX: What are the green circles in the image?
We will be making use of the `pyephem` package ⤴ package in the rest of this chapter to help us clarify and better understand some theoretical concepts. The two classes we will be using are the Observer and the Body class. The Observer class acts as a proxy for an array, while the Body class embodies a specific celestial object. In this section we will only make use of the Body class.
Earlier in this section I mentioned that the celestial objects do not move on the celestial sphere and therefore have fixed equatorial coordinates. This is not entirely true. Due to the precession (the change in the orientation of the earth's rotational axis) the location of the stars do in fact change minutely during the course of one generation. That is why we need to link the equatorial coordinates of a celestial object in a catalogue to a specific observational epoch (a specific instant in time). We can then easily compute the true coordinates as they would be today given the equatorial coordinates from a specific epoch as a starting point. There are two popular epochs that are often used, namely J2000 and B1950. Expressed in UT (Universal Time) ⤴:
The 'B' and the 'J' serve as a shorthand for the Besselian year and the Julian year respectively. They indicate the lenght of time used to measure one year while choosing the exact instant in time associated with J2000 and B1950. The Besselian year is based on the concept of a tropical year ⤴ and is not used anymore. The Julian year consists of 365.25 days. In the code snippet below we use pyephem to determine the J2000 and B1950 equatorial coordinates of Arcturus.
In [ ]:
arcturus = ephem.star('Arcturus')
arcturus.compute('2016/2/8',epoch=ephem.J2000)
print('J2000: RA:%s DEC:%s' % (arcturus.ra, arcturus.dec))
arcturus.compute('2016/2/8', epoch=ephem.B1950)
print('B1950: RA:%s DEC:%s' % (arcturus.a_ra, arcturus.a_dec))
To finish things off, let's make sure that given the concepts we have learned in this section we are able to interpret a radio skymap correctly. We will be plotting and inspecting the Haslam 408 MHz map ⤴. We load the Haslam map with read_map and view it with cartview. These two functions form part of the `healpy` package ⤴.
In [ ]:
haslam = hp.read_map('../data/fits/haslam/lambda_haslam408_nofilt.fits')
In [ ]:
matplotlib.rcParams.update({'font.size': 10})
proj_map = hp.cartview(haslam,coord=['G','C'], max=2e5, xsize=2000,return_projected_map=True,title="Haslam 408 MHz with no filtering",cbar=False)
hp.graticule()
The cartview function also produces a projected map as a byproduct (it takes the form of a 2D numpy array). We can now replot this projected map using matplotlib (see Fig. 3.1.3 ⤵ ). We do so in the code snippet that follows.
In [ ]:
fig = plt.figure()
ax = fig.add_subplot(111)
matplotlib.rcParams.update({'font.size': 22})
#replot the projected healpy map
ax.imshow(proj_map[::-1,:],vmax=2e5, extent=[12,-12,-90,90],aspect='auto')
names = np.array(["Vernal Equinox","Cassiopeia A","Sagitarius A","Cygnus A","Crab Nebula","Fornax A","Pictor A"])
ra = np.array([0,(23 + 23./60 + 24./3600)-24,(17 + 42./60 + 9./3600)-24,(19 + 59./60 + 28./3600)-24,5+34./60+32./3600,3+22./60+41.7/3600,5+19./60+49.7/3600])
dec = np.array([0,58+48./60+54./3600,-28-50./60,40+44./60+2./3600,22+52./3600,-37-12./60-30./3600,-45-46./60-44./3600])
#mark the positions of important radio sources
ax.plot(ra,dec,'ro',ms=20,mfc="None")
for k in xrange(len(names)):
ax.annotate(names[k], xy = (ra[k],dec[k]), xytext=(ra[k]+0.8, dec[k]+5))
#create userdefined axis labels and ticks
ax.set_xlim(12,-12)
ax.set_ylim(-90,90)
ticks = np.array([-90,-80,-70,-60,-50,-40,-30,-20,-10,0,10,20,30,40,50,60,70,80,90])
plt.yticks(ticks)
ticks = np.array([12,10,8,6,4,2,0,-2,-4,-8,-6,-10,-12])
plt.xticks(ticks)
plt.xlabel("Right Ascension [$h$]")
plt.ylabel("Declination [$^{\circ}$]")
plt.title("Haslam 408 MHz with no filtering")
#relabel the tick values
fig.canvas.draw()
labels = [item.get_text() for item in ax.get_xticklabels()]
labels = np.array(["12$^h$","10$^h$","8$^h$","6$^h$","4$^h$","2$^h$","0$^h$","22$^h$","20$^h$","18$^h$","16$^h$","14$^h$","12$^h$"])
ax.set_xticklabels(labels)
labels = [item.get_text() for item in ax.get_yticklabels()]
labels = np.array(["-90$^{\circ}$","-80$^{\circ}$","-70$^{\circ}$","-60$^{\circ}$","-50$^{\circ}$","-40$^{\circ}$","-30$^{\circ}$","-20$^{\circ}$","-10$^{\circ}$","0$^{\circ}$","10$^{\circ}$","20$^{\circ}$","30$^{\circ}$","40$^{\circ}$","50$^{\circ}$","60$^{\circ}$","70$^{\circ}$","80$^{\circ}$","90$^{\circ}$"])
ax.set_yticklabels(labels)
ax.grid('on')
To make sure that you are able to interpret Fig. 3.1.3 ⤵ correctly try to find the location of Pictor A from first principles. The equatorial coordinates of Pictor A is $\alpha = 5^\text{h}19^\text{m}49.7^\text{s}$ and $\delta = -45^\circ 46'44''$.