Add G09_MWAvg model

docs/dust_extinction/choose_model.rst

@@ -13,7 +13,8 @@ Simple Average Curves
 ---------------------
 
 These are straightforward averages of observed extinction curves.  They are the
-simpliest models and include models for the Mw (TBA), the LMC
+simpliest models and include models for the MW
+(:class:`~dust_extinction.dust_extinction.G09_MWAvg`), the LMC
 (:class:`~dust_extinction.dust_extinction.G03_LMCAvg`,
 :class:`~dust_extinction.dust_extinction.G03_LMC2`) and the SMC
 (:class:`~dust_extinction.dust_extinction.G03_SMCBar`).
@@ -25,9 +26,11 @@ Way, the usual average used is R(V) = 3.1.
 +------------+-------------+------------------+--------------+
 | Model      | x range     | wavelength range |       galaxy | 
 +============+=============+==================+==============+
+| G09_MWAvg  | 0.3 - 10.96 |     0.0912 - 3.3 |           MW |
++------------+-------------+------------------+--------------+
 | G03_LMCAvg |  0.3 - 10.0 |        0.1 - 3.3 |          LMC |
 +------------+-------------+------------------+--------------+
-| G03_LMC2   |  0.3 - 10.0 |        0.1 - 3.3 | 30 Dor (LMC) |
+| G03_LMC2   |  0.3 - 10.0 |        0.1 - 3.3 | LMC (30 Dor) |
 +------------+-------------+------------------+--------------+
 | G03_SMCBar |  0.3 - 10.0 |        0.1 - 3.3 |          SMC |
 +------------+-------------+------------------+--------------+

docs/dust_extinction/fit_extinction.rst

@@ -23,7 +23,7 @@ extinction curve for the LMC outside of the LMC2 supershell region
 
    from astropy.modeling.fitting import LevMarLSQFitter
 
-   from dust_extinction.dust_extinction import G03_LMCAvg, FM90
+   from dust_extinction.dust_extinction import (G03_LMCAvg, FM90)
 
    # get an observed extinction curve to fit
    g03_model = G03_LMCAvg()
@@ -64,51 +64,65 @@ Example: P92 Fit
 ================
 
 In this example, the P92 model is used to fit the observed average
-extinction curve for the MW as tabulted by Pei (1992). 
+extinction curve for the MW (G09_MWAvg ``dust_extinction`` model).
+The fit is done using the observed uncertainties that are passed
+as weights.  The weights assume the noise is Gaussian and not correlated
+between data points.
 
 .. plot::
    :include-source:
 
    import matplotlib.pyplot as plt
    import numpy as np
 
+   import warnings
+   from astropy.utils.exceptions import AstropyWarning
+
    from astropy.modeling.fitting import LevMarLSQFitter
 
-   from dust_extinction.dust_extinction import P92
-
-   # Milky Way observed extinction as tabulated by Pei (1992)
-   MW_x = [0.21, 0.29, 0.45, 0.61, 0.80, 1.11, 1.43, 1.82,
-           2.27, 2.50, 2.91, 3.65, 4.00, 4.17, 4.35, 4.57, 4.76,
-           5.00, 5.26, 5.56, 5.88, 6.25, 6.71, 7.18, 7.60,
-           8.00, 8.50, 9.00, 9.50, 10.00]
-   MW_x = np.array(MW_x)
-   MW_exvebv = [-3.02, -2.91, -2.76, -2.58, -2.23, -1.60, -0.78, 0.00,
-                1.00, 1.30, 1.80, 3.10, 4.19, 4.90, 5.77, 6.57, 6.23,
-                5.52, 4.90, 4.65, 4.60, 4.73, 4.99, 5.36, 5.91, 
-                6.55, 7.45, 8.45, 9.80, 11.30]
-   MW_exvebv = np.array(MW_exvebv)
-   Rv = 3.08
-   MW_axav = MW_exvebv/Rv + 1.0
+   from dust_extinction.dust_extinction import (P92, G09_MWAvg)
+
+   # get an observed extinction curve to fit
+   g09_model = G09_MWAvg()
 
    # get an observed extinction curve to fit
-   x = MW_x   
-   y = MW_axav
+   x = g09_model.obsdata_x
+   y = g09_model.obsdata_axav
+   y_unc = g09_model.obsdata_axav_unc
 
    # initialize the model
    p92_init = P92()
 
+   # fix a number of the parameters
+   #   mainly to avoid fitting parameters that are constrained at
+   #   wavelengths where the observed data for this case does not exist
+   p92_init.FUV_lambda.fixed = True
+   p92_init.SIL1_amp.fixed = True
+   p92_init.SIL1_lambda.fixed = True
+   p92_init.SIL1_b.fixed = True
+   p92_init.SIL2_amp.fixed = True
+   p92_init.SIL2_lambda.fixed = True
+   p92_init.SIL2_b.fixed = True
+   p92_init.FIR_amp.fixed = True
+   p92_init.FIR_lambda.fixed = True
+   p92_init.FIR_b.fixed = True
+
    # pick the fitter
    fit = LevMarLSQFitter()
 
+   # set to avoid the "fit may have been unsuccessful" warning
+   #   fit is fine, but this means the build of the docs fails
+   warnings.simplefilter('ignore', category=AstropyWarning)
+
    # fit the data to the P92 model using the fitter
    #   use the initialized model as the starting point
    #   accuracy set to avoid warning the fit may have failed
-   p92_fit = fit(p92_init, x, y, acc=1e-3)
+   p92_fit = fit(p92_init, x, y, weights=1.0/y_unc)
 
    # plot the observed data, initial guess, and final fit
    fig, ax = plt.subplots()
 
-   ax.plot(x, y, 'ko', label='Observed Curve')
+   ax.errorbar(x, y, yerr=y_unc, fmt='ko', label='Observed Curve')
    ax.plot(x, p92_init(x), label='Initial guess')
    ax.plot(x, p92_fit(x), label='Fitted model')
 

docs/dust_extinction/model_flavors.rst

@@ -10,18 +10,19 @@ Average models
 
    These models provide averages from the literature with the ability to
    interpolate between the observed data points.
-   For the Milky Way, one of the R(V) dependent models  with R(V) = 3.1
-   (see next section) are often used for the Milky Way 'average'.
-   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).
+   Models are provided for the Milky Way (Gordon, Cartlege, & Clayton 2009)
+   and the Magellanic Clouds (Gordon et al. 2003).
+
+   For the Milky Way, one of the R(V) dependent models with R(V) = 3.1
+   (see next section) can also be used for the Milky Way 'average'.
 
 .. plot::
 
    import numpy as np
    import matplotlib.pyplot as plt
    import astropy.units as u
 
-   from dust_extinction.dust_extinction import (F99,
+   from dust_extinction.dust_extinction import (G09_MWAvg,
                                                 G03_SMCBar,
                                                 G03_LMCAvg,
 					        G03_LMC2)
@@ -31,8 +32,8 @@ Average models
    # generate the curves and plot them
    x = np.arange(0.3,10.0,0.1)/u.micron
 
-   ext_model = F99()
-   ax.plot(x,ext_model(x),label='MW Average (F99 w/ $R(V)=3.1$)')
+   ext_model = G09_MWAvg()
+   ax.plot(x,ext_model(x),label='G09 MWAvg')
 
    ext_model = G03_SMCBar()
    ax.plot(x,ext_model(x),label='G03 SMCBar')

docs/index.rst

@@ -78,8 +78,6 @@ contributors page on Github
 
 Reference API
 =============
-.. toctree::
-   :maxdepth: 1
 
-   dust_extinction/reference_api.rst
+.. automodapi:: dust_extinction.dust_extinction