diff --git a/docs/_static/specreduce.css b/docs/_static/specreduce.css
index 229d8787..e69de29b 100644
--- a/docs/_static/specreduce.css
+++ b/docs/_static/specreduce.css
@@ -1,15 +0,0 @@
-@import url("bootstrap-astropy.css");
-
-div.topbar a.brand {
-  background: transparent url("logo_icon.png") no-repeat 8px 3px;
-  background-image: url("logo_icon.png"), none;
-  background-size: 32px 32px;
-}
-
-#logotext1 {
-  color: #519EA8;
-}
-
-#logotext2 {
-  color: #FF5000;
-}
diff --git a/docs/api.rst b/docs/api.rst
index e4b4de6b..a6fa6948 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -1,7 +1,7 @@
 .. _api_index:
 
-API Index
-=========
+API
+===
 
 .. automodapi:: specreduce
     :no-inheritance-diagram:
diff --git a/docs/background.rst b/docs/background.rst
new file mode 100644
index 00000000..28dc7c46
--- /dev/null
+++ b/docs/background.rst
@@ -0,0 +1,29 @@
+Background correction
+=====================
+
+The `specreduce.background` module generates and subtracts a background image from
+the input 2D spectral image.  The `~specreduce.background.Background` object is
+defined by one or more windows, and can be generated with:
+
+* `~specreduce.background.Background`
+* `Background.one_sided <specreduce.background.Background.one_sided>`
+* `Background.two_sided <specreduce.background.Background.two_sided>`
+
+The center of the window can either be passed as a float/integer or as a trace
+
+.. code-block:: python
+
+  bg = specreduce.background.Background.one_sided(image, trace, separation=5, width=2)
+
+or, equivalently
+
+.. code-block:: python
+
+  bg = specreduce.background.Background.one_sided(image, 15, separation=5, width=2)
+
+The background image can be accessed via `~specreduce.background.Background.bkg_image`
+and the background-subtracted image via `~specreduce.background.Background.sub_image`
+(or ``image - bg``).
+
+The background and trace steps can be done iteratively, to refine an automated
+trace using the background-subtracted image as input.
\ No newline at end of file
diff --git a/docs/conf.py b/docs/conf.py
index cf78023f..5079e8be 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -33,7 +33,7 @@
 from specreduce import __version__
 
 try:
-    from sphinx_astropy.conf.v1 import *  # noqa
+    from sphinx_astropy.conf.v2 import *  # noqa
 except ImportError:
     print('ERROR: the documentation requires the sphinx-astropy package to be installed')
     sys.exit(1)
@@ -52,7 +52,7 @@
 highlight_language = 'python3'
 
 # If your documentation needs a minimal Sphinx version, state it here.
-#needs_sphinx = '1.2'
+needs_sphinx = '3.0'
 
 # To perform a Sphinx version check that needs to be more specific than
 # major.minor, call `check_sphinx_version("x.y.z")` here.
@@ -87,40 +87,46 @@
 
 # -- Options for HTML output --------------------------------------------------
 
-# A NOTE ON HTML THEMES
-# The global astropy configuration uses a custom theme, 'bootstrap-astropy',
-# which is installed along with astropy. A different theme can be used or
-# the options for this theme can be modified by overriding some of the
-# variables set in the global configuration. The variables set in the
-# global configuration are listed below, commented out.
-
-
-# Add any paths that contain custom themes here, relative to this directory.
-# To use a different custom theme, add the directory containing the theme.
-#html_theme_path = []
-
-# The theme to use for HTML and HTML Help pages.  See the documentation for
-# a list of builtin themes. To override the custom theme, set this to the
-# name of a builtin theme or the name of a custom theme in html_theme_path.
-#html_theme = None
 html_static_path = ['_static']  # html_theme = None
 html_style = 'specreduce.css'
 
-
-html_theme_options = {
-    'logotext1': 'spec',  # white,  semi-bold
-    'logotext2': 'reduce',  # orange, light
-    'logotext3': ':docs'   # white,  light
+html_theme_options.update(
+    {
+        "github_url": "https://github.com/astropy/specreduce",
+        "use_edit_page_button": False,
+        "navigation_with_keys": False,
+        "logo": {
+            "text": f"{project}",
+            "image_light": "_static/logo_icon.png",
+            "image_dark": "_static/logo_icon.png",
+        },
     }
+)
+
+html_context = {
+    "default_mode": "light",
+    "version_slug": os.environ.get("READTHEDOCS_VERSION") or "",
+    "to_be_indexed": ["stable", "latest"],
+    "github_user": "astropy",
+    "github_repo": "specreduce",
+    "github_version": "main",
+    "doc_path": "docs",
+    "edit_page_url_template": "{{ astropy_custom_edit_url(github_user, github_repo, github_version, doc_path, file_name, default_edit_page_url_template) }}",
+    "default_edit_page_url_template": "https://github.com/{github_user}/{github_repo}/edit/{github_version}/{doc_path}{file_name}",
+    # Tell Jinja2 templates the build is running on Read the Docs
+    "READTHEDOCS": os.environ.get("READTHEDOCS", "") == "True",
+}
+
+#html_theme_options = {
+#    'logotext1': 'spec',  # white,  semi-bold
+#    'logotext2': 'reduce',  # orange, light
+#    'logotext3': ':docs'   # white,  light
+#    }
 
 # Custom sidebar templates, maps document names to template names.
 #html_sidebars = {}
-html_sidebars['**'] = ['localtoc.html']
-html_sidebars['index'] = ['globaltoc.html', 'localtoc.html']
-
-# The name of an image file (relative to this directory) to place at the top
-# of the sidebar.
-#html_logo = ''
+#html_sidebars['**'] = ['localtoc.html']
+#html_sidebars['index'] = ['globaltoc.html', 'localtoc.html']
 
 # The name of an image file (within the static path) to use as favicon of the
 # docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
diff --git a/docs/contributing.rst b/docs/contributing.rst
new file mode 100644
index 00000000..af621423
--- /dev/null
+++ b/docs/contributing.rst
@@ -0,0 +1,11 @@
+Contributing
+============
+
+Please feel free to contribute code and suggestions through `GitHub. <https://github
+.com/astropy/specreduce>`_
+
+.. toctree::
+    :maxdepth: 1
+
+    process/index
+    terms
\ No newline at end of file
diff --git a/docs/extraction.rst b/docs/extraction.rst
new file mode 100644
index 00000000..de1e6d50
--- /dev/null
+++ b/docs/extraction.rst
@@ -0,0 +1,110 @@
+Spectrum Extraction
+===================
+
+The `specreduce.extract` module extracts a 1D spectrum from an input 2D spectrum
+(likely a background-extracted spectrum from the previous step) and a defined
+window, using one of the following implemented methods:
+
+* `~specreduce.extract.BoxcarExtract`
+* `~specreduce.extract.HorneExtract`
+
+Each of these takes the input image and trace as inputs (see the :ref:`api_index` for
+other required and optional parameters)
+
+.. code-block:: python
+
+  extract = specreduce.extract.BoxcarExtract(image-bg, trace, width=3)
+
+or
+
+.. code-block:: python
+
+  extract = specreduce.extract.HorneExtract(image-bg, trace)
+
+For the Horne algorithm, the variance array is required. If the input image is
+an `~astropy.nddata.NDData` object with ``image.uncertainty`` provided,
+then this will be used. Otherwise, the ``variance`` parameter must be set.
+
+.. code-block:: python
+
+  extract = specreduce.extract.HorneExtract(image-bg, trace, variance=var_array)
+
+An optional mask array for the image may be supplied to HorneExtract as well.
+This follows the same convention and can either be attached to ``image`` if it
+is an `~astropy.nddata.NDData` object, or supplied as a keyword argument.
+
+The extraction methods automatically detect non-finite pixels in the input
+image and combine them with the user-supplied mask to prevent them from biasing the
+extraction. In the boxcar extraction, the treatment of these pixels is controlled by
+the ``mask_treatment`` option. When set to ``exclude`` (the default), non-finite
+pixels within the extraction window are excluded from the extraction, and the extracted
+flux is scaled according to the effective number of unmasked pixels. When using other
+options (``filter`` or ``omit``), the non-finite values may be propagated or treated
+differently as documented in the API.
+
+The previous examples in this section show how to initialize the BoxcarExtract
+or HorneExtract objects with their required parameters. To extract the 1D
+spectrum
+
+.. code-block:: python
+
+  spectrum = extract.spectrum
+
+The ``extract`` object contains all the set options.  The extracted 1D spectrum
+can be accessed via the ``spectrum`` property or by calling (e.g ``extract()``)
+the ``extract`` object (which also allows temporarily overriding any values)
+
+.. code-block:: python
+
+  spectrum2 = extract(width=6)
+
+or, for example to override the original ``trace_object``
+
+.. code-block:: python
+
+  spectrum2 = extract(trace_object=new_trace)
+
+Spatial profile options
+-----------------------
+The Horne algorithm provides two options for fitting the spatial profile to the
+cross dispersion direction of the source: a Gaussian fit (default),
+or an empirical ``interpolated_profile`` option.
+
+If the default Gaussian option is used, an optional background model may be
+supplied as well (default is a 2D Polynomial) to account
+for residual background in the spatial profile. This option is not supported for
+``interpolated_profile``.
+
+If  the ``interpolated_profile`` option is used, the image will be sampled in various
+wavelength bins (set by ``n_bins_interpolated_profile``), averaged in those bins, and
+samples are then interpolated between (linear by default, interpolation degree can
+be set with ``interp_degree_interpolated_profile``, which defaults to linear in
+x and y) to generate an empirical interpolated spatial profile. Since this option
+has two optional parameters to control the fit, the input can either be a string
+to indicate that ``interpolated_profile`` should be used for the spatial profile
+and to use the defaults for bins and interpolation degree, or to override these
+defaults a dictionary can be passed in.
+
+For example, to use the ``interpolated_profile`` option with default bins and
+interpolation degree
+
+.. code-block:: python
+
+    interp_profile_extraction = extract(spatial_profile='interpolated_profile')
+
+Or, to override the default of 10 samples and use 20 samples
+
+.. code-block:: python
+
+    interp_profile_extraction = extract(spatial_profile={'name': 'interpolated_profile',
+                                    'n_bins_interpolated_profile': 20)
+
+Or, to do a cubic interpolation instead of the default linear
+
+.. code-block:: python
+
+    interp_profile_extraction = extract(spatial_profile={'name': 'interpolated_profile',
+                                    'interp_degree_interpolated_profile': 3)
+
+As usual, parameters can either be set when instantiating the HorneExtraxt object,
+or supplied/overridden when calling the extraction method on that object.
\ No newline at end of file
diff --git a/docs/extraction_quickstart.rst b/docs/extraction_quickstart.rst
index febd4198..48158407 100644
--- a/docs/extraction_quickstart.rst
+++ b/docs/extraction_quickstart.rst
@@ -24,7 +24,9 @@ Supported trace types include:
 
 Each of these trace classes takes the 2D spectral image as input, as well as
 additional information needed to define or determine the trace (see the API docs
-above for required parameters for each of the available trace classes)::
+above for required parameters for each of the available trace classes)
+
+.. code-block:: python
 
   trace = specreduce.tracing.FlatTrace(image, 15)
 
@@ -45,16 +47,18 @@ defined by one or more windows, and can be generated with:
 * `Background.one_sided <specreduce.background.Background.one_sided>`
 * `Background.two_sided <specreduce.background.Background.two_sided>`
 
-The center of the window can either be passed as a float/integer or as a trace::
+The center of the window can either be passed as a float/integer or as a trace
+
+.. code-block:: python
 
   bg = specreduce.background.Background.one_sided(image, trace, separation=5, width=2)
 
+or, equivalently
 
-or, equivalently::
+.. code-block:: python
 
   bg = specreduce.background.Background.one_sided(image, 15, separation=5, width=2)
 
-
 The background image can be accessed via `~specreduce.background.Background.bkg_image`
 and the background-subtracted image via `~specreduce.background.Background.sub_image`
 (or ``image - bg``).
@@ -73,23 +77,29 @@ window, using one of the following implemented methods:
 * `~specreduce.extract.HorneExtract`
 
 Each of these takes the input image and trace as inputs (see the API above for
-other required and optional parameters)::
+other required and optional parameters)
+
+.. code-block:: python
 
   extract = specreduce.extract.BoxcarExtract(image-bg, trace, width=3)
 
-or::
+or
+
+.. code-block:: python
 
   extract = specreduce.extract.HorneExtract(image-bg, trace)
 
 For the Horne algorithm, the variance array is required. If the input image is
-an ``astropy.NDData`` object with ``image.uncertainty`` provided,
-then this will be used. Otherwise, the ``variance`` parameter must be set.::
+an `~astropy.nddata.NDData` object with ``image.uncertainty`` provided,
+then this will be used. Otherwise, the ``variance`` parameter must be set.
+
+.. code-block:: python
 
   extract = specreduce.extract.HorneExtract(image-bg, trace, variance=var_array)
 
 An optional mask array for the image may be supplied to HorneExtract as well. 
 This follows the same convention and can either be attached to ``image`` if it
-is an ``astropy.NDData`` object, or supplied as a keyword argument.
+is an ``~astropy.nddata.NDData`` object, or supplied as a keyword argument.
 
 The extraction methods automatically detect non-finite pixels in the input
 image and combine them with the user-supplied mask to prevent them from biasing the
@@ -102,17 +112,24 @@ differently as documented in the API.
 
 The previous examples in this section show how to initialize the BoxcarExtract
 or HorneExtract objects with their required parameters. To extract the 1D
-spectrum::
+spectrum
+
+.. code-block:: python
 
   spectrum = extract.spectrum
 
 The ``extract`` object contains all the set options.  The extracted 1D spectrum
 can be accessed via the ``spectrum`` property or by calling (e.g ``extract()``)
-the ``extract`` object (which also allows temporarily overriding any values)::
+the ``extract`` object (which also allows temporarily overriding any values)
+
+.. code-block:: python
 
   spectrum2 = extract(width=6)
 
-or, for example to override the original ``trace_object``::
+or, for example to override the original ``trace_object``
+
+.. code-block:: python
+
   spectrum2 = extract(trace_object=new_trace)
 
 Spatial profile options
@@ -138,16 +155,22 @@ and to use the defaults for bins and interpolation degree, or to override these
 defaults a dictionary can be passed in.
 
 For example, to use the ``interpolated_profile`` option with default bins and
-interpolation degree::
+interpolation degree
+
+.. code-block:: python
 
   interp_profile_extraction = extract(spatial_profile='interpolated_profile')
 
-Or, to override the default of 10 samples and use 20 samples::
+Or, to override the default of 10 samples and use 20 samples
+
+.. code-block:: python
 
   interp_profile_extraction = extract(spatial_profile={'name': 'interpolated_profile',
                                     'n_bins_interpolated_profile': 20)
 
-Or, to do a cubic interpolation instead of the default linear::
+Or, to do a cubic interpolation instead of the default linear
+
+.. code-block:: python
 
     interp_profile_extraction = extract(spatial_profile={'name': 'interpolated_profile',
                                     'interp_degree_interpolated_profile': 3)
@@ -163,7 +186,9 @@ wavelength units will be pixels. Wavelength and flux calibration steps are not
 included here.
 
 Putting all these steps together, a simple extraction process might look
-something like::
+something like
+
+.. code-block:: python
 
     from specreduce.tracing import FlatTrace
     from specreduce.background import Background
diff --git a/docs/index.rst b/docs/index.rst
index 40727b1a..5c96bdf1 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -1,20 +1,28 @@
+.. _docroot:
+
 ########################
 Specreduce Documentation
 ########################
 
-The `specreduce <https://specreduce.readthedocs.io/en/latest/index.html>`_ package
-aims to provide a data reduction toolkit for optical
+The :ref:`specreduce <docroot>` package aims to provide a data reduction toolkit for optical
 and infrared spectroscopy, on which applications such as pipeline processes for
-specific instruments can be built. The scope of its functionality is limited to
-basic spectroscopic reduction, currently encompassing the following three tasks:
+specific instruments can be built. The scope of its functionality is limited to basic spectroscopic
+reduction, currently encompassing the following tasks:
 
-#. Determining the trace of a spectrum dispersed in a 2D image, either by setting a flat
+#. **Spectrum trace determination:** Determining the trace of a spectrum dispersed in a 2D
+   image, either by setting a flat
    trace, providing a custom trace array, or fitting a spline, polynomial, or other model
    to the positions of the dispersed spectrum.
-#. Generating a background based on a region on one or both sides of this trace, and making
+#. **2D tilt correction:** Calculating a 2D tilt model from arc spectra and applying it to
+   correct the wavelength curvature along the cross-dispersion axis.
+#. **1D wavelength calibration:** Calculating a 1D wavelength solution from arc spectra and line
+   lists.
+#. **Background subtraction:** Generating a background based on a region on one or both sides of
+   this trace, and making
    available the background image, 1D spectrum of the background, and the
    background-subtracted image.
-#. Performing either a Horne (a.k.a. "optimal") or boxcar extraction on either the original
+#. **Spectrum extraction:** Performing either a Horne (a.k.a. "optimal") or boxcar extraction on
+   either the original
    or background-subtracted 2D spectrum, using the trace generated by the first task to
    generate a 1D spectrum.
 
@@ -26,39 +34,23 @@ and visualization by `matplotlib <https://matplotlib.org/>`_. A few
 examples will be provided that feature ``specreduce`` in conjunction with these
 complementary packages.
 
-.. note::
-
-    Specreduce is currently an incomplete work-in-progress and is liable to
-    change. Please feel free to contribute code and suggestions through github.
 
+.. toctree::
+    :maxdepth: 2
 
-.. _spectral-extraction:
-
-*******************
-Spectral Extraction
-*******************
+    quickstart
 
 .. toctree::
     :maxdepth: 2
 
-    extraction_quickstart.rst
+    user_guide.rst
 
-***********
-Calibration
-***********
 
 .. toctree::
-    :maxdepth: 1
-
-    wavelength_calibration.rst
-    extinction.rst
-    specphot_standards.rst
-    mask_treatment/mask_treatment.rst
+    :maxdepth: 2
 
+    contributing.rst
 
-*****
-Index
-*****
 
 .. toctree::
     :maxdepth: 1
@@ -68,15 +60,3 @@ Index
 * :ref:`genindex`
 * :ref:`modindex`
 * :ref:`search`
-
-
-****************
-Development Docs
-****************
-
-.. toctree::
-    :maxdepth: 1
-
-    process/index
-    terms
-
diff --git a/docs/installation.rst b/docs/installation.rst
new file mode 100644
index 00000000..dd531292
--- /dev/null
+++ b/docs/installation.rst
@@ -0,0 +1,24 @@
+.. _installation:
+
+Installation
+============
+
+Specreduce can be installed from PyPI using ``pip``
+
+.. code-block:: console
+
+    pip install specreduce
+
+from the conda-forge repository using ``conda``
+
+.. code-block:: console
+
+    conda install conda-forge::specreduce
+
+or by cloning the repository from `GitHub <https://github.com/astropy/specreduce>`_
+
+.. code-block:: console
+
+    git clone https://github.com/astropy/specreduce.git
+    cd specreduce
+    pip install .
\ No newline at end of file
diff --git a/docs/quickstart.rst b/docs/quickstart.rst
new file mode 100644
index 00000000..7f06cb2a
--- /dev/null
+++ b/docs/quickstart.rst
@@ -0,0 +1,10 @@
+.. _spectral-extraction:
+
+Quickstart
+==========
+
+.. toctree::
+    :maxdepth: 2
+
+    installation.rst
+    extraction_quickstart.rst
\ No newline at end of file
diff --git a/docs/terms.rst b/docs/terms.rst
index 815233e0..304bdfac 100644
--- a/docs/terms.rst
+++ b/docs/terms.rst
@@ -476,10 +476,9 @@ Processing Steps
       what redshift?)
 
 
-----
 
 Mentioned but not defined
--------------------------
+=========================
 
 - WCS & Database archive
 - Cloud archiving
diff --git a/docs/tilt_correction/tilt_correction.rst b/docs/tilt_correction/tilt_correction.rst
new file mode 100644
index 00000000..dcea1ff7
--- /dev/null
+++ b/docs/tilt_correction/tilt_correction.rst
@@ -0,0 +1,2 @@
+2D Tilt Correction
+==================
\ No newline at end of file
diff --git a/docs/trace.rst b/docs/trace.rst
new file mode 100644
index 00000000..7bb86063
--- /dev/null
+++ b/docs/trace.rst
@@ -0,0 +1,166 @@
+Tracing
+=======
+
+The `specreduce.tracing` module provides three ``Trace`` classes that are used to define the
+spatial position (trace) of a spectrum across a 2D detector image: `~specreduce.tracing.FlatTrace`,
+`~specreduce.tracing.FitTrace`, and `~specreduce.tracing.ArrayTrace`. Each trace class requires
+the 2D spectral image as input, along with trace-specific parameters that control how the trace
+is determined.
+
+Flat trace
+----------
+
+`~specreduce.tracing.FlatTrace` assumes that the spectrum follows a straight line across the
+detector, and is best for well-aligned spectrographs with minimal optical distortion. To
+initialize a `~specreduce.tracing.FlatTrace`, we need to specify the fixed cross-dispersion
+pixel position:
+
+.. code-block:: python
+
+    trace = specreduce.tracing.FlatTrace(image, trace_pos=6)
+
+
+.. plot::
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from scipy.stats import norm
+
+    fw = 10
+    nd, ncd = 31, 13
+    xd, xcd = np.arange(nd), np.arange(ncd)
+
+    spectrum = np.zeros((ncd, nd))
+    spectrum[:,:] = norm(6.0, 1.5).pdf(xcd)[:, None]
+
+    plt.rc('font', size=13)
+    fig, ax = plt.subplots(figsize=(fw, fw*(ncd/nd)), constrained_layout=True)
+    ax.imshow(spectrum, origin='lower')
+    ax.plot((0, nd-1), (6, 6), c='k')
+    ax.set_xticks(xd+0.5, minor=True)
+    ax.set_yticks(xcd+0.5, minor=True)
+    ax.grid(alpha=0.25, lw=1, which='minor')
+    plt.setp(ax, xlabel='Dispersion axis', ylabel='Cross-dispersion axis')
+    fig.show()
+
+FitTrace
+--------
+
+`~specreduce.tracing.FitTrace` fits a polynomial function to automatically detected spectrum
+positions, and is suitable for typical spectra with smooth, continuous trace profiles. The trace
+model can be chosen from `~astropy.modeling.polynomial.Chebyshev1D`,
+`~astropy.modeling.polynomial.Legendre1D`, `~astropy.modeling.polynomial.Polynomial1D`,
+or `~astropy.modeling.spline.Spline1D`, and the fitting can be optimized by binning the spectrum
+along the dispersion axis.
+
+.. code-block:: python
+
+    trace = specreduce.tracing.FitTrace(image, bins=10, trace_model=Polynomial1D(3))
+
+.. plot::
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from scipy.stats import norm
+    plt.rc('font', size=13)
+
+    fw = 10
+    nd, ncd = 31, 13
+    xd, xcd = np.arange(nd), np.arange(ncd)
+
+    tr = np.poly1d([-0.01, 0.2, 7.0])
+    spectrum = np.zeros((ncd, nd))
+    for i,x in enumerate(xd):
+        spectrum[:,i] = norm(tr(x), 1.0).pdf(xcd)
+
+    fig, ax = plt.subplots(figsize=(fw, fw*(ncd/nd)), constrained_layout=True)
+    ax.imshow(spectrum, origin='lower')
+    ax.plot(xd, tr(xd), 'k')
+    ax.set_xticks(xd+0.5, minor=True)
+    ax.set_yticks(xcd+0.5, minor=True)
+    ax.grid(alpha=0.25, lw=1, which='minor')
+    plt.setp(ax, xlabel='Dispersion axis', ylabel='Cross-dispersion axis')
+    fig.show()
+
+The method works by (optionally) binning the 2D spectrum along the dispersion axis, finding
+the PSF peak position along the cross-dispersion for each bin, and then fitting a 1D polynomial to
+the cross-dispersion and dispersion axis positions. Binning is optional, and the native image
+resolution is used by default. Binning can significantly increase the reliability of the fitting
+with low SNR spectra, and always increases the speed.
+
+The peak detection method can be chosen from ``max``, ``centroid``, and ``gaussian``. Of these
+methods, ``max`` is the fastest but yields an integer pixel precision.  Both ``centroid`` and
+``gaussian`` can be used when sub-pixel precision is required, and ``gaussian``, while being the
+slowest method of the three, is the best option if the data is significantly contaminated by
+non-finite values.
+
+.. plot::
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from scipy.stats import norm
+    from numpy.random import seed, normal
+    from astropy.modeling.models import Gaussian1D
+    from astropy.modeling.fitting import DogBoxLSQFitter
+    plt.rc('font', size=13)
+    seed(5)
+
+    fw = 10
+    nd, ncd = 31, 13
+    xd, xcd = np.arange(nd), np.arange(ncd)
+
+    psf = norm(5.4, 1.5).pdf(xcd) + normal(0, 0.01, ncd)
+    fitter = DogBoxLSQFitter()
+    m = fitter(Gaussian1D(), xcd, psf)
+
+    fig, ax = plt.subplots(figsize=(fw, fw*(ncd/nd)), constrained_layout=True)
+    ax.step(xcd, psf, where='mid', c='k')
+    ax.axvline(xcd[np.argmax(psf)], label='max')
+    ax.axvline(np.average(xcd, weights=psf), ls='--', label='centroid')
+    ax.axvline(m.mean.value, ls=':', label='gaussian')
+    ax.plot(xcd, m(xcd), ls=':')
+    ax.legend()
+    plt.setp(ax, yticks=[], ylabel='Flux', xlabel='Cross-dispersion axis [pix]', xlim=(0, ncd-1))
+    fig.show()
+
+ArrayTrace
+----------
+
+`~specreduce.tracing.ArrayTrace` uses a pre-defined array of positions for maximum flexibility,
+and is ideal for complex or unusual trace shapes that are difficult to model mathematically.
+`~specreduce.tracing.ArrayTrace` initialization requires an array of cross-dispersion pixel
+positions. The size of the array must match the number of dispersion-axis pixels in the image.
+
+.. code-block:: python
+
+    trace = specreduce.tracing.ArrayTrace(image, positions)
+
+.. plot::
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from scipy.stats import norm
+    plt.rc('font', size=13)
+
+    fw = 10
+    nd, ncd = 31, 13
+    xd, xcd = np.arange(nd), np.arange(ncd)
+
+    tr = np.full_like(xd, 6)
+    tr[:6] = 4
+    tr[15:23] = 8
+
+    spectrum = np.zeros((ncd, nd))
+
+    for i,x in enumerate(xd):
+        spectrum[:,i] = norm(tr[i], 1.0).pdf(xcd)
+
+    plt.rc('font', size=13)
+    fig, ax = plt.subplots(figsize=(fw, fw*(ncd/nd)), constrained_layout=True)
+    ax.imshow(spectrum, origin='lower')
+    ax.plot(xd, tr, 'k')
+    ax.set_xticks(xd+0.5, minor=True)
+    ax.set_yticks(xcd+0.5, minor=True)
+    ax.grid(alpha=0.25, lw=1, which='minor')
+    plt.setp(ax, xlabel='Dispersion axis', ylabel='Cross-dispersion axis')
+    fig.show()
diff --git a/docs/user_guide.rst b/docs/user_guide.rst
new file mode 100644
index 00000000..cea5dcf6
--- /dev/null
+++ b/docs/user_guide.rst
@@ -0,0 +1,21 @@
+User Guide
+==========
+
+.. toctree::
+    :maxdepth: 1
+    :caption: Core Functionality
+
+    tilt_correction/tilt_correction
+    wavelength_calibration
+    trace
+    background
+    extraction
+
+
+.. toctree::
+    :maxdepth: 1
+    :caption: Additional Topics
+
+    extinction
+    specphot_standards
+    mask_treatment/mask_treatment
\ No newline at end of file
diff --git a/docs/wavelength_calibration.rst b/docs/wavelength_calibration.rst
index f34d7dec..45921352 100644
--- a/docs/wavelength_calibration.rst
+++ b/docs/wavelength_calibration.rst
@@ -1,7 +1,7 @@
 .. _wavelength_calibration:
 
-Wavelength Calibration
-======================
+1D Wavelength Calibration
+=========================
 
 Wavelength calibration is currently supported for 1D spectra. Given a list of spectral
 lines with known wavelengths and estimated pixel positions on an input calibration
@@ -18,14 +18,16 @@ The `~specreduce.wavelength_calibration.WavelengthCalibration1D` class can be us
 to fit a dispersion model to a list of line positions and wavelengths. Future development
 will implement catalogs of known lamp spectra for use in matching observed lines. In the
 example below, the line positions (``pixel_centers``) have already been extracted from
-``lamp_spectrum``::
+``lamp_spectrum``
+
+.. code-block:: python
 
     import astropy.units as u
-	 from specreduce import WavelengthCalibration1D
-	 pixel_centers = [10, 22, 31, 43]
-	 wavelengths = [5340, 5410, 5476, 5543]*u.AA
-	 test_cal = WavelengthCalibration1D(lamp_spectrum, line_pixels=pixel_centers,
-										line_wavelengths=wavelengths)
+    from specreduce import WavelengthCalibration1D
+    pixel_centers = [10, 22, 31, 43]
+    wavelengths = [5340, 5410, 5476, 5543]*u.AA
+    test_cal = WavelengthCalibration1D(lamp_spectrum, line_pixels=pixel_centers,
+        line_wavelengths=wavelengths)
     calibrated_spectrum = test_cal.apply_to_spectrum(science_spectrum)
 
 The example above uses the default model (`~astropy.modeling.functional_models.Linear1D`)
@@ -34,7 +36,9 @@ spectrum (``science_spectrum``). Any other 1D ``astropy`` model can be provided
 input ``model`` parameter to the `~specreduce.wavelength_calibration.WavelengthCalibration1D`.
 In the above example, the model fit and WCS construction is all done as part of the
 ``apply_to_spectrum()`` call, but you could also access the `~gwcs.wcs.WCS` object itself
-by calling::
+by calling
+
+.. code-block:: python
 
     test_cal.wcs
 
@@ -44,7 +48,9 @@ fit.
 
 You can also provide the input pixel locations and wavelengths of the lines as an
 `~astropy.table.QTable` with (at minimum) columns ``pixel_center`` and ``wavelength``,
-using the ``matched_line_list`` input argument::
+using the ``matched_line_list`` input argument
+
+.. code-block:: python
 
     from astropy.table import QTable
     pixels = [10, 20, 30, 40]*u.pix
diff --git a/pyproject.toml b/pyproject.toml
index 52565e49..20f18f5e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -23,7 +23,7 @@ test = [
     "tox",
 ]
 docs = [
-    "sphinx-astropy",
+    "sphinx-astropy[confv2]",
     "matplotlib>=3.7",
     "photutils>=1.0",
     "synphot",