From 14a9fdb95707a55c9a519dfca13acd7ce62c59bb Mon Sep 17 00:00:00 2001
From: Eric Teunis de Boone <ericteunis@deboone.nl>
Date: Mon, 12 Dec 2022 19:52:34 +0100
Subject: [PATCH] ZH: show phase due to t_trace[0] from measurement

This is not required to be known when finding the phase
as lib.find_beacon_in_traces already accounts for it.
---
 .../ba_measure_beacon_phase.py                | 40 ++++++++++++++++---
 .../lib/tests/test_beacon_fourier.py          | 16 +++++++-
 2 files changed, 48 insertions(+), 8 deletions(-)

diff --git a/simulations/airshower_beacon_simulation/ba_measure_beacon_phase.py b/simulations/airshower_beacon_simulation/ba_measure_beacon_phase.py
index c56e4a8..ad16f2e 100755
--- a/simulations/airshower_beacon_simulation/ba_measure_beacon_phase.py
+++ b/simulations/airshower_beacon_simulation/ba_measure_beacon_phase.py
@@ -92,7 +92,7 @@ if __name__ == "__main__":
                 # TODO: refine masking
                 # use beacon but remove where E_AxB-Beacon != 0
                 if True:
-                    if True:
+                    if not True:
                         t_mask = np.isclose(h5ant['E_AxB'][1], h5ant['traces'][4], rtol=1e-3, atol=1e-3)
                     else:
                         t_mask = np.ones(len(t_trace), dtype=bool)
@@ -156,24 +156,52 @@ if __name__ == "__main__":
             amplitude = amps[idx]
             orientation = orients[idx]
 
+            # Correct for phase by t_trace[0]
+            corr_phase = lib.phase_mod(2*np.pi*f_beacon*t_trace[0])
+            if False:
+                # Subtract phase due to not starting at t=0
+                # This is already done in beacon_find_traces
+                phase = lib.phase_mod(phase + corr_phase)
+
             # for reporting using plots
             found_data[i] = frequency, phase, amplitude
 
-            if (show_plots or fig_dir) and (i == 60 or i == 72):
+            if (show_plots or fig_dir) and (i == 0 or i == 72 or i == 70):
+                p2t = lambda phase: phase/(2*np.pi*f_beacon)
+
                 fig, ax = plt.subplots()
                 ax.set_title(f"Beacon at antenna {h5ant.attrs['name']}\nF:{frequency:.2e}, P:{phase:.4f}, A:{amplitude:.1e}")
                 ax.set_xlabel("t [ns]")
                 ax.set_ylabel("Amplitude")
 
-                for i, trace in enumerate(test_traces):
-                    ax.plot(t_trace, test_traces[i], marker='.', label='trace '+orients[i])
+                if True:
+                    # let the trace start at t=0
+                    t_0 = min(t_trace)
+                    extra_phase = corr_phase
+                else:
+                    t_0 = 0
+                    extra_phase = -1*corr_phase
 
-                myt = np.linspace(min(t_trace), max(t_trace), 10*len(t_trace))
-                ax.plot(myt, lib.sine_beacon(frequency, myt, amplitude=amplitude, phase=phase), ls='dashed', label='simulated beacon')
+                for j, trace in enumerate(test_traces):
+                    ax.plot(t_trace - t_0, test_traces[j], marker='.', label='trace '+orients[j])
+
+                myt = np.linspace(min(t_trace), max(t_trace), 10*len(t_trace)) - t_0
+                ax.plot(myt, lib.sine_beacon(frequency, myt, amplitude=amplitude, t0=0, phase=phase+extra_phase), ls='dotted', label='simulated beacon')
+
+                ax.axvline( p2t(lib.phase_mod(-1*(phase+extra_phase), low=0)), color='r', ls='dashed', label='$t_\\varphi$')
+
+                ax.axvline(0,color='grey',alpha=0.5)
+                ax.axhline(0,color='grey',alpha=0.5)
 
                 ax.legend()
 
                 if fig_dir:
+                    old_xlims = ax.get_xlim()
+                    ax.set_xlim(min(t_trace)-t_0-10,min(t_trace)-t_0+40)
+
+                    fig.savefig(path.join(fig_dir, __file__ + f".A{h5ant.attrs['name']}.zoomed.pdf"))
+
+                    ax.set_xlim(*old_xlims)
                     fig.savefig(path.join(fig_dir, __file__ + f".A{h5ant.attrs['name']}.pdf"))
 
             # save to file
diff --git a/simulations/airshower_beacon_simulation/lib/tests/test_beacon_fourier.py b/simulations/airshower_beacon_simulation/lib/tests/test_beacon_fourier.py
index c749dc3..ffaca39 100755
--- a/simulations/airshower_beacon_simulation/lib/tests/test_beacon_fourier.py
+++ b/simulations/airshower_beacon_simulation/lib/tests/test_beacon_fourier.py
@@ -1,5 +1,12 @@
 #!/usr/bin/env python3
 
+"""
+Test the functions in lib concerning
+beacon generation and phase measuring
+work correctly together.
+"""
+
+
 import matplotlib.pyplot as plt
 import numpy as np
 
@@ -7,20 +14,25 @@ import lib
 
 seed = 12345
 dt = 1 # ns
-t = np.arange(0, 10*int(1e3), dt)
 frequency = 45e-3 # GHz
-
 N = 5e2
+
+t = np.arange(0, 10*int(1e3), dt, dtype=float)
 rng = np.random.default_rng(seed)
 
 phase_res = np.zeros(int(N))
 
+# Vary both the base time and the phase
 for i in range(int(N)):
+    t_extra = (2*rng.uniform(size=1) - 1) *1e3
+    t += t_extra
+
     phase = lib.phase_mod(np.pi*(2*rng.uniform(size=1) -1)) # rad
     beacon = lib.sine_beacon(frequency, t, t0=0, phase=phase)
 
     measured = lib.find_beacon_in_traces([beacon], t, frequency, frequency_fit=False)
 
+    t -= t_extra
     phase_res[i] = lib.phase_mod(measured[1][0] - phase)
 
 fig, ax = plt.subplots()