Measure the Milky Way with Stars

RR Lyrae Stars as Standard Candles

How can we measure astronomical distances? One technique is to use "standard candles", which are astronomical objects that have well-known brightnesses. By measuring how bright an object appears to us and knowing the actual brightness from some analysis, we can work out how far the object is from us. We employ the "inverse square law" for light to make this happen and do a bit of math.

Look for the variable stars. They seem to blink.

Need more on how we use variable stars to measure distance?

If you are new to the idea of variable stars as tools for measuring distances, start with these:

Inverse Square Law

When light leaves the surface of an astronomical object, like a star's photosphere, the resulting light spreads out radially in all directions at the speed of light. All the energy emitted at that same instant can be thought of as being spread out over surface of a sphere. As this imaginary sphere gets bigger, the light gets more spread out. This is why distance matters when looking at a light source. Note the equation for the surface area of a sphere: \begin{equation*} A=4\pi{R^2} \end{equation*} The closer you are, the more concentrated the light is because the sphere over which it is spread is smaller (smaller radius, R). The farther you get from the source, the bigger the sphere must be when it reaches you (larger radius, R) and consequently, the light is spread out more and so the source seems dimmer.

Illustration of inverse square law.

The physical properties of electromagnetic radiation eminating from a star or other astronomical object is best described using flux which is related to power which is related to energy. Another common way to describe light leaving an astronomical object is through the concept of luminosity. The Stefan-Boltzman relation connections the energy density of a star (how much energy leaves 1 square meter of the star's surfce) and the star's temperature.

\begin{equation*} L=\sigma{T^4}\>[\mathbf{W/m^2}] \end{equation*}

We can find the power emitted by the star if we include the entire spherical surface of the star

\begin{equation*} L=\sigma{T^4}(4\pi{R^2})\>[\mathbf{W}] \end{equation*}

Don't forget that astronomers can determine the surface temperature of a star using Wien's law. Once the peak wavelength of the light from the star is found, we can get the temperature using the relation:

\begin{equation*} \lambda_{peak}\cdot{T}=constant\>[\mathbf{m\cdot{K}}] \end{equation*}

Magnitude vs Luminosity

For historical and convenience reasons, astronomers use a system called magnitude dating back to Hipparcos which uses the star Vega as a reference star. There are many peculiarlities with magnitude. For example a negative magnitude represents a brighter object than does a positive magnitude. Another useful point is that just using light from an image, astronomers can quickly measure the magnitude as seen from earth, the apparent magnitude. That leaves some work to do to get turn the physical brightness or luminosity of the object into the absolute magnitude. The magnitude scale is clunky but still in use - for now.

Globular Clusters

A globular cluster is a very old grouping of stars usually found in the halo or the central region of a galaxy. NGC 3201 is one such globular cluster found in the Milky Way Galaxy. The images used for this activity were taken using the Skynet Robitic Telescope Network run by the University of North Carolina. Yes... they called a robotic telescope system Skynet. The globular cluster was imaged 5 times in a single night by the system. This way, any star that has a periodic change of a day or so will show up in the series of images.

Globular Cluster NGC 3201

For this activity, we are going to study images taken of a variable star in a globular cluster over a single night. This star is known as an RR Lyrae variable star which varies over a period of less than one day. All stars that fall into this class have approximately the same absolute luminosity. That means if all RR Lyrae stars were the same distance from the observer, they would all appear the same brightness.

But stars are, of course, spread out across a galaxy and ours is no exception. So RR Lyrae stars appear to be different brightnesses due to their distances.

Choosing a target

We will study one such star in a globular cluster visible from Earth's southern hemisphere.

In order to determine the magnitude of the star, we will need a known calibration star that is in the field and near the target star. Note that here we have zoomed in a lot on our target star. At this level in this image the only things visible are stars (no galaxies or such) and almost every star visible is a part of the globular cluster.

Zooming in on the target.

Using the Aladin Sky Atlas, we can overlay data from previous surveys onto the image we took. Note the cirlcles and labels. Some objects have known proper motion values and some are known to be variable stars.

Image loaded into Aladin Sky Atlas with with Simbad source data overlaid.

Note the types of objects seen:

  • inCl* - source is a star and part of the cluster.
  • BlueStraggler - source is likely result of smaller, redder stars merging
  • RGB* - source is a star that has left the main sequence and moved to the red giant branch of the HR diagram
  • RRLyr - source is a variable star similar to RR Lyrae with a period less than one day and a known absolute magnitude
  • Star - source is a star and not necessarily a member of the cluster
  • X - source is an X-ray emitter such as a pulsar (neutron star), or black hole.

List of targets visible in the field above.

In [1]:
# If you are using Microsoft Azure Notebook, uncomment the line below by removing the #
# We need this line to include the non-standard SEP aperture photometry package.
#!pip install --no-deps sep
In [2]:
# Let's import the packages and libraries to help us do our analysis.
# The NumPy package is extremely useful and very common in scientific Python coding.
import numpy as np
# The SEP package is rather specific to our task of performing aperture photometry
import sep
# We import the math package to make our calculations easier to code.
import math
In [3]:
# additional setup for reading the test image and displaying plots
from astropy.utils.data import download_file
from astropy.io import fits

# The MatplotLib package is extremely useful and very common for making plots and displaying images.
import matplotlib.pyplot as plt
from matplotlib import rcParams

# Let's make Jupyter and Matplotlib work together
%matplotlib inline

# This determines the size of the images we display. 
rcParams['figure.figsize'] = [8., 8.]
In [4]:
# The known apparent magnitude of calibration star Cl* NGC 3201 CWFD 3-106 
# in the visual or green band according to the Simbad database.
# http://simbad.u-strasbg.fr/simbad/sim-id?Ident=Cl*+NGC+3201+CWFD+3-106
m_calibration = 14.81

# The known absolute magnitude of target star Cl* NGC 3201 LS 358
# https://iopscience.iop.org/article/10.1086/344948
M_target = 0.58
In [5]:
# Convenience function for mirroring a image about X and Y
# One of the FITS images was flipped horizontally and vertically and this function helps us fix that.
def mirror_data(data):
    dataout = np.array(data)
    rows = len(data)
    cols = len(data[0])
    for r in range(rows):
        for c in range(cols):
            dataout[rows-1-r][cols-1-c] = data[r][c]
    return dataout

Distance Modulus

If the apparent magnitude and the absolute magnitude of a source are known, once can determine the distance to the source in parsecs using the distance modulus relation. Small m is the apparent magnitude of the source. Big M is the absolute magnitude of the source. The distance is then d. Note the base of the logarithm is 10 in this case. The distance modulus is logarithmic in nature but it is also based on historical uses of magnitudes. Hence a distance modulus of 5 corresponds to a distance of 100 pc. The absolute magnitude scale is based on objects all being located 100 pc away from the viewer. Another example: a magnitude 1 star is exactly 100 times brighter than a magnitude 6 star. \begin{equation*} m - M = -5 + 5\cdot log_{10}(d) \end{equation*}

Coding and Questions

If you need to install Python and Jupyter notebook, the best bet is to install Anaconda.

  • Follow the directions and answer the questions listed in the comments of a given block of code.
  • Don't be afraid to experiment! This is an interactive lesson using coding. There is always the undo button!
  • The code used here is Python 3 which is common for scientific computing. This isn't the only language used in science by any means. It just works well for our purposes here.
  • Answer questions directly in the code block using comments. That will be part of your finished product!
  • The # means just that line is a comment
  • The """ ... """ means everthing between the triple quotes is a comment and can span multiple lines.
  • Include your name, the date, and class or period information. Be sure not to leave out any group members!
In [6]:
# Example of a single line comment
""" 
    Example of multi-line comment
"""
Out[6]:
' \n    Example of multi-line comment\n'
In [7]:
""" 
    Function to find the distance of an astronomical object in parsecs 
    if given the apparent, m, and absolute, M, magnitudes of the object.

    Use some algebra and re-arrange the distance modulus equation to return 
    the distance in parsecs of the object
    Don't forget how exponents work: x^2 (x squared) would be written as x**2
"""
def distance_modulus(m,M):
    return 10**( (m-M)/5 +1)
    # Students see:
    # return ...
In [8]:
"""
    Testing the distance modulus: should be approximately 1.0 for the Sun
        m = -26.76 and M =4.81
    There are 648000 arcseconds in 180 degrees of arc. Why do we need that conversion here?
    Is your result close enough to be useful? How can you tell?
    Compare your result with another student.
"""
print(distance_modulus(-26.76,4.81)*648000/math.pi)
1.0009793845988386
In [9]:
""" 
    Function to find the apparent magnitude of an object if given the
    flux (or counts) found via aperture photometry.
    http://classic.sdss.org/dr7/algorithms/fluxcal.html
    This function depends on the filters used and the mininum about of 
    flux the filter allows. Likely you shouldn't change anything here.
"""
def app_mag(flux):
    return -2.5*math.log10(flux/(25.11*10**8))
In [10]:
""" 
    Testing the apparent magnitude function: 
    should be 14.81 for our calibration star according to literature.
    https://iopscience.iop.org/article/10.1086/344948
"""
print(app_mag(2967.8))
14.81853020554703
In [11]:
"""
    This functions consumes as input the filename as a string and a "flag"
    for whether or not this image should be mirrored. To mirror an image 
    is to flip the image in both the horizontal and vertical direction.
    
    The data is returned as a 2-D NumPy array. This array represents the
    light received by each pixel in the camera.

"""
def get_image_data(filename, mirror=False):
    # process filename as FITS file
    file = filename
    image_file = download_file(file, cache=True )

    data = fits.getdata(image_file)
    
    if mirror==True:
        data = mirror_data(data)
    
    # We need this line to allow the SEP package to work with our 
    # AstroPy formatted FITS file
    # See: https://sep.readthedocs.io/en/v1.0.x/tutorial.html#Finally-a-brief-word-on-byte-order
    data = data.byteswap().newbyteorder()
    
    return data
In [12]:
"""
    This functions consumes as input a 2D NumPy array from a FITS image.
    The data returned is the data with the background subtracted and also
    background data to use in other functions.
    
    Removing the background improves the signal to noise ratio so our sources
    can be extracted more easily.
"""
def subtract_background(data):
      # measure a spatially varying background on the image
    bkg = sep.Background(data)
    
    # evaluate background as 2-d array, same size as original image
    bkg_image = bkg.back()
    # bkg_image = np.array(bkg) # equivalent to above
    
    # evaluate the background noise as 2-d array, same size as original image
    bkg_rms = bkg.rms()
    data_sub = data - bkg
    
    return data_sub, bkg

Why subtract the background?

If we have an astronomical image, the things we care about, like stars, might be hidden in the excess light and noise of our image. By using the SEP package to find and remove the background, we improve the signal to noice ratio. We want more signal and less noise. Finding ways to improve the S/N is a very common practice is scientific analysis. The 1st image has the background, the 2nd shows the background and the 3rd image is missing the background.

In [13]:
data = get_image_data('http://jimmynewland.com/astronomy/ngc-3201/ngc3201_1.fits')
data_sub, bkg = subtract_background(data)
bkg_image = bkg.back()
In [14]:
m1, s1 = np.mean(data), np.std(data)
plt.imshow(data, interpolation='nearest', cmap='gray',vmin=m1-s1, vmax=m1+s1, origin='lower')
Out[14]:
<matplotlib.image.AxesImage at 0x10cf03b00>
In [15]:
# show the background
plt.imshow(bkg_image, interpolation='nearest', cmap='gray', origin='lower')
Out[15]:
<matplotlib.image.AxesImage at 0x121fe4780>
In [16]:
m2, s2 = np.mean(data_sub), np.std(data_sub)
plt.imshow(data_sub, interpolation='nearest', cmap='gray',vmin=m2-s2, vmax=m2+s2, origin='lower')
Out[16]:
<matplotlib.image.AxesImage at 0x1234e3f28>
In [17]:
"""
    This function consumes the subtracted data, the background, and 
    the x and y coordinates for the subsection of the image containing the target star.
    There are 2 XY pairs for the "bounding box" we want within the original image.
    We then plot ellipses around each detected source and finally return the list
    of sources extracted.
"""
def extract_sources(data_sub, bkg, x1, y1, x2, y2):
    # perform the extraction
    objects = sep.extract(data_sub, 1.5, err=bkg.globalrms)
    
    # display the image
    from matplotlib.patches import Ellipse

    # plot background-subtracted image
    fig, ax = plt.subplots()
    m, s = np.mean(data_sub), np.std(data_sub)
    im = ax.imshow(data_sub, interpolation='nearest', cmap='gray',
               vmin=m-s, vmax=m+s, origin='lower')
    plt.xlim(x1,x2)
    plt.ylim(y1,y2)

    # plot an ellipse for each object
    for i in range(len(objects)):
        e = Ellipse(xy=(objects['x'][i], objects['y'][i]),
                width=6*objects['a'][i],
                height=6*objects['b'][i],
                angle=objects['theta'][i] * 180. / np.pi)
        e.set_facecolor('none')
        e.set_edgecolor('red')
        ax.add_artist(e)
    plt.show()
    
    return objects
In [18]:
"""
    This function helps us chose the center of the given region
    so we can find the target closest to the center.
    
    We consume the 2 XY pairs and return the ordered pair for the 
    center of the image as a tuple.
"""
def select_target(x1, y1, x2, y2):
    # select our target star by x and y position in the image
    cen_x = x1 + (x2-x1)/2
    cen_y = y1 + (y2-y1)/2
    
    return cen_x, cen_y
In [19]:
"""
    This function consumes our list of extracted targets along with 
    the data after subtraction, the background, and the center x and
    center y values.
    
    We then go through the list and find the flux for the target closest
    to the center (within 5).
    
    The flux for the target star is then returned.
"""
def get_target_flux(objects, data_sub, bkg, cen_x, cen_y):
    target_flux = 0
    
    for i in range(len(objects)):
        curr_x = objects[i]['x']
        curr_y = objects[i]['y']
        
        # Get the flux for target star (center)
        if curr_x > cen_x-5 and curr_x < cen_x+5 and curr_y > cen_y-5 and curr_y < cen_y+5:
            # find the flux of object nearest the center
            flux, fluxerr, flag = sep.sum_circle(data_sub, objects[i]['x'], objects[i]['y'],
                                     3.0, err=bkg.globalrms, gain=1.0)
            target_flux = round(int(flux))
    return target_flux
In [20]:
""" 
    This function consumes our list of extracted targets along with 
    the data after subtraction, the background, and the center x and
    center y values.
    
    We then go through the list and find the flux for the calibration star.
    Note: this will work ONLY for these images. Changing the calibration star
    would require a tedious by-hand analysis to locate the new position.
    
    The flux for the calibration star is returned.
"""
def get_calib_flux(objects, data_sub, bkg, cen_x, cen_y):
    for i in range(len(objects)):
        curr_x = objects[i]['x']
        curr_y = objects[i]['y']
            
        # Get flux of calibration star (below)
        if curr_x > cen_x-5 and curr_x < cen_x+5 and curr_y > cen_y-81-15 and curr_y < cen_y-81+15:
            # find the flux of object nearest the center
            flux, fluxerr, flag = sep.sum_circle(data_sub, objects[i]['x'], objects[i]['y'],
                                     3.0, err=bkg.globalrms, gain=1.0)
            calibration_flux = round(int(flux))
    return calibration_flux
In [21]:
"""
    This functions consumes:
        a filename which can be local or a URL for a FITS file, 
        the starting and ending points marking the bounding box for the image, 
        and whether or not the image needs to be flipped (both vertical and horizontal)
    This function returns:
        the apparent magnitude of the target star
    The process image function is where we bring the steps together.
    1) Get the image data using the filename and mirror flag.
    2) Get the subtracted data and background using the data from step 1
    3) Extract a list of sources using the subtracted data, the background
        and the given x and y pairs (2 x values and 2 y values).
    4) Set the target as the center of the image.
    5) Get the flux for our target star.
    6) Get the flux for out calibration star.
    7) Determine the apparent magnitude of the calibration star.
    8) Calibrate the apparent magnitude of the target star.
    9) Print and return our target star magnitude.
"""
def process_image(filename, x1=0, y1=0, x2=2048, y2=2048, mirror=False):
    # Get image data
    data = get_image_data(filename, mirror)
    
    # Measure and subtract background
    data_sub, bkg = subtract_background(data)
    
    # Extract sources
    objects = extract_sources(data_sub, bkg, x1, y1, x2, y2)
    
    # Select target (center)
    cen_x, cen_y = select_target(x1, y1, x2, y2)
    
    # Find target star flux
    target_flux = get_target_flux(objects, data_sub, bkg, cen_x, cen_y)
    
    # Find the calibration star flux
    calibration_flux = get_calib_flux(objects, data_sub, bkg, cen_x, cen_y)
    
    # Convert the flux into an apparent magnitude
    m_cal_ap = app_mag(calibration_flux)
    
    # Calibrate the magnitudes for the target and calibration star
    f_cal = m_cal_ap/m_calibration # Calibration factor
    m_cal = m_cal_ap/f_cal
    
    # Convert the flux into an apparent magnitude
    m_tar_ap = app_mag(target_flux)
    
    # Calibrate the magnitude of the target
    m_tar = m_tar_ap/f_cal
        
    print("App. mag of target star m = : "+str(m_tar))
    print("App. mag of calib. star m = : "+str(m_cal))
    
    # Return the magnitude for the distance calculation
    return m_tar
In [22]:
# Empty list meant to hold calculated magnitudes for target star returned from process_image function
target_mags = []
In [23]:
"""
    Run process_image using the following parameters:
    filename='http://jimmynewland.com/astronomy/ngc-3201/ngc3201_1.fits'
    x1=860
    x2=1060
    y1=668
    y2=868
"""
m_out = process_image('http://jimmynewland.com/astronomy/ngc-3201/ngc3201_1.fits',x1=860,x2=1060,y1=668,y2=868)
# Students see
# m_out = process_image(..., ..., ..., ..., ...)
target_mags.append(m_out)
App. mag of target star m = : 14.703044003350854
App. mag of calib. star m = : 14.809999999999999
In [24]:
"""
    Run process_image using the following parameters:
    filename='http://jimmynewland.com/astronomy/ngc-3201/ngc3201_2.fits'
    x1=869
    x2=1069
    y1=673
    y2=873
"""
m_out = process_image('http://jimmynewland.com/astronomy/ngc-3201/ngc3201_2.fits',x1=869,x2=1069,y1=673,y2=873)
# Students see
# m_out = process_image(..., ..., ..., ..., ...)
target_mags.append(m_out)
App. mag of target star m = : 15.140060378737928
App. mag of calib. star m = : 14.81
In [25]:
"""
    Run process_image using the following parameters:
    filename='http://jimmynewland.com/astronomy/ngc-3201/ngc3201_3.fits'
    x1=205
    x2=435
    y1=190
    y2=420
"""
m_out = process_image('http://jimmynewland.com/astronomy/ngc-3201/ngc3201_3.fits',x1=205,x2=435,y1=190,y2=420)
# Students see
# m_out = process_image(..., ..., ..., ..., ...)
target_mags.append(m_out)