Astro 5830 Homework 1

This notebook illustrates how to retreive files from a web site, read them into an astropy table, and plot them interactively with matplotlib. This is most of the code necessary to complete the stellar classification exercise in the homework.

Step 1: Define Modules. Note that the '%matplotlib notebook' magic command enables interactive plots with matplotlib.

from astropy.io import ascii
%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
import urllib.request

Step 2: Define variables for files to retrieve. Change 'lastname' to your lastname.

url = "http://www.astronomy.ohio-state.edu/~martini/Astro5830/StellarSpectra/"
files = []
lastname = "martini" # Your last name here
for i in range(5):
  files.append(lastname+str(i)+".txt")

Step 3: Retrieve the files with urllib.request.

for file in files: 
  urllib.request.urlretrieve(url+file, file)

Step 4: Interactively plot the first object. The interactive mode is useful to zoom into spectra features to inspect their strength and wavelength. Hit the stop button when you are done with interactive mode.

index = 0
star = files[index]
spec = ascii.read(star)
spec['col2'] = 1e12*spec['col2'] # rescale 
ax = plt.axes()
ax.plot(spec['col1'], spec['col2'], 'k-', label=star)
ax.legend(loc="upper right", frameon=False)
ax.set_xlabel(r'Wavelength ($\AA$)')
ax.set_ylabel('Flux')
plt.show()

Step 5: Create a plot with labels. For this example I've created a dictionary with information for lines, written a routine to draw the line, and then I made the plot. Note there is an option to download the plot.

lines=dict()
lines['Hbeta']=dict()
lines['Hbeta']['lambda']=4861.
lines['Hbeta']['label']=r"H$\beta$"
lines['Hbeta']['index']=np.searchsorted(spec['col1'], lines['Hbeta']['lambda'], side="left")
lines['Hbeta']['fluxval']=np.max(spec['col2'][lines['Hbeta']['index']-15:lines['Hbeta']['index']+15])

lines['Halpha']=dict()
lines['Halpha']['lambda']=6563.
lines['Halpha']['label']=r"H$\alpha$"
lines['Halpha']['index']=np.searchsorted(spec['col1'], lines['Halpha']['lambda'], side="left")
lines['Halpha']['fluxval']=np.max(spec['col2'][lines['Halpha']['index']-15:lines['Halpha']['index']+15])

def plotline(line, dy):
    ax.plot([lines[line]['lambda'], lines[line]['lambda']], [lines[line]['fluxval']+0.1*dy, lines[line]['fluxval']+0.2*dy], 'k-')
    ax.text(lines[line]['lambda'], lines[line]['fluxval']+0.25*dy, lines[line]['label'], fontsize=14, ha='center')

spclass = ["F6V", "M5III", "K0III", "G6V", "B8V"] # my classifications
index = 0
star = files[index]
spec = ascii.read(star)
spec['col2'] = 1e12*spec['col2']
ax = plt.axes()
newlabel = files[index] + " (" + spclass[index] + ")"
ax.plot(spec['col1'], spec['col2'], 'k-', label=newlabel)
ax.legend(loc="upper right", frameon=False)
ax.set_xlabel(r'Wavelength ($\AA$)')
ax.set_ylabel('Flux')
[y1,y2] = ax.get_ylim()
dy = y2-y1
plotline("Hbeta", dy)
plotline("Halpha", dy)
plt.show()