diff --git a/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py b/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py
index 6e088be..e45c24c 100755
--- a/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py
+++ b/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py
@@ -62,7 +62,6 @@ if __name__ == "__main__":
     sigma_phase_matrix = sigma_phase_matrix - first_row[:,np.newaxis]
     sigma_phase_matrix = lib.phase_mod(sigma_phase_matrix)
 
-
     # Except for the first row, these are all separate measurements
     # Condense into phase offset per antenna
     if True: # do not use the first row
@@ -97,44 +96,57 @@ if __name__ == "__main__":
     ##############################
     # Compare actual time shifts #
     ##############################
-    antenna_time_shifts = { a.name: a.attrs['clock_offset'] for a in antennas }
+    antenna_time_shifts = { a.name: a.attrs['clock_offset'] for a in sorted(antennas, key=lambda a: int(a.name)) }
 
     if True:
-        # show means and std
-        fig, axs = plt.subplots(1,2, sharey=True)
-        fig.suptitle("Comparison Measured and Actual Phases")
-        axs[0].set_xlabel("Antenna no.")
-        axs[0].set_ylabel("Antenna Phase $\\Delta_\\varphi$")
-        axs[1].set_xlabel("#")
-        if True:
-            forward = lambda x: x/(2*np.pi*f_beacon)
-            inverse = lambda x: 2*np.pi*x*f_beacon
-            secax = axs[1].secondary_yaxis('right', functions=(forward, inverse))
-            secax.set_ylabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]')
+        actual_phase_shifts = [ -1*lib.phase_mod(2*np.pi*f_beacon*v) for k,v in antenna_time_shifts.items() ]
+        antenna_names = [int(k)-1 for k,v in antenna_time_shifts.items() ]
 
-        l = axs[0].errorbar(np.arange(N_ant), mean_sigma_phase, yerr=std_sigma_phase, marker='4', alpha=0.7, ls='none')
+        for i in range(2):
+            plot_residuals = i == 1
+            colors = ['blue', 'orange']
 
-        axs[1].hist(mean_sigma_phase, bins='sqrt', density=False, orientation='horizontal', color=l[0].get_color(), histtype='step')
-
-        # Actual time shifts
-        if True:
-            actual_phase_shifts = [ -1*lib.phase_mod(2*np.pi*f_beacon*v) for k,v in antenna_time_shifts.items() ]
-            antenna_names = [int(k)-1 for k,v in antenna_time_shifts.items() ]
-
-            # make sure to keep the same offset
-            if False:
-                phase_offset = mean_sigma_phase[0] - actual_phase_shifts[0]
-
-                actual_phase_shifts += phase_offset
-
-            l = axs[0].plot(antenna_names, actual_phase_shifts, ls='none', marker='3', alpha=0.8)
+            fig, axs = plt.subplots(1,2, sharey=True)
 
             if True:
-                axs[1].hist(actual_phase_shifts, bins='sqrt', density=False, orientation='horizontal', ls='dashed', color=l[0].get_color(), histtype='step')
+                forward = lambda x: x/(2*np.pi*f_beacon)
+                inverse = lambda x: 2*np.pi*x*f_beacon
+                secax = axs[-1].secondary_yaxis('right', functions=(forward, inverse))
+                secax.set_ylabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]')
 
-        fig.tight_layout()
-        if fig_dir:
-            fig.savefig(path.join(fig_dir, __file__ + f".residuals.pdf"))
+            if plot_residuals:
+                phase_residuals = lib.phase_mod(mean_sigma_phase - actual_phase_shifts)
+                fig.suptitle("Difference between Measured and Actual phases\n for Antenna $i$")
+                axs[0].set_ylabel("Antenna Phase Residual $\\Delta_\\varphi$")
+            else:
+                fig.suptitle("Comparison Measured and Actual phases\n for Antenna $i$")
+                axs[0].set_ylabel("Antenna Phase $\\Delta_\\varphi$")
+
+
+            i=0
+            axs[i].set_xlabel("Antenna no.")
+            if plot_residuals:
+                axs[i].plot(np.arange(N_ant), phase_residuals, alpha=0.6, ls='none', marker='x', color=colors[0])
+            else:
+                axs[i].errorbar(np.arange(N_ant), mean_sigma_phase, yerr=std_sigma_phase, marker='4', alpha=0.7, ls='none', color=colors[0], label='Measured')
+                axs[i].plot(antenna_names, actual_phase_shifts, ls='none', marker='3', alpha=0.8, color=colors[1], label='Actual')
+                axs[i].legend()
+
+
+            i=1
+            axs[i].set_xlabel("#")
+            if plot_residuals:
+                axs[i].hist(phase_residuals, bins='sqrt', alpha=0.8, color=colors[0], orientation='horizontal')
+            else:
+                axs[i].hist(mean_sigma_phase, bins='sqrt', density=False, orientation='horizontal', color=colors[0], histtype='step', label='Measured')
+                axs[i].hist(actual_phase_shifts, bins='sqrt', density=False, orientation='horizontal', ls='dashed', color=colors[1], histtype='step', label='Actual')
+
+            fig.tight_layout()
+            if fig_dir:
+                extra_name = "measured"
+                if plot_residuals:
+                    extra_name = "residuals"
+                fig.savefig(path.join(fig_dir, __file__ + f".phase.{extra_name}.pdf"))
 
     ##########################
     ##########################