diff --git a/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py b/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py
index 4295b28..b2efc53 100755
--- a/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py
+++ b/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py
@@ -8,6 +8,7 @@ import numpy as np
 
 import aa_generate_beacon as beacon
 import lib
+from lib import figlib
 
 if __name__ == "__main__":
     from os import path
@@ -155,42 +156,30 @@ if __name__ == "__main__":
                 ##
                 ## Histogram
                 ##
-                fig, axs = plt.subplots(2, 1, sharex=True, figsize=figsize)
-                colors = ['blue', 'orange']
-
-                if True:
-                    phase2time = lambda x: x/(2*np.pi*f_beacon)
-                    time2phase = lambda x: 2*np.pi*x*f_beacon
-                    secax = axs[0].secondary_xaxis('top', functions=(phase2time, time2phase))
-                    secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]')
+                fig = figlib.phase_comparison_figure(
+                        loc_c,
+                        actual_clock_phases,
+                        plot_residuals=plot_residuals,
+                        f_beacon=f_beacon,
+                        figsize=figsize
+                        )
 
                 if plot_residuals:
                     fig.suptitle("Difference between Measured and True Clock phases")
                 else:
                     fig.suptitle("Comparison Measured and True Clock Phases")
 
+                axs = fig.get_axes()
                 axs[-1].set_xlabel(f'Antenna {title} {color_label}')
 
                 #
-                hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step')
-
                 i=0
-                axs[i].set_ylabel("#")
-                axs[i].hist(loc_c, color=colors[0], label='Measured', ls='solid', **hist_kwargs)
-                if not plot_residuals: # also plot the true clock phases
-                    axs[i].hist(actual_clock_phases, color=colors[1], label='Actual', ls='dashed', **hist_kwargs)
-                    #axs[i].legend()
+                secax = axs[i].child_axes[0]
+                secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]')
 
                 #
-                plot_kwargs = dict(alpha=0.6, ls='none')
-
                 i=1
                 axs[i].set_ylabel("Antenna no.")
-                axs[i].plot(loc_c, np.arange(len(loc_c)), marker='x' if plot_residuals else '3', color=colors[0], label='Measured', **plot_kwargs)
-
-                if not plot_residuals: # also plot the true clock phases
-                    axs[i].plot(actual_clock_phases, np.arange(len(loc_c)), marker='4', color=colors[1], label='Actual', **plot_kwargs)
-                    axs[i].legend()
 
                 # Save figure
                 if fig_dir:
diff --git a/simulations/airshower_beacon_simulation/lib/figlib.py b/simulations/airshower_beacon_simulation/lib/figlib.py
new file mode 100644
index 0000000..bdd8bde
--- /dev/null
+++ b/simulations/airshower_beacon_simulation/lib/figlib.py
@@ -0,0 +1,65 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+def phase_comparison_figure(
+        measured_phases,
+        true_phases,
+        plot_residuals=True,
+        f_beacon=None,
+        hist_kwargs={},
+        sc_kwargs={},
+        colors=['blue', 'orange'],
+        legend_on_scatter=True,
+        **fig_kwargs
+        ):
+    """
+    Create a figure comparing measured_phase against true_phase
+    by both plotting the values, and the residuals.
+    """
+    default_fig_kwargs = dict(sharex=True)
+    fig_kwargs = {**default_fig_kwargs, **fig_kwargs}
+
+    do_hist_plot = hist_kwargs is not False
+    do_scatter_plot = sc_kwargs is not False
+
+    fig, axs = plt.subplots(0+do_hist_plot+do_scatter_plot, 1, **fig_kwargs)
+
+    if not hasattr(axs, '__len__'):
+        axs = [axs]
+
+    if f_beacon:
+        phase2time = lambda x: x/(2*np.pi*f_beacon)
+        time2phase = lambda x: 2*np.pi*x*f_beacon
+        secax = axs[0].secondary_xaxis('top', functions=(phase2time, time2phase))
+        secax.set_xlabel('Time $\\varphi/(2\\pi f_{beac})$ [ns]')
+
+    # Histogram
+    if do_hist_plot:
+        i=0
+        default_hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step')
+        hist_kwargs = {**default_hist_kwargs, **hist_kwargs}
+
+        axs[i].set_ylabel("#")
+
+        _counts, _bins, _patches = axs[i].hist(measured_phases, color=colors[0], label='Measured', ls='solid', **hist_kwargs)
+        if not plot_residuals: # also plot the true clock phases
+            axs[i].hist(true_phases, color=colors[1], label='Actual', ls='dashed', **{**hist_kwargs, **dict(bins=_bins)})
+
+    # Scatter plot
+    if do_scatter_plot:
+        i=1
+        default_sc_kwargs = dict(alpha=0.6, ls='none')
+        sc_kwargs = {**default_sc_kwargs, **sc_kwargs}
+
+        axs[i].set_ylabel("Antenna no.")
+        axs[i].plot(measured_phases, np.arange(len(measured_phases)), marker='x' if plot_residuals else '3', color=colors[0], label='Measured', **sc_kwargs)
+
+        if not plot_residuals: # also plot the true clock phases
+            axs[i].plot(true_phases, np.arange(len(true_phases)), marker='4', color=colors[1], label='Actual', **sc_kwargs)
+
+        if not plot_residuals and legend_on_scatter:
+            axs[i].legend()
+
+    fig.tight_layout()
+
+    return fig