This tutorial assumes you have gone through the Trace and Wavelength Calibration tutorials and have their results available.
Adam Ginsburg, Kelle Cruz, Lia Corrales, Jonathan Sick, Adrian Price-Whelan
Spectroscopy
This tutorial, the third in a series, shows how to apply a trace and a wavelength solution to science data. It then shows how to do basic analysis, i.e., line fitting.
from PIL import Image as PILImage
import numpy as np
import pylab as pl
pl.style.use('dark_background')
from astropy.modeling.models import Linear1D
from astropy import constants
from astropy import units as u
from astropy.visualization import quantity_support
quantity_support()
<astropy.visualization.units.quantity_support.<locals>.MplQuantityConverter at 0x7f46c6d11ac0>
image_array_2 = np.array(PILImage.open('deneb_3s_13.63g_1.bmp'))
from astropy.modeling.polynomial import Polynomial1D
from astropy.modeling.models import Gaussian1D
from astropy.modeling.fitting import LevMarLSQFitter, LinearLSQFitter
linfitter = LinearLSQFitter()
yaxis2 = np.repeat(np.arange(470, 520)[:,None], image_array_2.shape[1], axis=1)
xvals = np.arange(image_array_2.shape[1])
weighted_yaxis_values2 = np.average(yaxis2, axis=0, weights=image_array_2[470:520,:] - np.median(image_array_2))
polymodel2 = Polynomial1D(degree=3)
fitted_polymodel2 = linfitter(polymodel2, xvals, weighted_yaxis_values2)
trace_center2 = fitted_polymodel2(xvals)
npixels_to_cut = 15
trace_center = fitted_polymodel2(xvals)
cutouts = np.array([image_array_2[int(yval)-npixels_to_cut:int(yval)+npixels_to_cut, ii]
for yval, ii in zip(trace_center, xvals)])
cutouts.shape
mean_trace_profile = cutouts.mean(axis=0)
spectrum2 = np.array([np.average(image_array_2[int(yval)-npixels_to_cut:int(yval)+npixels_to_cut, ii],
weights=mean_trace_profile)
for yval, ii in zip(trace_center2, xvals)])
Next, we retrieve the wavelength solution derived in Part 2.
wlmodel = Linear1D(slope=-0.10213643, intercept=562.3862495)
wavelengths = wlmodel(xvals) * u.nm
pl.plot(wavelengths, spectrum2)
[<matplotlib.lines.Line2D at 0x7f46c4572790>]
Now, we go on to do some basic analysis on our fully extracted and wavelength-calibrated spectrum
We zoom in on the 4860 Angstrom line - H-Beta
pl.plot(wavelengths, spectrum2)
pl.axis([470,510,140,190])
(470.0, 510.0, 140.0, 190.0)
We can use astropy models to construct an absorption line model, consisting of a continuum level and a negative Gaussian to represent the absorption feature
absorption_model_guess = Linear1D(slope=0, intercept=175) + Gaussian1D(amplitude=-25, mean=486, stddev=3)
We can overplot our guessed model - it's not right, but it's in the right spot
pl.plot(wavelengths, spectrum2)
pl.plot(wavelengths, absorption_model_guess(wavelengths.value))
pl.axis([470,510,140,190])
(470.0, 510.0, 140.0, 190.0)
The Levenberg-Marquardt Least Squares fitter can be used to find the optimal fit to our data given the starting guess.
lmfitter = LevMarLSQFitter()
selection = (wavelengths > 470*u.nm) & (wavelengths < 500*u.nm)
fitted_absorption_model = lmfitter(model=absorption_model_guess, x=wavelengths.value[selection], y=spectrum2[selection])
pl.plot(wavelengths, spectrum2)
pl.plot(wavelengths, fitted_absorption_model(wavelengths.value))
pl.axis([470,510,140,190])
(470.0, 510.0, 140.0, 190.0)
We can now separate out the two components, the continuum and the absorption line model:
continuum_fit, absorption_fit = fitted_absorption_model
If we plot the data minus the absorption line model, we get a nice "line-free continuum":
pl.plot(wavelengths, spectrum2)
pl.plot(wavelengths, spectrum2 - absorption_fit(wavelengths.value))
pl.axis([470,510,140,190])
(470.0, 510.0, 140.0, 190.0)
The Equivalent Width of a spectral absorption line is defined to be the width of a feature that has the same integral as the absorption line, but goes all the way from the continuum level to zero.
We can compute this with our model, assuming our continuum is flat (has zero slope):
EQW = -absorption_fit(wavelengths.value[selection]).sum() / continuum_fit.intercept * u.nm
EQW
absorption_fit
<Gaussian1D(amplitude=-22.81529847, mean=486.34296327, stddev=0.4969129)>
We have identified the line as H-beta, so we can measure some of its properties now:
air_wavelength_hbeta = 486.135*u.nm # wikipedia https://en.wikipedia.org/wiki/Balmer_series
The Doppler Shift will tell us the velocity of the star. Note that this is the velocity in the topocentric frame, i.e., in the rest frame of the observatory. If we knew precisely when and where these observations were taken, we could convert this velocity to the heliocentric or LSR frames with astropy tools.
First, we do the calculation manually, following the optical definition $$v_{opt} = c \frac{\lambda-\lambda_0}{\lambda} $$
doppler_velocity = (absorption_fit.mean*u.nm - air_wavelength_hbeta) / (air_wavelength_hbeta) * constants.c
doppler_velocity.to(u.km/u.s)
We can do the same thing using astropy equivalencies:
doppler_velocity = (absorption_fit.mean*u.nm).to(u.km/u.s, u.doppler_optical(air_wavelength_hbeta))
doppler_velocity
We can also compute the line width from our fit:
linewidth_kms = (absorption_fit.stddev*u.nm) / air_wavelength_hbeta * constants.c.to(u.km/u.s)
linewidth_kms
That's it! You've extracted and wavelength-calibrated a spectrum.
The next tricky step is flux calibration, but we will do that in a subsequent tutorial because we need a different data set; this one didn't include the necessary data for flux calibration (though we could approximately calibrate the spectrum using a known magnitude over the observed band for Deneb)