G21 and G21_MWAvg models

docs/dust_extinction/choose_model.rst

@@ -21,7 +21,8 @@ simplest models and include models for the MW
 :class:`~dust_extinction.averages.CT06_MWLoc`,
 :class:`~dust_extinction.averages.CT06_MWGC`,
 :class:`~dust_extinction.averages.GCC09_MWAvg`,
-:class:`~dust_extinction.averages.F11_MWGC`;
+:class:`~dust_extinction.averages.F11_MWGC`,
+:class:`~dust_extinction.averages.G21_MWAvg`;
 Note the different valid wavelength ranges), the LMC
 (:class:`~dust_extinction.averages.G03_LMCAvg`,
 :class:`~dust_extinction.averages.G03_LMC2`) and the SMC
@@ -35,12 +36,14 @@ Way, the usual average used is R(V) = 3.1.
 | Model | x range | wavelength range | galaxy |
 | | [1/micron] | [micron] | |
 +==============+=============+==================+==============+
-| GCC09_MWAvg | 0.3 - 10.96 | 0.0912 - 3.3 | MW |
-+--------------+-------------+------------------+--------------+
 | B92_MWAvg | 1.3 - 2.9 | 0.34 - 0.78 | MW |
 +--------------+-------------+------------------+--------------+
 | I05_MWAvg | 0.13 - 0.8 | 1.24 - 7.76 | MW |
 +--------------+-------------+------------------+--------------+
+| GCC09_MWAvg | 0.3 - 10.96 | 0.0912 - 3.3 | MW |
++--------------+-------------+------------------+--------------+
+| G21_MWAvg | 0.3125 - 1 | 1 - 32 | MW |
++--------------+-------------+------------------+--------------+
 | CT06_MWLoc | 0.037 - 0.8 | 1.24 - 27.0 | MW (Local) |
 +--------------+-------------+------------------+--------------+
 | RL85_MWGC | 0.08 - 0.8 | 1.25 - 13.0 | MW (GCenter) |
@@ -167,7 +170,9 @@ 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.
+but only covers the UV wavelength range. The
+:class:`~dust_extinction.shapes.G21` model focuses on the NIR/MIR
+wavelength range from 1-40 micron.
 
 +------------+--------------+------------------+-------------------+
 | Model | x range | wavelength range | # of parameters |
@@ -177,3 +182,5 @@ but only covers the UV wavelength range.
 +------------+--------------+------------------+-------------------+
 | P92 | 0.001 - 1000 | 0.001 - 1000 | 19 (24 possible) |
 +------------+--------------+------------------+-------------------+
+| G21 | 0.025 - 1 | 1 - 40 | 10 |
++------------+--------------+------------------+-------------------+

docs/dust_extinction/model_flavors.rst

@@ -67,15 +67,16 @@ Average models
  I05_MWAvg,
  CT06_MWLoc,
  CT06_MWGC,
- F11_MWGC)
+ F11_MWGC,
+ G21_MWAvg)
 
  fig, ax = plt.subplots()
 
  # generate the curves and plot them
  x = 1.0 / (np.arange(1.0, 40.0 ,0.1) * u.micron)
 
  models = [RL85_MWGC, RRP89_MWGC, I05_MWAvg, CT06_MWLoc, CT06_MWGC,
- F11_MWGC]
+ F11_MWGC, G21_MWAvg]
 
  for cmodel in models:
  ext_model = cmodel()
@@ -435,6 +436,8 @@ Shape fitting models
  shape of the ultraviolet extinction.
  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
+ of the NIR/MIR 1-40 micron extinction.
 
 .. plot::
 
@@ -524,3 +527,41 @@ Shape fitting models
  ax.legend(loc='best')
  plt.tight_layout()
  plt.show()
+
+.. plot::
+
+ import numpy as np
+ import matplotlib.pyplot as plt
+ import astropy.units as u
+
+ from dust_extinction.shapes import G21
+
+ fig, ax = plt.subplots()
+
+ # generate the curves and plot them
+ lam = np.logspace(np.log10(1.01), np.log10(39.9), num=1000)
+ x = (1.0/lam)/u.micron
+
+ ext_model = G21()
+ ax.plot(1/x,ext_model(x),label='total')
+
+ ext_model = G21(sil1_amp=0.0, sil2_amp=0.0)
+ ax.plot(1./x,ext_model(x),label='power-law only')
+
+ ext_model = G21(sil2_amp=0.0)
+ ax.plot(1./x,ext_model(x),label='power-law+sil1 only')
+
+ ext_model = G21(sil1_amp=0.0)
+ ax.plot(1./x,ext_model(x),label='power-law+sil2 only')
+
+ ax.set_xscale('log')
+ ax.set_yscale('log')
+
+ ax.set_xlabel('$\lambda$ [$\mu$m]')
+ ax.set_ylabel('$A(x)/A(V)$')
+
+ ax.set_title('G21')
+
+ ax.legend(loc='best')
+ plt.tight_layout()
+ plt.show()

docs/dust_extinction/references.rst

@@ -43,9 +43,11 @@ G03: `Gordon et al. 2003, ApJ, 594, 279
 GCC09: `Gordon, Cartledge, & Clayton 2009, ApJ, 705, 1320
 <https://ui.adsabs.harvard.edu/abs/2009ApJ...705.1320G>`_
 
-G16: `Gordon et al. 2016, 826, 104
+G16: `Gordon et al. 2016, ApJ 826, 104
 <https://ui.adsabs.harvard.edu/abs/2016ApJ...826..104G>`_
 
+G21: Gordon et al. 2021, ApJ, submitted
+
 I05: `Indebetouw et al. 2005, ApJ, 619, 931
 <https://ui.adsabs.harvard.edu/abs/2005ApJ...619..931I>`_
 

dust_extinction/averages.py

@@ -8,7 +8,7 @@
 
 from .helpers import _get_x_in_wavenumbers, _test_valid_x_range
 from .baseclasses import BaseExtModel
-from .shapes import P92, _curve_F99_method
+from .shapes import P92, G21, _curve_F99_method
 
 __all__ = [
  "RL85_MWGC",
@@ -22,6 +22,7 @@
  "CT06_MWLoc",
  "GCC09_MWAvg",
  "F11_MWGC",
+ "G21_MWAvg",
 ]
 
 
@@ -425,7 +426,9 @@ def evaluate(self, in_x):
  xo = 4.6
  gamma = 1.0
 
- optnir_axav_x = 1.0 / np.array([2.198, 1.65, 1.25, 0.81, 0.65, 0.55, 0.44, 0.37])
+ optnir_axav_x = 1.0 / 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]
@@ -842,7 +845,9 @@ def __init__(self, **kwargs):
  # get the tabulated information
  data_path = pkg_resources.resource_filename("dust_extinction", "data/")
 
- a = Table.read(data_path + "CT06_pixiedust.dat", format="ascii.commented_header")
+ a = Table.read(
+ data_path + "CT06_pixiedust.dat", format="ascii.commented_header"
+ )
 
  self.obsdata_x = 1.0 / a["wave"].data
  # ext is A(lambda)/A(K)
@@ -942,7 +947,9 @@ def __init__(self, **kwargs):
  # get the tabulated information
  data_path = pkg_resources.resource_filename("dust_extinction", "data/")
 
- a = Table.read(data_path + "CT06_pixiedust.dat", format="ascii.commented_header")
+ a = Table.read(
+ data_path + "CT06_pixiedust.dat", format="ascii.commented_header"
+ )
 
  self.obsdata_x = 1.0 / a["wave"].data
  # ext is A(lambda)/A(K)
@@ -1080,7 +1087,11 @@ def __init__(self, **kwargs):
  (self.obsdata_axav_fuse, self.obsdata_axav_iue, self.obsdata_axav_bands)
  )
  self.obsdata_axav_unc = np.concatenate(
- (self.obsdata_axav_unc_fuse, self.obsdata_axav_unc_iue, self.obsdata_axav_unc_bands,)
+ (
+ self.obsdata_axav_unc_fuse,
+ self.obsdata_axav_unc_iue,
+ self.obsdata_axav_unc_bands,
+ )
  )
 
  # accuracy of the observed data based on published table
@@ -1202,7 +1213,9 @@ def __init__(self, **kwargs):
  # get the tabulated information
  data_path = pkg_resources.resource_filename("dust_extinction", "data/")
 
- a = Table.read(data_path + "fritz11_galcenter.dat", format="ascii.commented_header")
+ a = Table.read(
+ data_path + "fritz11_galcenter.dat", format="ascii.commented_header"
+ )
 
  self.obsdata_x = 1.0 / a["wave"].data
  # ext is total extinction to GalCenter
@@ -1247,3 +1260,142 @@ def evaluate(self, in_x):
  f = interp1d(self.obsdata_x, self.obsdata_axav)
 
  return f(x)
+
+
+class G21_MWAvg(BaseExtModel):
+ r"""
+ Gordon et al. (2021) Milky Way Average Extinction Curve
+
+ Parameters
+ ----------
+ None
+
+ Raises
+ ------
+ None
+
+ Notes
+ -----
+ From Gordon et al. (2021, ApJ, submitted)
+
+ Example showing the average curve
+
+ .. plot::
+ :include-source:
+
+ import numpy as np
+ import matplotlib.pyplot as plt
+ import astropy.units as u
+
+ from dust_extinction.averages import G21_MWAvg
+
+ fig, ax = plt.subplots()
+
+ # generate the curves and plot them
+ lam = np.logspace(np.log10(1.01), np.log10(31.9), num=1000)
+ x = (1.0/lam)/u.micron
+
+ # define the extinction model
+ ext_model = G21_MWAvg()
+
+ ax.plot(1.0/x,ext_model(x),label='G21_MWAvg')
+ ax.errorbar(1.0/ext_model.obsdata_x_irs, ext_model.obsdata_axav_irs,
+ yerr=ext_model.obsdata_axav_unc_irs,
+ fmt='ko', label='obsdata (IRS)')
+ ax.errorbar(1.0/ext_model.obsdata_x_bands, ext_model.obsdata_axav_bands,
+ yerr=ext_model.obsdata_axav_unc_bands,
+ fmt='g^', label='obsdata (photometry)')
+
+ ax.set_xlabel(r'$\lambda$ [$\mu m$]')
+ ax.set_ylabel(r'$A(x)/A(V)$')
+
+ ax.set_xscale('log')
+ ax.set_yscale('log')
+
+ ax.legend(loc='best')
+ plt.show()
+ """
+
+ x_range = [1.0 / 32.0, 1.0]
+
+ Rv = 3.17
+
+ def __init__(self, **kwargs):
+
+ # get the tabulated information
+ data_path = pkg_resources.resource_filename("dust_extinction", "data/")
+
+ # GCC09 sigma clipped average of 75 sightlines
+ a = Table.read(data_path + "G21_IRS.dat", format="ascii.commented_header")
+ b = Table.read(data_path + "G21_PHOT.dat", format="ascii.commented_header")
+
+ # IRS range
+ self.obsdata_x_irs = 1.0 / a["wave"].data
+ self.obsdata_axav_irs = a["ext"].data
+ self.obsdata_axav_unc_irs = a["unc"].data
+
+ # Opt/NIR range
+ self.obsdata_x_bands = 1.0 / b["wave"].data
+ self.obsdata_axav_bands = b["ext"].data
+ self.obsdata_axav_unc_bands = b["unc"].data
+
+ # put them together
+ self.obsdata_x = np.concatenate((self.obsdata_x_irs, self.obsdata_x_bands))
+ self.obsdata_axav = np.concatenate(
+ (self.obsdata_axav_irs, self.obsdata_axav_bands)
+ )
+ self.obsdata_axav_unc = np.concatenate(
+ (
+ self.obsdata_axav_unc_irs,
+ self.obsdata_axav_unc_bands,
+ )
+ )
+
+ # accuracy of the observed data based on published table
+ self.obsdata_tolerance = 5e-1
+
+ super().__init__(**kwargs)
+
+ def evaluate(self, in_x):
+ """
+ G21_MWAvg 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
+ """
+ x = _get_x_in_wavenumbers(in_x)
+
+ # check that the wavenumbers are within the defined range
+ _test_valid_x_range(x, self.x_range, self.__class__.__name__)
+
+ # G21 parameters fit to the data using uncs as weights
+ g21_fit = G21(
+ scale=0.366,
+ alpha=1.480,
+ sil1_amp=0.06893,
+ sil1_center=9.865,
+ sil1_fwhm=2.507,
+ sil1_asym=-0.232,
+ sil2_amp=0.02684,
+ sil2_center=19.973,
+ sil2_fwhm=16.989,
+ sil2_asym=-0.273,
+ )
+
+ # return A(x)/A(V)
+ return g21_fit(in_x)

dust_extinction/data/G21_IRS.dat

@@ -0,0 +1,46 @@
+# wave ext unc
+5.312549407981746 0.029308363795280457 0.003634442170167764
+5.530947850652181 0.027905302743117016 0.0026571168785127328
+5.758324634246696 0.029355501135190327 0.0025729641062889382
+5.995048857577364 0.02604266752799352 0.002619889850274867
+6.241504793076228 0.02616392821073532 0.002530134239443791
+6.498092510581479 0.02660254140694936 0.002367237532485079
+6.765228526767456 0.024048514664173126 0.0024828063332462684
+7.043346481272642 0.026687875390052795 0.0022931202003653086
+7.3328978406232235 0.02662523587544759 0.002213158210630499
+7.63435263109489 0.027789589017629623 0.002682915162113982
+7.948200201702518 0.03038164973258972 0.0024963033812044305
+8.274950018556295 0.03427134454250336 0.002269854126973733
+8.615132491873744 0.0381233811378479 0.002383284960132179
+8.969299836990164 0.053234229485193886 0.002141578808573921
+9.3380269707651 0.07119244833787282 0.00222967406133923
+9.721912444840044 0.08155429859956105 0.0022924791267714653
+10.121579417262236 0.07836014529069264 0.0024396674774119958
+10.53767666405189 0.06589555740356445 0.0021636736265298524
+10.970879632354798 0.05675037418093 0.0018643859620950042
+11.421891536889973 0.04748629033565521 0.00221778032892276
+11.891444501472208 0.039633107682069145 0.0021834892266003837
+12.380300747462408 0.029393449425697327 0.002193555483636868
+12.88925383107517 0.028628597656885784 0.0018674402916312946
+13.419129931551787 0.028001102308432262 0.0018286532585305023
+13.970789192290137 0.02380336821079254 0.0022981451320860996
+14.545127117108091 0.026327942808469135 0.0022807936379891686
+15.143076023907396 0.025927623112996418 0.002345864135562102
+15.765606558097339 0.02731793870528539 0.0024300971072726773
+16.413729268235368 0.03348533312479655 0.002600925103881555
+17.088496246441892 0.032762023309866585 0.0022261574621071064
+17.791002836252677 0.0352509468793869 0.0025123482124195283
+18.522389410680596 0.03661670039097468 0.0034635569054111705
+19.283843223373673 0.03391921520233154 0.003084574761091035
+20.07660033587375 0.03694053739309311 0.003986164119362912
+20.90194762410488 0.03770768344402313 0.0029403961077071575
+21.76122486734801 0.037564915418624875 0.0036820972099157277
+22.65582692309351 0.03561693131923675 0.004824016622403211
+23.587205991301417 0.035626145700613655 0.0029854581006330014
+24.556873971745496 0.036531063417593636 0.003518307169855367
+25.566404918267235 0.0370990534623464 0.0030506717447946743
+26.6174375939243 0.0243232399225235 0.0025511351982575527
+27.711678131180584 0.02782232811053594 0.0029772651129264575
+28.850902801456833 0.024926339586575825 0.0027219564658656926
+30.03696089853681 0.01981469492117564 0.0023255953486483487
+31.27177774051052 0.023677619795004528 0.0035135418754853108