ZH: integer period plotting

simulations/airshower_beacon_simulation/bc_beacon_periods.py

@@ -30,64 +30,83 @@ if __name__ == "__main__":
         ref_ant = antennas[0]
         baselines = list(zip_longest([], antennas, fillvalue=ref_ant))
 
-    integer_periods = None
-    # read integer ks from file if possible
-    # and save beacon_phase_true
-    with h5py.File(antennas_fname, 'a') as fp:
-        for i, ant in enumerate(antennas):
-            name = ant.name
-            # set true beacon_phase
-            fp['antennas'][name].attrs['beacon_phase_true'] = true_phases[i]
+    freq_names = antennas[0].beacon_info.keys()
+    if len(freq_names) > 1:
+        raise NotImplementedError
 
-        # read integer period from file
-        if True and 'beacon_ks' in fp:
-            integer_periods = np.array(fp['beacon_ks'])
+    freq_name = next(iter(freq_names))
 
 
     # Determine integer multiple of periods to shift
-    if integer_periods is None:
     integer_periods = np.empty( (len(baselines), 3) )
     for i, base in enumerate(baselines):
-            # Delta between first timestamp from both antennas
-            delta_t_a = base[0].t[0] - base[1].t[0]
-            # + phase difference
-            delta_t_p = np.diff([ant.attrs['beacon_phase_true'] for ant in base])[0]/(2*np.pi*f_beacon)
-
-            sampling_dt = (base[1].t[1] - base[1].t[0])
-
-            print("DT(A,P)", delta_t_a, delta_t_p, 1/f_beacon)
-
-            # which traces to keep track of
-            traces = [ base[0].Ex, base[1].Ex ]
-
-            # how many samples to shift
-            ks, maxima = lib.coherence_sum_maxima(-1*traces[0], -1*traces[1])
-            max_idx = np.argmax(maxima)
-            delta_t_c = sampling_dt*ks[max_idx] # ns
-            print("K", ks[max_idx], sampling_dt, '=', delta_t_c)
-
-            k, rest = np.divmod(delta_t_c, f_beacon)
-            integer_periods[i] = [int(base[0].name), int(base[1].name), k]
-
-
-            print(k, rest*f_beacon, delta_t_p)
-
-            # Only continue for two random combinations
-            if i not in [ 50, 51 ]:
+        if i not in [98, 99]:
             continue
 
+        # which traces to keep track of
+        traces = [ base[0].E_AxB, base[1].E_AxB ]
+
+        # how many samples do we need to shift
+        sampling_dt = (base[1].t[1] - base[1].t[0]) # ns
+        ks, maxima = lib.coherence_sum_maxima(traces[0], traces[1])
+        max_idx = np.argmax(maxima)
+        best_k = ks[max_idx]
+        delta_t_coherence = sampling_dt*best_k # ns
+
+        print("K", best_k, sampling_dt, '=', delta_t_coherence)
+
+
+        # get the amount of periods to move
+        f_beacon = base[0].beacon_info[freq_name]['freq']
+        k_period, rest = np.divmod(delta_t_coherence, 1/f_beacon)
+
+        # always keep the reference before traces[1]
+        if rest < 0:
+            k_period -= 1
+
+        # save k_period with antenna names
+        integer_periods[i] = [int(base[0].name), int(base[1].name), k_period]
+
+        if i in [ 98, 99 ]:
+            print('i',i,'k[T]',k_period, 'rest[ns]',rest, 'T[ns]',1/f_beacon)
+
+            # Show correlation maxima plot
+            if True:
                 fig, ax = plt.subplots()
+                ax.set_title(f"Correlation Maxima {i}")
                 ax.set_xlabel("k")
                 ax.set_ylabel("Maximum correlation")
                 ax.plot(ks, maxima)
-            ax.plot(ks[max_idx], maxima[max_idx], marker='X')
+                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
+            try:
+                true_phases = np.array([ant.beacon_info[freq_name]['true_phase'] for ant in base])
+                delta_true_phases = lib.phase_mod(true_phases[0] - true_phases[1])
+                delta_t_beacon = delta_true_phases/(2*np.pi*f_beacon)
+            except e:
+                # freq_name not found
+                # simply continue and set it them 0
+                print("No beacon")
+                delta_true_phases = 0
+                delta_t_beacon = 0
+
+            print("t0[ns]", delta_t_antennas, "t_beacon[ns]", delta_t_beacon, "phase", delta_true_phases)
             fig, ax = plt.subplots()
-            dt = base[1].t[1] - base[1].t[0]
             ax.set_xlabel('t')
-            ax.plot(base[0].t, traces[0], label='Reference')
-            ax.plot(base[1].t, traces[1], label='Original', alpha=0.4)
-            ax.plot(base[1].t + delta_t_a + delta_t_c, traces[1], label='Coherence', alpha=0.6)
+            ax.plot(base[0].t, traces[0], label=f'Reference {base[0].name}', alpha=0.5)
+            # plot vertical lines indicating f_beacon
+            min_t, max_t = base[0].t[0], base[0].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.5, label=None if i!=0 else 'P_beacon')
+
+            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 + 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='Beacon only + Periods', alpha=0.6)
 
             ax.legend()