From a4fa874b54754f214dec23113f4a89ac87378bc8 Mon Sep 17 00:00:00 2001
From: Eric Teunis de Boone <ericteunis@deboone.nl>
Date: Thu, 15 Dec 2022 13:35:33 +0100
Subject: [PATCH] ZH: test script for true_phase/calculated phase at tx

---
 .../lib/tests/test_beacon_fourier.py          | 11 +++-
 .../test_calculated_phase_measurement.py      | 63 +++++++++++++++++++
 2 files changed, 73 insertions(+), 1 deletion(-)
 create mode 100755 simulations/airshower_beacon_simulation/lib/tests/test_calculated_phase_measurement.py

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 ffaca39..4e86c34 100755
--- a/simulations/airshower_beacon_simulation/lib/tests/test_beacon_fourier.py
+++ b/simulations/airshower_beacon_simulation/lib/tests/test_beacon_fourier.py
@@ -23,19 +23,28 @@ rng = np.random.default_rng(seed)
 phase_res = np.zeros(int(N))
 
 # Vary both the base time and the phase
+t_extra = 0
 for i in range(int(N)):
+
+    # Change timebase
+    t -= t_extra
     t_extra = (2*rng.uniform(size=1) - 1) *1e3
     t += t_extra
 
+    # Randomly phased beacon
     phase = lib.phase_mod(np.pi*(2*rng.uniform(size=1) -1)) # rad
     beacon = lib.sine_beacon(frequency, t, t0=0, phase=phase)
 
+    if True: # blank part of the beacon
+        blank_low, blank_high = 2*int(1e3), 4*int(1e3)
+        beacon[blank_low:blank_high] = 0
+
     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()
+ax.set_title("Sine beacon phase determination\nwith random time shifts")
 ax.set_xlabel("$\\varphi_{meas} - \\varphi_{true}$ [rad]")
 ax.set_ylabel("#")
 ax.hist(phase_res, bins='sqrt')
diff --git a/simulations/airshower_beacon_simulation/lib/tests/test_calculated_phase_measurement.py b/simulations/airshower_beacon_simulation/lib/tests/test_calculated_phase_measurement.py
new file mode 100755
index 0000000..6a19b57
--- /dev/null
+++ b/simulations/airshower_beacon_simulation/lib/tests/test_calculated_phase_measurement.py
@@ -0,0 +1,63 @@
+#!/usr/bin/env python3
+
+"""
+Test the functions in lib concerning
+beacon generation and removing the geometrical phase
+work correctly together.
+"""
+
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from earsim import Antenna
+import lib
+
+seed = 12345
+dt = 1 # ns
+frequency = 45e-3 # GHz
+N = 5e2
+c_light = 3e8*1e-9
+
+t = np.arange(0, 10*int(1e3), dt, dtype=float)
+rng = np.random.default_rng(seed)
+phase_in = lib.phase_mod(np.pi*(2*rng.uniform(size=1) -1)) # rad
+
+tx = Antenna(x=0,y=0,z=0,name='tx')
+rx = Antenna(x=30,y=40,z=120,name='rx')
+
+# Vary both the base time and the phase
+t_extra = 0
+phase_res = np.zeros(int(N))
+for i in range(int(N)):
+    # Change timebase
+    t -= t_extra
+    t_extra = (2*rng.uniform(size=1) - 1) *1e3
+    t += t_extra
+
+    # Randomise Antenna Location
+    if True:
+        rx.x, rx.y, rx.z = (2*rng.uniform(size=3) -1) * 1e4
+
+    # Randomly phased beacon
+    # at Antenna
+    phase = lib.phase_mod(np.pi*(2*rng.uniform(size=1) -1)) # rad
+    beacon = lib.beacon_from(tx, rx, frequency, t, radiate_rsq=False, phase=phase, c_light=c_light)
+
+    if True: # blank part of the beacon
+        blank_low, blank_high = 2*int(1e3), 4*int(1e3)
+        beacon[blank_low:blank_high] = 0
+
+    measured = lib.find_beacon_in_traces([beacon], t, frequency, frequency_fit=False)
+
+    calculated_phase = lib.remove_antenna_geometry_phase(tx, rx, frequency, measured[1][0], c_light=c_light)
+
+    phase_res[i] = lib.phase_mod(calculated_phase - phase)
+
+fig, ax = plt.subplots()
+ax.set_title("Measured phase at Antenna - geometrical phase")
+ax.set_xlabel("$\\varphi_{meas} - \\varphi_{true}$ [rad]")
+ax.set_ylabel("#")
+ax.hist(phase_res, bins='sqrt')
+plt.show()
+