From 429e2fff1d9cca93580ef0a661965038460020e6 Mon Sep 17 00:00:00 2001
From: Eric Teunis de Boone <ericteunis@deboone.nl>
Date: Tue, 20 Dec 2022 15:26:10 +0100
Subject: [PATCH] ZH: plot showing (measured-actual) residuals for true_phases

---
 .../bb_measure_true_phase.py                  | 46 +++++++++++++++++--
 1 file changed, 43 insertions(+), 3 deletions(-)

diff --git a/simulations/airshower_beacon_simulation/bb_measure_true_phase.py b/simulations/airshower_beacon_simulation/bb_measure_true_phase.py
index fc8fd46..f766e1b 100755
--- a/simulations/airshower_beacon_simulation/bb_measure_true_phase.py
+++ b/simulations/airshower_beacon_simulation/bb_measure_true_phase.py
@@ -80,7 +80,7 @@ if __name__ == "__main__":
         # Plot True Phases at their locations
         if show_plots or fig_dir:
             fig, ax = plt.subplots()
-            spatial_unit=None
+            spatial_unit='m'
             fig.suptitle('Clock phases\nf_beacon= {:2.0f}MHz'.format(f_beacon*1e3))
 
             antenna_locs = list(zip(*[(ant.x, ant.y) for ant in antennas]))
@@ -88,16 +88,56 @@ if __name__ == "__main__":
             ax.set_ylabel('y' if spatial_unit is None else 'y [{}]'.format(spatial_unit))
             scatter_kwargs = {}
             scatter_kwargs['cmap'] = 'inferno'
-            scatter_kwargs['vmin'] = -np.pi
-            scatter_kwargs['vmax'] = +np.pi
+            #scatter_kwargs['vmin'] = -np.pi
+            #scatter_kwargs['vmax'] = +np.pi
             color_label='$\\varphi(\\sigma_t)$ [rad]'
 
             sc = ax.scatter(*antenna_locs, c=true_phases, **scatter_kwargs)
             fig.colorbar(sc, ax=ax, label=color_label)
 
+            if False:
+                for i, ant in enumerate(antennas):
+                    ax.text(ant.x, ant.y, ant.name)
+
+            if not True:
+                ax.plot(tx.x, tx.y, 'X', color='k', markeredgecolor='white')
+
             if fig_dir:
                 fig.savefig(path.join(fig_dir, __file__ + f".F{freq_name}.pdf"))
 
+        # Plot True Phases - Actual True Phases at their location
+        if show_plots or fig_dir:
+            fig, ax = plt.subplots()
+            fig.suptitle('Clock phase Residuals\nf_beacon={:2.0f}MHz'.format(f_beacon*1e3))
+
+            actual_true_phases = np.array([ -2*np.pi*a.attrs['clock_offset']*f_beacon for a in antennas ])
+
+            # Modify actual_true_phases, the same way as true_phases
+            # was modified
+            if remove_absolute_phase_offset_first_antenna or remove_absolute_phase_offset_minimum:
+                if remove_absolute_phase_offset_first_antenna: # just take the first phase
+                    minimum_phase = actual_true_phases[0]
+                else: # take the minimum
+                    minimum_phase = np.min(actual_true_phases, axis=-1)
+
+                actual_true_phases -= minimum_phase
+                actual_true_phases = lib.phase_mod(actual_true_phases)
+
+            true_phase_residuals = lib.phase_mod(true_phases - actual_true_phases)
+
+            antenna_locs = list(zip(*[(ant.x, ant.y) for ant in antennas]))
+            ax.set_xlabel('x' if spatial_unit is None else 'x [{}]'.format(spatial_unit))
+            ax.set_ylabel('y' if spatial_unit is None else 'y [{}]'.format(spatial_unit))
+            scatter_kwargs = {}
+            scatter_kwargs['cmap'] = 'inferno'
+            color_label='$\\Delta\\varphi(\\sigma_t) = \\varphi_{meas} - \\varphi_{true}$ [rad]'
+
+            sc = ax.scatter(*antenna_locs, c=true_phase_residuals, **scatter_kwargs)
+            fig.colorbar(sc, ax=ax, label=color_label)
+
+            if fig_dir:
+                fig.savefig(path.join(fig_dir, __file__ + f".residual.F{freq_name}.pdf"))
+
     print(f"True phases written to", antennas_fname)
 
     if show_plots: