ZH: script saving k_period to file

simulations/airshower_beacon_simulation/bc_beacon_periods.py

@@ -15,6 +15,9 @@ if __name__ == "__main__":
 
     fname = "ZH_airshower/mysim.sry"
 
+    show_plots = True
+    ref_ant_id = None # leave None for all baselines
+
     ####
     fname_dir = path.dirname(fname)
     antennas_fname = path.join(fname_dir, beacon.antennas_fname)
@@ -23,11 +26,12 @@ if __name__ == "__main__":
     f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
 
     # run over all baselines
-    if True:
+    if ref_ant_id is None:
+        print("Doing all baselines")
         baselines = list(combinations(antennas,2))
     # use ref_ant
     else:
-        ref_ant = antennas[0]
+        ref_ant = antennas[ref_ant_idx]
         baselines = list(zip_longest([], antennas, fillvalue=ref_ant))
 
     freq_names = antennas[0].beacon_info.keys()
@@ -36,22 +40,21 @@ if __name__ == "__main__":
 
     freq_name = next(iter(freq_names))
 
-
     # Determine integer multiple of periods to shift
-    integer_periods = np.empty( (len(baselines), 3) )
+    # and True phase differences
+    time_diffs = np.empty( (len(baselines), 3) )
     for i, base in enumerate(baselines):
         # 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
+        # how many samples do we need to shift
         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)
-
+        print('A1:', base[0].name, 'A2:', base[1].name, "K:", best_k, '= [ns]', delta_t_coherence)
 
         # get the amount of periods to move
         f_beacon = base[0].beacon_info[freq_name]['freq']
@@ -61,14 +64,24 @@ if __name__ == "__main__":
         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]
+        # 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])
+        except IndexError:
+            # freq_name not in beacon_info
+            # or true_phase not determined yet
+            true_phases_diff = np.nan
 
-        if i in [ 98, 99 ]:
+        # save k_period with antenna names
+        time_diffs[i] = [true_phases_diff, k_period, f_beacon]
+
+        # Plotting for one or two iterations
+        if show_plots and i in [ 0, 1 ]:
             print('i',i,'k[T]',k_period, 'rest[ns]',rest, 'T[ns]',1/f_beacon)
 
             # Show correlation maxima plot
-            if True:
+            if not True:
                 fig, ax = plt.subplots()
                 ax.set_title(f"Correlation Maxima {i}")
                 ax.set_xlabel("k")
@@ -80,18 +93,11 @@ if __name__ == "__main__":
             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
+            if true_phases_diff is np.nan:
+                true_phases_diff = 0
+            delta_t_beacon = true_phases_diff/(2*np.pi*f_beacon)
 
-            print("t0[ns]", delta_t_antennas, "t_beacon[ns]", delta_t_beacon, "phase", delta_true_phases)
+            print("t0[ns]", delta_t_antennas, "t_beacon[ns]", delta_t_beacon, "phase", true_phases_diff)
             fig, ax = plt.subplots()
             ax.set_xlabel('t')
             ax.plot(base[0].t, traces[0], label=f'Reference {base[0].name}', alpha=0.5)
@@ -108,13 +114,8 @@ if __name__ == "__main__":
             ax.legend()
 
     # Save integer periods to antennas
-    with h5py.File(antennas_fname, 'a') as fp:
-        group_name = 'beacon_ks'
-        if group_name in fp:
-            del fp[group_name]
+    beacon.write_baseline_time_diffs_hdf5(antennas_fname, baselines, time_diffs[:,0], time_diffs[:,1], time_diffs[:,2])
 
-        fp.create_dataset(group_name, data=integer_periods)
-
-    plt.show()
     # Report back to CLI
-    print("Period Multiples resolved in", antennas_fname)
+    print("Period Multiples resolved and written to ", antennas_fname)
+    plt.show()