diff --git a/anvil/config.py b/anvil/config.py
index 7eeb68e..44e5d54 100644
--- a/anvil/config.py
+++ b/anvil/config.py
@@ -185,15 +185,36 @@ def parse_bootstrap_sample_size(self, n_boot: int):
         log.warning(f"Using user specified bootstrap sample size: {n_boot}")
         return n_boot
 
-    def produce_bootstrap_seed(
-        self, manual_bootstrap_seed: (int, type(None)) = None):
+    def produce_bootstrap_seed(self, manual_bootstrap_seed: (int, type(None)) = None):
         if manual_bootstrap_seed is None:
             return randint(0, maxsize)
         # numpy is actually this strict but let's keep it sensible.
-        if (manual_bootstrap_seed < 0) or (manual_bootstrap_seed > 2**32):
+        if (manual_bootstrap_seed < 0) or (manual_bootstrap_seed > 2 ** 32):
             raise ConfigError("Seed is outside of appropriate range: [0, 2 ** 32]")
         return manual_bootstrap_seed
 
+    def parse_cosh_fit_min_separation(self, n: int, training_geometry):
+        """The smallest lattice separation to include in when fitting a cosh function
+        to the correlator, so as to the extract the correlation length.
+
+        See also: ``produce_cosh_fit_window``.
+        """
+        if n > training_geometry.length // 2 - 2:
+            raise ConfigError("Not enough points to for a three-parameter fit.")
+        return n
+
+    def produce_cosh_fit_window(self, training_geometry, cosh_fit_min_separation=None):
+        """Window of values corresponding to lattice separations, within which to fit
+        a cosh function to the correlator, so as to the extract the correlation length.
+        """
+        if cosh_fit_min_separation is None:
+            return slice(1, None)  # include all but (0, 0) separation by default
+
+        return slice(
+            cosh_fit_min_separation,
+            training_geometry.length - cosh_fit_min_separation + 1,
+        )
+
     @element_of("windows")
     def parse_window(self, window: float):
         """A numerical factor featuring in the calculation of the optimal 'window'
@@ -215,4 +236,4 @@ def produce_use_multiprocessing(self):
         """Don't use mp on MacOS"""
         if platform.system() == "Darwin":
             return False
-        return True
\ No newline at end of file
+        return True
diff --git a/anvil/observables.py b/anvil/observables.py
index 02930ba..c9657cf 100644
--- a/anvil/observables.py
+++ b/anvil/observables.py
@@ -5,47 +5,97 @@
 """
 import numpy as np
 from scipy.signal import correlate
-from math import ceil, pi, sin
+import scipy.optimize as optim
 import logging
 
 from anvil.utils import bootstrap_sample, Multiprocessing
 
-import scipy.optimize as optim
 
 log = logging.getLogger(__name__)
 
+
 def cosh_shift(x, xi, A, c):
     return A * np.cosh(-x / xi) + c
 
 
-def fit_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
-    # TODO should I bootstrap this whole process...?
-
-    T = training_geometry.length
-    # TODO: would be good to specify this in runcard
-    t0 = T // 4
-    window = slice(t0, T - t0 + 1)
+def fit_zero_momentum_correlator(
+    zero_momentum_correlator, training_geometry, cosh_fit_window
+):
+    r"""Uses scipy.optimize.curve_fit to fit a cosh function (i.e. exponential decay
+    with periodicity) to each correlator in the bootrap ensemble.
+
+    The correlator decays as a pure exponential in the limit of large separations,
+    and the characteristic scale of this decay is the correlation length, whose
+    reciprocal is a.k.a the (effective) pole mass.
+
+    Parameters
+    ----------
+    zero_momentum_correlator
+        The two point correlation function at zero spatial momentum, i.e. the
+        correlation between 1-d 'slices'.
+    training_geometry
+        The anvil.geometry object defining the lattice.
+    cosh_fit_window: slice object
+        A slice object which selects the points (i.e. separations) to include in the
+        fit. In general the signal at short separations will be contaminated by
+        shorter modes and should not be included in the fit.
+
+    Returns
+    -------
+    xi: list
+        List of optimal correlation lengths for each member of the bootstrap ensemble
+        for whom the fitting process converged successfully.
+    A: list
+        Same as above, but for the amplitude of the cosh function.
+    c: list
+        Same as above, but for the global shift in the fit (which should correspond
+        to the absolute value of the magnetization, squared.
+
+    See also
+    --------
+    :py:func:`anvil.observables.cosh_shift` : the function being fit to the data.
+    """
+    t = np.arange(training_geometry.length) - training_geometry.length // 2
+
+    # fit for each correlation func in the bootstrap ensemble
+    xi, A, c = [], [], []
+    for correlator in zero_momentum_correlator.transpose():
+        try:
+            popt, pcov = optim.curve_fit(
+                cosh_shift,
+                xdata=t[cosh_fit_window],
+                ydata=correlator[cosh_fit_window],
+            )
+            xi.append(popt[0])
+            A.append(popt[1])
+            c.append(popt[2])
+        except RuntimeError:
+            pass
+
+    n_boot = zero_momentum_correlator.shape[-1]
+    n_fits = len(xi)
+    log.info(
+        f"Cosh fit succeeded for {n_fits}/{n_boot} members of the bootstrap ensemble."
+    )
+    return xi, A, c
 
-    t = np.arange(T)
-    y = zero_momentum_correlator.mean(axis=-1)
-    yerr = zero_momentum_correlator.std(axis=-1)
 
-    try:
-        popt, pcov = optim.curve_fit(
-            cosh_shift,
-            xdata=t[window] - T // 2,
-            ydata=y[window],
-            sigma=yerr[window],
-        )
-        return (popt, pcov, t0)
-    except RuntimeError:
-        log.warning("Failed to fit cosh to correlation function.")
-        return None
+def correlation_length_from_fit(fit_zero_momentum_correlator):
+    """Returns numpy array containing a value for the  correlation length for each member
+    of the bootstrap ensemble for whom :py:func:`fit_zero_momentum_correlator` successfully
+    converged.
+    """
+    xi, _, _ = fit_zero_momentum_correlator
+    return np.array(xi)
 
 
-def correlation_length_from_fit(fit_zero_momentum_correlator):
-    popt, pcov, _ = fit_zero_momentum_correlator
-    return popt[0], np.sqrt(pcov[0, 0])
+def abs_magnetization_sq_from_fit(fit_zero_momentum_correlator):
+    """Returns numpy array containing a value for the absolute magnetization squared
+    for each member of the bootstrap ensemble for whom :py:func:`fit_zero_momentum_correlator`
+    successfully converged.
+    """
+    _, _, c = fit_zero_momentum_correlator
+    return np.array(c)
 
 
 def autocorrelation(chain):
@@ -102,16 +152,16 @@ def magnetization(configs, bootstrap_sample_size, bootstrap_seed):
     )
 
 
-def abs_magnetization_squared(magnetization):
+def abs_magnetization_sq(magnetization):
     return np.abs(magnetization).mean(axis=-1) ** 2  # <|m|>^2
 
 
-def magnetic_susceptibility(magnetization, abs_magnetization_squared):
-    return (magnetization ** 2).mean(axis=-1) - abs_magnetization_squared
+def magnetic_susceptibility(magnetization, abs_magnetization_sq):
+    return (magnetization ** 2).mean(axis=-1) - abs_magnetization_sq
 
 
 def magnetization_series(configs):
-    return configs.sum(axis=1)
+    return configs.sum(axis=1).numpy()
 
 
 def magnetization_autocorr(magnetization_series):
@@ -144,7 +194,9 @@ def __two_point_correlator(
             axis=-1  # sample average
         )
 
-    return correlator.reshape((training_geometry.length, training_geometry.length, -1))
+    return correlator.reshape(
+        (training_geometry.length, training_geometry.length, -1)
+    ).numpy()
 
 
 def two_point_correlator(
@@ -173,8 +225,8 @@ def two_point_correlator(
     return correlator.reshape((training_geometry.length, training_geometry.length, -1))
 
 
-def two_point_connected_correlator(two_point_correlator, abs_magnetization_squared):
-    return two_point_correlator - abs_magnetization_squared.view(1, 1, -1)
+def two_point_connected_correlator(two_point_correlator, abs_magnetization_sq):
+    return two_point_correlator - abs_magnetization_sq.view(1, 1, -1)
 
 
 def zero_momentum_correlator(two_point_correlator):
@@ -210,16 +262,6 @@ def ising_energy(two_point_correlator):
     return (two_point_correlator[1, 0] + two_point_correlator[0, 1]) / 2
 
 
-def inverse_pole_mass(effective_pole_mass, training_geometry):
-    T = training_geometry.length
-    t0 = T // 4
-    window = slice(t0, T - t0 + 1)
-
-    xi = np.reciprocal(effective_pole_mass)[window]
-
-    return np.nanmean(xi, axis=0)  # average over "large" t points
-
-
 def second_moment_correlation_length(two_point_correlator, susceptibility):
     """Second moment correlation length, defined as the normalised second
     moment of the two point correlator."""
@@ -239,11 +281,11 @@ def second_moment_correlation_length(two_point_correlator, susceptibility):
 def low_momentum_correlation_length(two_point_correlator, susceptibility):
     """A low-momentum estimate for the correlation length."""
     L = two_point_correlator.shape[0]
-    kernel = np.cos(2 * pi / L * np.arange(L)).reshape(L, 1, 1)
+    kernel = np.cos(2 * np.pi / L * np.arange(L)).reshape(L, 1, 1)
 
     g_tilde_00 = susceptibility
     g_tilde_10 = (kernel * two_point_correlator).sum(axis=(0, 1))
 
-    xi_sq = (g_tilde_00 / g_tilde_10 - 1) / (4 * sin(pi / L) ** 2)
+    xi_sq = (g_tilde_00 / g_tilde_10 - 1) / (4 * np.sin(np.pi / L) ** 2)
 
     return np.sqrt(xi_sq)
diff --git a/anvil/plot.py b/anvil/plot.py
index bea99ab..f6d1294 100644
--- a/anvil/plot.py
+++ b/anvil/plot.py
@@ -9,7 +9,7 @@
 import torch
 import numpy as np
 import matplotlib.pyplot as plt
-from mpl_toolkits.axes_grid1 import make_axes_locatable
+import matplotlib as mpl
 from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
 from matplotlib.font_manager import FontProperties
 
@@ -19,126 +19,80 @@
 from anvil.observables import cosh_shift
 
 
-def field_component(i, x_base, phi_model, base_neg, model_neg):
-    fig, ax = plt.subplots()
-
-    ax.hist(x_base, bins=50, density=True, histtype="step", label="base")
-    ax.hist(phi_model, bins=50, density=True, histtype="step", label="model, full")
-    ax.hist(
-        base_neg, bins=50, density=True, histtype="step", label="model, $M_{base} < 0$"
-    )
-    ax.hist(
-        model_neg, bins=50, density=True, histtype="step", label="model, $M_{mod} < 0$"
-    )
-    ax.set_title(f"Coordinate {i}")
-    ax.legend()
-    fig.tight_layout()
-    return fig
-
-
-def field_components(loaded_model, base_dist, lattice_size):
-    """Plot the distributions of base coordinates 'x' and output coordinates 'phi' and,
-    if known, plot the pdf of the target distribution."""
-    sample_size = 10000
-
-    # Generate a large sample from the base distribution and pass it through the trained model
-    with torch.no_grad():
-        x_base, _ = base_dist(sample_size)
-        sign = x_base.sum(dim=1).sign()
-        neg = (sign < 0).nonzero().squeeze()
-        phi_model, model_log_density = loaded_model(x_base, 0, neg)
-
-    base_neg = phi_model[neg]
-
-    sign = phi_model.sum(dim=1).sign()
-    neg = (sign < 0).nonzero().squeeze()
-    model_neg = phi_model[neg]
-
-    # Convert to shape (n_coords, sample_size * lattice_size)
-    # NOTE: this is all pointless for the 1-component scalar
-    x_base = x_base.reshape(sample_size * lattice_size, -1).transpose(0, 1)
-    phi_model = phi_model.reshape(sample_size * lattice_size, -1).transpose(0, 1)
-
-    base_neg = base_neg.reshape(1, -1)
-    model_neg = model_neg.reshape(1, -1)
-
-    for i in range(x_base.shape[0]):
-        yield field_component(i, x_base[i], phi_model[i], base_neg[i], model_neg[i])
-
-
-_plot_field_components = collect("field_components", ("training_context",))
-
-
-def example_configs(loaded_model, base_dist, training_geometry):
-    sample_size = 10
-
-    # Generate a large sample from the base distribution and pass it through the trained model
-    with torch.no_grad():
-        x_base, _ = base_dist(sample_size)
-        sign = x_base.sum(dim=1).sign()
-        neg = (sign < 0).nonzero().squeeze()
-        phi_model, model_log_density = loaded_model(x_base, 0, neg)
-
-    L = int(np.sqrt(phi_model.shape[1]))
-
-    phi_true = np.zeros((4, L, L))
-    phi_true[:, training_geometry.checkerboard] = phi_model[:4, : L ** 2 // 2]
-    phi_true[:, ~training_geometry.checkerboard] = phi_model[:4, L ** 2 // 2 :]
-
-    fig, axes = plt.subplots(2, 2, sharex=True, sharey=True)
-    for i, ax in enumerate(axes.flatten()):
-        conf = ax.imshow(phi_true[i])
-        fig.colorbar(conf, ax=ax)
-
-    fig.suptitle("Example configurations")
-
-    return fig
-
-
-_plot_example_configs = collect("example_configs", ("training_context",))
-
-
-@figure
-def plot_example_configs(_plot_example_configs):
-    return _plot_example_configs[0]
-
-
-@figuregen
-def plot_field_components(_plot_field_components):
-    yield from _plot_field_components[0]
-
-
 @figure
 def plot_zero_momentum_correlator(
     zero_momentum_correlator,
     training_geometry,
     fit_zero_momentum_correlator,
+    cosh_fit_window: type(slice),
+    plot_cosh_fit: bool = True,
 ):
-    """Plot zero_momentum_2pf as a function of t. Points are means across bootstrap
-    sample and errorbars are standard deviations across boostrap samples
+    r"""Plots the correlation function for pairs of one-dimensional 'slices', otherwise
+    referred to as the two point correlator at zero spatial momentum, as a function of
+    time.
+
+    Points and errorbars are means and standard deviations across a boostrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
+
+    Optionally plots a :math:`1\sigma` confidence interval for a pure-exponential (cosh)
+    fit performed for each member of the bootstrap sample in
+    :py:func:`fit_zero_momentum_correlator`.
+
+    Parameters
+    ---------
+    zero_momentum_correlator
+        Array containing bootstrapped correlation function
+    training_geometry
+        Geometry object defining the lattice
+    fit_zero_momentum_correlator
+        The parameters resulting from a least-squares fit of a cosh function to the
+        correlator
+    cosh_fit_window
+        Slice object which indexes the lattice separations that were used to perform
+        the cosh fit to the correlation function
+    plot_cosh_fit
+        If False, only plot the correlation function, and not the result of the fit
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See Also
+    --------
+    :py:func:`anvil.table.table_zero_momentum_correlator`
     """
-    L = training_geometry.length
-
     fig, ax = plt.subplots()
 
     ax.errorbar(
-        x=np.arange(L),
+        x=np.arange(training_geometry.length),
         y=zero_momentum_correlator.mean(axis=-1),
         yerr=zero_momentum_correlator.std(axis=-1),
         linestyle="",
+        zorder=2,
         label="sample statistics",
     )
 
-    if fit_zero_momentum_correlator is not None:
-        popt, pcov, x0 = fit_zero_momentum_correlator
-
-        x_2 = np.linspace(x0, L - x0, 100)
-        ax.plot(
-            x_2,
-            cosh_shift(x_2 - L // 2, *popt),
-            marker="",
-            label=r"fit: $A \cosh(-(x_2 - L/2) / \xi) + c$",
+    if plot_cosh_fit:
+        t = np.arange(training_geometry.length)[cosh_fit_window]
+        fit = []
+        for xi, A, c in zip(*fit_zero_momentum_correlator):
+            fit.append(
+                cosh_shift(t - training_geometry.length // 2, xi, A, c),
+            )
+        fit = np.array(fit).T  # (n_points, n_boot)
+
+        ax.fill_between(
+            t,
+            fit.mean(axis=1) - fit.std(axis=1),
+            fit.mean(axis=1) + fit.std(axis=1),
+            color="orange",
+            alpha=0.3,
+            zorder=1,
+            label=r"fit: $A \cosh(-(x_2 - L/2) / \xi) + c$"
+            + "\n"
+            + r"($1\sigma$ confidence)",
         )
+
     ax.set_yscale("log")
     ax.set_ylabel(r"$\sum_{x_1=0}^{L-1} G(x_1, x_2)$")
     ax.set_xlabel("$x_2$")
@@ -149,9 +103,26 @@ def plot_zero_momentum_correlator(
 
 @figure
 def plot_effective_pole_mass(training_geometry, effective_pole_mass):
-    """Plot effective pole mass as a function of x_2. The points are means
-    across bootstrap samples and the errorbars are standard deviation across
-    bootstrap.
+    r"""Plots the (effective) pole mass as a function of 'time' separation.
+    
+    Points and errorbars are means and standard deviations across a boostrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
+
+    Parameters
+    ----------
+    training_geometry
+        Geometry object defining the lattice.
+    effective_pole_mass
+        Array containing bootstrap ensemble of effective pole mass, for each
+        separation :math:`t = 1, \ldots, T - 1`
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See Also
+    --------
+    :py:func:`anvil.table.table_effective_pole_mass`
     """
     fig, ax = plt.subplots()
     ax.errorbar(
@@ -159,96 +130,305 @@ def plot_effective_pole_mass(training_geometry, effective_pole_mass):
         y=effective_pole_mass.mean(axis=-1),
         yerr=effective_pole_mass.std(axis=-1),
     )
-    ax.set_ylabel("$m_p^\mathrm{eff}$")
+    ax.set_ylabel(r"$m_p^\mathrm{eff}$")
     ax.set_xlabel("$x_2$")
     ax.set_title("Effective pole mass")
     return fig
 
 
 @figure
-def plot_two_point_correlator(two_point_correlator):
-    """Represent the two point function and it's error in x and t as heatmaps
-    of the respective matrices. Returns a figure with two plots, the left plot
-    is the mean two point function across bootstrap and the right plot is the
-    standard deviation divide by the mean (fractional error)
+def plot_correlation_length(
+    effective_pole_mass,
+    low_momentum_correlation_length,
+    correlation_length_from_fit,
+):
+    r"""Plots three estimates of correlation length.
+    
+    These are:
+        1. Estimate from fitting a cosh function to the correlation between
+            1-dimensional slices, using py:func:`correlation_length_from_fit`
+        2. Reciprocal of the effective pole mass estimator, using
+            :py:func:`effective_pole_mass` (evaluated at each separation, :math:`x_2`. 
+        3. Low momentum estimate, using :py:func:`low_momentum_correlation_length`
+
+    Points and errorbars are means and standard deviations across a boostrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
+    
+    Parameters
+    ----------
+    effective_pole_mass
+        Array containing estimate of the effective pole mass, for each separation
+        and each member of the bootstrap ensemble
+    low_momentum_correlation_length
+        Array containing a low-momentum estimate of the correlation length for
+        each member of the bootstrap ensemble.
+    correlation_length_from_fit
+        Array containing an estimate of the correlation length from a cosh fit
+        to the correlation function, for each member of the bootstrap
+        ensemble.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See Also
+    --------
+    :py:func:`anvil.observables.fit_zero_momentum_correlator`
+    :py:func:`anvil.table.table_correlation_length`
+    """
+    xi_arcosh = np.reciprocal(effective_pole_mass)
 
+    fig, ax = plt.subplots()
+
+    arcosh_points = ax.errorbar(
+        x=range(1, xi_arcosh.shape[0] + 1),
+        y=xi_arcosh.mean(axis=-1),
+        yerr=xi_arcosh.std(axis=-1),
+        zorder=3,
+    )
+
+    xi_lm = low_momentum_correlation_length.mean()
+    e_lm = low_momentum_correlation_length.std()
+    lm_hline = ax.axhline(xi_lm, linestyle="-", marker="", color="grey", zorder=2)
+    lm_fill = ax.fill_between(
+        ax.get_xlim(),
+        xi_lm + e_lm,
+        xi_lm - e_lm,
+        color="grey",
+        alpha=0.3,
+        zorder=1,
+    )
+
+    xi_fit = correlation_length_from_fit.mean()
+    e_fit = correlation_length_from_fit.std()
+    fit_hline = ax.axhline(xi_fit, linestyle="-", marker="", color="orange", zorder=2)
+    fit_fill = ax.fill_between(
+        ax.get_xlim(),
+        xi_fit + e_fit,
+        xi_fit - e_fit,
+        color="orange",
+        alpha=0.3,
+        zorder=1,
+    )
+
+    ax.legend(
+        handles=[arcosh_points, (lm_fill, lm_hline), (fit_fill, fit_hline)],
+        labels=["Estimate using arcosh", "Low momentum estimate", "Estimate from fit"],
+    )
+    return fig
+
+
+@figure
+def plot_two_point_correlator(two_point_correlator):
+    r"""Represents the two point correlator as a heatmap.
+
+    The data shown is the mean of a bootstrap sample of correlation functions, and is
+    normalised so that :math:`G(0, 0) = 1`. The colour axis is scaled using a symmetric
+    log scale, with a linear region spanning :math:`[-0.01, 0.01]`.
+    
+    The bootstrap dimension is assumed to be the last (``-1``) dimension of input arrays.
+    
+    Parameters
+    ----------
+    two_point_correlator
+        Array containing two point correlation function for each two-dimensional
+        separation :math:`(x_1, x_2)`, for each member of a bootstrap ensemble.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See Also
+    --------
+    :py:func:`anvil.plot.plot_two_point_correlator_error`
+    :py:func:`anvil.table.table_two_point_correlator`
     """
     corr = two_point_correlator.mean(axis=-1)
-    error = two_point_correlator.std(axis=-1)
-    norm = corr[0, 0]
+    corr /= corr[0, 0]
 
     L = corr.shape[0]
     corr = np.roll(corr, (-L // 2 - 1, -L // 2 - 1), (0, 1))
-    error = np.roll(error, (-L // 2 - 1, -L // 2 - 1), (0, 1))
 
-    fig, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)
-    ax1.set_title("$G(x, t)$")
-    ax2.set_title(r"$| \sigma_G / G |$")
-    ax1.set_xlabel("$x_1$")
-    ax2.set_xlabel("$x_1$")
-    ax1.set_ylabel("$x_2$")
+    fig, ax = plt.subplots()
+    ax.set_title("$G(x)$")
+    ax.set_xlabel("$x_1$")
+    ax.set_ylabel("$x_2$")
 
     tick_positions = [0, L // 2 - 1, L - 1]
     tick_labels = [r"$-\frac{L + 1}{2}$", 0, r"$\frac{L}{2}$"]
-    ax1.set_xticks(tick_positions)
-    ax1.set_yticks(tick_positions)
-    ax1.set_xticklabels(tick_labels)
-    ax1.set_yticklabels(tick_labels)
-    ax2.tick_params(axis="y", width=0)
+    ax.set_xticks(tick_positions)
+    ax.set_yticks(tick_positions)
+    ax.set_xticklabels(tick_labels)
+    ax.set_yticklabels(tick_labels)
 
-    im1 = ax1.imshow(corr / norm)
-    im2 = ax2.imshow(np.abs(error / corr))
+    norm = mpl.colors.SymLogNorm(linthresh=0.01, base=10)
+    im = ax.imshow(corr, norm=norm)
+    fig.colorbar(im, ax=ax, pad=0.01)
 
-    div1 = make_axes_locatable(ax1)
-    cax1 = div1.append_axes("right", size="5%", pad=0.05)
-    fig.colorbar(im1, cax=cax1)
+    return fig
 
-    div2 = make_axes_locatable(ax2)
-    cax2 = div2.append_axes("right", size="7%", pad=0.05)
-    fig.colorbar(im2, cax=cax2)
 
-    fig.tight_layout()
+@figure
+def plot_two_point_correlator_error(two_point_correlator):
+    r"""Heatmap of the error in the two point correlator for each separation.
+
+    The error is computed as the standard deviation over the bootstrap sample. The
+    data shown is this error divided by the mean of the bootstrap sample, i.e. the
+    fractional error.
+    
+    The bootstrap dimension is assumed to be the last (``-1``) dimension of input arrays.
+    
+    Parameters
+    ----------
+    two_point_correlator
+        Array containing two point correlation function for each two-dimensional
+        separation :math:`(x_1, x_2)`, for each member of a bootstrap ensemble.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See Also
+    --------
+    :py:func:`anvil.plot.plot_two_point_correlator`
+    :py:func:`anvil.table.table_two_point_correlator`
+    """
+    corr = two_point_correlator.mean(axis=-1)
+    error = two_point_correlator.std(axis=-1)
+
+    L = corr.shape[0]
+    corr = np.roll(corr, (-L // 2 - 1, -L // 2 - 1), (0, 1))
+    error = np.roll(error, (-L // 2 - 1, -L // 2 - 1), (0, 1))
+
+    fig, ax = plt.subplots()
+    ax.set_title(r"$| \sigma_G(x) / G(x) |$")
+    ax.set_xlabel("$x_1$")
+    ax.set_ylabel("$x_2$")
+
+    tick_positions = [0, L // 2 - 1, L - 1]
+    tick_labels = [r"$-\frac{L + 1}{2}$", 0, r"$\frac{L}{2}$"]
+    ax.set_xticks(tick_positions)
+    ax.set_yticks(tick_positions)
+    ax.set_xticklabels(tick_labels)
+    ax.set_yticklabels(tick_labels)
+
+    im = ax.imshow(np.abs(error / corr))
+    fig.colorbar(im, ax=ax, pad=0.01)
 
     return fig
 
 
 @figure
 def plot_magnetization(magnetization_series):
+    r"""Plots a histogram of the magnetization of each configuration in the Markov
+    chain resulting from the Metropolis-Hastings sampling phase.
+
+    Parameters
+    ----------
+    magnetization_series
+        Array containing the magnetization for each configuration in the output
+        sample from the Metropolis-Hastings sampling phase.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See also
+    --------
+    :py:func:`anvil.plot.plot_magnetization_series`
+    :py:func:`anvil.table.table_magnetization`
+    """
     fig, ax = plt.subplots()
     ax.set_title("Magnetization")
     ax.set_ylabel("Frequency")
     ax.set_xlabel("$M(t)$")
-    ax.hist(magnetization_series.numpy(), histtype="stepfilled", edgecolor="black")
+    ax.hist(magnetization_series, histtype="stepfilled", edgecolor="black")
     return fig
 
 
 @figure
 def plot_magnetization_series(magnetization_series, sample_interval):
-    chain_indices = np.arange(magnetization_series.shape[-1]) * sample_interval
-    fig, ax = plt.subplots()
-    ax.set_title("Magnetization")
-    ax.set_ylabel("$M(t)$")
-    ax.set_xlabel("$t$")
-    ax.plot(
-        chain_indices,
-        magnetization_series,
-        linestyle="-",
-        marker="",
-    )
+    r"""Plots the magnetization of each configuration in the Markov chain over the
+    course of the Metropolis-Hastings sampling phase.
+    
+    Parameters
+    ----------
+    magnetization_series
+        Array containing the magnetization for each configuration in the output
+        sample from the Metropolis-Hastings sampling phase.
+    sample_interval
+        The number of Metropolis updates which were discarded between each
+        configuration appearing in the input series.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See also
+    --------
+    :py:func:`anvil.plot.plot_magnetization`.
+    :py:func:`anvil.table.table_magnetization`
+    """
+    n_rows = 5
+    if magnetization_series.size % n_rows != 0:
+        magnetization_series = np.pad(
+            magnetization_series,
+            (0, magnetization_series.size % n_rows),
+            "empty",
+        )
+    magnetization_series = magnetization_series.reshape(n_rows, -1)
+    t = (np.arange(magnetization_series.size) * sample_interval).reshape(n_rows, -1)
+
+    fig, axes = plt.subplots(n_rows, 1, sharey=True)
+    for ax, x, y in zip(axes, t, magnetization_series):
+        ax.plot(x, y, linestyle="-", marker="")
+        ax.margins(0, 0)
+
+    axes[0].set_title("Magnetization")
+    axes[n_rows // 2].set_ylabel("$M(t)$")
+    axes[-1].set_xlabel("$t$")
+    fig.tight_layout()
     return fig
 
 
 @figure
 def plot_magnetization_autocorr(
-    magnetization_autocorr, magnetization_optimal_window, sample_interval
+        magnetization_autocorr, magnetization_optimal_window: int, sample_interval: int
 ):
+    r"""Plots the autocorrelation function for the magnetization of the sequence of
+    configurations generated in the Metropolis-Hastings sampling phase.
+
+    The x-axis corresponds to a number of steps separating pairs of configurations
+    in the sequence.
+    
+    Parameters
+    ----------
+    magnetization_autocorr
+        Array containing the autocorrelation function of the magnetization for each
+        configuration in the output sample from the Metropolis-Hastings sampling phase.
+    magnetization_optimal_window
+        The size of the window in which the integrated autocorrelation time is to be
+        computed such that the total error is minimized.
+    sample_interval
+        The number of Metropolis updates which were discarded between each
+        configuration appearing in the input series.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See also
+    --------
+    :py:func:`anvil.observables.optimal_window`
+    :py:func:`anvil.plot.plot_magnetization_integrated_autocorr`.
+    """
     cut = max(10, 2 * magnetization_optimal_window)
     chain_indices = np.arange(cut) * sample_interval
 
     fig, ax = plt.subplots()
     ax.set_title("Autocorrelation of magnetization")
     ax.set_ylabel(r"$\Gamma_M(\delta t)$")
-    ax.set_xlabel("$\delta t$")
+    ax.set_xlabel(r"$\delta t$")
 
     ax.plot(chain_indices, magnetization_autocorr[:cut])
 
@@ -279,9 +459,37 @@ def plot_magnetization_autocorr(
 @figure
 def plot_magnetization_integrated_autocorr(
     magnetization_integrated_autocorr,
-    magnetization_optimal_window,
-    sample_interval,
+    magnetization_optimal_window: int,
+    sample_interval: int,
 ):
+    r"""Plots the integrated autocorrelation function for the magnetization of the
+    sequence of configurations generated in the Metropolis-Hastings sampling phase.
+    
+    The x axis represents the size of the 'window' in which the summation is performed,
+    i.e. the point at which the autocorrelation function is truncated.
+    
+    Parameters
+    ----------
+    magnetization_integrated_autocorr
+        Array containing the cumulative sum of the autocorrelation function of the
+        magnetization for each configuration in the output sample from the Metropolis-
+        Hastings sampling phase.
+    magnetization_optimal_window
+        The size of the window in which the integrated autocorrelation time is to be
+        computed such that the total error is minimized.
+    sample_interval
+        The number of Metropolis updates which were discarded between each
+        configuration appearing in the input series.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+
+    See also
+    --------
+    :py:func:`anvil.observables.optimal_window`
+    :py:func:`anvil.plot.plot_magnetization_autocorr`.
+    """
     cut = max(10, 2 * np.max(magnetization_optimal_window))
     chain_indices = np.arange(cut) * sample_interval
     tau = magnetization_integrated_autocorr[magnetization_optimal_window]
@@ -300,11 +508,11 @@ def plot_magnetization_integrated_autocorr(
         ax.transData,
         sample_interval,
         f"sample interval: {sample_interval}",
-        "upper left",
+        "center right",
         pad=0.6,
         frameon=False,
         sep=4,
-        label_top=True,
+        label_top=False,
         fontproperties=FontProperties(size="x-large"),
     )
     ax.add_artist(scalebar)
diff --git a/anvil/table.py b/anvil/table.py
index 3346aa4..2e873dc 100644
--- a/anvil/table.py
+++ b/anvil/table.py
@@ -15,41 +15,104 @@
 @table
 def table_autocorrelation(
     magnetization_integrated_autocorr,
-    magnetization_optimal_window,
-    tau_chain,
-    acceptance,
+    magnetization_optimal_window: int,
+    tau_chain: float,
+    acceptance: float,
 ):
+    r"""
+    Tabulate some information related to the statistical efficiency of the Metropolis-
+    Hastings sampling phase.
+
+    Parameters
+    ----------
+    magnetization_integrated_autocorr
+        Array containing the cumulative sum of the autocorrelation function of the
+        magnetization for each configuration in the sample output by the Metropolis-
+        Hastings sampling phase.
+    magnetization_optimal_window
+        Integer corresponding to a window size in which the autocorrelation function
+        should be summed, such that the resulting estimate of the integrated
+        autocorrelation has the smallest possible total error.
+    tau_chain
+        Estimate of the integrated autocorrelation using the accept-reject statistics
+        of the sampling phase.
+    acceptance
+        Fraction of proposals which were accepted in the sampling phase.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.observables.optimal_window`
+    :py:func:`anvil.sample.calc_tau_chain`
+    """
     tau_mag = magnetization_integrated_autocorr[magnetization_optimal_window]
 
     df = pd.DataFrame(
         [acceptance, tau_mag, tau_chain],
-        index=["acceptance", "tau_mag", "tau_chain"],
+        index=["acceptance", "tau_from_magnetization", "tau_from_chain"],
         columns=["value"],
     )
     return df
 
 
 @table
-def table_fit(fit_zero_momentum_correlator, training_geometry):
-    if fit_zero_momentum_correlator is not None:
-        popt, pcov, t0 = fit_zero_momentum_correlator
-
-        res = [
-            [popt[0], np.sqrt(pcov[0, 0])],
-            [popt[2], np.sqrt(pcov[2, 2])],
-        ]
-        df = pd.DataFrame(
-            res,
-            columns=["Mean", "Standard deviation"],
-            index=["xi_fit", "m_fit"],
-        )
-        return df
+def table_fit(correlation_length_from_fit, abs_magnetization_sq_from_fit):
+    r"""Tabulate the correlation length and magnetization estimates resulting from the
+    fitting of a cosh to the correlation function.
+    
+    Values and errors are means and standard deviations over a bootstrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
+
+    Parameters
+    ----------
+    correlation_length_from_fit
+    abs_magnetization_sq_from_fit
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.observables.fit_zero_momentum_correlator`
+    """
+    res = [
+        [correlation_length_from_fit.mean(), correlation_length_from_fit.std()],
+        [abs_magnetization_sq_from_fit.mean(), abs_magnetization_sq_from_fit.std()],
+    ]
+    df = pd.DataFrame(
+        res,
+        columns=["mean", "error"],
+        index=["xi_from_fit", "abs_magnetization_sq_from_fit"],
+    )
+    return df
 
 
 @table
 def table_two_point_scalars(ising_energy, susceptibility):
-    """Table of the ising observables, with mean and standard deviation taken
-    across boostrap samples
+    r"""Table of scalar observables derived from the two point correlation function.
+    
+    Values and errors are means and standard deviations over a bootstrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
+
+    Parameters
+    ----------
+    ising_energy
+        Nearest-neighbour iteraction energy.
+    susceptibility
+        Magnetic susceptibility defined by the sum of the correlation function.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.observables.ising_energy`
+    :py:func:`anvil.observables.susceptibility`
     """
     res = [
         [ising_energy.mean(), ising_energy.std()],
@@ -57,42 +120,103 @@ def table_two_point_scalars(ising_energy, susceptibility):
     ]
     df = pd.DataFrame(
         res,
-        columns=["Mean", "Standard deviation"],
-        index=["Ising energy", "susceptibility"],
+        columns=["mean", "error"],
+        index=["ising_energy", "susceptibility"],
     )
     return df
 
 
 @table
-def table_magnetization(abs_magnetization_squared, magnetic_susceptibility):
+def table_magnetization(abs_magnetization_sq, magnetic_susceptibility):
+    r"""Table containing quantities derived from the sample-averaged magnetization.
+
+    Values and errors are means and standard deviations over a bootstrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
+
+    Parameters
+    ----------
+    abs_magnetization_sq
+        Array containing the sample mean of the absolute magnetization, squared, for each
+        member of the bootstrap ensemble.
+    magnetic_susceptibility
+        Array containing the susceptibility for each member of the bootstrap ensemble.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.tables.table_two_point_scalars`
+    :py:func:`anvil.tables.table_fit`
+    """
+
     res = [
-        [abs_magnetization_squared.mean(), abs_magnetization_squared.std()],
+        [abs_magnetization_sq.mean(), abs_magnetization_sq.std()],
         [magnetic_susceptibility.mean(), magnetic_susceptibility.std()],
     ]
     df = pd.DataFrame(
         res,
-        columns=["Mean", "Standard deviation"],
-        index=["<|m|>^2", "<m^2> - <|m|>^2"],
+        columns=["mean", "error"],
+        index=["abs_magnetization_sq", "magnetic_susceptibility"],
     )
     return df
 
 
 @table
 def table_correlation_length(
-    inverse_pole_mass,
+    effective_pole_mass,
     second_moment_correlation_length,
     low_momentum_correlation_length,
     correlation_length_from_fit,
     training_geometry,
 ):
-    """Tabulate four estimators of correlation length, with values and errors
-    taken as the mean and standard deviation of the bootstrap sample.
+    r"""Table containing four estimates of correlation length.
+
+    Values and errors are means and standard deviations over a bootstrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
+
+    Also displays the number of correlation lengths that can fit on the lattice, i.e.
+    :math:`\xi / L` where :math:`\xi` is the correlation length and :math:`L` is the
+    linear extent of the lattice.
+
+    Parameters
+    ----------
+    effective_pole_mass
+        Array containing estimate of the effective pole mass, for each separation
+        and each member of the bootstrap ensemble
+    second_moment_correlation_length
+        Estimate of the correlation length based on the second moment of the
+        two point correlation function, for each member of the bootstrap ensemble.
+    low_momentum_correlation_length
+        Array containing a low-momentum estimate of the correlation length for each
+        member of the bootstrap ensemble.
+    correlation_length_from_fit
+        Array containing an estimate of the correlation length from a cosh fit to
+        the correlation function, for each member of the bootstrap ensemble.
+    training_geometry
+        Geometry object defining the lattice.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.observables.fit_zero_momentum_correlator`
+    :py:func:`anvil.plot.plot_correlation_length`
     """
-    # TODO could add column with inverse (i.e. effective pole mass)
+    # Take the mean of the arcosh estimator over "large" separations
+    x0 = training_geometry.length // 4
+    window = slice(x0, training_geometry.length - x0 + 1)
+    xi_arcosh = np.nanmean(
+        np.reciprocal(effective_pole_mass)[window],
+        axis=0,
+    )
 
     res = [
-        list(correlation_length_from_fit),
-        [inverse_pole_mass.mean(), inverse_pole_mass.std()],
+        [correlation_length_from_fit.mean(), correlation_length_from_fit.std()],
+        [xi_arcosh.mean(), xi_arcosh.std()],
         [
             second_moment_correlation_length.mean(),
             second_moment_correlation_length.std(),
@@ -102,22 +226,41 @@ def table_correlation_length(
 
     df = pd.DataFrame(
         res,
-        columns=["Mean", "Standard deviation"],
+        columns=["mean", "error"],
         index=[
-            "Estimate from fit",
-            "Estimate using arcosh",
-            "Second moment estimate",
-            "Low momentum estimate",
+            "xi_from_fit",
+            "xi_from_arcosh",
+            "xi_from_second_moment",
+            "xi_from_low_momentum",
         ],
     )
-    df["No. correlation lengths"] = training_geometry.length / df["Mean"]
+    df["n_correlation_lengths"] = training_geometry.length / df["mean"]
     return df
 
 
 @table
 def table_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
-    """Table of zero_momentum_correlator, with mean and standard deviation
-    from bootstrap
+    r"""Table containing values of the two point correlation function in time-momentum
+    representation at zero momentum, for each separation.
+
+    Values and errors are means and standard deviations over a bootstrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
+
+    Parameters
+    ----------
+    zero_momentum_correlator
+        Array containing the correlation function for each 1-d separation, for each
+        member of the bootstrap ensemble.
+    training_geometry
+        Geometry object defining the lattice.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.plot.plot_zero_momentum_correlator`
     """
     means = zero_momentum_correlator.mean(axis=-1)[:, np.newaxis]
     stds = zero_momentum_correlator.std(axis=-1)[:, np.newaxis]
@@ -126,7 +269,7 @@ def table_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
 
     df = pd.DataFrame(
         data,
-        columns=["Mean", "Standard deviation"],
+        columns=["mean", "error"],
         index=range(training_geometry.length),
     )
     return df
@@ -134,8 +277,26 @@ def table_zero_momentum_correlator(zero_momentum_correlator, training_geometry):
 
 @table
 def table_effective_pole_mass(effective_pole_mass, training_geometry):
-    """Table of effective_pole_mass, with mean and standard deviation
-    from bootstrap
+    r"""Table containing values of the effective pole mass for each separation.
+
+    Values and errors are means and standard deviations over a bootstrap ensemble,
+    which is assumed to be the last (``-1``) dimension of input arrays.
+
+    Parameters
+    ----------
+    effective_pole_mass
+        Array containing the effective pole mass for each separation, for each
+        member of the bootstrap ensemble.
+    training_geometry
+        Geometry object defining the lattice.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.plot.plot_effective_pole_mass`
     """
     means = effective_pole_mass.mean(axis=-1)[:, np.newaxis]
     stds = effective_pole_mass.std(axis=-1)[:, np.newaxis]
@@ -143,16 +304,36 @@ def table_effective_pole_mass(effective_pole_mass, training_geometry):
     data = np.concatenate((means, stds), axis=1)
     df = pd.DataFrame(
         data,
-        columns=["Mean", "Standard deviation"],
+        columns=["mean", "error"],
         index=range(1, training_geometry.length - 1),
     )
     return df
 
 
 @table
-def table_two_point_correlator(training_geometry, two_point_correlator):
-    """For each x and t, tabulate the mean and standard deviation of the two
-    point function, estimated from bootstrap sample
+def table_two_point_correlator(two_point_correlator, training_geometry):
+    r"""Table containing values of the two point correlation function for each
+    two-dimensional separation.
+
+    Values and errors are means and standard deviations over a bootstrap ensemble
+    which is assumed to be the last (``-1``) dimension of input arrays.
+
+    Parameters
+    ----------
+    two_point_correlator
+        Array containing the correlation function for each 2-d separation, for each
+        member of the bootstrap ensemble.
+    training_geometry
+        Geometry object defining the lattice.
+
+    Returns
+    -------
+    pandas.core.frame.DataFrame
+
+    See Also
+    --------
+    :py:func:`anvil.plot.plot_two_point_correlator`
+    :py:func:`anvil.plot.plot_two_point_correlator_error`
     """
     corr = []
     index = []
@@ -163,5 +344,5 @@ def table_two_point_correlator(training_geometry, two_point_correlator):
         for j in range(training_geometry.length):
             corr.append([float(means[i, j]), float(stds[i, j])])
             index.append((i, j))
-    df = pd.DataFrame(corr, columns=["Mean", "Standard deviation"], index=index)
+    df = pd.DataFrame(corr, columns=["mean", "error"], index=index)
     return df
diff --git a/examples/runcards/report.md b/examples/runcards/report.md
index 9d71141..06f2d26 100644
--- a/examples/runcards/report.md
+++ b/examples/runcards/report.md
@@ -5,15 +5,18 @@
 {@table_two_point_scalars@}
 ## Two Point Correlator
 {@plot_two_point_correlator@}
+{@plot_two_point_correlator_error@}
 ## Zero momentum two point correlator
 {@table_fit@}
 {@plot_zero_momentum_correlator@}
 ## Correlation length
 {@plot_effective_pole_mass@}
 {@table_correlation_length@}
+{@plot_correlation_length@}
 ## Magnetisation
 {@table_magnetization@}
 {@plot_magnetization@}
 {@plot_magnetization_series@}
+## Autocorrelation
 {@plot_magnetization_autocorr@}
 {@plot_magnetization_integrated_autocorr@}