From 4c834ad8e7ceb98b01c2390216c1ffebfa4e3295 Mon Sep 17 00:00:00 2001
From: Eric Teunis de Boone <ericteunis@deboone.nl>
Date: Mon, 12 Dec 2022 19:01:02 +0100
Subject: [PATCH] ZH: use correct DTFT convention

Only affects phase determination.
Introduces a minus sign for the s_k terms and changes arctan2 parameters
---
 .../airshower_beacon_simulation/lib/lib.py    |  8 ++---
 .../airshower_beacon_simulation/lib/tests/lib |  1 +
 .../lib/tests/test_beacon_fourier.py          | 30 +++++++++++++++++++
 3 files changed, 35 insertions(+), 4 deletions(-)
 create mode 120000 simulations/airshower_beacon_simulation/lib/tests/lib
 create mode 100755 simulations/airshower_beacon_simulation/lib/tests/test_beacon_fourier.py

diff --git a/simulations/airshower_beacon_simulation/lib/lib.py b/simulations/airshower_beacon_simulation/lib/lib.py
index 3a3f2c1..70a67ce 100644
--- a/simulations/airshower_beacon_simulation/lib/lib.py
+++ b/simulations/airshower_beacon_simulation/lib/lib.py
@@ -12,7 +12,7 @@ def sine_beacon(f, t, t0=0, amplitude=1, baseline=0, phase=0):
     """
     Return a sine appropriate as a beacon
     """
-    return amplitude * np.sin(2*np.pi*f*(t+t0) + phase) + baseline
+    return amplitude * np.cos(2*np.pi*f*(t+t0) + phase) + baseline
 
 def phase_mod(phase, low=np.pi):
     """
@@ -62,8 +62,8 @@ def ft_corr_vectors(freqs, time):
     """
     freqtime = np.outer(freqs, time)
 
-    c_k = np.cos(2*np.pi*freqtime)
-    s_k = np.sin(2*np.pi*freqtime)
+    c_k =    np.cos(2*np.pi*freqtime)
+    s_k = -1*np.sin(2*np.pi*freqtime)
 
     return c_k, s_k
 
@@ -197,7 +197,7 @@ def find_beacon_in_traces(
 
         real, imag = direct_fourier_transform(freq, t_trace, traces[i])
 
-        phases[i] = np.arctan2(real, imag)
+        phases[i] = np.arctan2(imag, real)
         amplitudes[i] = 2/n_samples * (real**2 + imag**2)**0.5
 
     return frequencies, phases, amplitudes
diff --git a/simulations/airshower_beacon_simulation/lib/tests/lib b/simulations/airshower_beacon_simulation/lib/tests/lib
new file mode 120000
index 0000000..b870225
--- /dev/null
+++ b/simulations/airshower_beacon_simulation/lib/tests/lib
@@ -0,0 +1 @@
+../
\ No newline at end of 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
new file mode 100755
index 0000000..c749dc3
--- /dev/null
+++ b/simulations/airshower_beacon_simulation/lib/tests/test_beacon_fourier.py
@@ -0,0 +1,30 @@
+#!/usr/bin/env python3
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import lib
+
+seed = 12345
+dt = 1 # ns
+t = np.arange(0, 10*int(1e3), dt)
+frequency = 45e-3 # GHz
+
+N = 5e2
+rng = np.random.default_rng(seed)
+
+phase_res = np.zeros(int(N))
+
+for i in range(int(N)):
+    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)
+
+    phase_res[i] = lib.phase_mod(measured[1][0] - phase)
+
+fig, ax = plt.subplots()
+ax.set_xlabel("$\\varphi_{meas} - \\varphi_{true}$ [rad]")
+ax.set_ylabel("#")
+ax.hist(phase_res, bins='sqrt')
+plt.show()