Over the course of this semester, I have been utilizing the program Sherpa, written in Python, to model Potassium (K) absorption lines and make three specific measurements. Specifically, these measurements are the equivalent width, full width half maximum, and line-to-continuum flux ratio. In this entry, I will show the commands I utilized to make these measurements.
When I use Sherpa, I import it as a python module as such:
>>> from sherpa.astro.ui import * # All Sherpa functions are found here
If you followed my post on modeling with Sherpa, then I have used 3 separate model components to create my model. These include the power law (pl), the lorentzian (l1) and the gaussian (g1). The power law describes the continuum, and the three together represent the data.
To measure the equivalent width of the K line, I use Sherpa’s eqwidth() function. The syntax is:
>>>eqwidth(continuum, continuum+line, lo=lower wl limit, hi=upper wl limit, id=data ID #)
Using my data for a K line whose line goes from 1.252µ to 1.253µ, I measure the eq. width like so:
>>> eqwidth(pl, pl+l1+g1, lo=1.252, hi=1.253, id=1)
The output will be the equivalent width value, in microns. Nice and simple. The hard part was constructing the models and eyeballing the lower/upper limits for the eq. width measurement. All of this I did manually by eye.
The line-to-continuum flux ratio and full width half maximum measurements are harder. They require locating the actual arrays of values for the various models. In order to access the model’s x and y values (continuum+lorentzian+gaussian), you do:
>>> modelArrays = get_model_plot(1) # The 1 is the ID number
# modelArrays.x –> array of x values
# modelArrays.y –> array of y values
To access the y values for the continuum only (pl), you have to do:
>>> key = pl._cache.keys()
>>> cont_yvals = pl._cache[key]
# This may not work if some of the values in pl.__dict__ are not set properly. To be honest, # I’m not fully sure which ones may screw this up, but for me, I have:
>>> pl._NoNewAttributesAfterInit__initialized = True
>>> pl._use_caching = True
# For me, since I use a power law with an exponent of 0, all of the y values are the same. So I can access this value instead by doing:
>>> cont_yval = pl.ampl # The continuum is a constant value: the amplitude
Armed with the model’s x and y value arrays, and the continuum’s y value array (or value), we can make the measurements. To calculate the line-to-continuum flux ratio, we must first find the minimum (maximum for emission lines) value of the model, and then compare it with the corresponding continuum flux value:
>>> F_model = numpy.min(modelArrays.y) # Minimum of the model
>>> min_ind = numpy.where(modelArrays.y == F_model) # min’s array index
>>> F_cont = cont_yvals[min_ind] # Continuum value
>>> depth = 1.0 – (F_model/F_cont) # The line-to-continuum flux ratio!
The full width half maximum was a little tougher. The algorithm I used is as follows:
>>> F_halfmax = (F_cont + F_max) / 2.0 # Flux at the half maximum
>>> indices = numpy.where(modelArrays.y <= F_halfmax)
>>> wls = modelArrays.x[indices]
>>> fwhm = wls[-1] – wls
The problem with this algorithm is that sometimes there are pockets of spectrum that are equal to or lower than the half max, which get into the wls array. For those, I had to go through them manually and fix the wls array.
That’s how I measured the equivalent width, line-to-continuum flux ratio, and full width half maximum using Sherpa. If you have any questions or know of a better (or different) way of calculating these values, feel free to leave a comment or email me at firstname.lastname@example.org
Thank you for another good post. A few suggestions:
To get the actual x/y data of a single component of your fitted model, sherpa has a command called get_model_component_plot(pl). Like always, the different arrays can be accessed by:
cont = get_model_component_plot(pl)
cont_yvals = a.y
Similar goes for a.x, a.xerr, a.yerr, etc., and the x/y values of the line components can of course be accessed by calling get_model_component_plot() with the relevant component as argument. This way, you don’t have to mess around with pl.__dict__ and other hairy stuff.
Also, when you do your FWHM-stuff, there is a nice little trick to avoid contamination by other lines, which is simply to only consider a wavelength interval centered around the line, i.e. something like:
init_width = [Width of a wl interval around the line]
line_center = [wavelength of wherever the line center is]
indices = scipy.where((modelArrays.y = line_center-init_width) &
(modelArrays.x <= line_center+init_width))
wls = modelArrays.x[indices] (etc.)
This is possible because the size of the x and y arrays are the same, so a statement like " a = y[where(x <= some_wavelength)] " is completely valid.
Thanks for the tips! I will definitely try that instead.
Oooops… Of course, the first piece of code should read:
cont = get_model_component_plot(pl)
cont_yvals = cont.y
Sorry ’bout that!
I am having a bit of trouble finding a reasonable equivalent width. I am focusing on a set of data from .fits files around a sodium absorption line and I was able to get a good fit with set_source(1, powlaw1d.pl+gauss1d.g1+lorentz1d.l1). I set g1.pos and l1.pos = 6360 which is where the center of the absorption is located and tried finding the equivalent width with the method you’ve mentioned above: eqwidth(pl, pl+l1+g1, lo=6340., hi=6385., id=1)
however, I am getting some outrageous number (10^39). I’ve also tried converting my lo and hi values to microns in case that was the problem. Do you have any suggestions as to what I might be doing incorrectly?