ZH: resolve integer multiples for each combination of antennas

simulations/airshower_beacon_simulation/bb_beacon_periods.py

@@ -2,6 +2,7 @@
 # vim: fdm=indent ts=4
 
 import h5py
+from itertools import combinations, zip_longest
 import matplotlib.pyplot as plt
 import numpy as np
 
@@ -40,9 +41,8 @@ if __name__ == "__main__":
 
             ant.attrs['beacon_phase_true'] = true_phases[i]
 
-    # True phase (without geometry determined)
+    # Plot True Phases
-
+    if True:
-    if True: # show antenna phases
         fig, ax = plt.subplots()
         spatial_unit=None
         fig.suptitle('f= {:2.0f}MHz'.format(f_beacon*1e3))
@@ -59,6 +59,83 @@ if __name__ == "__main__":
         sc = ax.scatter(*antenna_locs, c=true_phases, **scatter_kwargs)
         fig.colorbar(sc, ax=ax, label=color_label)
 
+    # run over all baselines
+    if True:
+        baselines = list(combinations(antennas,2))
+    # use ref_ant
+    else:
+        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]
+
+        # read integer period from file
+        if True and 'beacon_ks' in fp:
+            integer_periods = np.array(fp['beacon_ks'])
+
+
+    # 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 ]:
+                continue
+
+            fig, ax = plt.subplots()
+            ax.set_xlabel("k")
+            ax.set_ylabel("Maximum correlation")
+            ax.plot(ks, maxima)
+            ax.plot(ks[max_idx], maxima[max_idx], marker='X')
+
+            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.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]
+
+        fp.create_dataset(group_name, data=integer_periods)
+
     plt.show()
     # Report back to CLI
-    #print("Period Multiples resolved in", antenna_fname)
+    print("Period Multiples resolved in", antennas_fname)

simulations/airshower_beacon_simulation/lib.py

@@ -34,6 +34,12 @@ def distance(x1, x2):
 
     return np.sqrt( np.sum( (x1-x2)**2 ) )
 
+def geometry_time(dist, x2=None, c_light=3e8):
+    if x2 is not None:
+        dist = distance(dist, x2)
+
+    return dist/c_light
+
 def beacon_from(tx, rx, f, t=0, t0=0, c_light=3e8, radiate_rsq=True, amplitude=1,**kwargs):
     dist = distance(tx,rx)
     t0 = t0 + dist/c_light
@@ -81,6 +87,8 @@ def direct_fourier_transform(freqs, time, samplesets_iterable):
     return np.dot(c_k, samplesets_iterable), np.dot(s_k, samplesets_iterable)
 
 def phase_field_from_tx(x, y, tx, f_beacon, c_light=3e8, t0=0, wrap_phase=True, return_meshgrid=True):
+    """
+    """
 
     assert type(tx) in [Antenna]
 
@@ -191,3 +199,19 @@ def find_beacon_in_traces(
         amplitudes[i] = 2/n_samples * (real**2 + imag**2)**0.5
 
     return frequencies, phases, amplitudes
+
+def coherence_sum_maxima(ref_x, y, k_step=1):
+    """
+    Use the maximum of a coherent sum to determine
+    the best number of samples to move
+    """
+    max_k = int( len(ref_x) )
+    ks = np.arange(0, max_k, step=k_step)
+
+    maxima = np.empty(len(ks))
+
+    for i,k in enumerate(ks, 0):
+        augmented_y = np.roll(y, k)
+        maxima[i] = max(ref_x + augmented_y)
+
+    return ks, maxima