Skip to content

Commit

Permalink
Start of FM90 as a fittable model
Browse files Browse the repository at this point in the history
  • Loading branch information
@karllark
karllark committed May 25, 2017
1 parent 9111ac8 commit 60b4e07
Showing 1 changed file with 71 additions and 40 deletions.
111 changes: 71 additions & 40 deletions dust_extinction/dust_extinction.py
View file
@@ -1,4 +1,3 @@

from __future__ import (absolute_import, print_function, division)

# STDLIB
Expand All @@ -8,7 +7,8 @@
from scipy import interpolate

import astropy.units as u
from astropy.modeling import Model, Parameter, InputParameterError
from astropy.modeling import (Model, Fittable1DModel,
Parameter, InputParameterError)

__all__ = ['BaseExtModel','BaseExtRvModel',
'CCM89', 'FM90', 'F99']
Expand All @@ -23,51 +23,51 @@ class BaseExtModel(Model):
"""
inputs = ('x',)
outputs = ('axav',)

def extinguish(self, x, Av=None, Ebv=None):
"""
Calculate the extinction as a fraction
Parameters
----------
x: float
expects either x in units of wavelengths or frequency
or assumes wavelengths in wavenumbers [1/micron]
expects either x in units of wavelengths or frequency
or assumes wavelengths in wavenumbers [1/micron]
internally wavenumbers are used
internally wavenumbers are used
Av: float
A(V) value of dust column
Av or Ebv must be set
A(V) value of dust column
Av or Ebv must be set
Ebv: float
E(B-V) value of dust column
Av or Ebv must be set
E(B-V) value of dust column
Av or Ebv must be set
Returns
-------
frac_ext: np array (float)
fractional extinction as a function of x
fractional extinction as a function of x
"""
# get the extinction curve
axav = self(x)

# check that av or ebv is set
if (Av is None) and (Ebv is None):
raise InputParameterError('neither Av or Ebv passed, one required')

# if Av is not set and Ebv set, convert to Av
if Av is None:
Av = self.Rv*Ebv

# return fractional extinction
return np.power(10.0,-0.4*axav*Av)

class BaseExtRvModel(BaseExtModel):
"""
Base Extinction R(V)-dependent Model. Do not use.
"""

Rv = Parameter(description="R(V) = A(V)/E(B-V) = " \
+ "total-to-selective extinction",
default=3.1)
Expand All @@ -92,7 +92,7 @@ def Rv(self, value):
+ str(self.Rv_range[0])
+ " and "
+ str(self.Rv_range[1]))

class CCM89(BaseExtRvModel):
"""
CCM89 extinction model calculation
Expand Down Expand Up @@ -128,15 +128,15 @@ class CCM89(BaseExtRvModel):
# generate the curves and plot them
x = np.arange(0.5,10.0,0.1)/u.micron
Rvs = ['2.0','3.0','4.0','5.0','6.0']
for cur_Rv in Rvs:
ext_model = CCM89(Rv=cur_Rv)
ax.plot(x,ext_model(x),label='R(V) = ' + str(cur_Rv))
ax.set_xlabel('$x$ [$\mu m^{-1}$]')
ax.set_ylabel('$A(x)/A(V)$')
ax.legend(loc='best')
plt.show()
"""
Expand Down Expand Up @@ -170,7 +170,7 @@ def evaluate(in_x, Rv):
# 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)
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
Expand All @@ -185,7 +185,7 @@ def evaluate(in_x, Rv):
+ ' <= x <= '
+ str(x_range_CCM89[1])
+ ', x has units 1/micron]')

# setup the a & b coefficient vectors
n_x = len(x)
a = np.zeros(n_x)
Expand All @@ -197,31 +197,31 @@ def evaluate(in_x, Rv):
nuv_indxs = np.where(np.logical_and(3.3 <= x,x <= 8.0))
fnuv_indxs = np.where(np.logical_and(5.9 <= x,x <= 8))
fuv_indxs = np.where(np.logical_and(8 < x,x <= 10))

# Infrared
y = x[ir_indxs]**1.61
a[ir_indxs] = .574*y
b[ir_indxs] = -0.527*y

# NIR/optical
y = x[opt_indxs] - 1.82
a[opt_indxs] = np.polyval((.32999, -.7753, .01979, .72085, -.02427,
-.50447, .17699, 1), y)
b[opt_indxs] = np.polyval((-2.09002, 5.3026, -.62251, -5.38434,
1.07233, 2.28305, 1.41338, 0), y)

# NUV
a[nuv_indxs] = 1.752-.316*x[nuv_indxs] \
- 0.104/((x[nuv_indxs] - 4.67)**2 + .341)
b[nuv_indxs] = -3.09 + \
1.825*x[nuv_indxs] \
+ 1.206/((x[nuv_indxs] - 4.62)**2 + .263)

# far-NUV
y = x[fnuv_indxs] - 5.9
a[fnuv_indxs] += -.04473*(y**2) - .009779*(y**3)
b[fnuv_indxs] += .2130*(y**2) + .1207*(y**3)
b[fnuv_indxs] += .2130*(y**2) + .1207*(y**3)

# FUV
y = x[fuv_indxs] - 8.0
a[fuv_indxs] = np.polyval((-.070, .137, -.628, -1.073), y)
Expand All @@ -230,10 +230,10 @@ def evaluate(in_x, Rv):
# return A(x)/A(V)
return a + b/Rv

class FM90(Model):
class FM90(Fittable1DModel):
"""
FM90 extinction model calculation
Parameters
----------
C1: float
Expand All @@ -246,11 +246,11 @@ class FM90(Model):
amplitude of "2175 A" bump
C4: float
amplitude of FUV rise
amplitude of FUV rise
xo: float
centroid of "2175 A" bump
gamma: float
width of "2175 A" bump
Expand All @@ -277,7 +277,7 @@ class FM90(Model):
# generate the curves and plot them
x = np.arange(3.8,8.6,0.1)/u.micron
ext_model = FM90()
ax.plot(x,ext_model(x),label='total')
Expand All @@ -292,13 +292,13 @@ class FM90(Model):
ax.set_xlabel('$x$ [$\mu m^{-1}$]')
ax.set_ylabel('$E(\lambda - V)/E(B - V)$')
ax.legend(loc='best')
plt.show()
"""
inputs = ('x',)
outputs = ('exvebv',)

C1 = Parameter(description="linear term: y-intercept",
default=0.10)
C2 = Parameter(description="linear term: slope",
Expand Down Expand Up @@ -340,7 +340,7 @@ def evaluate(in_x, C1, C2, C3, C4, xo, gamma):
# 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)
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
Expand Down Expand Up @@ -372,6 +372,38 @@ def evaluate(in_x, C1, C2, C3, C4, xo, gamma):
# return E(x-V)/E(B-V)
return exvebv

@staticmethod
def fit_deriv(in_x, C1, C2, C3, C4, xo, gamma):
"""
Derivatives of the FM90 function with respect to the parameters
"""

# useful quantitites
x2 = x**2
xo2 = xo**2
g2 = gamma**2
x2mxo2_2 = (x2 - xo2)**2
denom = (x2mxo2_2 - x2*g2)**2

# derivatives

d_C1 = 1.
d_C2 = x

d_C3 = (x2/(x2mxo2_2 + x2*g2))

d_xo = (4.*C2*x2*xo*(x2 - xo2))/denom

d_gamma = (2.*C2*(x2**2)*gamma)/denom

d_C4 = np.zeros((len(x)))
fnuv_indxs = np.where(x >= 5.9)
if len(fnuv_indxs) > 0:
y = x[fnuv_indxs] - 5.9
d_C4[fuv_indxs] = (0.5392*(y**2) + 0.05644*(y**3))

return [d_C1, d_C2, d_C3, d_C4, d_xo, d_gamm]

class F99(BaseExtRvModel):
"""
F99 extinction model calculation
Expand Down Expand Up @@ -410,15 +442,15 @@ class F99(BaseExtRvModel):
# generate the curves and plot them
x = np.arange(0.2,8.7,0.1)/u.micron
Rvs = ['2.0','3.0','4.0','5.0','6.0']
for cur_Rv in Rvs:
ext_model = F99(Rv=cur_Rv)
ax.plot(x,ext_model(x),label='R(V) = ' + str(cur_Rv))
ax.set_xlabel('$x$ [$\mu m^{-1}$]')
ax.set_ylabel('$A(x)/A(V)$')
ax.legend(loc='best')
plt.show()
"""
Expand Down Expand Up @@ -452,7 +484,7 @@ def evaluate(in_x, Rv):
# 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)
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
Expand All @@ -467,7 +499,7 @@ def evaluate(in_x, Rv):
+ ' <= x <= '
+ str(x_range_F99[1])
+ ', x has units 1/micron]')

# ensure Rv is a single element, not numpy array
Rv = Rv[0]

Expand Down Expand Up @@ -514,7 +546,7 @@ def evaluate(in_x, Rv):
# ingore the spline points
if len(indxs_uv) > 0:
axebv[indxs_uv] = axebv_fm90[2:]

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

Expand All @@ -533,7 +565,7 @@ def evaluate(in_x, Rv):
-0.050 + 1.0016*Rv,
0.701 + 1.016*Rv,
1.208 + 1.0032*Rv - 0.00033*(Rv**2)])
y_splineval_ir = np.array([0.0,0.265,0.829])*Rv/3.1
y_splineval_ir = np.array([0.0,0.265,0.829])*Rv/3.1
y_splineval_optir = np.concatenate([y_splineval_ir,y_splineval_opt])

spline_x = np.concatenate([x_splineval_optir, x_splineval_uv])
Expand All @@ -544,4 +576,3 @@ def evaluate(in_x, Rv):

# return A(x)/A(V)
return axebv/Rv

0 comments on commit 60b4e07

Please sign in to comment.