Skip to content

Commit

Permalink
Merge pull request #20 from karllark/move_aves
Browse files Browse the repository at this point in the history 
Move aves
  • Loading branch information
@karllark
karllark committed Apr 26, 2017
2 parents 3b58bc9 + 2c280cd commit 9111ac8
Show file tree
Hide file tree
Showing 4 changed files with 641 additions and 340 deletions.
63 changes: 62 additions & 1 deletion docs/dust_extinction/index.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,65 @@ Interstellar dust extinction curves implemnted as astropy models

Uses the astropy affiliated package template

Flavors of Models
=================

There are three differnet types of models (to be completed).

1. Average models

These models provide averages from the literature with the ability to
interpolate between the observed data points.
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).

.. plot::

import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u

from dust_extinction.dust_extinction_averages import (G03_SMCBar,
G03_LMCAvg,
G03_LMC2)

fig, ax = plt.subplots()

# generate the curves and plot them
x = np.arange(0.3,10.0,0.1)/u.micron

ext_model = G03_SMCBar()
ax.plot(x,ext_model(x),label='G03 SMCBar')

ext_model = G03_LMCAvg()
ax.plot(x,ext_model(x),label='G03 LMCAvg')

ext_model = G03_LMC2()
ax.plot(x,ext_model(x),label='G03 LMC2')

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

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

2. Shape fitting models

These models are used to fit the detailed shape of dust extinction curves.

- FM90
- others needed (P92)

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

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

- CCM89 [function of R(V)]
- F99 [function of R(V)]
- others needed (GCC09, G16, etc)

Repository
==========

Expand All @@ -15,8 +74,10 @@ Reference API
.. toctree::
:maxdepth: 1

.. automodapi:: dust_extinction.dust_extinction_averages

.. automodapi:: dust_extinction.dust_extinction
Indices and tables
==================

Expand Down
231 changes: 1 addition & 230 deletions dust_extinction/dust_extinction.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,12 +11,11 @@
from astropy.modeling import Model, Parameter, InputParameterError

__all__ = ['BaseExtModel','BaseExtRvModel',
'CCM89', 'FM90', 'F99', 'G03']
'CCM89', 'FM90', 'F99']

x_range_CCM89 = [0.3,10.0]
x_range_FM90 = [1.0/0.32,1.0/0.0912]
x_range_F99 = [1.0/6.0,8.7]
x_range_G03 = [0.3,10.0]

class BaseExtModel(Model):
"""
Expand Down Expand Up @@ -546,231 +545,3 @@ def evaluate(in_x, Rv):
# return A(x)/A(V)
return axebv/Rv

class G03(BaseExtModel):
"""
G03 average extinction curves
Parameters
----------
AveNameIndex: int
index of allowed name
[allowed numbers = 1, 2, 3]
[allowed names = 'SMCBar', 'LMCave', 'LMC2']
Raises
------
InputParameterError
AveNameIndex not in allowed list
Notes
-----
G03 average extinction curve
From Gordon et al. (2003, ApJ, 594, 279)
Example showing the G03 average curves
.. plot::
:include-source:
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from dust_extinction.dust_extinction import G03
fig, ax = plt.subplots()
# generate the curves and plot them
x = np.arange(0.3,10.0,0.1)/u.micron
AveNames = ['SMCBar', 'LMCAve', 'LMC2']
AveName_Indexs = [1, 2, 3]
for cur_ave_idx in AveName_Indexs:
ext_model = G03(AveNameIndex=cur_ave_idx)
ax.plot(x,ext_model(x),label=str(cur_ave_idx) +
": " + AveNames[cur_ave_idx-1])
ax.set_xlabel('$x$ [$\mu m^{-1}$]')
ax.set_ylabel('$A(x)/A(V)$')
ax.legend(loc='best')
plt.show()
"""

AveNameIndex = Parameter(description="G03 Average Curve Name Index",
default=1)

x_range = x_range_G03

@AveNameIndex.validator
def AveNameIndex(self, value):
"""
Check that AveName has a valid value
Parameters
----------
value: int
index of allowed name
Raises
------
InputParameterError
Input AveName number not valid
"""
allowed_names = ['SMCBar', 'LMCAve', 'LMC2']
allowed_numbers = [1, 2, 3]
if value not in allowed_numbers:
raise InputParameterError("parameter AveNameIndex must be in "
+ " ".join(str(an) for an in allowed_numbers))

@staticmethod
def evaluate(in_x, AveNameIndex):
"""
G03 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
with u.add_enabled_equivalencies(u.spectral()):
x_quant = u.Quantity(in_x, 1.0/u.micron, dtype=np.float64)

# strip the quantity to avoid needing to add units to all the
# polynomical coefficients
x = x_quant.value

# check that the wavenumbers are within the defined range
if np.logical_or(np.any(x < x_range_G03[0]),
np.any(x > x_range_G03[1])):
raise ValueError('Input x outside of range defined for G03' \
+ ' ['
+ str(x_range_G03[0])
+ ' <= x <= '
+ str(x_range_G03[1])
+ ', x has units 1/micron]')

# initialize extinction curve storage
axav = np.zeros(len(x))

# get the parameters for the requested AveName
# From Tables 2, 3, and 4 of Gordon et al. (2003)

if AveNameIndex == 2: # LMCAve
Rv = 3.41
C1 = -0.890
C2 = 0.998
C3 = 2.719
C4 = 0.400
xo = 4.579
gamma = 0.934

optnir_axav_x = 1./np.array([2.198, 1.65, 1.25,
0.55, 0.44, 0.37])
# value at 2.198 changed to provide smooth interpolation
# as noted in Gordon et al. (2016, ApJ, 826, 104) for SMCBar
optnir_axav_y = [0.10, 0.186, 0.257,
1.000, 1.293, 1.518]
elif AveNameIndex == 3: # LMC2
Rv = 2.76
C1 = -1.475
C2 = 1.132
C3 = 1.463
C4 = 0.294
xo = 4.558
gamma = 0.945

optnir_axav_x = 1./np.array([2.198, 1.65, 1.25,
0.55, 0.44, 0.37])
# value at 1.65 changed to provide smooth interpolation
# as noted in Gordon et al. (2016, ApJ, 826, 104) for SMCBar
optnir_axav_y = [0.101, 0.15, 0.299,
1.000, 1.349, 1.665]
else: # AveNameIndex = 1 # SMCBar
Rv = 2.74
C1 = -4.959
C2 = 2.264
C3 = 0.389
C4 = 0.461
xo = 4.6
gamma = 1.0

optnir_axav_x = 1./np.array([2.198, 1.65, 1.25, 0.81, 0.65,
0.55, 0.44, 0.37])

# values at 2.198 and 1.25 changed to provide smooth interpolation
# as noted in Gordon et al. (2016, ApJ, 826, 104)
optnir_axav_y = [0.11, 0.169, 0.25, 0.567, 0.801,
1.00, 1.374, 1.672]

# x value above which FM90 parametrization used
x_cutval_uv = 10000.0/2700.0

# required UV points for spline interpolation
x_splineval_uv = 10000.0/np.array([2700.0,2600.0])

# UV points in input x
indxs_uv, = np.where(x >= x_cutval_uv)

# add in required spline points, otherwise just spline points
if len(indxs_uv) > 0:
xuv = np.concatenate([x_splineval_uv, x[indxs_uv]])
else:
xuv = x_splineval_uv

# FM90 model and values
fm90_model = FM90(C1=C1, C2=C2, C3=C3, C4=C4, xo=xo, gamma=gamma)
# evaluate model and get results in A(x)/A(V)
axav_fm90 = fm90_model(xuv)/Rv + 1.0

# save spline points
y_splineval_uv = axav_fm90[0:2]

# ingore the spline points
if len(indxs_uv) > 0:
axav[indxs_uv] = axav_fm90[2:]

# **Optical Portion**
# using cubic spline anchored in UV, optical, and IR

# optical/NIR points in input x
indxs_opir, = np.where(x < x_cutval_uv)

if len(indxs_opir) > 0:
# spline points
x_splineval_optir = optnir_axav_x

# determine optical/IR values at spline points
y_splineval_optir = optnir_axav_y

# add in zero extinction at infinite wavelength
x_splineval_optir = np.insert(x_splineval_optir, 0, 0.0)
y_splineval_optir = np.insert(y_splineval_optir, 0, 0.0)

spline_x = np.concatenate([x_splineval_optir, x_splineval_uv])
spline_y = np.concatenate([y_splineval_optir, y_splineval_uv])
spline_rep = interpolate.splrep(spline_x, spline_y)
axav[indxs_opir] = interpolate.splev(x[indxs_opir],
spline_rep, der=0)

# return A(x)/A(V)
return axav

Loading

0 comments on commit 9111ac8

Please sign in to comment.