diff --git a/routines/__init__.py b/routines/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..298e3378860c408b896f3b1ec70273338053fe26
--- /dev/null
+++ b/routines/__init__.py
@@ -0,0 +1,23 @@
+"""Routines for analysing post-processed data that don't follow the parameters-processes structure.
+
+Routines for analysing post-processed data that don't follow the
+parameters-processes structure.
+
+Currently, these consist of plotting routines.  There is one module for each
+plotting routine.  Each module consists of two components and is structured as
+follows:
+1. An internal class.
+    The class defines the parameters and defines additional class attributes to
+    help with plotting.  The class also has one method (`plot`) that takes a
+    `matplotlib.Axes` object as an argument.  This method draws the plot on the
+    `Axes` object.
+2. A plotting function.
+    The user accesses this function.  This function defines the default
+    parameters in its arguments.  Within the function, a 'plotter' object is
+    defined using the internal class and then the function draws the plot on a
+    specified `matplotlib.Axes` object.
+
+This structure follows that of plotting routines in `seaborn`
+(https://github.com/mwaskom/seaborn), a Python visualisation library that is
+based on `matplotlib`.
+"""
diff --git a/routines/boxplot.py b/routines/boxplot.py
new file mode 100644
index 0000000000000000000000000000000000000000..d189859ec1591128c9b5342b1dafa2c0206224f1
--- /dev/null
+++ b/routines/boxplot.py
@@ -0,0 +1,113 @@
+#!/usr/bin/env python3
+
+import matplotlib.pyplot as plt
+import matplotlib.ticker as ticker
+import seaborn as sns
+
+from plottingabc import BasePlotter
+
+
+class _BoxplotPlotter(BasePlotter):
+    """Draw boxplots over time"""
+
+    def __init__(
+        self,
+        trace_df,
+        trace_name,
+        unit_scaling,
+        box_color,
+        xtick_step,
+        xlabel,
+        plot_title,
+    ):
+        super().__init__(trace_name, unit_scaling, xlabel, plot_title)
+        # Define attributes from arguments
+        self.trace_df = trace_df
+        self.box_color = box_color
+        self.xtick_step = xtick_step
+
+        # Define some labels
+        self.ylabel = "Normalised " + self.trace_name + " fluorescence (AU)"
+
+        # Define horizontal axis ticks and labels
+        # hacky! -- redefine column names
+        trace_df.columns = trace_df.columns * self.unit_scaling
+        self.fmt = ticker.FuncFormatter(
+            lambda x, pos: "{0:g}".format(
+                x / (self.xtick_step / self.unit_scaling)
+            )
+        )
+
+    def plot(self, ax):
+        """Draw the heatmap on the provided Axes."""
+        super().plot(ax)
+        ax.xaxis.set_major_formatter(self.fmt)
+        sns.boxplot(
+            data=self.trace_df,
+            color=self.box_color,
+            linewidth=1,
+            ax=ax,
+        )
+        ax.xaxis.set_major_locator(
+            ticker.MultipleLocator(self.xtick_step / self.unit_scaling)
+        )
+
+
+def boxplot(
+    trace_df,
+    trace_name,
+    unit_scaling=1,
+    box_color="b",
+    xtick_step=60,
+    xlabel="Time (min)",
+    plot_title="",
+    ax=None,
+):
+    """Draw series of boxplots from an array of time series of traces
+
+    Draw series of boxplots from an array of time series of traces, showing the
+    distribution of values at each time point over time.
+
+    Parameters
+    ----------
+    trace_df : pandas.DataFrame
+        Time series of traces (rows = cells, columns = time points).
+    trace_name : string
+        Name of trace being plotted, e.g. 'flavin'.
+    unit_scaling : int or float
+        Unit scaling factor, e.g. 1/60 to convert minutes to hours.
+    box_color : string
+        matplolib colour string, specifies colour of boxes in boxplot
+    xtick_step : int or float
+        Interval length, in unit time, to draw x axis ticks.
+    xlabel : string
+        x axis label.
+    plot_title : string
+        Plot title.
+    ax : matplotlib Axes
+        Axes in which to draw the plot, otherwise use the currently active Axes.
+
+    Returns
+    -------
+    ax : matplotlib Axes
+        Axes object with the heatmap.
+
+    Examples
+    --------
+    FIXME: Add docs.
+
+    """
+
+    plotter = _BoxplotPlotter(
+        trace_df,
+        trace_name,
+        unit_scaling,
+        box_color,
+        xtick_step,
+        xlabel,
+        plot_title,
+    )
+    if ax is None:
+        ax = plt.gca()
+    plotter.plot(ax)
+    return ax
diff --git a/routines/heatmap.py b/routines/heatmap.py
new file mode 100644
index 0000000000000000000000000000000000000000..e0f7b629cbd97d8b94758a7b5f640703e65823d2
--- /dev/null
+++ b/routines/heatmap.py
@@ -0,0 +1,181 @@
+#!/usr/bin/env python3
+
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib import cm, ticker
+
+from postprocessor.core.processes.standardscaler import standardscaler
+from plottingabc import BasePlotter
+
+
+class _HeatmapPlotter(BasePlotter):
+    """Draw heatmap"""
+
+    def __init__(
+        self,
+        trace_df,
+        trace_name,
+        buddings_df,
+        cmap,
+        unit_scaling,
+        xtick_step,
+        scale,
+        robust,
+        xlabel,
+        ylabel,
+        cbarlabel,
+        plot_title,
+    ):
+        super().__init__(trace_name, unit_scaling, xlabel, plot_title)
+        # Define attributes from arguments
+        self.trace_df = trace_df
+        self.buddings_df = buddings_df
+        self.cmap = cmap
+        self.xtick_step = xtick_step
+        self.scale = scale
+        self.robust = robust
+
+        # Define some labels
+        self.cbarlabel = cbarlabel
+        self.ylabel = ylabel
+
+        # Scale
+        if self.scale:
+            self.trace_scaled = standardscaler.as_function(self.trace_df)
+        else:
+            self.trace_scaled = self.trace_df
+
+        # If robust, redefine colormap scale to remove outliers
+        if self.robust:
+            self.vmin = np.nanpercentile(self.trace_scaled, 2)
+            self.vmax = np.nanpercentile(self.trace_scaled, 98)
+            # Make axes even
+            if self.scale:
+                if np.abs(self.vmin) > np.abs(self.vmax):
+                    self.vmax = -self.vmin
+                else:
+                    self.vmin = -self.vmax
+        else:
+            self.vmin = None
+            self.vmax = None
+
+        # Define horizontal axis ticks and labels
+        # hacky! -- redefine column names
+        trace_df.columns = trace_df.columns * self.unit_scaling
+        self.fmt = ticker.FuncFormatter(
+            lambda x, pos: "{0:g}".format(x * self.unit_scaling)
+        )
+
+    def plot(self, ax, cax):
+        """Draw the heatmap on the provided Axes."""
+        super().plot(ax)
+        ax.xaxis.set_major_formatter(self.fmt)
+        # Draw trace heatmap
+        trace_heatmap = ax.imshow(
+            self.trace_scaled,
+            cmap=self.cmap,
+            interpolation="none",
+            vmin=self.vmin,
+            vmax=self.vmax,
+        )
+        # Horizontal axis labels as multiples of xtick_step, taking
+        # into account unit scaling
+        ax.xaxis.set_major_locator(
+            ticker.MultipleLocator(self.xtick_step / self.unit_scaling)
+        )
+        # Overlay buddings, if present
+        if self.buddings_df is not None:
+            # Must be masked array for transparency
+            buddings_array = self.buddings_df.to_numpy()
+            buddings_heatmap_mask = np.ma.masked_where(
+                buddings_array == 0, buddings_array
+            )
+            # Overlay
+            ax.imshow(
+                buddings_heatmap_mask,
+                interpolation="none",
+            )
+        # Draw colour bar
+        ax.figure.colorbar(mappable=trace_heatmap, cax=cax, ax=ax, label=self.cbarlabel)
+
+
+def heatmap(
+    trace_df,
+    trace_name,
+    buddings_df=None,
+    cmap=cm.RdBu,
+    unit_scaling=1,
+    xtick_step=60,
+    scale=True,
+    robust=True,
+    xlabel="Time (min)",
+    ylabel="Cell",
+    cbarlabel="Normalised fluorescence (AU)",
+    plot_title="",
+    ax=None,
+    cbar_ax=None,
+):
+    """Draw heatmap from an array of time series of traces
+
+    Parameters
+    ----------
+    trace_df : pandas.DataFrame
+        Time series of traces (rows = cells, columns = time points).
+    trace_name : string
+        Name of trace being plotted, e.g. 'flavin'.
+    buddings_df : pandas.DataFrame
+        Birth mask (rows = cells, columns = time points).  Elements should be
+        0 or 1.
+    cmap : matplotlib ColorMap
+        Colour map for heatmap.
+    unit_scaling : int or float
+        Unit scaling factor, e.g. 1/60 to convert minutes to hours.
+    xtick_step : int or float
+        Interval length, in unit time, to draw x axis ticks.
+    scale : bool
+        Whether to use standard scaler to scale the trace time series.
+    robust : bool
+        If True, the colour map range is computed with robust quantiles instead
+        of the extreme values.
+    xlabel : string
+        x axis label.
+    ylabel : string
+        y axis label.
+    cbarlabel : string
+        Colour bar label.
+    plot_title : string
+        Plot title.
+    ax : matplotlib Axes
+        Axes in which to draw the plot, otherwise use the currently active Axes.
+    cbar_ax : matplotlib Axes
+        Axes in which to draw the colour bar, otherwise take space from the main
+        Axes.
+
+    Returns
+    -------
+    ax : matplotlib Axes
+        Axes object with the heatmap.
+
+    Examples
+    --------
+    FIXME: Add docs.
+
+    """
+    plotter = _HeatmapPlotter(
+        trace_df,
+        trace_name,
+        buddings_df,
+        cmap,
+        unit_scaling,
+        xtick_step,
+        scale,
+        robust,
+        xlabel,
+        ylabel,
+        cbarlabel,
+        plot_title,
+    )
+    if ax is None:
+        ax = plt.gca()
+    plotter.plot(ax, cbar_ax)
+    return ax
diff --git a/routines/histogram.py b/routines/histogram.py
new file mode 100644
index 0000000000000000000000000000000000000000..eff7f8d720b707c694ec655da515ab53f3e0d539
--- /dev/null
+++ b/routines/histogram.py
@@ -0,0 +1,137 @@
+#!/usr/bin/env python3
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+class _HistogramPlotter:
+    """Draw histogram"""
+
+    def __init__(
+        self,
+        values,
+        label,
+        color,
+        binsize,
+        lognormal,
+        lognormal_base,
+        xlabel,
+        ylabel,
+        plot_title,
+    ):
+        # Define attributes from arguments
+        self.values = values
+        self.label = label
+        self.color = color
+        self.binsize = binsize
+        self.lognormal = lognormal
+        self.lognormal_base = lognormal_base
+        self.xlabel = xlabel
+        self.ylabel = ylabel
+        self.plot_title = plot_title
+
+        # Define median
+        self.median = np.median(self.values)
+
+        # Define bins
+        if self.lognormal:
+            self.bins = np.logspace(
+                0,
+                np.ceil(
+                    np.log(np.nanmax(values)) / np.log(self.lognormal_base)
+                ),
+                base=self.lognormal_base,
+            )  # number of bins will be 50 by default, as it's the default in np.logspace
+        else:
+            self.bins = np.arange(
+                self.binsize * (np.nanmin(values) // self.binsize - 2),
+                self.binsize * (np.nanmax(values) // self.binsize + 2),
+                self.binsize,
+            )
+
+    def plot(self, ax):
+        """Plot histogram onto specified Axes."""
+        ax.set_ylabel(self.ylabel)
+        ax.set_xlabel(self.xlabel)
+        ax.set_title(self.plot_title)
+
+        if self.lognormal:
+            ax.set_xscale("log")
+        ax.hist(
+            self.values,
+            self.bins,
+            alpha=0.5,
+            color=self.color,
+            label=self.label,
+        )
+        ax.axvline(
+            self.median,
+            color=self.color,
+            alpha=0.75,
+            label="median " + self.label,
+        )
+        ax.legend(loc="upper right")
+
+
+def histogram(
+    values,
+    label,
+    color="b",
+    binsize=5,
+    lognormal=False,
+    lognormal_base=10,
+    xlabel="Time (min)",
+    ylabel="Number of occurences",
+    plot_title="Distribution",
+    ax=None,
+):
+    """Draw histogram with median indicated
+
+    Parameters
+    ----------
+    values : array_like
+        Input values for histogram
+    label : string
+        Name of value being plotting, e.g. cell division cycle length.
+    color : string
+        Colour of bars.
+    binsize : float
+        Bin size.
+    lognormal : bool
+        Whether to use a log scale for the horizontal axis.
+    lognormal_base : float
+        Base of the log scale, if lognormal is True.
+    xlabel : string
+        x axis label.
+    ylabel : string
+        y axis label.
+    plot_title : string
+        Plot title.
+    ax : matplotlib Axes
+        Axes in which to draw the plot, otherwise use the currently active Axes.
+
+    Returns
+    -------
+    ax : matplotlib Axes
+        Axes object with the histogram.
+
+    Examples
+    --------
+    FIXME: Add docs.
+
+    """
+    plotter = _HistogramPlotter(
+        values,
+        label,
+        color,
+        binsize,
+        lognormal,
+        lognormal_base,
+        xlabel,
+        ylabel,
+        plot_title,
+    )
+    if ax is None:
+        ax = plt.gca()
+    plotter.plot(ax)
+    return ax
diff --git a/routines/mean_plot.py b/routines/mean_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..25959ef100d75fc366176cb4db973596331f7970
--- /dev/null
+++ b/routines/mean_plot.py
@@ -0,0 +1,132 @@
+#!/usr/bin/env python3
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from plottingabc import BasePlotter
+
+
+class _MeanPlotter(BasePlotter):
+    """Draw mean time series plus standard error."""
+
+    def __init__(
+        self,
+        trace_df,
+        trace_name,
+        unit_scaling,
+        label,
+        mean_color,
+        error_color,
+        mean_linestyle,
+        mean_marker,
+        xlabel,
+        ylabel,
+        plot_title,
+    ):
+        super().__init__(trace_name, unit_scaling, xlabel, plot_title)
+        # Define attributes from arguments
+        self.trace_df = trace_df
+        self.label = label
+        self.mean_color = mean_color
+        self.error_color = error_color
+        self.mean_linestyle = mean_linestyle
+        self.mean_marker = mean_marker
+
+        # Define some labels
+        self.ylabel = ylabel
+
+        # Mean and standard error
+        self.trace_time = (
+            np.array(self.trace_df.columns, dtype=float) * self.unit_scaling
+        )
+        self.mean_ts = self.trace_df.mean(axis=0)
+        self.stderr = self.trace_df.std(axis=0) / np.sqrt(len(self.trace_df))
+
+    def plot(self, ax):
+        """Draw lines and shading on provided Axes."""
+        super().plot(ax)
+        ax.plot(
+            self.trace_time,
+            self.mean_ts,
+            color=self.mean_color,
+            alpha=0.75,
+            linestyle=self.mean_linestyle,
+            marker=self.mean_marker,
+            label="Mean, " + self.label,
+        )
+        ax.fill_between(
+            self.trace_time,
+            self.mean_ts - self.stderr,
+            self.mean_ts + self.stderr,
+            color=self.error_color,
+            alpha=0.5,
+            label="Standard error, " + self.label,
+        )
+        ax.legend(loc="upper right")
+
+
+def mean_plot(
+    trace_df,
+    trace_name="flavin",
+    unit_scaling=1,
+    label="wild type",
+    mean_color="b",
+    error_color="lightblue",
+    mean_linestyle="-",
+    mean_marker="",
+    xlabel="Time (min)",
+    ylabel="Normalised flavin fluorescence (AU)",
+    plot_title="",
+    ax=None,
+):
+    """Plot mean time series of a DataFrame, with standard error shading.
+
+    Parameters
+    ----------
+    trace_df : pandas.DataFrame
+        Time series of traces (rows = cells, columns = time points).
+    trace_name : string
+        Name of trace being plotted, e.g. 'flavin'.
+    unit_scaling : int or float
+        Unit scaling factor, e.g. 1/60 to convert minutes to hours.
+    label : string
+        Name of group being plotted, e.g. a strain name.
+    mean_color : string
+        matplotlib colour string for the mean trace.
+    error_color : string
+        matplotlib colour string for the standard error shading.
+    mean_linestyle : string
+        matplotlib linestyle argument for the mean trace.
+    mean_marker : string
+        matplotlib marker argument for the mean trace.
+    xlabel : string
+        x axis label.
+    ylabel : string
+        y axis label.
+    plot_title : string
+        Plot title.
+    ax : matplotlib Axes
+        Axes in which to draw the plot, otherwise use the currently active Axes.
+
+    Examples
+    --------
+    FIXME: Add docs.
+
+    """
+    plotter = _MeanPlotter(
+        trace_df,
+        trace_name,
+        unit_scaling,
+        label,
+        mean_color,
+        error_color,
+        mean_linestyle,
+        mean_marker,
+        xlabel,
+        ylabel,
+        plot_title,
+    )
+    if ax is None:
+        ax = plt.gca()
+    plotter.plot(ax)
+    return ax
diff --git a/routines/median_plot.py b/routines/median_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..cab5fc3753bd064bbc44c6fc3aa039f76869b1f5
--- /dev/null
+++ b/routines/median_plot.py
@@ -0,0 +1,133 @@
+#!/usr/bin/env python3
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from plottingabc import BasePlotter
+
+
+class _MedianPlotter(BasePlotter):
+    """Draw median time series plus interquartile range."""
+
+    def __init__(
+        self,
+        trace_df,
+        trace_name,
+        unit_scaling,
+        label,
+        median_color,
+        error_color,
+        median_linestyle,
+        median_marker,
+        xlabel,
+        ylabel,
+        plot_title,
+    ):
+        super().__init__(trace_name, unit_scaling, xlabel, plot_title)
+        # Define attributes from arguments
+        self.trace_df = trace_df
+        self.label = label
+        self.median_color = median_color
+        self.error_color = error_color
+        self.median_linestyle = median_linestyle
+        self.median_marker = median_marker
+
+        # Define some labels
+        self.ylabel = ylabel
+
+        # Median and interquartile range
+        self.trace_time = (
+            np.array(self.trace_df.columns, dtype=float) * self.unit_scaling
+        )
+        self.median_ts = self.trace_df.median(axis=0)
+        self.quartile1_ts = self.trace_df.quantile(0.25)
+        self.quartile3_ts = self.trace_df.quantile(0.75)
+
+    def plot(self, ax):
+        """Draw lines and shading on provided Axes."""
+        super().plot(ax)
+        ax.plot(
+            self.trace_time,
+            self.median_ts,
+            color=self.median_color,
+            alpha=0.75,
+            linestyle=self.median_linestyle,
+            marker=self.median_marker,
+            label="Median, " + self.label,
+        )
+        ax.fill_between(
+            self.trace_time,
+            self.quartile1_ts,
+            self.quartile3_ts,
+            color=self.error_color,
+            alpha=0.5,
+            label="Interquartile range, " + self.label,
+        )
+        ax.legend(loc="upper right")
+
+
+def median_plot(
+    trace_df,
+    trace_name="flavin",
+    unit_scaling=1,
+    label="wild type",
+    median_color="b",
+    error_color="lightblue",
+    median_linestyle="-",
+    median_marker="",
+    xlabel="Time (min)",
+    ylabel="Normalised flavin fluorescence (AU)",
+    plot_title="",
+    ax=None,
+):
+    """Plot median time series of a DataFrame, with interquartile range
+    shading.
+
+    Parameters
+    ----------
+    trace_df : pandas.DataFrame
+        Time series of traces (rows = cells, columns = time points).
+    trace_name : string
+        Name of trace being plotted, e.g. 'flavin'.
+    unit_scaling : int or float
+        Unit scaling factor, e.g. 1/60 to convert minutes to hours.
+    label : string
+        Name of group being plotted, e.g. a strain name.
+    median_color : string
+        matplotlib colour string for the median trace.
+    error_color : string
+        matplotlib colour string for the interquartile range shading.
+    median_linestyle : string
+        matplotlib linestyle argument for the median trace.
+    median_marker : string
+        matplotlib marker argument for the median trace.
+    xlabel : string
+        x axis label.
+    ylabel : string
+        y axis label.
+    plot_title : string
+        Plot title.
+    ax : matplotlib Axes
+        Axes in which to draw the plot, otherwise use the currently active Axes.
+
+    Examples
+    --------
+    FIXME: Add docs.
+    """
+    plotter = _MedianPlotter(
+        trace_df,
+        trace_name,
+        unit_scaling,
+        label,
+        median_color,
+        error_color,
+        median_linestyle,
+        median_marker,
+        xlabel,
+        ylabel,
+        plot_title,
+    )
+    if ax is None:
+        ax = plt.gca()
+    plotter.plot(ax)
+    return ax
diff --git a/routines/plot_utils.py b/routines/plot_utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..cb33009995d1235805c6780b6db1f6183bb5eb87
--- /dev/null
+++ b/routines/plot_utils.py
@@ -0,0 +1,16 @@
+#!/usr/bin/env python3
+
+import numpy as np
+from matplotlib import cm, colors
+
+
+def generate_palette_map(df):
+    """Create a palette map based on the strains in a dataframe"""
+    strain_list = np.unique(df.index.get_level_values("strain"))
+    palette_cm = cm.get_cmap("Set1", len(strain_list) + 1)
+    palette_rgb = [
+        colors.rgb2hex(palette_cm(index / len(strain_list))[:3])
+        for index, _ in enumerate(strain_list)
+    ]
+    palette_map = dict(zip(strain_list, palette_rgb))
+    return palette_map
diff --git a/routines/plottingabc.py b/routines/plottingabc.py
new file mode 100644
index 0000000000000000000000000000000000000000..97b89aa7e95e01e5b937003233f8ba3d37da7f52
--- /dev/null
+++ b/routines/plottingabc.py
@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+
+from abc import ABC
+
+
+class BasePlotter(ABC):
+    """Base class for plotting handler classes"""
+
+    def __init__(self, trace_name, unit_scaling, xlabel, plot_title):
+        """Common attributes"""
+        self.trace_name = trace_name
+        self.unit_scaling = unit_scaling
+
+        self.xlabel = xlabel
+        self.ylabel = None
+        self.plot_title = plot_title
+
+    def plot(self, ax):
+        """Template for drawing on provided Axes"""
+        ax.set_ylabel(self.ylabel)
+        ax.set_xlabel(self.xlabel)
+        ax.set_title(self.plot_title)
+        # Derived classes extends this with plotting functions
+
+
+# TODO: something about the plotting functions at the end of the modules.
+# Decorator?
diff --git a/routines/single_birth_plot.py b/routines/single_birth_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..b571d15376467c5bbf0667b1820644a8202f1dd5
--- /dev/null
+++ b/routines/single_birth_plot.py
@@ -0,0 +1,136 @@
+#!/usr/bin/env python3
+
+import matplotlib.pyplot as plt
+
+from single_plot import _SinglePlotter
+
+
+class _SingleBirthPlotter(_SinglePlotter):
+    """Draw a line plot of a single time series, but with buddings overlaid"""
+
+    def __init__(
+        self,
+        trace_timepoints,
+        trace_values,
+        trace_name,
+        birth_mask,
+        unit_scaling,
+        trace_color,
+        birth_color,
+        trace_linestyle,
+        birth_linestyle,
+        xlabel,
+        ylabel,
+        plot_title,
+    ):
+        # Define attributes from arguments
+        super().__init__(
+            trace_timepoints,
+            trace_values,
+            trace_name,
+            unit_scaling,
+            trace_color,
+            trace_linestyle,
+            xlabel,
+            ylabel,
+            plot_title,
+        )
+        # Add some more attributes useful for buddings
+        self.birth_mask = birth_mask
+        self.birth_color = birth_color
+        self.birth_linestyle = birth_linestyle
+
+    def plot(self, ax):
+        """Draw the line plots on the provided Axes."""
+        trace_time = self.trace_timepoints * self.unit_scaling
+        super().plot(ax)
+        birth_mask_bool = self.birth_mask.astype(bool)
+        for occurence, birth_time in enumerate(trace_time[birth_mask_bool]):
+            if occurence == 0:
+                label = "birth event"
+            else:
+                label = None
+            ax.axvline(
+                birth_time,
+                color=self.birth_color,
+                linestyle=self.birth_linestyle,
+                label=label,
+            )
+        ax.legend()
+
+
+def single_birth_plot(
+    trace_timepoints,
+    trace_values,
+    trace_name="flavin",
+    birth_mask=None,
+    unit_scaling=1,
+    trace_color="b",
+    birth_color="k",
+    trace_linestyle="-",
+    birth_linestyle="--",
+    xlabel="Time (min)",
+    ylabel="Normalised flavin fluorescence (AU)",
+    plot_title="",
+    ax=None,
+):
+    """Plot time series of trace, overlaid with buddings
+
+    Parameters
+    ----------
+    trace_timepoints : array_like
+        Time points (as opposed to the actual times in time units)
+    trace_values : array_like
+        Trace to plot
+    trace_name : string
+        Name of trace being plotted, e.g. 'flavin'.
+    birth_mask : array_like
+        Mask to indicate where buddings are. Expect values of '0' and '1' or
+        'False' and 'True' in the elements.
+    unit_scaling : int or float
+        Unit scaling factor, e.g. 1/60 to convert minutes to hours.
+    trace_color : string
+        matplotlib colour string for the trace
+    birth_color : string
+        matplotlib colour string for the vertical lines indicating buddings
+    trace_linestyle : string
+        matplotlib linestyle argument for the trace
+    birth_linestyle : string
+        matplotlib linestyle argument for the vertical lines indicating buddings
+    xlabel : string
+        x axis label.
+    ylabel : string
+        y axis label.
+    plot_title : string
+        Plot title.
+    ax : matplotlib Axes
+        Axes in which to draw the plot, otherwise use the currently active Axes.
+
+    Returns
+    -------
+    ax : matplotlib Axes
+        Axes object with the plot.
+
+    Examples
+    --------
+    FIXME: Add docs.
+
+    """
+    plotter = _SingleBirthPlotter(
+        trace_timepoints,
+        trace_values,
+        trace_name,
+        birth_mask,
+        unit_scaling,
+        trace_color,
+        birth_color,
+        trace_linestyle,
+        birth_linestyle,
+        xlabel,
+        ylabel,
+        plot_title,
+    )
+    if ax is None:
+        ax = plt.gca()
+    plotter.plot(ax)
+    return ax
diff --git a/routines/single_plot.py b/routines/single_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..4e6952481f396a6d54f22135bbb77eaa35226dbf
--- /dev/null
+++ b/routines/single_plot.py
@@ -0,0 +1,106 @@
+#!/usr/bin/env python3
+
+import matplotlib.pyplot as plt
+
+from plottingabc import BasePlotter
+
+
+class _SinglePlotter(BasePlotter):
+    """Draw a line plot of a single time series."""
+
+    def __init__(
+        self,
+        trace_timepoints,
+        trace_values,
+        trace_name,
+        unit_scaling,
+        trace_color,
+        trace_linestyle,
+        xlabel,
+        ylabel,
+        plot_title,
+    ):
+        super().__init__(trace_name, unit_scaling, xlabel, plot_title)
+        # Define attributes from arguments
+        self.trace_timepoints = trace_timepoints
+        self.trace_values = trace_values
+        self.trace_color = trace_color
+        self.trace_linestyle = trace_linestyle
+
+        # Define some labels
+        self.ylabel = ylabel
+
+    def plot(self, ax):
+        """Draw the line plot on the provided Axes."""
+        super().plot(ax)
+        ax.plot(
+            self.trace_timepoints * self.unit_scaling,
+            self.trace_values,
+            color=self.trace_color,
+            linestyle=self.trace_linestyle,
+            label=self.trace_name,
+        )
+
+
+def single_plot(
+    trace_timepoints,
+    trace_values,
+    trace_name="flavin",
+    unit_scaling=1,
+    trace_color="b",
+    trace_linestyle="-",
+    xlabel="Time (min)",
+    ylabel="Normalised flavin fluorescence (AU)",
+    plot_title="",
+    ax=None,
+):
+    """Plot time series of trace.
+
+    Parameters
+    ----------
+    trace_timepoints : array_like
+        Time points (as opposed to the actual times in time units).
+    trace_values : array_like
+        Trace to plot.
+    trace_name : string
+        Name of trace being plotted, e.g. 'flavin'.
+    unit_scaling : int or float
+        Unit scaling factor, e.g. 1/60 to convert minutes to hours.
+    trace_color : string
+        matplotlib colour string, specifies colour of line plot.
+    trace_linestyle : string
+        matplotlib linestyle argument.
+    xlabel : string
+        x axis label.
+    ylabel : string
+        y axis label.
+    plot_title : string
+        Plot title.
+    ax : matplotlib Axes
+        Axes in which to draw the plot, otherwise use the currently active Axes.
+
+    Returns
+    -------
+    ax : matplotlib Axes
+        Axes object with the plot.
+
+    Examples
+    --------
+    FIXME: Add docs.
+
+    """
+    plotter = _SinglePlotter(
+        trace_timepoints,
+        trace_values,
+        trace_name,
+        unit_scaling,
+        trace_color,
+        trace_linestyle,
+        xlabel,
+        ylabel,
+        plot_title,
+    )
+    if ax is None:
+        ax = plt.gca()
+    plotter.plot(ax)
+    return ax