We are
bounded in a nutshell of Infinite space: Blog Post #14, Worksheet # 5.2: Close star, Distant
Star, Near Star….that just became a Far Star….
As
we previously learned, Cepheid variables are a special class of stars that
radially pulsate in a predictable way. In 1908, Henrietta Swan Leavitt
discovered that there is a distinct relationship between a Cepheid’s luminosity
and pulsation period by examining many stars in the Magellanic Clouds.
Henrietta was a member of “Harvard’s computers,” a group of women hired by
Edward Pickering to analyze stellar spectra and light curves. In this
worksheet, we will use Henrietta’s original data set to find our own
Period-Luminosity relation for Cepheid variables.
1.
The data file, “Cepheid variables.csv,” contains data for 25 Cepheid variables
located in the Small Magellanic Cloud (SMC). Each line contains a specific
Cepheids: (1) ID number, (2) Maximum apparent magnitude, (3) Minimum apparent
magnitude and (4) Period. Calculate the mean apparent magnitude for each Cepheid.
2. The distance to the SMC is about 60 kpc,
where kpc = 1000 pc. Convert your mean apparent magnitudes into mean absolute
magnitudes. Plot the Cepheid mean absolute magnitudes as a function of period.
This plot should look exponential.
3.
It is often handy to plot exponential (or power-law) functions with one or more
logarithmic axes, which “straightens out” the data. Magnitudes are already
exponential, so we don’t need to adjust that axis. Plot the Cepheid mean
absolute magnitudes as a function of \(\log_{10}
Period\) . Verify that the plot now
looks linear.
4.
Now that the data look linear, we can estimate the parameters of a linear
relation, \(M_V (P) = A\log_{10}(Period) +B\) . A
and B are “free parameters” that
allow the function to match the data.
5.
BONUS: You can determine the precise values of A and B by minimizing the
difference between the observed points and the model using the metric: where \[\chi^2 = \sum\limits_{i=1}^{25} (O_i – C_i )^2\] \(O_i\) is the observed value and \(C_i\) is the predicted (model) value.
We are now working with how astronomers
measure extreme distances once tools like parallax become impossible to use.
One way larger distances are measured is by finding clear examples of
astronomical markers, objects in the sky whose characteristics are constant
throughout the universe, whose traits can be analyzed and a correlation to
distance and intrinsic properties can be found. Cepheid variables are a type of
star that pulsate with light in a very predictable and measurable manner, and
when many of these variables are compared and plotted, a distance formula
emerges, as discovered by Henrietta Levitt.
In the case of this particular problem we
have been presented with, we are going to go through the entire worksheet in
one fail swoop, working it as individual components of a more complete system. Furthermore,
some programming experience would definitely help with plotting and fitting the
functions we will be seeing in this problem, but in any case we will be putting
up all the code used, and brief explanations for what they do. So for this
problem, we will be using some of the first data ever used to measure these
Cepheid Variable Stars, which you can access through
and with it have data as prescribed by the instructions
of the problem.
Using the magnitude-distance equation: \[M_{absolute}
= M_{apparent} – 5 \times (\log_{10} (Distance )-1 ) ,\] where the \(M_{apparent}\) is \[\frac{M_{max} +
M_{min} }{2} , \] we can plot the values given in the file and see how these
line up in a graph.
Up to know, we’ve had to use a lot of
different measurements and formulas for calculating what leads to the solution
to this problem, so let’s see the code that can help us find the answer: (%%
indicates a comment, not part of the code)
import numpy as np
import math
import matplotlib
import matplotlib.pyplot as plt
%%The
basic Python introductions to the code, which allows for the use of mathematical
and number systems, as well as a plotting element
c=np.loadtxt(fname='Cepheid_variables.csv',
delimiter=',')
id, max, min, period = c.T
%%Presents
the file (stored in the used computer) to the code and we give each column a
name that will be used in later
equations.
mean_apparent= (max+min)/2.0
mean_absolute = mean_apparent-
5.0*(math.log10(60000)-1)
%%The
two equations we were seeing earlier, which establish the values used to plot
the figure, with 60 times 1000 parsecs as the distance in the equation.
plt.figure(1)
%%This
establishes that this is the first of other figures to be presented at the same
time.
plt.plot(period, mean_absolute,
"bo")
plt.xlabel('Period (Days)')
plt.ylabel('Mean Absolute
Magnitude')
%%Now
we use these commands to name the axes of our graph, give some color to
definition to the set of points, and tell it to plot the exact values we have
been finding with the previous commands, and finally plot it:
Next,
we will begin analyzing the plot, and change it as necessary to find what we
wish to extract form it. Basically, this means we have to adjust the period
axis to it follows a logarithmic increase, and this convert what looks like an exponential
or logarithmic function into a linear function. This is done in Python by
writing:
plt.figure(2)
LogPeriod = np.log10(period) %% This changes the period
into a logarithmic quantity
plt.plot(LogPeriod,
mean_absolute, "bo")
plt.xlabel('Period (Days)')
plt.ylabel('Mean Absolute
Magnitude')
%%And
now we plot the figure again and get:
Now
finally we can start interpreting the data by placing a linear fit to the data
points and extracting the coefficients and constants A and B to find the complete
description of Cepheid periods and magnitudes. We use the script:
A,B = np.polyfit(LogPeriod,mean_absolute,1) %%The linear fit of the
equation defined…
plt.plot(LogPeriod, A*LogPeriod+B, "r") %%...and placed
plt.show()
%%The
final command which plots all outlined figures and presents them as
illustrations of the set of points
print(A,B)
%%Tells
the script to also print the values of A and B, the main question of the final problem.
And with this, plots the final plot and values we have been looking for:
(A,B) = (-2.0331864794721257, -2.7275593585821594)
Altogether,
the code looks like:
import numpy as np
import math
import matplotlib
import matplotlib.pyplot as plt
c=np.loadtxt(fname='Cepheid_variables.csv',
delimiter=',')
id, max, min, period = c.T
mean_apparent= (max+min)/2.0
mean_absolute = mean_apparent- 5.0*(math.log10(60000)-1)
plt.figure(1)
plt.plot(period, mean_absolute, "bo")
plt.xlabel('Period (Days)')
plt.ylabel('Mean Absolute Magnitude')
plt.figure(2)
LogPeriod = np.log10(period)
plt.plot(LogPeriod, mean_absolute, "bo")
plt.xlabel('Period (Days)')
plt.ylabel('Mean Absolute Magnitude')
A,B = np.polyfit(LogPeriod,mean_absolute,1)
plt.plot(LogPeriod, A*LogPeriod+B, "r")
plt.show()
print(A,B)
