Skip to content

Commit

Permalink
Merge 07a81f7 into 5f2c511
Browse files Browse the repository at this point in the history
  • Loading branch information
@karllark
karllark committed Apr 19, 2019
2 parents 5f2c511 + 07a81f7 commit c958d1b
Show file tree
Hide file tree
Showing 4 changed files with 220 additions and 44 deletions.
72 changes: 30 additions & 42 deletions docs/dust_extinction/model_flavors.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,8 @@ R(V) (+ other variables) dependent prediction models
import matplotlib.pyplot as plt
import astropy.units as u

from dust_extinction.parameter_averages import (CCM89, O94, F99, F04, M14)
from dust_extinction.parameter_averages import (CCM89, O94, F99, F04,
GCC09, M14)

fig, ax = plt.subplots()

Expand All @@ -86,20 +87,15 @@ R(V) (+ other variables) dependent prediction models

Rv = 3.1

ext_model = CCM89(Rv=Rv)
ax.plot(x,ext_model(x),label='CCM89')
models = [CCM89, O94, F99, F04, GCC09, M14]

ext_model = O94(Rv=Rv)
ax.plot(x,ext_model(x),label='O94')

ext_model = F99(Rv=Rv)
ax.plot(x,ext_model(x),label='F99')

ext_model = F04(Rv=Rv)
ax.plot(x,ext_model(x),label='F04')

ext_model = M14(Rv=Rv)
ax.plot(x[x<3.3/u.micron],ext_model(x[x<3.3/u.micron]),label='M14')
for cmodel in models:
ext_model = cmodel(Rv=Rv)
indxs, = np.where(np.logical_and(
x.value >= ext_model.x_range[0],
x.value <= ext_model.x_range[1]))
yvals = ext_model(x[indxs])
ax.plot(x[indxs], yvals, label=ext_model.__class__.__name__)

ax.set_xlabel('$x$ [$\mu m^{-1}$]')
ax.set_ylabel('$A(x)/A(V)$')
Expand All @@ -116,7 +112,8 @@ R(V) (+ other variables) dependent prediction models
import matplotlib.pyplot as plt
import astropy.units as u

from dust_extinction.parameter_averages import (CCM89, O94, F99, F04, M14)
from dust_extinction.parameter_averages import (CCM89, O94, F99, F04,
GCC09, M14)

fig, ax = plt.subplots()

Expand All @@ -125,20 +122,15 @@ R(V) (+ other variables) dependent prediction models

Rv = 2.0

ext_model = CCM89(Rv=Rv)
ax.plot(x,ext_model(x),label='CCM89')

ext_model = O94(Rv=Rv)
ax.plot(x,ext_model(x),label='O94')

ext_model = F99(Rv=Rv)
ax.plot(x,ext_model(x),label='F99')
models = [CCM89, O94, F99, F04, GCC09, M14]

ext_model = F04(Rv=Rv)
ax.plot(x,ext_model(x),label='F04')

ext_model = M14(Rv=Rv)
ax.plot(x[x<3.3/u.micron],ext_model(x[x<3.3/u.micron]),label='M14')
for cmodel in models:
ext_model = cmodel(Rv=Rv)
indxs, = np.where(np.logical_and(
x.value >= ext_model.x_range[0],
x.value <= ext_model.x_range[1]))
yvals = ext_model(x[indxs])
ax.plot(x[indxs], yvals, label=ext_model.__class__.__name__)

ax.set_xlabel('$x$ [$\mu m^{-1}$]')
ax.set_ylabel('$A(x)/A(V)$')
Expand All @@ -156,7 +148,8 @@ R(V) (+ other variables) dependent prediction models
import matplotlib.pyplot as plt
import astropy.units as u

from dust_extinction.parameter_averages import (CCM89, O94, F99, F04, M14)
from dust_extinction.parameter_averages import (CCM89, O94, F99, F04,
GCC09, M14)

fig, ax = plt.subplots()

Expand All @@ -165,20 +158,15 @@ R(V) (+ other variables) dependent prediction models

Rv = 5.5

ext_model = CCM89(Rv=Rv)
ax.plot(x,ext_model(x),label='CCM89')

ext_model = O94(Rv=Rv)
ax.plot(x,ext_model(x),label='O94')

ext_model = F99(Rv=Rv)
ax.plot(x,ext_model(x),label='F99')

ext_model = F04(Rv=Rv)
ax.plot(x,ext_model(x),label='F04')
models = [CCM89, O94, F99, F04, GCC09, M14]

ext_model = M14(Rv=Rv)
ax.plot(x[x<3.3/u.micron],ext_model(x[x<3.3/u.micron]),label='M14')
for cmodel in models:
ext_model = cmodel(Rv=Rv)
indxs, = np.where(np.logical_and(
x.value >= ext_model.x_range[0],
x.value <= ext_model.x_range[1]))
yvals = ext_model(x[indxs])
ax.plot(x[indxs], yvals, label=ext_model.__class__.__name__)

ax.set_xlabel('$x$ [$\mu m^{-1}$]')
ax.set_ylabel('$A(x)/A(V)$')
Expand Down
114 changes: 113 additions & 1 deletion dust_extinction/parameter_averages.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,12 +8,15 @@
from .averages import G03_SMCBar
from .shapes import _curve_F99_method

__all__ = ['CCM89', 'O94', 'F99', 'F04', 'M14', 'G16']
__all__ = ['CCM89',
'O94', 'F99', 'F04', 'GCC09',
'M14', 'G16']

x_range_CCM89 = [0.3, 10.0]
x_range_O94 = x_range_CCM89
x_range_F99 = [0.3, 10.0]
x_range_F04 = [0.3, 10.0]
x_range_GCC09 = [3.3, 11.0]
x_range_M14 = [0.3, 3.3]
x_range_G16 = [0.3, 10.0]

Expand Down Expand Up @@ -512,6 +515,115 @@ def evaluate(self, in_x, Rv):
self.x_range, 'F04')


class GCC09(BaseExtRvModel):
"""
GCC09 extinction model calculation
Parameters
----------
Rv: float
R(V) = A(V)/E(B-V) = total-to-selective extinction
Raises
------
InputParameterError
Input Rv values outside of defined range
Notes
-----
GCC09 Milky Way R(V) dependent extinction model
From Gordon, Cartledge, & Clayton (2009, ApJ, 705, 1320)
Including erratum: 2014, ApJ, 781, 128
Example showing GCC09 curves for a range of R(V) values.
.. plot::
:include-source:
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from dust_extinction.parameter_averages import GCC09
fig, ax = plt.subplots()
# generate the curves and plot them
x = np.arange(3.3, 11, 0.1)/u.micron
Rvs = ['2.0','3.0','4.0','5.0','6.0']
for cur_Rv in Rvs:
ext_model = GCC09(Rv=cur_Rv)
ax.plot(x,ext_model(x),label='R(V) = ' + str(cur_Rv))
ax.set_xlabel(r'$x$ [$\mu m^{-1}$]')
ax.set_ylabel(r'$A(x)/A(V)$')
ax.legend(loc='best')
plt.show()
"""
Rv_range = [2.0, 6.0]
x_range = x_range_GCC09

@staticmethod
def evaluate(in_x, Rv):
"""
GCC09 function
Parameters
----------
in_x: float
expects either x in units of wavelengths or frequency
or assumes wavelengths in wavenumbers [1/micron]
internally wavenumbers are used
Returns
-------
axav: np array (float)
A(x)/A(V) extinction curve [mag]
Raises
------
ValueError
Input x values outside of defined range
"""
# convert to wavenumbers (1/micron) if x input in units
# otherwise, assume x in appropriate wavenumber units
x = _get_x_in_wavenumbers(in_x)

# check that the wavenumbers are within the defined range
_test_valid_x_range(x, x_range_GCC09, 'GCC09')

# check that the wavenumbers are within the defined range
_test_valid_x_range(x, x_range_GCC09, 'GCC09')

# setup the a & b coefficient vectors
n_x = len(x)
a = np.zeros(n_x)
b = np.zeros(n_x)

# define the ranges
nuv_indxs = np.where(np.logical_and(3.3 <= x, x <= 11.0))
fnuv_indxs = np.where(np.logical_and(5.9 <= x, x <= 11.0))

# NUV
a[nuv_indxs] = (1.894 - .373*x[nuv_indxs]
- 0.0101/((x[nuv_indxs] - 4.57)**2 + .0384))
b[nuv_indxs] = (-3.490
+ 2.057*x[nuv_indxs]
+ 0.706/((x[nuv_indxs] - 4.59)**2 + .169))

# far-NUV
y = x[fnuv_indxs] - 5.9
a[fnuv_indxs] += -.110*(y**2) - .0100*(y**3)
b[fnuv_indxs] += .531*(y**2) + .0544*(y**3)

# return A(x)/A(V)
return a + b/Rv


class M14(BaseExtRvModel):

"""
Expand Down
2 changes: 1 addition & 1 deletion dust_extinction/tests/test_g03.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ def test_invalid_micron(x_invalid_micron, tmodel):

@pytest.mark.parametrize("x_invalid_angstrom", u.angstrom*1e4/x_bad)
@pytest.mark.parametrize("tmodel", models)
def test_invalid_anstrom(x_invalid_angstrom, tmodel):
def test_invalid_angstrom(x_invalid_angstrom, tmodel):
_invalid_x_range(x_invalid_angstrom, tmodel, 'G03')


Expand Down
76 changes: 76 additions & 0 deletions dust_extinction/tests/test_gcc09.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
import numpy as np
import pytest

import astropy.units as u

from ..parameter_averages import GCC09
from .helpers import _invalid_x_range


x_bad = [-1.0, 0.1, 12.0, 100.]


@pytest.mark.parametrize("x_invalid", x_bad)
def test_invalid_wavenumbers(x_invalid):
_invalid_x_range(x_invalid, GCC09(), 'GCC09')


@pytest.mark.parametrize("x_invalid_wavenumber", x_bad/u.micron)
def test_invalid_wavenumbers_imicron(x_invalid_wavenumber):
_invalid_x_range(x_invalid_wavenumber, GCC09(), 'GCC09')


@pytest.mark.parametrize("x_invalid_micron", u.micron/x_bad)
def test_invalid_micron(x_invalid_micron):
_invalid_x_range(x_invalid_micron, GCC09(), 'GCC09')


@pytest.mark.parametrize("x_invalid_angstrom", u.angstrom*1e4/x_bad)
def test_invalid_angstrom(x_invalid_angstrom):
_invalid_x_range(x_invalid_angstrom, GCC09(), 'GCC09')


def get_axav_cor_vals(Rv):
# testing wavenumbers
x = np.array([10., 9., 8., 7.,
6., 5., 4.6, 4.,
3.4])

# add units
x = x/u.micron

# correct values
if Rv == 3.1:
cor_vals = np.array([5.23161, 4.20810, 3.45123, 2.92264, 2.61283,
2.85130, 3.19451, 2.34301, 1.89256])
elif Rv == 2.0:
cor_vals = np.array([10.5150, 8.07274, 6.26711, 5.00591, 4.24237,
4.42844, 4.99482, 3.42585, 2.59322])
elif Rv == 3.0:
cor_vals = np.array([5.55181, 4.44232, 3.62189, 3.04890, 2.71159,
2.94688, 3.30362, 2.40863, 1.93502])
elif Rv == 4.0:
cor_vals = np.array([3.07020, 2.62711, 2.29927, 2.07040, 1.94621,
2.20610, 2.45801, 1.90003, 1.60592])
elif Rv == 5.0:
cor_vals = np.array([1.58123, 1.53798, 1.50571, 1.48330, 1.48697,
1.76164, 1.95065, 1.59486, 1.40846])
elif Rv == 6.0:
cor_vals = np.array([0.588581, 0.811898, 0.976660, 1.09190, 1.18082,
1.46533, 1.61241, 1.39142, 1.27682])
else:
cor_vals = np.array([0.0])

return (x, cor_vals)


@pytest.mark.parametrize("Rv", [2.0, 3.0, 3.1, 4.0, 5.0, 6.0])
def test_extinction_GCC09_values(Rv):
# get the correct values
x, cor_vals = get_axav_cor_vals(Rv)

# initialize extinction model
tmodel = GCC09(Rv=Rv)

# test
np.testing.assert_allclose(tmodel(x), cor_vals, rtol=1e-5)

0 comments on commit c958d1b

Please sign in to comment.