diff --git a/docs/source/conf.py b/docs/source/conf.py
index 4add61c4c3..06f8a39a34 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -23,7 +23,7 @@ sys.path.insert(0, os.path.abspath('../sphinxext'))
 
 # Add any Sphinx extension module names here, as strings. They can be extensions
 # coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
-extensions = ['sphinx.ext.pngmath', 'sphinxcontrib.tikz']
+extensions = ['sphinx.ext.pngmath', 'sphinxcontrib.tikz', 'sphinx.ext.numfig']
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
@@ -196,11 +196,12 @@ latex_elements = {
 \setlistdepth{9}
 \usepackage{tikz}
 \usetikzlibrary{shapes,snakes,shadows,arrows,calc,decorations.markings,patterns,fit,matrix,spy}
+\usepackage{fixltx2e}
 \hypersetup{bookmarksdepth=3}
-%\renewcommand{\thechapter}{\Roman{chapter}}
 \setcounter{tocdepth}{2}
 \numberwithin{equation}{section}
-"""
+""",
+'printindex': r""
 }
 
 # The name of an image file (relative to this directory) to place at the top of
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 83eac22e0a..7566429403 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -18,9 +18,11 @@ free to send a message to the User's Group `mailing list`_.
 .. _Massachusetts Institute of Technology: http://web.mit.edu
 .. _mailing list: https://groups.google.com/forum/?fromgroups=#!forum/openmc-users
 
---------
-Contents
---------
+.. only:: html
+
+    --------
+    Contents
+    --------
 
 .. toctree::
     :maxdepth: 1
diff --git a/docs/source/methods/cmfd.rst b/docs/source/methods/cmfd.rst
index 098f0de204..49cb86c674 100644
--- a/docs/source/methods/cmfd.rst
+++ b/docs/source/methods/cmfd.rst
@@ -439,7 +439,7 @@ by the user. Once all diffusion parameters are calculated, CMFD matrices are
 formed where energy groups are the inner most iteration index. In OpenMC,
 compressed row storage sparse matrices are used due to the sparsity of CMFD
 operators. An example of this sparsity is shown for the 3-D BEAVRS model in
-figures :ref:`fig_loss` and :ref:`fig_prod` [BEAVRS]_. These matrices represent
+figures :num:`fig-loss` and :num:`fig-prod` [BEAVRS]_. These matrices represent
 an assembly radial mesh, 24 cell mesh in the axial direction and two energy
 groups. The loss matrix is 99.92% sparse and the production matrix is 99.99%
 sparse. Although the loss matrix looks like it is tridiagonal, it is really a
@@ -473,14 +473,14 @@ no fission neutrons appear with energies in the thermal group.
    | \ :math:`\left\langle\overline{J}^{u,g}_{l\pm 1/2,m,n}\Delta_m^v\Delta_n^w\right\rangle`   | current        | mesh, energy              |
    +--------------------------------------------------------------------------------------------+----------------+---------------------------+
 
-.. _fig_loss:
+.. _fig-loss:
 
 .. figure:: ../_images/loss.png
    :scale: 50
 
    Sparsity of Neutron Loss Operator
 
-.. _fig_prod:
+.. _fig-prod:
 
 .. figure:: ../_images/prod.png
    :scale: 50
@@ -525,9 +525,9 @@ Once weight adjustment factors were calculated, each neutron's statistical
 weight in the source bank was modified according to its location and energy.
 Examples of CMFD simulations using OpenMC can be found in [HermanThesis]_.
 
-----------
-References
-----------
+.. only:: html
+
+   .. rubric:: References
 
 .. [BEAVRS] Nick Horelik, Bryan Herman. *Benchmark for Evaluation And Verification of Reactor
             Simulations*. Massachusetts Institute of Technology, http://crpg.mit.edu/pub/beavrs
diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst
index db2e5156ee..5126252cd7 100644
--- a/docs/source/methods/cross_sections.rst
+++ b/docs/source/methods/cross_sections.rst
@@ -63,9 +63,9 @@ Other Methods
 A good survey of other energy grid techniques, including unionized energy grids,
 can be found in a paper by Leppanen_.
 
-----------
-References
-----------
+.. only:: html
+
+   .. rubric:: References
 
 .. [Brown] Forrest B. Brown, "New Hash-based Energy Lookup Algorithm for Monte
            Carlo codes," LA-UR-14-24530, Los Alamos National Laboratory (2014).
diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst
index 8dda10e013..9eece1e919 100644
--- a/docs/source/methods/geometry.rst
+++ b/docs/source/methods/geometry.rst
@@ -42,9 +42,11 @@ One can confirm that any point inside this sphere will correspond to
 In OpenMC, every surface defined by the user is assigned an integer to uniquely
 identify it. We can then refer to either of the two half-spaces created by a
 surface by a combination of the unique ID of the surface and a positive/negative
-sign. The following illustration shows an example of an ellipse with unique ID 1
+sign.  Figure :num:`fig-halfspace` shows an example of an ellipse with unique ID 1
 dividing space into two half-spaces.
 
+.. _fig-halfspace:
+
 .. figure:: ../_images/halfspace.*
    :align: center
    :figclass: align-center
@@ -57,9 +59,11 @@ to be defined by intersections, unions, and differences or half-spaces, OpenMC
 is currently limited to cells defined only as intersections of
 half-spaces. Thus, the specification of the cell must include a list of
 half-space references whose intersection defines the region. The region is then
-assigned a material defined elsewhere. The following illustration shows an
+assigned a material defined elsewhere. Figure :num:`fig-union` shows an
 example of a cell defined as the intersection of an ellipse and two planes.
-   
+
+.. _fig-union:
+
 .. figure:: ../_images/union.*
    :align: center
    :figclass: align-center
@@ -399,10 +403,12 @@ Rectilinear Lattice Indexing
 ----------------------------
 
 Indices are assigned to tiles in a rectilinear lattice based on the tile's
-position along the :math:`x`, :math:`y`, and :math:`z` axes.  The figure below
-maps the indices for a 2D lattice.  The indices, (1, 1), map to the
-lower-left tile.  (5, 1) and (5, 5) map to the lower-right and upper-right
-tiles, respectively.
+position along the :math:`x`, :math:`y`, and :math:`z` axes.  Figure
+:num:`fig-rect-lat` maps the indices for a 2D lattice.  The indices, (1, 1),
+map to the lower-left tile.  (5, 1) and (5, 5) map to the lower-right and
+upper-right tiles, respectively.
+
+.. _fig-rect-lat:
 
 .. figure:: ../_images/rect_lat.*
    :align: center
@@ -431,12 +437,15 @@ corner of the lattice, and :math:`p_0, p_1, p_2` are the pitches along the
 Hexagonal Lattice Indexing
 --------------------------
 
-A skewed coordinate system is used for indexing hexagonal lattice tiles.  Rather
-than a :math:`y`-axis, another axis is used that is rotated 30 degrees
+A skewed coordinate system is used for indexing hexagonal lattice tiles.
+Rather than a :math:`y`-axis, another axis is used that is rotated 30 degrees
 counter-clockwise from the :math:`y`-axis.  This axis is referred to as the
-:math:`\alpha`-axis.  The figure below shows how 2D hexagonal tiles are mapped
-with the :math:`(x, \alpha)` basis.  In this system, (0, 0) maps to the center
-tile, (0, 2) to the top tile, and (2, -1) to the middle tile on the right side.
+:math:`\alpha`-axis.  Figure :num:`fig-hex-lat` shows how 2D hexagonal tiles
+are mapped with the :math:`(x, \alpha)` basis.  In this system, (0, 0) maps to
+the center tile, (0, 2) to the top tile, and (2, -1) to the middle tile on the
+right side.
+
+.. _fig-hex-lat:
 
 .. figure:: ../_images/hex_lat.*
    :align: center
diff --git a/docs/source/methods/parallelization.rst b/docs/source/methods/parallelization.rst
index a6ec2b885e..63c4a9d212 100644
--- a/docs/source/methods/parallelization.rst
+++ b/docs/source/methods/parallelization.rst
@@ -599,9 +599,9 @@ is actually independent of the number of nodes:
     E \left [ \Lambda_{j_{\text{max}}} \right ] = \sqrt{ \frac{N\sigma^2}{2\pi
     k^2}}.
 
-----------
-References
-----------
+.. only:: html
+
+   .. rubric:: References
 
 .. [Troubetzkoy] E. Troubetzkoy, H. Steinberg, and M. Kalos, "Monte Carlo
    Radiation Penetration Calculations on a Parallel Computer,"
diff --git a/docs/source/methods/physics.rst b/docs/source/methods/physics.rst
index db6cfd89e9..ed10b32343 100644
--- a/docs/source/methods/physics.rst
+++ b/docs/source/methods/physics.rst
@@ -1552,9 +1552,9 @@ default, the cutoff weight in OpenMC is :math:`w_c = 0.25` and the survival
 weight is :math:`w_s = 1.0`. These parameters vary from one Monte Carlo code to
 another.
 
-----------
-References
-----------
+.. only:: html
+
+   .. rubric:: References
 
 .. [Doyas] Richard J. Doyas and Sterrett T. Perkins, "Interpolation of Tabular
    Secondary Neutron and Photon Energy Distributions," *Nucl. Sci. Eng.*,
diff --git a/docs/source/methods/random_numbers.rst b/docs/source/methods/random_numbers.rst
index 88c4d65d52..3ea61719c5 100644
--- a/docs/source/methods/random_numbers.rst
+++ b/docs/source/methods/random_numbers.rst
@@ -62,9 +62,9 @@ Note that :eq:`lcg-skipahead` has the same general form as \eqref{eq:lcg}, so
 the idea is to determine the new multiplicative and additive constants in
 :math:`O(\log_2 N)` operations.
 
-----------
-References
-----------
+.. only:: html
+
+   .. rubric:: References
 
 .. [LEcuyer] P. L’Ecuyer, "Tables of Linear Congruential Generators of
    Different Sizes and Good Lattice Structures," *Math. Comput.*, **68**, 249
diff --git a/docs/source/methods/tallies.rst b/docs/source/methods/tallies.rst
index 6c16f08136..fd0812245a 100644
--- a/docs/source/methods/tallies.rst
+++ b/docs/source/methods/tallies.rst
@@ -475,9 +475,9 @@ normal distribution, we use an `unpublished rational approximation`_. After
 using the rational approximation, one iteration of Newton's method is applied to
 improve the estimate of the percentile.
 
-----------
-References
-----------
+.. only:: html
+
+   .. rubric:: References
 
 .. [George] E. E. Olusegun George and Meenakshi Sivaram, "A modification of the
    Fisher-Cornish approximation for the student t percentiles," Communication