From c15f4e08af40342d7d1968698ef449530ca3cef6 Mon Sep 17 00:00:00 2001
From: Eric Teunis de Boone <ericteunis@deboone.nl>
Date: Thu, 24 Nov 2022 17:54:48 +0100
Subject: [PATCH] ZH: simple time diff inspection script

---
 .../bd_time_diffs.py                          | 58 +++++++++++++++++++
 1 file changed, 58 insertions(+)
 create mode 100755 simulations/airshower_beacon_simulation/bd_time_diffs.py

diff --git a/simulations/airshower_beacon_simulation/bd_time_diffs.py b/simulations/airshower_beacon_simulation/bd_time_diffs.py
new file mode 100755
index 0000000..ce7e1c9
--- /dev/null
+++ b/simulations/airshower_beacon_simulation/bd_time_diffs.py
@@ -0,0 +1,58 @@
+#!/usr/bin/env python3
+# vim: fdm=indent ts=4
+
+import h5py
+from itertools import combinations, zip_longest
+import matplotlib.pyplot as plt
+import numpy as np
+
+import aa_generate_beacon as beacon
+import lib
+
+
+if __name__ == "__main__":
+    from os import path
+    import sys
+
+    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)
+
+    basenames, time_diffs, f_beacon, true_phase_diffs, k_periods = beacon.read_baseline_time_diffs_hdf5(antennas_fname)
+
+    f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
+
+    antenna_time_shifts = { a.name: a.attrs['clock_offset'] for a in antennas }
+
+
+    my_time_shifts = []
+    for i,b in enumerate(basenames):
+        actual_time_shift = antenna_time_shifts[b[0]] - antenna_time_shifts[b[1]]
+        my_time_shifts.append(actual_time_shift)
+
+        print(
+                f'B({b[0]}, {b[1]}):',
+                time_diffs[i],
+                k_periods[i],
+                'A', actual_time_shift
+            )
+
+    # Make a plot
+    N = len(basenames)
+    fig, ax = plt.subplots()
+    ax.set_xlabel("Baseline combo")
+    ax.set_ylabel("[ns]")
+    ax.plot(np.arange(N), my_time_shifts, marker='+', label='actual time shifts')
+    ax.plot(np.arange(N), time_diffs, marker='x', label='calculated')
+
+    ax.legend()
+
+    plt.show()
+
+
+