Practicing plotting: Cepheid light curves

Let's put together some of the skills you've learned! You have a file called data/cepheid.txt, which contains flux measurements from a variable star over time. In this notebook, your task will be to measure the approximate distance to this star. Bold text makes clear what your task is in each cell.

Read in and plot the data

We have used the function ascii.read to read in data before (see 06-plotting.ipynb for example). Use ascii.read to load the file data/cepheid.txt into a variable data_table:



In [ ]:

This data table has a column called Time which tells you the time of each brightness measurement (units of days) and a column called Red and Blue which tells you the flux of the star through two filters.

Plot the red fluxes as a function of time and reproduce this plot using matplotlib. Don't forget to import matplotlib and to use the %matplotlib inline magic function within the notebook:



In [ ]:

This type of star is called a Cepheid variable star. These stars are useful to astronomers because their flux varies in a predictable sinusoidal pattern. It turns out that the period of the flux oscillation is directly related to how intrinsically bright the star is – so if you know how long the period is, and you measure how dim the star appears to be, you can estimate how far away it is because you know how bright it really is.

A very rough "model"

In the cell below, roughly measure the period of this Cepheid by plotting a sinusoidal function over the data. The function that you should play with is:

$$\textrm{flux} = \textrm{amplitude} \times \sin \left( \frac{2\pi}{\textrm{period}} \left( \textrm{time} - \textrm{offset} \right) \right)$$

Hint: start with an offset of zero and a period somewhere between 24 and 25 days, and remember that the sin function is within numpy.



In [ ]:

Period-luminosity relation

The brightness of this Cepheid variable star is related to its period by the period-luminosity relation:

$$ M_{v}=-2.43 \left(\log _{{10}}( \textrm{period})-1\right) - 4.05 $$

where $M_v$ is the absolute magnitude of the star. Note that brighter stars have more negative absolute magnitudes. Write a function called period_to_absolute_magnitude which takes the period as the argument, and returns the absolute magnitude.

Create a range of periods and plot the corresponding absolute magnitudes using your function, in order to reproduce this plot:



In [ ]:

Use your function from above to calculate the absolute magnitude of the star in the data, and add a point for the star to the plot:



In [ ]:

Measure the distance

To get the distance to the star, we need to use one more equation, which calculates the distance $d$ in parsecs to a star given its absolute magnitude $M$ and its apparent magnitude $m$:

$$ M = m - 5 (\log_{10}{d} - 1) $$

$$ d = 100 ^{\left(\frac{1 + \frac{m - M}{5}}{2}\right)}$$

If the apparent magnitude of this star is $m = 4$, calcuate the distance to the star using your absolute magnitude calculated above.



In [ ]:

Phase-folded light curve

When you find the period of some periodic light curve, we often want to line up each cycle of the observations and line them up with one another, so you can see the pattern in detail. This is called phase folding. You can phase fold the light curve replacing the times in your plotting commands with the modulus of the times in the light curve and the period:

$$\textrm{folded times} = \textrm{times} \,\%\, \textrm{period}$$



In [ ]:

Making a periodogram

Astronomers are often looking for periodic signals in various forms of data. One of the diagnostic plots that astronomers often make when looking for periodicities is the Lomb-Scargle periodogram, which shows how periodic a signal is over a range of possible periods, allowing you to pick out the approximate period.

The function lomb_scargle_periodogram takes the times and fluxes from your light curve, and returns the strongest period in your light curve and makes a plot of the periodogram. Run this function on both the red and the blue data, and identify the period of the light curve. How close was your guess from earlier?



In [23]:

    
from gatspy.periodic import LombScargle
import warnings

def lomb_scargle_periodogram(times, fluxes):
    """
    Calculate the best period for a light curve using
    a Lomb-Scargle periodogram, and plot the periodogram.
    """
    # Use gatspy to calculate the LS-periodogram 
    model = LombScargle(fit_period=True)
    model.optimizer.period_range = (20, 30)
    model.optimizer.quiet = True
    model.fit(times, fluxes)
    
    periods = np.linspace(10, 40, 10000)
    
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        scores = model.score(periods)
        
    # Plot the periodogram
    plt.figure()
    plt.plot(periods, scores)
    plt.xlabel('Period [days]')
    plt.ylabel('Power')
    plt.title("Best period: {0:.4f} days".format(model.best_period))
    
    # Return the best period
    return model.best_period



In [ ]: