Friday, June 10, 2016 - Week Recap
This week my team’s advisor Jeff Bary was away at the Cool Stars 19 meeting in Upsala, Sweden. In his absence, we had the freedom to explore our new project -- something I feel I took full advantage of.
One of the main goals Jeff gave us for the week was to get used to using TripleSpecTool, a GUI constructed to do reductions of images taken with the TriplesSpec spectrograph on the ARC 3.5m telescope at APO. On Friday of last week, I finished getting IDL set up on the lab computers so we were ready to work this week.
I quickly felt comfortable with the GUI, which did many of the harder parts of spectrum reduction (from my experience with IRAF) automatically. I reduced the data for one night of prior TripleSpec data, and felt comfortable enough not to continue doing the same thing all week.
Instead, I transferred some of the reduced spectra to my laptop and began looking at how we should be analyzing the spectra in the coming weeks. We are investigating the reddening laws in dense cloud cores by taking spectra of bright stars behind the clouds. When the light from these stars passes through the dense cloud, the light is extincted, and reddened. To measure that reddening, we have to know what spectral type (OBAFGKM) and luminosity class (I-V) these background stars are. Thus, spectral classification is going to be a large part of this project.
After searching the literature, I found several papers focusing on stellar classification in the near infrared regime we will be studying, (Wallace 1997, Meyer 1998, Coelho 2005, Winge 2009) and noticed that they all report either continuum subtracted or continuum normalized spectra. I figured a good place to start was to look at determination of the stellar continuum (the overall shape of the spectral energy distribution in the region we are studying).
I read that polynomial fitting was common, so I coded it up in Python and met with less than stellar results. There seemed to be several regions of the spectrum that were not well fit, even with high (<15) degree polynomials. So I looked at using other methods, and thought to use one of my favorites: the spline.
At the heart, a spline is an interpolation technique for fitting the regions in between data points. But unlike polynomial interpolation, where the interpolating function is valid globally, spline is a piecewise function with a different polynomial fit between each of its “knots,” in a pure interpolation spline, each knot is a data point.
However, there are also smoothing splines, where the number of knots is less than the number of data points. When N(knots) < N(data), the spline interpolation is rougher than it would be if all data points are used as knots: thus, the smoothing spline, as its name suggests, can produce a smoothed fit. I used a cubic spline, which is continuous up to its second derivative, with a constrained error of 0.8 Jansky. The algorithm adds knots up until the total error of the spline is less than the constrained error.
The results of this analysis are summarized in the plot below. I chose to subtract the continuum, but another valid technique is to divide the original spectrum by the continuum, which results in a normalised spectrum (average value 1, not 0).
I also note the wavelength ranges where other studies overlap with ours in the box to the lower left.
After completing this, I still felt frustrated about the quality of the continuum fit. It still seems to me that there is a slight trend in the region 1.4-1.8 microns and in the CO overtone bands where the continuum was not fit perfectly, because of the strength of the absorption lines in these regions. Thus I looked on for another solution.
I understood that IRAF has splot, which does things similar to what I’m trying to do, but since updating the lab computers to OS X 10.11, it seems the IRAF installations broke, and after looking on iraf.net, it seems many others experienced this as well. Without admin access to fix the installation, I had to try something else.
I found a simple GUI program written in Python at http://python4esac.github.io/plotting/specnorm.html for spectrum normalisation, and I thought it was worthwhile to investigate.
The core functionality is simple but robust: a user marks several continuum spots over the spectrum, then a cubic spline is fit to those points and the spectrum is divided by that spline and replotted.
After looking at the code, I realized that the original author had set up the interactive part of the code for me, so that adding more functionality would be as easy as picking a key and writing code. So I got to it.
I altered the code to read in .fits files correctly, converted some units from angstroms to microns that caused bugs with our spectra, and, most importantly, implemented equivalent width fitting functionality.
As evidenced by a thread on astrobetter in which professional astronomers express frustration at the current tools available for analysis of spectral features, I decided none of the suggested tools looked particularly like what I wanted. So I started tinkering in the specnorm.py code and implemented spectral feature fitting of my own.
One of the quantities of these spectral features we are most interested is their equivalent width (see https://en.wikipedia.org/wiki/Equivalent_width ) which requires a robust integral of the spectral feature and an estimate of the continuum around the feature.
I currently have equivalent width fitting implemented in two ways: after the user selects a point on either side of the feature, both a cubic spline and a gaussian profile are fit to the data in the user-defined region. The fits are integrated and the equivalent widths and the full width at half maximum of the gaussian are printed to the terminal.
I also experimented with MCMC in python (emcee package) this week, as I thought perhaps one way we might determine spectral types would be by fitting spectral classification templates to our data, and the best way to do that would probably be MCMC.
Overall I think it was a good week. I might see about improving the gaussian fitting of the code, but we’ll see what Jeff thinks of it next week. I might also want to look into fitting Voigt profiles -- but that would have to be for another post.
