G23 added with minimal tests and some docs

CHANGES.rst

@@ -1,7 +1,7 @@
 1.2 (unreleased)
 ================
 
-- No changes yet
+- Added G23 parameter average model
 
 1.1 (2022-03-30)
 ================

docs/dust_extinction/model_flavors.rst

@@ -24,6 +24,7 @@ Average models
 
  import numpy as np
  import matplotlib.pyplot as plt
+ from matplotlib.ticker import ScalarFormatter
  import astropy.units as u
 
  from dust_extinction.averages import (GCC09_MWAvg,
@@ -45,10 +46,13 @@ Average models
  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.plot(1./x[indxs], yvals, label=ext_model.__class__.__name__)
 
- ax.set_xlabel('$x$ [$\mu m^{-1}$]')
- ax.set_ylabel('$A(x)/A(V)$')
+ ax.set_xscale('log')
+ ax.xaxis.set_major_formatter(ScalarFormatter())
+
+ ax.set_xlabel(r'$\lambda$ [$\mu$m]')
+ ax.set_ylabel(r'$A(x)/A(V)$')
  ax.set_title('Ultraviolet to Near-Infrared Models')
 
  ax.legend(loc='best')
@@ -60,6 +64,7 @@ Average models
 
  import numpy as np
  import matplotlib.pyplot as plt
+ from matplotlib.ticker import ScalarFormatter
  import astropy.units as u
 
  from dust_extinction.averages import (RL85_MWGC,
@@ -87,8 +92,10 @@ Average models
  yvals = ext_model(x[indxs])
  ax.plot(1.0 / x[indxs], yvals, label=ext_model.__class__.__name__)
 
- ax.set_yscale("log")
- ax.set_xlabel(r'$\lambda$ [$\mu m$]')
+ ax.set_xscale('log')
+ ax.xaxis.set_major_formatter(ScalarFormatter())
+
+ ax.set_xlabel(r'$\lambda$ [$\mu$m]')
  ax.set_ylabel(r'$A(\lambda)/A(V)$')
  ax.set_title('Near- to Mid-Infrared Models')
 
@@ -121,10 +128,12 @@ R(V) (+ other variables) dependent prediction models
 
  import numpy as np
  import matplotlib.pyplot as plt
+ from matplotlib.ticker import ScalarFormatter
  import astropy.units as u
 
  from dust_extinction.parameter_averages import (CCM89, O94, F99, F04,
- VCG04, GCC09, M14, F19, D22)
+ VCG04, GCC09, M14, F19, D22,
+ G23)
 
  fig, ax = plt.subplots()
 
@@ -133,17 +142,20 @@ R(V) (+ other variables) dependent prediction models
 
  Rv = 3.1
 
- models = [CCM89, O94, F99, F04, VCG04, GCC09, M14, F19, D22]
+ models = [CCM89, O94, F99, F04, VCG04, GCC09, M14, F19, D22, G23]
 
  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.plot(1./x[indxs], yvals, label=ext_model.__class__.__name__)
+
+ ax.set_xscale('log')
+ ax.xaxis.set_major_formatter(ScalarFormatter())
 
- ax.set_xlabel('$x$ [$\mu m^{-1}$]')
+ ax.set_xlabel(r'$\lambda$ [$\mu$m]')
  ax.set_ylabel('$A(x)/A(V)$')
 
  ax.set_title('R(V) = 3.1')
@@ -157,10 +169,12 @@ R(V) (+ other variables) dependent prediction models
 
  import numpy as np
  import matplotlib.pyplot as plt
+ from matplotlib.ticker import ScalarFormatter
  import astropy.units as u
 
  from dust_extinction.parameter_averages import (CCM89, O94, F99, F04,
- VCG04, GCC09, M14, F19, D22)
+ VCG04, GCC09, M14, F19, D22,
+ G23)
 
  fig, ax = plt.subplots()
 
@@ -169,17 +183,20 @@ R(V) (+ other variables) dependent prediction models
 
  Rv = 2.5
 
- models = [CCM89, O94, F99, F04, VCG04, GCC09, M14, F19, D22]
+ models = [CCM89, O94, F99, F04, VCG04, GCC09, M14, F19, D22, G23]
 
  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.plot(1./x[indxs], yvals, label=ext_model.__class__.__name__)
+
+ ax.set_xscale('log')
+ ax.xaxis.set_major_formatter(ScalarFormatter())
 
- ax.set_xlabel('$x$ [$\mu m^{-1}$]')
+ ax.set_xlabel(r'$\lambda$ [$\mu$m]')
  ax.set_ylabel('$A(x)/A(V)$')
 
  ax.set_title('R(V) = 2.5')
@@ -193,10 +210,12 @@ R(V) (+ other variables) dependent prediction models
 
  import numpy as np
  import matplotlib.pyplot as plt
+ from matplotlib.ticker import ScalarFormatter
  import astropy.units as u
 
  from dust_extinction.parameter_averages import (CCM89, O94, F99, F04,
- VCG04, GCC09, M14, F19, D22)
+ VCG04, GCC09, M14, F19, D22,
+ G23)
 
  fig, ax = plt.subplots()
 
@@ -205,17 +224,20 @@ R(V) (+ other variables) dependent prediction models
 
  Rv = 5.5
 
- models = [CCM89, O94, F99, F04, VCG04, GCC09, M14, F19, D22]
+ models = [CCM89, O94, F99, F04, VCG04, GCC09, M14, F19, D22, G23]
 
  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.plot(1./x[indxs], yvals, label=ext_model.__class__.__name__)
+
+ ax.set_xscale('log')
+ ax.xaxis.set_major_formatter(ScalarFormatter())
 
- ax.set_xlabel('$x$ [$\mu m^{-1}$]')
+ ax.set_xlabel(r'$\lambda$ [$\mu$m]')
  ax.set_ylabel('$A(x)/A(V)$')
 
  ax.set_title('R(V) = 5.5')
@@ -228,6 +250,7 @@ R(V) (+ other variables) dependent prediction models
 
  import numpy as np
  import matplotlib.pyplot as plt
+ from matplotlib.ticker import ScalarFormatter
  import astropy.units as u
 
  from dust_extinction.parameter_averages import G16
@@ -243,9 +266,12 @@ R(V) (+ other variables) dependent prediction models
  Rvs = [2.0, 3.0, 4.0, 5.0, 6.0]
  for cur_Rv in Rvs:
  ext_model = G16(RvA=cur_Rv, fA=1.0)
- ax.plot(x,ext_model(x),label=r'$R_A(V) = ' + str(cur_Rv) + '$')
+ ax.plot(1./x,ext_model(x),label=r'$R_A(V) = ' + str(cur_Rv) + '$')
+
+ ax.set_xscale('log')
+ ax.xaxis.set_major_formatter(ScalarFormatter())
 
- ax.set_xlabel('$x$ [$\mu m^{-1}$]')
+ ax.set_xlabel(r'$\lambda$ [$\mu$m]')
  ax.set_ylabel('$A(x)/A(V)$')
 
  ax.set_title('G16; $f_A = 1.0$; $R(V)_A$ variable')
@@ -258,6 +284,7 @@ R(V) (+ other variables) dependent prediction models
 
  import numpy as np
  import matplotlib.pyplot as plt
+ from matplotlib.ticker import ScalarFormatter
  import astropy.units as u
 
  from dust_extinction.parameter_averages import G16
@@ -273,9 +300,12 @@ R(V) (+ other variables) dependent prediction models
  fAs = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
  for cur_fA in fAs:
  ext_model = G16(RvA=3.1, fA=cur_fA)
- ax.plot(x,ext_model(x),label=r'$f_A = ' + str(cur_fA) + '$')
+ ax.plot(1./x,ext_model(x),label=r'$f_A = ' + str(cur_fA) + '$')
+
+ ax.set_xscale('log')
+ ax.xaxis.set_major_formatter(ScalarFormatter())
 
- ax.set_xlabel('$x$ [$\mu m^{-1}$]')
+ ax.set_xlabel(r'$\lambda$ [$\mu$m]')
  ax.set_ylabel('$A(x)/A(V)$')
 
  ax.set_title('G16; $f_A$ variable; $R(V)_A = 3.1$')
@@ -346,6 +376,7 @@ Grain models
 
  import numpy as np
  import matplotlib.pyplot as plt
+ from matplotlib.ticker import ScalarFormatter
  import astropy.units as u
 
  from dust_extinction.grain_models import DBP90, WD01, D03, ZDA04, C11, J13
@@ -381,6 +412,7 @@ Grain models
  ax.set_title('Grain Models')
 
  ax.set_xscale('log')
+ ax.xaxis.set_major_formatter(ScalarFormatter())
  ax.set_yscale('log')
 
  ax.set_title('Milky Way - Ultraviolet to Mid-Infrared')
@@ -444,6 +476,7 @@ Shape fitting models
 
  import numpy as np
  import matplotlib.pyplot as plt
+ from matplotlib.ticker import ScalarFormatter
  import astropy.units as u
 
  from dust_extinction.shapes import FM90
@@ -454,18 +487,21 @@ Shape fitting models
  x = np.arange(3.8,11.0,0.1)/u.micron
 
  ext_model = FM90()
- ax.plot(x,ext_model(x),label='total')
+ ax.plot(1./x,ext_model(x),label='total')
 
  ext_model = FM90(C3=0.0, C4=0.0)
- ax.plot(x,ext_model(x),label='linear term')
+ ax.plot(1./x,ext_model(x),label='linear term')
 
  ext_model = FM90(C1=0.0, C2=0.0, C4=0.0)
- ax.plot(x,ext_model(x),label='bump term')
+ ax.plot(1./x,ext_model(x),label='bump term')
 
  ext_model = FM90(C1=0.0, C2=0.0, C3=0.0)
- ax.plot(x,ext_model(x),label='FUV rise term')
+ ax.plot(1./x,ext_model(x),label='FUV rise term')
 
- ax.set_xlabel('$x$ [$\mu m^{-1}$]')
+ ax.set_xscale('log')
+ ax.xaxis.set_major_formatter(ScalarFormatter())
+ ax.xaxis.set_minor_formatter(ScalarFormatter())
+ ax.set_xlabel(r'$\lambda$ [$\mu$m]')
  ax.set_ylabel('$E(\lambda - V)/E(B - V)$')
 
  ax.set_title('FM90')
@@ -533,6 +569,7 @@ Shape fitting models
 
  import numpy as np
  import matplotlib.pyplot as plt
+ from matplotlib.ticker import ScalarFormatter
  import astropy.units as u
 
  from dust_extinction.shapes import G21
@@ -556,6 +593,7 @@ Shape fitting models
  ax.plot(1./x,ext_model(x),label='power-law+sil2 only')
 
  ax.set_xscale('log')
+ ax.xaxis.set_major_formatter(ScalarFormatter())
  ax.set_yscale('log')
 
  ax.set_xlabel('$\lambda$ [$\mu$m]')

dust_extinction/helpers.py

@@ -1,11 +1,10 @@
-from __future__ import absolute_import, print_function, division
-
 import warnings
 
 import numpy as np
+from scipy.special import comb
 import astropy.units as u
 
-__all__ = ["_get_x_in_wavenumbers", "_test_valid_x_range"]
+__all__ = ["_get_x_in_wavenumbers", "_test_valid_x_range", "_smoothstep"]
 
 
 def _get_x_in_wavenumbers(in_x):
@@ -67,3 +66,15 @@ def _test_valid_x_range(x, x_range, outname):
  + str(x_range[1])
  + ", x has units 1/micron]"
  )
+
+
+def _smoothstep(x, x_min=0, x_max=1, N=1):
+ x = np.clip((x - x_min) / (x_max - x_min), 0, 1)
+
+ result = 0
+ for n in range(0, N + 1):
+ result += comb(N + n, n) * comb(2 * N + 1, N - n) * (-x) ** n
+
+ result *= x ** (N + 1)
+
+ return result