ZH: beacon phase finding working

simulations/airshower_beacon_simulation/ba_beacon_phases.py

@@ -12,11 +12,13 @@ if __name__ == "__main__":
     from os import path
     import sys
 
-    
+
     f_beacon_band = (49e-3,55e-3) #GHz
-    allow_frequency_fitting = True
+    allow_frequency_fitting = False
     read_frequency_from_file = True
 
+    show_plots = True
+
     fname = "ZH_airshower/mysim.sry"
 
     ####
@@ -65,40 +67,72 @@ if __name__ == "__main__":
                 sys.exit(1)
 
             traces = ant_group['traces']
+            t_trace = traces[0]
 
-            freqs, phases, amps = find_beacon_in_traces( 
-                    traces[1:-1], traces[0],
-                    f_beacon_estimate=f_beacon,
-                    frequency_fit=allow_frequency_fitting,
-                    f_beacon_estimate_band=f_beacon_estimate_band
-                    )
+            if not True:
+                # only take the Beacon trace
+                test_traces = traces[4]
+                orients = ['B']
+            else:
+                test_traces = traces[1:]
+                orients = ['Ex', 'Ey', 'Ez', 'B']
 
-            # only take Ex for now
-            frequency = freqs[-1]
-            phase = phases[-1]
-            amplitude = amps[-1]
+            if True:
+                freqs, phases, amps = lib.find_beacon_in_traces(
+                        test_traces, t_trace,
+                        f_beacon_estimate=f_beacon,
+                        frequency_fit=allow_frequency_fitting,
+                        f_beacon_estimate_band=f_beacon_estimate_band
+                        )
+            else:
+                # Testing
+                freqs = [f_beacon]
+                t0 = ant_group.attrs['t0']
 
-            print(frequency, phase, amplitude)
+                phases = [ 2*np.pi*t0*f_beacon ]
+                amps = [ 3e-7 ]
+
+            # choose highest amp
+            #idx = np.argmax(amps, axis=1)
+            idx = 0
+
+            frequency = freqs[idx]
+            phase = phases[idx]
+            amplitude = amps[idx]
+            orientation = orients[idx]
+
+            if show_plots and (i == 60 or i == 72):
+                fig, ax = plt.subplots()
+                trace_amp = max(traces[-1]) - min(traces[-1])
+
+                myt = np.linspace(min(traces[0]), max(traces[0]), 10*len(traces[0]))
+                ax.plot(t_trace, traces[-1], marker='.', label='trace')
+                ax.plot(myt, lib.sine_beacon(frequency, myt, amplitude=amplitude, phase=phase), ls='dashed', label='simulated')
+                ax.set_title(f"Beacon at antenna {ant_group}\nF:{frequency}, P:{phase}, A:{amplitude}")
+                ax.legend()
 
             ant_group.attrs['beacon_freq'] = frequency
             ant_group.attrs['beacon_phase'] = phase
             ant_group.attrs['beacon_amplitude'] = amplitude
-            ant_group.attrs['beacon_orientation'] = 'Ex'
+            ant_group.attrs['beacon_orientation'] = orientation
 
             found_data[i] = frequency, phase, amplitude
 
+    print("Beacon Phases, Amplitudes and Frequencies written to", antennas_fname)
+
     # show histogram of found frequencies
-    if True:
+    if show_plots:
         if True or allow_frequency_fitting:
             fig, ax = plt.subplots()
             ax.set_xlabel("Frequency")
             ax.set_ylabel("Counts")
-            ax.hist(found_data[:,0], bins='auto', density=False)
+            ax.axvline(f_beacon, ls='dashed', color='g')
+            ax.hist(found_data[:,0], bins='sqrt', density=False)
 
         if True:
             fig, ax = plt.subplots()
             ax.set_xlabel("Amplitudes")
             ax.set_ylabel("Counts")
-            ax.hist(found_data[:,2], bins='auto', density=False)
+            ax.hist(found_data[:,2], bins='sqrt', density=False)
 
         plt.show()

simulations/airshower_beacon_simulation/lib.py

@@ -11,7 +11,7 @@ def sine_beacon(f, t, t0=0, amplitude=1, baseline=0, phase=0):
     """
     Return a sine appropriate as a beacon
     """
-    return amplitude * np.cos(2*np.pi*f*(t+t0) + phase) + baseline
+    return amplitude * np.sin(2*np.pi*f*(t+t0) + phase) + baseline
 
 
 def distance(x1, x2):
@@ -74,6 +74,9 @@ 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]
+
     xs, ys = np.meshgrid(x, y, sparse=True)
 
     x_distances = (tx.x - xs)**2
@@ -128,8 +131,8 @@ def find_beacon_in_traces(
             real, imag = ft_amp
             amps = 1/n_samples * ( real**2 + imag**2)**0.5
 
-            # find frequency peak and surrounding
-            # bins valid for parabola fitting
+            # find frequency peak and surrounding bins
+            # that are valid for the parabola fit
             max_amp_idx = np.argmax(amps)
             max_amp = amps[max_amp_idx]
 
@@ -170,9 +173,12 @@ def find_beacon_in_traces(
             freq = deriv.roots()[0]
 
             frequencies[i] = freq
-    else:
+
+    else: # no frequency finding
         frequencies[:] = f_beacon_estimate
 
+
+    n_samples = len(t_trace)
     # evaluate fourier transform at freq for each trace
     for i, freq in enumerate(frequencies):
         if freq is np.nan:
@@ -183,6 +189,6 @@ def find_beacon_in_traces(
         real, imag = direct_fourier_transform(freq, t_trace, traces[i])
 
         phases[i] = np.arctan2(real, imag)
-        amplitudes[i] = 2/len(t_trace) * (real**2 + imag**2)**0.5
+        amplitudes[i] = 2/n_samples * (real**2 + imag**2)**0.5
 
     return frequencies, phases, amplitudes