histograms the phase slices such as to build a histogram of the position
    distribution at each phase value.

    Parameters
    ----------
    phase_slices : ndarray
        2d array containing slices from many oscillations at each phase
    phase : ndarray
        1d array of phases corresponding to slices
    histbins : int, optional (default=200)
        number of bins to use in histogramming data
    show_plot : bool, optional (default=False)
        if true plots and shows the heatmap of the
        phase against the positon distribution

    Returns
    -------
    counts_array : ndarray
        2d array containing the number of counts varying with
        phase and position.
    bin_edges : ndarray
        positions of bin edges

def histogram_phase(phase_slices, phase, histbins=200, show_plot=False):
    """
    histograms the phase slices such as to build a histogram of the position
    distribution at each phase value.

    Parameters
    ----------
    phase_slices : ndarray
        2d array containing slices from many oscillations at each phase
    phase : ndarray
        1d array of phases corresponding to slices
    histbins : int, optional (default=200)
        number of bins to use in histogramming data
    show_plot : bool, optional (default=False)
        if true plots and shows the heatmap of the
        phase against the positon distribution

    Returns
    -------
    counts_array : ndarray
        2d array containing the number of counts varying with
        phase and position.
    bin_edges : ndarray
        positions of bin edges

    """
    counts_array = _np.zeros([len(phase), histbins])

    histedges = [phase_slices.min(), phase_slices.max()]
    for i, phase_slice in enumerate(phase_slices): # for each value of phase
        counts, bin_edges = _np.histogram(phase_slice, bins=histbins, range=histedges) # histogram the position distribution at that phase
        counts_array[i] = counts
    counts_array = _np.array(counts_array)
    counts_array_transposed = _np.transpose(counts_array).astype(float)

    if show_plot == True:
        fig = _plt.figure(figsize=(12, 6))
        ax = fig.add_subplot(111)
        ax.set_title('Phase Distribution')
        ax.set_xlabel("phase (°)")
        ax.set_ylabel("x")
        _plt.imshow(counts_array_transposed, cmap='hot', interpolation='nearest', extent=[phase[0], phase[-1], histedges[0], histedges[1]])
        ax.set_aspect('auto')
        _plt.show()

    return counts_array_transposed, bin_edges