Exploring UV Extinction Curves#

Learning Goals#

By the end of this tutorial, you will be able to:

  • Perform data queries on the MAST archives using astroquery

  • Narrow down query results to the desired spectrum

  • Clean, analyze, and plot spectral data

  • Fit extinction curves to a parameterized model


In this tutorial, you will learn how to download IUE data using the MAST Archive API, astroquery. We’ll obtain the spectra of two stars, one reddened and one not. We will use these spectra to construct and plot UV extinction curves. An extinction curve demonstrates the wavelength dependence of dust extinction. It compares the dust-free Spectral Energy Distribution (SED) of a star to the observed SED. Normally these curves are created by plotting $k(\lambda-V)$ versus $1/\lambda$, with $\lambda$ being the wavelength (see ‘Useful Equations’ below).

Extinction is relevant in a diverse range of scenarios. Sometimes, dust is found very near the observed object, like in stars with disks or proto-stellar clouds. However, the target does not need to be intrinsically “dusty”; dust in a galaxy that is light-years from a target star might still block the line-of-sight view and impact the observation.

In this lesson, our data comes from the International Ultraviolet Explorer (IUE). IUE performed spectrophotometry in the wavelength range from 1150 Å to 3200 Å. MAST hosts these observations, which include more than 100,000 spectra.

Defining some terms:

  • Color index: difference between magnitude of a star in two different passbands, typically B and V. Symbol: $(B-V)$.

  • Extinction: absoption and scattering of light by dust and gas between an object and the observer. It is a measure of the interstellar reddening quantified by the difference between the magnitude of the band when observed through dust versus a dust-free environment. Symbol: $A(\lambda)$.

  • Magellanic Clouds: Irregular dwarf galaxies that orbit the Milky Way galaxy and are located in the southern celestial hemisphere. Two distinct groups can be differentiated: the Large Magellanic Cloud (LMC) and the Small Magellanic Cloud (SMC).

  • Spectral type: stellar classification from hotter (O stars) to cooler (M stars). Temperature defines a star’s “color” and surface brightness.

The “extinction equation” we’ll use to create our plots:

  • $k(\lambda-V) = \frac{m(\lambda-V)-m(\lambda-V)_o}{(B-V)-(B-V)_o} = \frac{A(\lambda)-A(V)}{A(B)-A(V)}$

_Note: the $x_o$ terms refer to the star that is nearly unaffected by dust, i.e. $(B-V)$ corresponds to the observed color index and $(B-V)o$ to the color index if there were no extinction due to dust. The stars should have the same spectral type in order to perform this comparison.

Table of Contents#

  • Imports

  • Searching by Target Name and Downloading Spectra

  • Exploring the downloaded data files

  • Data access and manipulation

  • Saving values into your local computer

  • Exercises

  • Additional Resources


The first step will be to import the libraries we will be using throughout this tutorial:

  • Observations from astroquery.mast to query the MAST Archive.

  • fits from astropy.io for accessing FITS files

  • Simbad from astroquery.simbad to query the SIMBAD database.

  • numpy for array manipulations

  • matplotlib for plotting

from astroquery.mast import Observations
from astropy.io import fits
from astroquery.simbad import Simbad

import numpy as np

import matplotlib.pyplot as plt

Warning: If you have not installed the astroquery package on your computer, you will need to. Information about astroquery can be found on the readthedocs.

Background and Targets#

The first step is to find the data files we want to use. We will using the Observations Class from astroquery.mast both to search for and download UV data from the IUE.

We’ll be targeting two stars: AzV 456 and AzV 70. This former is reddened by dust, while the latter is not. This particular example is from Gordon et al. (2003), but you can use any pair of reddened and unreddened stars with the same spectral type and luminosity class.

Searching and Downloading in astroquery.mast#

We’ll make our targets flexible here by creating a “dusty” and “no dust” option. If you want to explore other stars, change the targets below!

target_dusty = "Azv456"
target_nodust = "Azv70"

The query_criteria function allows you to filter by mission, as well as object name. For a full list of fields you can query, as well as their associated metadata names, see the CAOM fields description page. For this search, we’ll filter for IUE data taken on the “no dust” target.

Note: We’ll use a radius of zero to exactly match our target; depending on your particular query, this may result in some missing observations. Use with caution!

# We use a radius of zero to find exact matches
obs_table_nodust = Observations.query_criteria(obs_collection="IUE",
                                           radius = "0m")
# Display the matching observations, if you want
Table length=2
0scienceIUE--LWP--LOW DISPUVSK 35----lwp1238712.577660313500019-72.6341409203OBJEFFitzpatrickspectrum247157.5873747157.604031438.556185118000000.0334760000000.0Energy Distributions of B Supergiants in the Small Magellanic Cloudnan--PUBLIC--2822185090965090960.0
1scienceIUE--SWP--LOW DISPUVSK 35----swp1883012.577660313500019-72.6341409203OD90BWalbornspectrum245323.6822645323.703081798.573115059000000.0197870000000.0SNC OB Supergiantsnan--PUBLIC--3159885428665428660.0

There are two matching observations. One covers a short wavelength domain and the other covers the long domain (short and long are relative: both are UV observations). Let’s use both so that we’ll have a more complete picture of the extinction behavior over different wavelengths.

Let’s find the data products associated with these observations. We’ll use the get_products function to do that.

data_products_nodust = Observations.get_product_list(obs_table_nodust)
# Display the products, if you want
Table length=18
7282218IUEspectrumlwp12387MXLOSmast:IUE/url/pub/iue/data/lwp/12000/lwp12387.mxlo.gzSCIENCEMinimum Recommended Products--------OBJEFlwp12387.mxlo.gz18343282218PUBLIC2
8282218IUEspectrumlwp12387(extracted spectra/vo spectral container/SSAP) fits fileSmast:IUE/url/pub/vospectra/iue2/lwp12387mxlo_vo.fitsSCIENCEMinimum Recommended Products(extracted spectra/vo spectral container/SSAP) fits file------OBJEFlwp12387mxlo_vo.fits51840282218PUBLIC2
16315988IUEspectrumswp18830MXLOSmast:IUE/url/pub/iue/data/swp/18000/swp18830.mxlo.gzSCIENCEMinimum Recommended Products--------OD90Bswp18830.mxlo.gz15975315988PUBLIC2
17315988IUEspectrumswp18830(extracted spectra/vo spectral container/SSAP) fits fileSmast:IUE/url/pub/vospectra/iue2/swp18830mxlo_vo.fitsSCIENCEMinimum Recommended Products(extracted spectra/vo spectral container/SSAP) fits file------OD90Bswp18830mxlo_vo.fits51840315988PUBLIC2

This looks good, but a fair number of these results have the productType ‘AUXILIARY’; they were orginially intended for calibration or engineering purposes. We want only the the products marked ‘SCIENCE’; more specifically, we’re looking for extracted spectra in .fits files. Let’s add some filters to match those results.

filtered_products_nodust = Observations.filter_products(data_products_nodust,
                                                     extension = '.fits')
# Display the results, if you want
Table length=2
0282218IUEspectrumlwp12387(extracted spectra/vo spectral container/SSAP) fits fileSmast:IUE/url/pub/vospectra/iue2/lwp12387mxlo_vo.fitsSCIENCEMinimum Recommended Products(extracted spectra/vo spectral container/SSAP) fits file------OBJEFlwp12387mxlo_vo.fits51840282218PUBLIC2
1315988IUEspectrumswp18830(extracted spectra/vo spectral container/SSAP) fits fileSmast:IUE/url/pub/vospectra/iue2/swp18830mxlo_vo.fitsSCIENCEMinimum Recommended Products(extracted spectra/vo spectral container/SSAP) fits file------OD90Bswp18830mxlo_vo.fits51840315988PUBLIC2

Great! All that remains is to download the files. You can do this by passing the entire product table to the download_products function.

manifest_nodust = Observations.download_products(filtered_products_nodust)
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:IUE/url/pub/vospectra/iue2/lwp12387mxlo_vo.fits to ./mastDownload/IUE/lwp12387/lwp12387mxlo_vo.fits ... [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:IUE/url/pub/vospectra/iue2/swp18830mxlo_vo.fits to ./mastDownload/IUE/swp18830/swp18830mxlo_vo.fits ... [Done]

Rinse and Repeat for Target Two#

Now, let’s repeat this process for our second, dusty target. We’ll condense the code into one cell for convenience.

obs_table_dusty = Observations.query_criteria(obs_collection="IUE",
                                            radius = "0m",
                                            proposal_pi="Prevot") # We specify the PI since this search returns more results

data_products_dusty = Observations.get_product_list(obs_table_dusty)

# Note that you can skip the 'filter products' step
# Instead, you can pass the filters directly to 'download_products'
manifest_dusty = Observations.download_products(data_products_dusty,
                                             extension = ".fits")
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:IUE/url/pub/vospectra/iue2/lwr12347mxlo_vo.fits to ./mastDownload/IUE/lwr12347/lwr12347mxlo_vo.fits ... [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:IUE/url/pub/vospectra/iue2/swp16051mxlo_vo.fits to ./mastDownload/IUE/swp16051/swp16051mxlo_vo.fits ... [Done]

Exploring the Downloaded Data Files#

Now, let’s explore the FITS files we just downloaded. The filepaths were saved in our manifest variable when we used download_products.

# Pull the filepaths from the manifests
filenames_dusty = manifest_dusty['Local Path']
filenames_nodust = manifest_nodust['Local Path']

# The first path in the list is for the 'long wavelength' data
lw_dusty = filenames_dusty[0]

# We'll create the other files now, for convenience
sw_dusty = filenames_dusty[1]
lw_nodust = filenames_nodust[0]
sw_nodust = filenames_nodust[1]

# Print the information from the first file only
Filename: ./mastDownload/IUE/lwr12347/lwr12347mxlo_vo.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     359   ()      
  1  Spectral Container    1 BinTableHDU    141   1R x 4C   [562E, 562E, 562E, 562I]   
  • No. 0 (PRIMARY): This HDU contains metadata related to the entire file.

  • No. 1 (Spectral Container): This HDU contains the spectral profile of the target as a function of wavelength.

The header of the file contains additional information about the data.

Let’s take a look at some of the columns that are more relevant to us:

with fits.open(lw_dusty) as hdulist: 
    header_lw_dusty = hdulist[1].header

# Remove to the [9:15] index to see the entire header
COMMENT  *** Column names ***                                                   
TTYPE1  = 'WAVE    '           /                                                
TTYPE2  = 'FLUX    '           /                                                
TTYPE3  = 'SIGMA   '           /                                                
TTYPE4  = 'QUALITY '           /                                                

This tell us that first column corresponds to the wavelengths, and the second column to the fluxes. We can also see the units:

COMMENT   *** Column units ***                                                  
TUNIT1  = 'angstrom'           / wavelength unit is Angstrom                    
TUNIT2  = 'erg/cm**2/s/angstrom' / flux unit is ergs/cm2/sec/A                  
TUNIT3  = 'erg/cm**2/s/angstrom' / sigma unit is ergs/cm2/sec/A                 

Data Access and Validation#

The data contained in this fits file can be accessed using io.fits and the .data attribute. This can be tedious to do more than once, so let’s create a helper function to do it for us.

# This helper function will take a filename and return the wavelength/flux data
def extractData(filename):
    # Open the file and read the data
    with fits.open(filename) as hdulist:
        spectrum = hdulist[1].data
    # Extract our desired data from the corresponding columns
    wav = spectrum[0][0]   # wavelength, angstrom, A
    flux = spectrum [0][1] # flux, ergs/cm2/sec/A

    return wav, flux

Let’s make sure we get a sensible result when we extract the data. We can make a quick plot to see if anything is wrong with the results.

# Let's run our helper function on our data
wav_lw_dusty, flux_lw_dusty = extractData(lw_dusty)
wav_sw_dusty, flux_sw_dusty = extractData(sw_dusty)
wav_lw_nodust, flux_lw_nodust = extractData(lw_nodust)
wav_sw_nodust, flux_sw_nodust = extractData(sw_nodust)

# Let's also make a rough plot of our data, just so we can see what we're working with
plt.plot(wav_lw_nodust, np.log10(flux_lw_nodust))
plt.plot(wav_lw_dusty, np.log10(flux_lw_dusty))
[<matplotlib.lines.Line2D at 0x7fb61fa404d0>]

This looks good on inspection, but we should check that both targets have consistent data. We’d like to compare one star to another and our later analysis requires that corresponding datasets have the same number of points.

# Inspect the lengths by eye
print(len(wav_lw_nodust), len(flux_lw_nodust), len(wav_lw_dusty), len(flux_lw_dusty)) # check long wavelength data
print(len(wav_sw_nodust), len(flux_sw_nodust), len(wav_sw_dusty), len(flux_sw_dusty)) # check short
563 563 562 562
495 495 495 495

Whoops! Looks like the the long wavelength data is not consistent. Let’s adjust and check again.

# Shorten the lw_nodust data so it matches the dusty data
# This may not work for your targets! Always check your data
wav_lw_nodust = wav_lw_nodust[:len(wav_lw_dusty)]
flux_lw_nodust = flux_lw_nodust[:len(flux_lw_dusty)]
print(len(wav_lw_nodust), len(flux_lw_nodust), len(wav_lw_dusty), len(flux_lw_dusty))
print(len(wav_sw_nodust), len(flux_sw_nodust), len(wav_sw_dusty), len(flux_sw_dusty))
562 562 562 562
495 495 495 495

Data Plotting#

After searching for, downloading, and cleaning our data, we’re finally ready to display the plot. We’ll plot the flux against the inverse wavelength; this is customary for this type of study and follows the example in Gordon et al. (2003).

# Set up the plot size, labels, units
fig = plt.figure()
ax = plt.subplot(111)
ax.set_xlabel('1/$\lambda$ ($\mu m^{-1}$)')
ax.set_ylabel(r'log(Flux (ergs $cm^{-2}$ $s^{-1}$ $\AA^{-1}$))')
ax.set_ylim([-13.6, -11.5])

# Create a helper function to plot the data. We want to plot inverse wavelength against log of flux
def Plot(wav, flux, style, lbl):
    plt.plot(10**4/wav, np.log10(flux), style, label=lbl)

# Call the helper function four times: for long/short wavelengths of both targets
Plot(wav_lw_dusty, flux_lw_dusty, 'r:', 'lw_dusty')
Plot(wav_sw_dusty, flux_sw_dusty, 'r', 'sw_dusty')
Plot(wav_lw_nodust, flux_lw_nodust, 'b:', 'lw_nodust')
Plot(wav_sw_nodust, flux_sw_nodust, 'b', 'sw_nodust')

# Add some labels to our plot to clarify which spectrum is which
plt.text(6.7, -13.4, target_dusty, fontsize = 11, color='r')
plt.text(6.1, -12.6, target_nodust, fontsize = 11, color='b')

# Create the legend and show our plot

What a wonderful example of extinction! Both stars have the same spectral features, but we can see a significant dimming of AzV456, especially in the shorter, bluer wavelengths; it is safe to say this star is reddened.

Extinction Curve#

Let’s now use the SIMBAD database to look for the fluxes in the B and V bands for both stars. We can do a simple query using the identifier of both stars. The magnitudes can be found under the 8th subgroup presented below the name of the stars, called ‘Fluxes’, since SIMBAD can provide you with either the flux or the magnitude of the star in those bands.

# Request the B and V magnitudes

# Query for dusty star, in our case AvZ 456
table_dusty = Simbad.query_object(target_dusty)
# Examine the data, if you want
Table length=1
0SK 14301 10 55.7567-72 42 56.22314140.4390.42590AO2020yCat.1350....0G12.9812.891

Great. Let’s get the data for AzV70 as well.

table_nodust = Simbad.query_object(target_nodust)

From these values we can directly calculate $E(B-V) = (B-V)-(B-V)_o$. We can do this since AvZ 70 is essentially an unreddened version of AvZ 456.

# Extract these values from the table
V_dusty = table_dusty['FLUX_V'][0]
B_dusty = table_dusty['FLUX_B'][0]

V_nodust = table_nodust['FLUX_V'][0]
B_nodust = table_nodust['FLUX_B'][0]

# Plug into the formula
EBV = (B_dusty-V_dusty)-(B_nodust-V_nodust)
print(f"The value of E(B-V) is equal to {EBV:.3f}")
The value of E(B-V) is equal to 1.530

Now we can fully plot the extinction curve! Recall the equation from the start of this Notebook.

$$k(\lambda-V) = \frac{m(\lambda-V)-m(\lambda-V)_o}{(B-V)-(B-V)_o} = \frac{A(\lambda)-A(V)}{E(B-V)} $$

Let’s expand this out and be explicit: which values do we need to use to calculate extinction?

$$\frac{A(\lambda)-A(V)}{E(B-V)}=\frac{log(V_ratio)-log(flux_ratio)}{E(B-V)} = \frac{log (V_ratio \div flux_ratio)}{E(B-V)}$$

V_rat = V_dusty/V_nodust

# Another helper function to create our plots, and save the output data for further analysis
def Extinction(wav, flux1, flux2):
    flux_rat = flux1/flux2
    # Calcuating extinction from the formula 
    ext = np.log(V_rat/flux_rat)/EBV
    inv_wav= 10**4/wav
    plt.plot(inv_wav, ext, 'o', markersize=2)
    return inv_wav, ext 

# Create the figure, then use the helper function to create the plots and save the data
s_inv_wav, s_ext = Extinction(wav_sw_dusty, flux_sw_dusty, flux_sw_nodust)
l_inv_wav, l_ext = Extinction(wav_lw_dusty, flux_lw_dusty, flux_lw_nodust)

# Add a title, axes labels, and a legend
plt.legend(['short wavelength', 'long wavelength'])
plt.xlabel(r'$1/\lambda$ $[\mu m^{-1}]$')

The values at the edge of the detectors don’t look quite right. Particularly the right edge of the short wavelength; would we really expect the reddened star to be brighter in the short wavelengths region Let’s excise the data from the edges and try again.

# How many points do we want to remove from either end of the data?
s_crop = 50
l_crop = -50

# Use helper function with cropped data
s_inv_wav, s_ext = Extinction(wav_sw_dusty[s_crop:], flux_sw_dusty[s_crop:], flux_sw_nodust[s_crop:])
l_inv_wav, l_ext = Extinction(wav_lw_dusty[:l_crop], flux_lw_dusty[:l_crop], flux_lw_nodust[:l_crop])

# Add a title, axes labels, and a legend
plt.legend(['short wavelength', 'long wavelength'])
plt.xlabel(r'$1/\lambda$ $[\mu m^{-1}]$')

Much better.

This is the typical shape encountered for extinction curves corresponding to the Small Magellanic Cloud (SMC), as can be seen in Gordon et al. (2003). Note in particular the “bump” at the transition between long and short wavelengths. This is the mysterious 2175 $\mathring A$ bump; we don’t fully understand its physical origin!


If you’re interested in exploring further, go back to the beginning of this Notebook and select a new pair of stars (Gordon et al. have some examples!). Can you find a star that doesn’t have the 2175 $\mathring A$ bump?

For convenience, you might want to add a more sophisticated method of error checking. Recall that above, we printed out the lengths of the data to ensure they matched. Can you think of a way to check for this automatically and adjust as needed?

Additional Resources#

For more information about the MAST archive and details about mission data:

For more information about extinction curves and their parametrization:

About this Notebook#

Authors: Clara Puerto Sánchez, Claire Murray, Thomas Dutkiewicz
Keyword(s): UV, reddening, extinction-curve
Last Updated: Apr 2023

For support, please contact the Archive HelpDesk at archive@stsci.edu.

Top of Page Space Telescope Logo