Merge pull request #201 from karllark/add_FM90_B3

Add version of FM90 using B3 instead of C3

docs/dust_extinction/choose_model.rst

@@ -183,8 +183,11 @@ used to fit measured extinction curves. These models allow features in the
 extinction curve to be measured (e.g., 2175 A bump or 10 micron silicate
 feature). The :class:`~dust_extinction.shapes.P92` is the most
 general as it covers the a very broad wavelength range. The
-:class:`~dust_extinction.shapes.FM90` model has been extensively used,
-but only covers the UV wavelength range. The
+:class:`~dust_extinction.shapes.FM90` model covers the UV wavelength range
+and has been extensively shown to fit all known UV extinction curves. 
+The :class:`~dust_extinction.shapes.FM90_B3` model provides a variant
+of the FM90 model that uses B3 instead of C3 as B3 = explicit 2175 A 
+bump height = C3/gamma^2.
 :class:`~dust_extinction.shapes.G21` model focuses on the NIR/MIR
 wavelength range from 1-40 micron.
 
@@ -194,6 +197,8 @@ wavelength range from 1-40 micron.
 +============+==============+==================+===================+
 | FM90 | 3.13 - 11.0 | 0.0912 - 0.32 | 6 |
 +------------+--------------+------------------+-------------------+
+| FM90_B3 | 3.13 - 11.0 | 0.0912 - 0.32 | 6 |
++------------+--------------+------------------+-------------------+
 | P92 | 0.001 - 1000 | 0.001 - 1000 | 19 (24 possible) |
 +------------+--------------+------------------+-------------------+
 | G21 | 0.025 - 1 | 1 - 40 | 10 |

docs/dust_extinction/model_flavors.rst

@@ -492,6 +492,9 @@ 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.
+ Note there are two forms of the FM90 model, FM90 that implements the model
+ as published and FM90_B3 that B3 = C3/gamma^2 as the explicit amplitude of
+ the 2175 A bump (easier to interpret).
  The P92 (Pei 1992) uses 19 parameters to fit the shape of the X-ray to
  far-infrared extinction.
  The G21 (Gordon et al. 2021) models uses 10 parameters to fit the shape

dust_extinction/__init__.py

@@ -2,5 +2,5 @@
 
 try:
  __version__ = _version(__name__)
-except PackageNotFoundError:
+except PackageNotFoundError: # pragma: no cover
  pass

dust_extinction/shapes.py

@@ -6,7 +6,7 @@
 
 from .helpers import _get_x_in_wavenumbers, _test_valid_x_range
 
-__all__ = ["FM90", "P92", "G21"]
+__all__ = ["FM90", "FM90_B3", "P92", "G21"]
 
 x_range_FM90 = [1.0 / 0.35, 1.0 / 0.09]
 x_range_P92 = [1.0 / 1e3, 1.0 / 1e-3]
@@ -179,7 +179,7 @@ class FM90(Fittable1DModel):
  slope of liner term
 
  C3: float
- amplitude of "2175 A" bump
+ strength of "2175 A" bump (true amplitude is C3/gamma^2)
 
  C4: float
  amplitude of FUV rise
@@ -332,6 +332,128 @@ def fit_deriv(in_x, C1, C2, C3, C4, xo, gamma):
  # 'C4': outputs_unit[self.outputs[0]]}
 
 
+class FM90_B3(Fittable1DModel):
+ r"""
+ Fitzpatrick & Massa (1990) 6 parameter ultraviolet shape model
+ Version with bump amplitude B3 = C3/gamma^2
+
+ Parameters
+ ----------
+ C1: float
+ y-intercept of linear term
+
+ C2: float
+ slope of liner term
+
+ B3: float
+ amplitude of "2175 A" bump
+
+ C4: float
+ amplitude of FUV rise
+
+ xo: float
+ centroid of "2175 A" bump
+
+ gamma: float
+ width of "2175 A" bump
+
+ Notes
+ -----
+ From Fitzpatrick & Massa (1990, ApJS, 72, 163)
+
+ Only applicable at UV wavelengths
+
+ Example showing a FM90 curve with components identified.
+
+ .. plot::
+ :include-source:
+
+ import numpy as np
+ import matplotlib.pyplot as plt
+ import astropy.units as u
+
+ from dust_extinction.shapes import FM90_B3
+
+ fig, ax = plt.subplots()
+
+ # generate the curves and plot them
+ x = np.arange(3.8,8.6,0.1)/u.micron
+
+ ext_model = FM90_B3()
+ ax.plot(x,ext_model(x),label='total')
+
+ ext_model = FM90_B3(B3=0.0, C4=0.0)
+ ax.plot(x,ext_model(x),label='linear term')
+
+ ext_model = FM90_B3(C1=0.0, C2=0.0, C4=0.0)
+ ax.plot(x,ext_model(x),label='bump term')
+
+ ext_model = FM90_B3(C1=0.0, C2=0.0, B3=0.0)
+ ax.plot(x,ext_model(x),label='FUV rise term')
+
+ ax.set_xlabel(r'$x$ [$\mu m^{-1}$]')
+ ax.set_ylabel(r'$E(\lambda - V)/E(B - V)$')
+
+ ax.legend(loc='best')
+ plt.show()
+ """
+ n_inputs = 1
+ n_outputs = 1
+
+ C1 = Parameter(description="linear term: y-intercept", default=0.10)
+ C2 = Parameter(description="linear term: slope", default=0.70)
+ B3 = Parameter(description="bump: amplitude", default=3.23)
+ C4 = Parameter(description="FUV rise: amplitude", default=0.41)
+ xo = Parameter(description="bump: centroid", default=4.60, min=0.0)
+ gamma = Parameter(description="bump: width", default=0.99, min=0.0)
+
+ x_range = x_range_FM90
+
+ @staticmethod
+ def evaluate(in_x, C1, C2, B3, C4, xo, gamma):
+ """
+ FM90 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
+ -------
+ exvebv: np array (float)
+ E(x-V)/E(B-V) extinction curve [mag]
+
+ Raises
+ ------
+ ValueError
+ Input x values outside of defined range
+ """
+ x = _get_x_in_wavenumbers(in_x)
+
+ # check that the wavenumbers are within the defined range
+ _test_valid_x_range(x, x_range_FM90, "FM90_B3")
+
+ # linear term
+ exvebv = C1 + C2 * x
+
+ # bump term
+ x2 = x ** 2
+ exvebv += (B3 * gamma ** 2) * (x2 / ((x2 - xo ** 2) ** 2 + x2 * (gamma ** 2)))
+
+ # FUV rise term
+ fnuv_indxs = np.where(x >= 5.9)
+ if len(fnuv_indxs) > 0:
+ y = x[fnuv_indxs] - 5.9
+ exvebv[fnuv_indxs] += C4 * (0.5392 * (y ** 2) + 0.05644 * (y ** 3))
+
+ # return E(x-V)/E(B-V)
+ return exvebv
+
+
 class P92(Fittable1DModel):
  r"""
  Pei (1992) 24 parameter shape model

dust_extinction/tests/helpers.py

@@ -13,7 +13,7 @@
  D22,
  G23,
 )
-from dust_extinction.shapes import FM90, P92, G21
+from dust_extinction.shapes import FM90, FM90_B3, P92, G21
 from dust_extinction.averages import (
  RL85_MWGC,
  RRP89_MWGC,
@@ -34,7 +34,7 @@
 param_ave_models_Rv = [CCM89, O94, F99, F04, VCG04, GCC09, M14, F19, D22, G23]
 param_ave_models_Rv_fA = [G16]
 param_ave_models = param_ave_models_Rv + param_ave_models_Rv_fA
-shape_models = [FM90, P92, G21]
+shape_models = [FM90, FM90_B3, P92, G21]
 ave_models = [
  RL85_MWGC,
  RRP89_MWGC,