Stellar Temperatures from Spectra#
Learning Objectives
learn to read in tabular data to Python
be able to calculate a chi-squared goodness of fit statistic
practice for loops and indexing
learn data visualization skills
Understand how stellar spectra change with effective temperature & how to quantify how well two datasets match.
Introduction#
In this activity, we’re going to look at spectra of some real stars and of computer models, compare stars of different spectral types, and finding the temperatures of our mystery stars.
The files for this activity contain the model spectra for different effective temperatures (“models.csv”) and the observed spectra (“Star#.csv”) for 3 real stars. The files have:
Wavelength – units of nm
Normalized Flux – the fluxes have been “continuum normalized”, meaning that the overall blackbody shape has been divided out. So the flux values vary from 0-1 in value, rather than normal erg/cm2/s/nm units.
You can download them here: models.csv | Star 1 | Star 2 | Star 3
Hint
Be sure to have your data files in the same directory as your code where you’re loading the data!
Part 1 - Model spectra#
Start by importing the python packages we’ll need by running the cell below.
# start by importing useful packages
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import ascii
%matplotlib inline
from astropy.table import Table
Our model spectra correspond to effective temperatures of 4000, 6000, 8000 and 10,000 K. We will use the custom function below (“get_model”) to read in the model for a given temperature, then return the wavelength and flux of the model spectrum in datatable format (“model”).
Run the cell below so we can use this function.
def get_model(temperature):
# input = effective temperature
templist = [4000, 6000, 8000, 10000]
if temperature not in templist:
print(' Please choose a new temperature:')
print(' possible values =', templist)
pass
else:
data = ascii.read("models.csv")
model = Table([data["wavelength"], data["flux_"+str(temperature)]])
model.rename_column("flux_"+str(temperature), 'flux')
return(model)
Exercise 1
Let’s practice reading in a model spectrum! Choose a temperature in our model range and write it in the empty parentheses in the cell below. Then run the cell to print out the resulting model table values. An example scaffold is provided.
model = get_model( )
print(model)
Exercise 2
Write code to plot the model spectrum. You can pull columns from astropy tables (like "wavelength" above) by indexing the object with just the column name. Include reasonable axis labels and change the plot ranges as needed.
Pro tip - Add minor ticks using after the plt.plot() command:
plt.minorticks_on()
Exercise 3
Let’s look at the data for the other models; to do this, we’re going to make a loop that runs over effective temperature to (1) read in each model and (2) plot the spectrum on our graph.
Make a “teff” numpy array that holds the effective temperature values of our models (4000-10000 in steps of 2000 K).
Make a For-loop that iterates over the effective temperature. In each iteration, read in the model for this temperature, and then make a plot of the spectrum. Bonus - Add a title with the model’s temperature so we can tell them apart.
As a reminder, a basic For Loop looks like:
for i in range(0, 10): print('iteration', i)
Exercise 4
How does the spectrum change as the model’s effective temperature changes? Qualitatively describe your intuition from examining the model plots.
Part 2 - observed spectra#
Now, we’re going to compare the observed spectra of a real star to the models from Part I and determine this star’s spectral type (and surface temperature).
Exercise 5
Read in the spectrum for Star 1 using the ascii.read() package. (Like in Step 1).
Plot the spectrum. What do you notice about the spectrum? Take some time to modify your plot boundaries/axes, labels, and aspect ratio to best see the relevant features in the spectrum
One way to determine the spectral type of Star 1 would be “by eye”, where you visually choose which model looks most similar to Star 1. A better way would be comparing the spectra using a chi-squared (\(\chi^2\)) goodness-of-fit statistic:
This statistic quantifies how different (or similar) two data sets are, where the “best fitting” model is going to have the lowest \(\chi^2\) value. Notice that this metric is computed as the sum over the normalized difference between every equivalent point in the model to the data; individual points with very large discrepancies can thus drive the \(\chi^2\) metric up considerably.
Exercise 6
Calculate the \(\chi^2\) value for Star 1 compared to the 4000 K model, then print out the result. You should get a \(\chi^2\) value around 400, which isn’t very good. We want a low number, so let’s try other spectral types…
Hint - You can use np.sum() to add up all the values in an array; this calculation does not require a for-loop
We’re going to calculate \(\chi^2\) value for each of the temperatures, and save these values in an array called “chi2”. This array will have the same size as the effective temperature array you made earlier.
First, we need to make an empty numpy array to store the \(\chi^2\) values. You can make a blank numpy array using:
x = np.zeros(length)
where “length” is the number of elements in the array.
Exercise 7
Write the code to make a “chi2” array with the same number of elements as your “teff” array.
Hint – You can find the length of an array using “len(array)”.
Next, create a for-loop that iterates over temperature (as in Part I), calculating the \(\chi^2\) for each model with Star 1, and store the output in the chi2 array.
Now we can find the minimum \(\chi^2\) value and find the best-fit temperature for Star 1.
Exercise 8
Make a plot of effective temperature versus chi2, including axis labels and whatever colors/symbols you like! Which model temperature best fits the data? Put another way, using our modeling, what is our best estimate for the temperature of Star 1?
Exercise 9
Make a plot of Star 1’s spectrum with the best fitting model (determined from your \(\chi^2\) calculation) plotted on top to see how well they match. (They should match pretty well!)
Part 3 - New Mystery Stars#
Now we’ll move on to the other stars, repeating the process from Part II.
Exercise 10
Repeat the process in Part 2 for the other mystery stars (Stars 2 & 3) to find their spectral types. For each star,
Read in the mystery star’s spectrum.
Calculate the \(\chi^2\) values for all of the model spectra. Make a chi2 vs temperature plot to find the best fitting spectral type.
Plot the observed spectrum and the best-fit model spectrum to confirm that they match.
Determine the effective temperature for the two stars.
Concluding thoughts#
Question: How well did the model spectra fit the real stars’ spectra? Do you think our fitting code worked successfully?
Question: What’s one situation where our code would NOT be able to find the temperature of the star? Is there a way we could update our code / analysis to solve the problem?
Write out a discussion for these two questions
Bonus/Extra Credit: Array Manipulation#
In the exercises above, we used a for-loop to iterate over all the models in our library, comparing each to an individual stellar spectrum. Each model has the same length as the real spectrum (which is why the model-data array-wise subtraction worked). This means that it is in fact possible to compute our chi2 array for all models and a given star in one step, without using a for-loop.
You may have learned that to compute the point-wise, e.g., sum, multiplication, etc., between two arrays, they have to be the same shape. That’s not strictly true, however, due to the ability for numpy arrays to be broadcast. You’ve probably used this technique before actually! When we compute
a = np.array([1,2,3,4,5])
b = a - 2
print(b)
array([-1,0,1,2,3])
the actual “array” subtraction is array([1,2,3,4,5]) - array([2,2,2,2,2]). Numpy is smart enough to extend our integer (2) to match the length of a to facilitate the subtraction. Sometimes, in more complex scenarios (like our spectra), we can be a bit more explicit and tell numpy how to broadcast our arrays to facilitate a single, vectorized computation.
Exercise 11
Doing some googling if neccessary, determine how to use array broadcasting to turn the computation of the chi2 array for a star and all our models into a single line calculation using numpy arrays.
For an even trickier challenge, determine how to find all the chi2 arrays for an arbitrary number of mystery stars with all the models in a single numpy operation. For this one, you’ll need to stack up the mystery stars into one array, and the models into another.
Note: You may need to adjust how you load in the data to get all temperature models into a single array.