karllark · karllark · Nov 9, 2017 · Oct 27, 2017 · Oct 27, 2017 · Nov 7, 2017
diff --git a/.travis.yml b/.travis.yml
@@ -106,6 +106,8 @@ matrix:
           env: PYTHON_VERSION=3.5 NUMPY_VERSION=1.11
         - os: linux
           env: NUMPY_VERSION=1.12
+        - os: linux
+          env: NUMPY_VERSION=1.13
 
         # Try numpy pre-release
         - os: linux
@@ -117,6 +119,14 @@ matrix:
           env: MAIN_CMD='pycodestyle dust_extinction --count' SETUP_CMD=''
 
     allow_failures:
+        - os: linux
+          env: ASTROPY_VERSION=development
+               EVENT_TYPE='pull_request push cron'
+
+        - os: linux
+          env: NUMPY_VERSION=prerelease
+               EVENT_TYPE='pull_request push cron'
+
         # Do a PEP8 test with pycodestyle
         # (allow to fail unless your code completely compliant)
         - os: linux

diff --git a/docs/conf.py b/docs/conf.py
@@ -105,8 +105,8 @@
 
 # Please update these texts to match the name of your package.
 html_theme_options = {
-    'logotext1': 'package',  # white,  semi-bold
-    'logotext2': '-template',  # orange, light
+    'logotext1': 'dust_',  # white,  semi-bold
+    'logotext2': 'extinction',  # orange, light
     'logotext3': ':docs'   # white,  light
     }
 
@@ -131,6 +131,10 @@
 # pixels large.
 #html_favicon = ''
 
+#  example from imexam
+#from os.path import join
+#html_favicon = join('_static', 'imexam.ico')
+
 # If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
 # using the given strftime format.
 #html_last_updated_fmt = ''
@@ -142,6 +146,13 @@
 # Output file base name for HTML help builder.
 htmlhelp_basename = project + 'doc'
 
+# more examples from imexam
+#html_logo = '_static/imexam_logo_trans.png'
+
+# Static files to copy after template files
+#html_static_path = ['_static']
+#html_style = 'imexam.css'
+
 
 # -- Options for LaTeX output -------------------------------------------------
 

diff --git a/docs/dust_extinction/extinguish.rst b/docs/dust_extinction/extinguish.rst
@@ -0,0 +1,66 @@
+###############################
+Extinguish or Unextinguish Data
+###############################
+
+Two of the three flavors of models include a function to calculate the 
+factor to multiple (extinguish) or divide (unextinguish) a spectrum by
+to add or remove the effects of dust, respectively.  
+
+Extinguish is also often called reddening.  Extinguishing a spectrum often
+reddens the flux, but not always (e.g, on the short wavelength side of the 
+2175 A bump.  So extinguish is the more generic term.
+
+Example: Extinguish a Blackbody
+===============================
+
+.. plot::
+   :include-source:
+
+   import matplotlib.pyplot as plt
+   import numpy as np
+
+   import astropy.units as u
+   from astropy.modeling.blackbody import blackbody_lambda
+
+   from dust_extinction.dust_extinction import CCM89
+
+   # generate wavelengths between 0.1 and 3 microns
+   #    within the valid range for the CCM89 R(V) dependent relationship
+   lam = np.logspace(np.log10(0.1), np.log10(3.0), num=1000)
+
+   # setup the inputs for the blackbody function
+   wavelengths = lam*1e4*u.AA
+   temperature = 10000*u.K
+
+   # get the blackbody flux
+   flux = blackbody_lambda(wavelengths, temperature)
+
+   # initialize the model
+   ext = CCM89(Rv=3.1)
+
+   # get the extinguished blackbody flux for different amounts of dust
+   flux_ext_av05 = flux*ext.extinguish(wavelengths, Av=0.5)
+   flux_ext_av15 = flux*ext.extinguish(wavelengths, Av=1.5)
+   flux_ext_ebv10 = flux*ext.extinguish(wavelengths, Ebv=1.0)
+
+   # plot the intrinsic and extinguished fluxes
+   fig, ax = plt.subplots()
+
+   ax.plot(wavelengths, flux, label='Intrinsic')
+   ax.plot(wavelengths, flux_ext_av05, label='$A(V) = 0.5$')
+   ax.plot(wavelengths, flux_ext_av15, label='$A(V) = 1.5$')
+   ax.plot(wavelengths, flux_ext_ebv10, label='$E(B-V) = 1.0$')
+
+   ax.set_xlabel('$\lambda$ [$\AA$]')
+   ax.set_ylabel('$Flux$')
+
+   ax.set_xscale('log')
+   ax.set_yscale('log')
+
+   ax.set_title('Example extinguishing a blackbody')
+
+   ax.legend(loc='best')
+   plt.tight_layout()
+   plt.show()
+
+
diff --git a/docs/dust_extinction/fit_extinction.rst b/docs/dust_extinction/fit_extinction.rst
@@ -0,0 +1,122 @@
+#####################
+Fit Extinction Curves
+#####################
+
+The ``dust_extinction`` package is built on the `astropy.modeling
+<http://docs.astropy.org/en/stable/modeling/>`_ package.  Fitting is
+done in the standard way for this package where the model is initialized
+with a starting point (either the default or user input), the fitter
+is choosen, and the fit performed.  
+
+Example: FM90 Fit
+=================
+
+In this example, the FM90 model is used to fit the observed average
+extinction curve for the LMC outside of the LMC2 supershell region 
+(G03_LMCAvg ``dust_extinction`` model). 
+
+.. plot::
+   :include-source:
+
+   import matplotlib.pyplot as plt
+   import numpy as np
+
+   from astropy.modeling.fitting import LevMarLSQFitter
+
+   from dust_extinction.dust_extinction import G03_LMCAvg, FM90
+
+   # get an observed extinction curve to fit
+   g03_model = G03_LMCAvg()
+
+   x = g03_model.obsdata_x
+   # convert to E(x-V)/E(B0V)
+   y = (g03_model.obsdata_axav - 1.0)*g03_model.Rv
+   # only fit the UV portion (FM90 only valid in UV)
+   gindxs, = np.where(x > 3.125)
+
+   # initialize the model
+   fm90_init = FM90()
+
+   # pick the fitter
+   fit = LevMarLSQFitter()
+
+   # fit the data to the FM90 model using the fitter
+   #   use the initialized model as the starting point
+   g03_fit = fit(fm90_init, x[gindxs], y[gindxs])
+
+   # plot the observed data, initial guess, and final fit
+   fig, ax = plt.subplots()
+
+   ax.plot(x, y, 'ko', label='Observed Curve')
+   ax.plot(x[gindxs], fm90_init(x[gindxs]), label='Initial guess')
+   ax.plot(x[gindxs], g03_fit(x[gindxs]), label='Fitted model')
+
+   ax.set_xlabel('$x$ [$\mu m^{-1}$]')
+   ax.set_ylabel('$E(x-V)/E(B-V)$')
+
+   ax.set_title('Example FM90 Fit to G03_LMCAvg curve')
+
+   ax.legend(loc='best')
+   plt.tight_layout()
+   plt.show()
+
+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). 
+
+.. plot::
+   :include-source:
+
+   import matplotlib.pyplot as plt
+   import numpy as np
+
+   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
+
+   # get an observed extinction curve to fit
+   x = MW_x   
+   y = MW_axav
+
+   # initialize the model
+   p92_init = P92()
+
+   # pick the fitter
+   fit = LevMarLSQFitter()
+
+   # 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)
+
+   # plot the observed data, initial guess, and final fit
+   fig, ax = plt.subplots()
+
+   ax.plot(x, y, 'ko', label='Observed Curve')
+   ax.plot(x, p92_init(x), label='Initial guess')
+   ax.plot(x, p92_fit(x), label='Fitted model')
+
+   ax.set_xlabel('$x$ [$\mu m^{-1}$]')
+   ax.set_ylabel('$A(x)/A(V)$')
+
+   ax.set_title('Example P92 Fit to MW average curve')
+
+   ax.legend(loc='best')
+   plt.tight_layout()
+   plt.show()
diff --git a/docs/dust_extinction/index.rst b/docs/dust_extinction/index.rst
diff --git a/docs/dust_extinction/model_flavors.rst b/docs/dust_extinction/model_flavors.rst
@@ -0,0 +1,280 @@
+#############
+Model Flavors
+#############
+
+There are three differnet types of models: average, R(V)+ dependent prediction,
+and shape fitting.
+
+Average models
+==============
+
+   These models provide averages from the literature with the ability to
+   interpolate between the observed data points.
+   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).
+
+.. plot::
+
+   import numpy as np
+   import matplotlib.pyplot as plt
+   import astropy.units as u
+
+   from dust_extinction.dust_extinction import (G03_SMCBar,
+                                                G03_LMCAvg,
+					        G03_LMC2)
+
+   fig, ax = plt.subplots()
+
+   # generate the curves and plot them
+   x = np.arange(0.3,10.0,0.1)/u.micron
+
+   ext_model = G03_SMCBar()
+   ax.plot(x,ext_model(x),label='G03 SMCBar')
+
+   ext_model = G03_LMCAvg()
+   ax.plot(x,ext_model(x),label='G03 LMCAvg')
+
+   ext_model = G03_LMC2()
+   ax.plot(x,ext_model(x),label='G03 LMC2')
+
+   ax.set_xlabel('$x$ [$\mu m^{-1}$]')
+   ax.set_ylabel('$A(x)/A(V)$')
+
+   ax.legend(loc='best')
+   plt.tight_layout()
+   plt.show()
+
+R(V) (+ other variables) dependent prediction models
+====================================================
+
+   These models provide predictions of the shape of the dust extinction
+   given input variable(s).  
+
+   These include CCM89 the original R(V) dependent model from 
+   Cardelli, Clayton, and Mathis (1989) and updated versions F99 
+   (Fitzpatrick 1999).  These models are based on the average
+   behavior of extinction in the Milky Way.
+
+   In addition, the (R(V), f_A) two parameter relationship from 
+   Gordon et al. (2016) is included.  This model is based on the average
+   behavior of extinction in the Milky Way, Large Magellanic Cloud, and 
+   Small Magellanic Cloud.
+
+.. plot::
+
+   import numpy as np
+   import matplotlib.pyplot as plt
+   import astropy.units as u
+
+   from dust_extinction.dust_extinction import CCM89
+
+   fig, ax = plt.subplots()
+
+   # 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.set_title('CCM89')
+
+   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.dust_extinction import F99
+
+   fig, ax = plt.subplots()
+
+   # temp model to get the correct x range
+   text_model = F99()
+
+   # generate the curves and plot them
+   x = np.arange(text_model.x_range[0], text_model.x_range[1],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.set_title('F99')
+
+   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.dust_extinction import G16
+
+   fig, ax = plt.subplots()
+
+   # temp model to get the correct x range
+   text_model = G16()
+
+   # generate the curves and plot them
+   x = np.arange(text_model.x_range[0], text_model.x_range[1],0.1)/u.micron
+
+   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.set_xlabel('$x$ [$\mu m^{-1}$]')
+   ax.set_ylabel('$A(x)/A(V)$')
+
+   ax.set_title('G16; $f_A = 1.0$; $R(V)_A$ variable')
+
+   ax.legend(loc='best', title=r'$f_A = 1.0$')
+   plt.tight_layout()
+   plt.show()
+
+.. plot::
+
+   import numpy as np
+   import matplotlib.pyplot as plt
+   import astropy.units as u
+
+   from dust_extinction.dust_extinction import G16
+
+   fig, ax = plt.subplots()
+
+   # temp model to get the correct x range
+   text_model = G16()
+
+   # generate the curves and plot them
+   x = np.arange(text_model.x_range[0], text_model.x_range[1],0.1)/u.micron
+
+   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.set_xlabel('$x$ [$\mu m^{-1}$]')
+   ax.set_ylabel('$A(x)/A(V)$')
+
+   ax.set_title('G16; $f_A$ variable; $R(V)_A = 3.1$')
+
+   ax.legend(loc='best', title=r'$R_A(V) = 3.1$')
+   plt.tight_layout()
+   plt.show()
+
+
+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.
+   The P92 (Pei 1992) uses 19 parameters to fit the shape of the X-ray to 
+   far-infrared extinction.
+
+.. plot::
+
+   import numpy as np
+   import matplotlib.pyplot as plt
+   import astropy.units as u
+
+   from dust_extinction.dust_extinction import FM90
+
+   fig, ax = plt.subplots()
+
+   # 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')
+
+   ext_model = FM90(C3=0.0, C4=0.0)
+   ax.plot(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')
+
+   ext_model = FM90(C1=0.0, C2=0.0, C3=0.0)
+   ax.plot(x,ext_model(x),label='FUV rise term')
+
+   ax.set_xlabel('$x$ [$\mu m^{-1}$]')
+   ax.set_ylabel('$E(\lambda - V)/E(B - V)$')
+
+   ax.set_title('FM90')
+
+   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.dust_extinction import P92
+
+   fig, ax = plt.subplots()
+
+   # generate the curves and plot them
+   lam = np.logspace(-3.0, 3.0, num=1000)
+   x = (1.0/lam)/u.micron
+
+   ext_model = P92()
+   ax.plot(1/x,ext_model(x),label='total')
+
+   ext_model = P92(FUV_amp=0., NUV_amp=0.0, 
+                   SIL1_amp=0.0, SIL2_amp=0.0, FIR_amp=0.0)
+   ax.plot(1./x,ext_model(x),label='BKG only')
+
+   ext_model = P92(NUV_amp=0.0, 
+                   SIL1_amp=0.0, SIL2_amp=0.0, FIR_amp=0.0)
+   ax.plot(1./x,ext_model(x),label='BKG+FUV only')
+
+   ext_model = P92(FUV_amp=0.,
+                   SIL1_amp=0.0, SIL2_amp=0.0, FIR_amp=0.0)
+   ax.plot(1./x,ext_model(x),label='BKG+NUV only')
+
+   ext_model = P92(FUV_amp=0., NUV_amp=0.0, 
+                   SIL2_amp=0.0)
+   ax.plot(1./x,ext_model(x),label='BKG+FIR+SIL1 only')
+
+   ext_model = P92(FUV_amp=0., NUV_amp=0.0, 
+                   SIL1_amp=0.0)
+   ax.plot(1./x,ext_model(x),label='BKG+FIR+SIL2 only')
+
+   ext_model = P92(FUV_amp=0., NUV_amp=0.0, 
+                   SIL1_amp=0.0, SIL2_amp=0.0)
+   ax.plot(1./x,ext_model(x),label='BKG+FIR only')
+
+   ax.set_xscale('log')
+   ax.set_yscale('log')
+
+   ax.set_ylim(1e-3,10.)
+
+   ax.set_xlabel('$\lambda$ [$\mu$m]')
+   ax.set_ylabel('$A(x)/A(V)$')
+
+   ax.set_title('P92')
+
+   ax.legend(loc='best')
+   plt.tight_layout()
+   plt.show()
+
+
+
diff --git a/docs/dust_extinction/reference_api.rst b/docs/dust_extinction/reference_api.rst
@@ -0,0 +1,2 @@
+.. automodapi:: dust_extinction.dust_extinction
+
diff --git a/docs/index.rst b/docs/index.rst
@@ -1,17 +1,82 @@
-Documentation
-=============
+############################
+Interstellar Dust Extinction
+############################
+
+``dust_extinction`` is a python package to provide interstellar dust extinction
+curves.
+
+While there are other python packages that provide some of the extinction
+curves contained here, the explicit motivation for this package is to provide
+extinction curves to those using them to model/correct their data and those
+studying extinction curves directly to better undertand interstellar dust.
+
+This package is developed in the 
+`astropy affiliated package <http://www.astropy.org/affiliated/>`_ 
+template and uses the 
+`astropy.modeling <http://docs.astropy.org/en/stable/modeling/>`_ 
+framework.
+
+Installation
+============
+
+To be added
+
+User Documenation
+=================
+
+.. toctree::
+   :maxdepth: 2
+
+   Flavors of Models <dust_extinction/model_flavors.rst>
+   Extinguish (or unextinguish) data <dust_extinction/extinguish.rst>
+   Fitting extinction curves <dust_extinction/fit_extinction.rst>
+
+
+Reporting Issues
+================
+
+If you have found a bug in ``dust_extinction`` please report it by creating a
+new issue on the ``dust_extinction`` `GitHub issue tracker
+<https://github.com/karllark/dust_extinction/issues>`_.
 
-This is an affiliated package for the AstroPy package. The documentation for
-this package is here:
+Please include an example that demonstrates the issue sufficiently so that
+the developers can reproduce and fix the problem. You may also be asked to
+provide information about your operating system and a full Python
+stack trace.  The developers will walk you through obtaining a stack
+trace if it is necessary.
 
+Contributing
+============
+
+Like the `Astropy`_ project, ``dust_extinction`` is made both by and for its
+users.  We accept contributions at all levels, spanning the gamut from
+fixing a typo in the documentation to developing a major new feature.
+We welcome contributors who will abide by the `Python Software
+Foundation Code of Conduct
+<https://www.python.org/psf/codeofconduct/>`_.
+
+``dust_extinction`` follows the same workflow and coding guidelines as
+`Astropy`_.  The following pages will help you get started with
+contributing fixes, code, or documentation (no git or GitHub
+experience necessary):
+
+* `How to make a code contribution <http://astropy.readthedocs.io/en/stable/development/workflow/development_workflow.html>`_
+
+* `Coding Guidelines <http://docs.astropy.io/en/latest/development/codeguide.html>`_
+
+* `Try the development version <http://astropy.readthedocs.io/en/stable/development/workflow/get_devel_version.html>`_
+
+* `Developer Documentation <http://docs.astropy.org/en/latest/#developer-documentation>`_
+
+
+For the complete list of contributors please see the `dust_extinction
+contributors page on Github
+<https://github.com/karllark/dust_extinction/graphs/contributors>`_.
+
+Reference API
+=============
 .. toctree::
-  :maxdepth: 2
+   :maxdepth: 1
 
-  dust_extinction/index.rst
+   dust_extinction/reference_api.rst
 
-.. note:: The layout of this directory is simply a suggestion.  To follow
-          traditional practice, do *not* edit this page, but instead place
-          all documentation for the affiliated package inside ``packagename/``.
-          The traditional practice was intended to allow the affiliated
-          package to eventually be merged into the main astropy package.
-          You can follow this practice or choose your own layout.
diff --git a/dust_extinction/tests/test_p92.py b/dust_extinction/tests/test_p92.py
@@ -99,17 +99,20 @@ def test_P92_fitting():
 
     x = x_quant.value
 
+    # change from defaults to make the best fit harder to find
     p92_init = P92()
+
     fit = LevMarLSQFitter()
-    p92_fit = fit(p92_init, x, y)
+    # accuracy set to avoid warning the fit may have failed
+    p92_fit = fit(p92_init, x, y, acc=1e-3)
 
     fit_vals = p92_fit._parameters
 
-    good_vals = [222.840030742, 0.0468548209216, 87.5860789967, 2.0,
-                 17.5852152025, 0.0754660056762, 2.08487806377, 6.13359917075,
-                 0.0496273802215, 0.218586483035, -1.95320679992, 2.0,
-                 0.175702371244, 13.0, 8.82541783881, 2.0,
-                 0.0, 15.0, -136.539895555, 2.0,
-                 0.0, 20.0, -2174.83099919, 2.0]
-    
+    good_vals = [222.820720554, 0.0468116978742, 87.5683270158, 2.0, 
+                 18.1622695074, 0.0770527204445, 2.83090815471, 6.40496987921, 
+                 0.0518566585953, 0.218641802256, -1.95087797308, 2.0, 
+                 0.00600017900912, 13.0, 109.369888618, 2.0, 
+                 0.09441454528, 15.0, -251.33735749, 2.0, 
+                 0.189219789029, 20.0, -1130.90875587, 2.0]
+
     np.testing.assert_allclose(good_vals, fit_vals)
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1,2 @@
		.. automodapi:: dust_extinction.dust_extinction