Close star, Distant Star, Near Star….that just became a Far Star….

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

https://drive.google.com/file/d/0B-EDLU7A4RFSTllFckJRdnlVa2M/view?usp=sharing

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)