ZH: script to determine true phases (and periods)

Periods still need to be done

simulations/airshower_beacon_simulation/bb_beacon_periods.py

@@ -0,0 +1,64 @@
+#!/usr/bin/env python3
+# vim: fdm=indent ts=4
+
+import h5py
+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"
+
+    ####
+    fname_dir = path.dirname(fname)
+    antennas_fname = path.join(fname_dir, beacon.antennas_fname)
+
+    if not path.isfile(antennas_fname):
+        print("Antenna file cannot be found, did you try generating a beacon?")
+        sys.exit(1)
+
+    # Read in antennas from file
+    f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
+
+    if True and 'beacon_phase_true' in antennas[0].attrs:
+        true_phases = np.array([a.attrs['beacon_phase_true'] for a in antennas])
+    else:
+        true_phases = np.empty( (len(antennas)) )
+
+        for i, ant in enumerate(antennas):
+            measured_phase = ant.attrs['beacon_phase_measured']
+
+            geom_time = lib.geometry_time(tx, ant, c_light=3e8*1e-9)
+            geom_phase = geom_time * 2*np.pi*f_beacon
+
+            true_phases[i] = lib.phase_mod(measured_phase) - lib.phase_mod(geom_phase)
+
+            ant.attrs['beacon_phase_true'] = true_phases[i]
+
+    # True phase (without geometry determined)
+
+    if True: # show antenna phases
+        fig, ax = plt.subplots()
+        spatial_unit=None
+        fig.suptitle('f= {:2.0f}MHz'.format(f_beacon*1e3))
+
+        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'] = 'Spectral_r'
+        scatter_kwargs['vmin'] = -np.pi
+        scatter_kwargs['vmax'] = +np.pi
+        color_label='$\\varphi$'
+
+        sc = ax.scatter(*antenna_locs, c=true_phases, **scatter_kwargs)
+        fig.colorbar(sc, ax=ax, label=color_label)
+
+    plt.show()
+    # Report back to CLI
+    #print("Period Multiples resolved in", antenna_fname)