Interstellar Dust Extinction

dust_extinction is a python package to provide interstellar dust extinction curves.

While there are other python packages that provide some of the extinction curves contained here, the explicit motivation for this package is to provide extinction curves to those using them to model/correct their data and those studying extinction curves directly to better undertand interstellar dust.

This package is developed in the astropy affiliated package template and uses the astropy.modeling framework.

Installation

To be added

User Documenation

Model Flavors

There are three differnet types of models: average, R(V)+ dependent prediction, and shape fitting.

Average models

These models provide averages from the literature with the ability to interpolate between the observed data points. For the Milky Way, one of the R(V) dependent models with R(V) = 3.1 (see next section) are often used for the Milky Way ‘average’. Models are provided for the Magellanic Clouds from Gordon et al. (2003). Models for the Milky Way still to be added (both UV/optical/NIR and IR).

(Source code, png, hires.png, pdf)

_images/model_flavors-1.png

R(V) (+ other variables) dependent prediction models

These models provide predictions of the shape of the dust extinction given input variable(s).

These include CCM89 the original R(V) dependent model from Cardelli, Clayton, and Mathis (1989) and updated versions F99 (Fitzpatrick 1999). These models are based on the average behavior of extinction in the Milky Way.

In addition, the (R(V), f_A) two parameter relationship from Gordon et al. (2016) is included. This model is based on the average behavior of extinction in the Milky Way, Large Magellanic Cloud, and Small Magellanic Cloud.

(Source code, png, hires.png, pdf)

_images/model_flavors-2.png

(Source code, png, hires.png, pdf)

_images/model_flavors-3.png

(Source code, png, hires.png, pdf)

_images/model_flavors-4.png

(Source code, png, hires.png, pdf)

_images/model_flavors-5.png

Shape fitting models

These models are used to fit the detailed shape of dust extinction curves. The FM90 (Fitzpatrick & Mass 1990) model uses 6 parameters to fit the shape of the ultraviolet extinction. The P92 (Pei 1992) uses 19 parameters to fit the shape of the X-ray to far-infrared extinction.

(Source code, png, hires.png, pdf)

_images/model_flavors-6.png

(Source code, png, hires.png, pdf)

_images/model_flavors-7.png

Extinguish or Unextinguish Data

Two of the three flavors of models include a function to calculate the factor to multiple (extinguish) or divide (unextinguish) a spectrum by to add or remove the effects of dust, respectively.

Extinguish is also often called reddening. Extinguishing a spectrum often reddens the flux, but sometimes ‘bluens’ the flux (e.g, on the short wavelength side of the 2175 A bump). So extinguish is the more generic term.

Example: Extinguish a Blackbody

import matplotlib.pyplot as plt
import numpy as np

import astropy.units as u
from astropy.modeling.blackbody import blackbody_lambda

from dust_extinction.dust_extinction import CCM89

# generate wavelengths between 0.1 and 3 microns
#    within the valid range for the CCM89 R(V) dependent relationship
lam = np.logspace(np.log10(0.1), np.log10(3.0), num=1000)

# setup the inputs for the blackbody function
wavelengths = lam*1e4*u.AA
temperature = 10000*u.K

# get the blackbody flux
flux = blackbody_lambda(wavelengths, temperature)

# initialize the model
ext = CCM89(Rv=3.1)

# get the extinguished blackbody flux for different amounts of dust
flux_ext_av05 = flux*ext.extinguish(wavelengths, Av=0.5)
flux_ext_av15 = flux*ext.extinguish(wavelengths, Av=1.5)
flux_ext_ebv10 = flux*ext.extinguish(wavelengths, Ebv=1.0)

# plot the intrinsic and extinguished fluxes
fig, ax = plt.subplots()

ax.plot(wavelengths, flux, label='Intrinsic')
ax.plot(wavelengths, flux_ext_av05, label='$A(V) = 0.5$')
ax.plot(wavelengths, flux_ext_av15, label='$A(V) = 1.5$')
ax.plot(wavelengths, flux_ext_ebv10, label='$E(B-V) = 1.0$')

ax.set_xlabel('$\lambda$ [$\AA$]')
ax.set_ylabel('$Flux$')

ax.set_xscale('log')
ax.set_yscale('log')

ax.set_title('Example extinguishing a blackbody')

ax.legend(loc='best')
plt.tight_layout()
plt.show()

(Source code, png, hires.png, pdf)

_images/extinguish-1.png

Fit Extinction Curves

The dust_extinction package is built on the astropy.modeling package. Fitting is done in the standard way for this package where the model is initialized with a starting point (either the default or user input), the fitter is choosen, and the fit performed.

Example: FM90 Fit

In this example, the FM90 model is used to fit the observed average extinction curve for the LMC outside of the LMC2 supershell region (G03_LMCAvg dust_extinction model).

import matplotlib.pyplot as plt
import numpy as np

from astropy.modeling.fitting import LevMarLSQFitter

from dust_extinction.dust_extinction import G03_LMCAvg, FM90

# get an observed extinction curve to fit
g03_model = G03_LMCAvg()

x = g03_model.obsdata_x
# convert to E(x-V)/E(B0V)
y = (g03_model.obsdata_axav - 1.0)*g03_model.Rv
# only fit the UV portion (FM90 only valid in UV)
gindxs, = np.where(x > 3.125)

# initialize the model
fm90_init = FM90()

# pick the fitter
fit = LevMarLSQFitter()

# fit the data to the FM90 model using the fitter
#   use the initialized model as the starting point
g03_fit = fit(fm90_init, x[gindxs], y[gindxs])

# plot the observed data, initial guess, and final fit
fig, ax = plt.subplots()

ax.plot(x, y, 'ko', label='Observed Curve')
ax.plot(x[gindxs], fm90_init(x[gindxs]), label='Initial guess')
ax.plot(x[gindxs], g03_fit(x[gindxs]), label='Fitted model')

ax.set_xlabel('$x$ [$\mu m^{-1}$]')
ax.set_ylabel('$E(x-V)/E(B-V)$')

ax.set_title('Example FM90 Fit to G03_LMCAvg curve')

ax.legend(loc='best')
plt.tight_layout()
plt.show()

(Source code, png, hires.png, pdf)

_images/fit_extinction-1.png

Example: P92 Fit

In this example, the P92 model is used to fit the observed average extinction curve for the MW as tabulted by Pei (1992).

import matplotlib.pyplot as plt
import numpy as np

from astropy.modeling.fitting import LevMarLSQFitter

from dust_extinction.dust_extinction import P92

# Milky Way observed extinction as tabulated by Pei (1992)
MW_x = [0.21, 0.29, 0.45, 0.61, 0.80, 1.11, 1.43, 1.82,
        2.27, 2.50, 2.91, 3.65, 4.00, 4.17, 4.35, 4.57, 4.76,
        5.00, 5.26, 5.56, 5.88, 6.25, 6.71, 7.18, 7.60,
        8.00, 8.50, 9.00, 9.50, 10.00]
MW_x = np.array(MW_x)
MW_exvebv = [-3.02, -2.91, -2.76, -2.58, -2.23, -1.60, -0.78, 0.00,
             1.00, 1.30, 1.80, 3.10, 4.19, 4.90, 5.77, 6.57, 6.23,
             5.52, 4.90, 4.65, 4.60, 4.73, 4.99, 5.36, 5.91,
             6.55, 7.45, 8.45, 9.80, 11.30]
MW_exvebv = np.array(MW_exvebv)
Rv = 3.08
MW_axav = MW_exvebv/Rv + 1.0

# get an observed extinction curve to fit
x = MW_x
y = MW_axav

# initialize the model
p92_init = P92()

# pick the fitter
fit = LevMarLSQFitter()

# fit the data to the P92 model using the fitter
#   use the initialized model as the starting point
#   accuracy set to avoid warning the fit may have failed
p92_fit = fit(p92_init, x, y, acc=1e-3)

# plot the observed data, initial guess, and final fit
fig, ax = plt.subplots()

ax.plot(x, y, 'ko', label='Observed Curve')
ax.plot(x, p92_init(x), label='Initial guess')
ax.plot(x, p92_fit(x), label='Fitted model')

ax.set_xlabel('$x$ [$\mu m^{-1}$]')
ax.set_ylabel('$A(x)/A(V)$')

ax.set_title('Example P92 Fit to MW average curve')

ax.legend(loc='best')
plt.tight_layout()
plt.show()

(Source code, png, hires.png, pdf)

_images/fit_extinction-2.png

Reporting Issues

If you have found a bug in dust_extinction please report it by creating a new issue on the dust_extinction GitHub issue tracker.

Please include an example that demonstrates the issue sufficiently so that the developers can reproduce and fix the problem. You may also be asked to provide information about your operating system and a full Python stack trace. The developers will walk you through obtaining a stack trace if it is necessary.

Contributing

Like the Astropy project, dust_extinction is made both by and for its users. We accept contributions at all levels, spanning the gamut from fixing a typo in the documentation to developing a major new feature. We welcome contributors who will abide by the Python Software Foundation Code of Conduct.

dust_extinction follows the same workflow and coding guidelines as Astropy. The following pages will help you get started with contributing fixes, code, or documentation (no git or GitHub experience necessary):

For the complete list of contributors please see the dust_extinction contributors page on Github.

Reference API

dust_extinction.dust_extinction Module

Classes

BaseExtModel(*args, **kwargs) Base Extinction Model.
BaseExtRvModel(*args, **kwargs) Base Extinction R(V)-dependent Model.
BaseExtAve(*args, **kwargs) Base Extinction Average.
CCM89([Rv]) CCM89 extinction model calculation
FM90([C1, C2, C3, C4, xo, gamma]) FM90 extinction model calculation
P92([BKG_amp, BKG_lambda, BKG_b, BKG_n, …]) P92 extinction model calculation
F99([Rv]) F99 extinction model calculation
G03_SMCBar(*args, **kwargs) G03 SMCBar Average Extinction Curve
G03_LMCAvg(*args, **kwargs) G03 LMCAvg Average Extinction Curve
G03_LMC2(*args, **kwargs) G03 LMC2 Average Extinction Curve
G16([RvA, fA]) G16 extinction model calculation

Class Inheritance Diagram

Inheritance diagram of dust_extinction.dust_extinction.BaseExtModel, dust_extinction.dust_extinction.BaseExtRvModel, dust_extinction.dust_extinction.BaseExtAve, dust_extinction.dust_extinction.CCM89, dust_extinction.dust_extinction.FM90, dust_extinction.dust_extinction.P92, dust_extinction.dust_extinction.F99, dust_extinction.dust_extinction.G03_SMCBar, dust_extinction.dust_extinction.G03_LMCAvg, dust_extinction.dust_extinction.G03_LMC2, dust_extinction.dust_extinction.G16