ZH: improve plot with beacon delay + k_period

simulations/airshower_beacon_simulation/bc_beacon_periods.py

@@ -70,7 +70,7 @@ if __name__ == "__main__":
         # Get true phase diffs
         try:
             true_phases = np.array([ant.beacon_info[freq_name]['true_phase'] for ant in base])
-            true_phases_diff = lib.phase_mod(true_phases[0] - true_phases[1])
+            true_phases_diff = lib.phase_mod(lib.phase_mod(true_phases[0]) - lib.phase_mod(true_phases[1]))
         except IndexError:
             # true_phase not determined yet
             print(f"Missing true_phases for {freq_name} in baseline {base[0].name},{base[1].name}")
@@ -80,7 +80,7 @@ if __name__ == "__main__":
         time_diffs[i] = [true_phases_diff, k_period, f_beacon]
 
         # Plotting for one or two iterations
-        if show_plots and (i in [ 0, 1 ] or k_period > 3):
+        if show_plots and (i in [ 1, 57 ] or k_period > 3):
             # More than three periods is quite much so report it
             print('i',i,'k[T]',k_period, 'rest[ns]',t_rest, 'T[ns]',1/f_beacon, 'dT_coher[ns]', delta_t_coherence)
 
@@ -93,37 +93,72 @@ if __name__ == "__main__":
                 ax.plot(ks, maxima)
                 ax.plot(best_k, maxima[max_idx], marker='X')
 
-            # Delta between first timestamp from both antennas
-            delta_t_antennas = base[0].t[0] - base[1].t[0]
-
             # Delta t due to the beacon
-            if true_phases_diff is np.nan:
-                true_phases_diff = 0
-            delta_t_beacon = true_phases_diff/(2*np.pi*f_beacon)
+            # Note that we want to show some overlapping waveforms
+            # Therefore we use the phase from the original waveforms
+            # and not the true_phases (we lost t_d information there)
+            phases = np.array([ant.beacon_info[freq_name]['phase'] for ant in base])
+            phases_diff = lib.phase_mod(lib.phase_mod(phases[0]) - lib.phase_mod(phases[1]))
+            delta_t_beacon = phases_diff/(2*np.pi*f_beacon)
 
-            fig, ax = plt.subplots()
-            ax.set_title(
+
+            # Do we make it shared plot with both the
+            # signal and the beacon?
+            beacons = None
+            if True:
+                try:
+                    beacons = [ base[0].beacon, base[1].beacon ]
+                except:
+                    # No beacon waveforms available..
+                    pass
+
+            # Start the figure
+            fig, axs = plt.subplots(1+(beacons is not None), 1, sharex=True)
+            if beacons is None:
+                axs = [axs]
+
+            fig.suptitle(
                 ", ".join([
-                    f"$\\Delta$t0 [ns] : {delta_t_antennas:.2f}",
                     f"$\\Delta$t_beacon [ns]: {delta_t_beacon:.2f}",
+                    f"$\\Delta\\varphi$: {phases_diff:.4f}",
                     f"$\\Delta\\sigma_\\varphi$: {true_phases_diff:.4f}",
                     f"",
                 ])
             )
-            ax.set_xlabel('Sampling t [ns]')
-            ax.set_ylabel('Amplitude [a.u.]')
-            ax.plot(base[0].t, traces[0], label=f'Reference: {base[0].name}', alpha=0.5)
+
+
+            axs[-1].set_xlabel('t [ns]')
+            axs[0].set_ylabel('Amplitude [a.u.]')
+
             # plot vertical lines indicating f_beacon
-            min_t, max_t = base[0].t[0], base[0].t[-1]
+            min_t, max_t = min(base[0].t[0], base[1].t[0]), max(base[0].t[-1], base[1].t[-1])
             N_lines = int( (max_t - min_t)*f_beacon) +1
             for i, t in enumerate(np.arange(N_lines)/f_beacon):
-                ax.axvline( min_t + t, color='k', alpha=0.3)
+                for ax in axs:
+                    ax.axvline( min_t + t, color='k', alpha=0.3)
 
-            ax.plot(base[1].t + delta_t_antennas, traces[1], label=f'Original: {base[1].name} (t0 removed)', alpha=0.4, marker='+', ms=5)
-            ax.plot(base[1].t + delta_t_antennas + k_period/f_beacon + t_rest, traces[1], label='Coherence', alpha=0.3, marker='x', ms=5)
-            ax.plot(base[1].t + delta_t_antennas + k_period/f_beacon + delta_t_beacon, traces[1], label=f'$\\Delta t_\\varphi$ + $k={k_period:.0f}$ Periods', alpha=0.6)
+            # Plot traces
+            l1 = axs[0].plot(base[0].t, traces[0], label=f'Ref:  {base[0].name}', alpha=0.8)
+            l2 = axs[0].plot(base[1].t, traces[1], label=f'Orig: {base[1].name}', alpha=0.3, marker='+', ms=5)
+            axs[0].plot(base[0].t + k_period/f_beacon + t_rest, traces[1], label='Coherence', alpha=0.3, marker='x', ms=5)
+            l3 = axs[0].plot(base[1].t - delta_t_beacon +k_period/f_beacon, traces[1], label=f'$\\Delta t_{{\\sigma\\varphi}}$ + ($k={k_period:.0f}$)T', alpha=0.8)
 
-            ax.legend(fancybox=True, framealpha=0.5)
+            axs[0].legend(fancybox=True, framealpha=0.5)
+
+            # Plot beacon if available
+            if beacons is not None:
+                ax = axs[1]
+                ax.set_title("Original Beacons")
+                ax.plot(base[0].t, beacons[0], label=f'Ref: {base[0].name}', alpha=0.8, color=l1[0].get_color())
+                ax.plot(base[1].t, beacons[1], label=f'Orig: {base[1].name}', alpha=0.3, marker='+', ms=5, color=l2[0].get_color())
+                ax.plot(base[1].t -delta_t_beacon +k_period/f_beacon, beacons[1], label=f'$\\Delta t_{{\\sigma\\varphi}}$ + ($k={k_period:.0f}$)T', alpha=0.8, color=l3[0].get_color())
+
+
+            if True:
+                fig.savefig(__file__ + f"_i{i}_k{k_period}_zoomed_out.pdf")
+                ax.set_xlim(base[0].t[0]-1/f_beacon, base[0].t[0] + 5/f_beacon)
+
+            fig.savefig(__file__ + f"_i{i}_k{k_period}.pdf")
 
     # Save integer periods to antennas
     beacon.write_baseline_time_diffs_hdf5(antennas_fname, baselines, time_diffs[:,0], time_diffs[:,1], time_diffs[:,2])