Thesis: final PDF

documents/thesis/chapters/colophon.tex

@@ -8,4 +8,19 @@
 }
 
 \begin{document}
+\vspace*{\fill}
+A digital copy of this document can be found at\\
+\url{https://etdeboone.nl/masters-thesis/thesis.pdf}\\[0.2cm]
+
+
+The \LaTeX\ source, including figures, can be found at\\
+\url{https://gitlab.science.ru.nl/mthesis-edeboone/m.internship-documentation}\\
+--- or ---\\
+\url{https://etdeboone.nl/masters-thesis/documentation}\\[0.2cm].
+
+Code for generating figures, as well as the beacon synchronising pipeline can be found at\\
+\url{https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction/}\\
+--- or ---\\
+\url{https://etdeboone.nl/masters-thesis/airshower_beacon_simulation}.
+
 \end{document}

documents/thesis/chapters/conclusion.tex

@@ -40,7 +40,7 @@ Unfortunately, with the \gls{GNSS} timing accuracy estimated in the same order o
 
 % Combined sine beacon + air shower %<<<
 Recording an air shower, in addition to such a narrow-band beacon, might provide a method to determine the correct beacon period.
-Radio interferometeric analysis of the air shower depends on the coherence of the received signals.
+Radio interferometric analysis of the air shower depends on the coherence of the received signals.
 Any synchronicity problems in the radio antennas decrease the coherence and thus the power mapping used to derive properties of the air shower.
 With a limited set of periods to test, this power can be maximised while simultaneously inferring the correct beacon period.
 \\

documents/thesis/chapters/grand_characterisation.tex

@@ -39,8 +39,8 @@ This chapter starts an investigation into these systematic delays within \gls{GR
 % ADC
 At the base of every single antenna, a \gls{DU} is mounted.
 Its protective encasing has three inputs to which the different polarisations of the antenna are connected.
-These inputs are connected to their respective filterchains, leaving a fourth filterchain as spare.
-Each filterchain bandpasses the signal between $30\MHz$ and $200\MHz$.
+These inputs are connected to their respective filter chains, leaving a fourth filter chain as spare.
+Each filter chain band-passes the signal between $30\MHz$ and $200\MHz$.
 Finally, the signals are digitised by a four channel 14-bit \gls{ADC} sampling at $500\MHz$.
 %The input voltage ranges from $-900\mV$ to $+900\mV$.
 In our setup, the channels are read out after one of two internal ``monitoring'' triggers fire with the ten-second trigger (TD) linked to the 1~\acrlong{PPS} of the \gls{GNSS} chip and the other (MD) a variable randomising trigger.
@@ -63,7 +63,7 @@ In our setup, the channels are read out after one of two internal ``monitoring''
 
 % >>>
 %\section{Filterchain Relative Time Delays}% <<<
-Both the \gls{ADC} and the filterchains introduce systematic delays.
+Both the \gls{ADC} and the filter chains introduce systematic delays.
 Since each channel corresponds to a polarisation, it is important that relative systematic delays between the channels can be accounted for.
 \\
 
@@ -76,7 +76,7 @@ Since each channel corresponds to a polarisation, it is important that relative
 	}
 	\label{fig:channel-delay-setup}
 \end{figure}
-Figure~\ref{fig:channel-delay-setup} illustrates a setup to measure the relative time delays of the filterchain and \gls{ADC}.
+Figure~\ref{fig:channel-delay-setup} illustrates a setup to measure the relative time delays of the filter chain and \gls{ADC}.
 Two \gls{DU}-channels receive the same signal from a signal generator where one of the channels takes an extra time delay $\Delta t_\mathrm{cable}$ due to extra cable length.
 In this ``forward'' setup, both channels are read out at the same time, and a time delay is derived from the channels' traces.
 Afterwards, the cables are interchanged and a second (``backward'') time delay is measured.
@@ -154,7 +154,7 @@ Figures~\ref{fig:grand:phaseshift:measurements} and~\ref{fig:grand:phaseshift} s
 
 % Conclusion
 Figure~\ref{fig:channel-delays} shows the measured total time delays and the resulting signal chain time delays between both channels 1 and 2, and channels 2 and 4.
-Apart from two exceptional time delays upto $0.2\ns$, the signal chain time delays are in general below $0.05\ns$.
+Apart from two exceptional time delays up to $0.2\ns$, the signal chain time delays are in general below $0.05\ns$.
 \\
 Note that the reported signal chain time delays must be taken to be indications due to systematic behaviours (see below).
 \\
@@ -168,13 +168,13 @@ Nevertheless, the result at these frequencies must be interpreted with some caut
 
 % Discussion
 The time delays for both TD- and MD-triggered events in Figure~\ref{fig:grand:phaseshift:measurements} show a systematic behaviour of increasing total time delays for the forward setup.
-However, in the backward setup, this is not as noticable.
+However, in the backward setup, this is not as noticeable.
 \\
 This skewing of the channel time delays in one of the setups is also found at other frequencies (see Figures~\ref{fig:grand:phaseshift:ch1ch2} and~\ref{fig:grand:phaseshift:ch2ch4}), raising questions on the stability of the setup.
 Unfortunately, it is primarily visible in the larger datasets which correspond to measurements over larger timescales.
 As the number of these large datasets is limited, further investigation with the current datasets is prohibited.
 \\
-The skewing might also be an artifact of the short waveforms ($N\sim500\;\mathrm{samples}$) the data acquisition system was able to retrieve at the time of measurement.
+The skewing might also be an artefact of the short waveforms ($N\sim500\;\mathrm{samples}$) the data acquisition system was able to retrieve at the time of measurement.
 Since the data acquisition system is now able to retrieve the maximum size waveforms, this systematic behaviour can be investigated in a further experiment.
 \\
 

documents/thesis/chapters/grand_characterisation.tex.bak

@@ -0,0 +1,207 @@
+% vim: fdm=marker fmr=<<<,>>>
+\documentclass[../thesis.tex]{subfiles}
+
+\graphicspath{
+	{.}
+	{../../figures/}
+	{../../../figures/}
+}
+
+\begin{document}
+\chapter{GRAND signal chain characterisation}
+\label{sec:gnss_accuracy}
+
+% systematic delays important to obtain the best synchronisation
+The beacon synchronisation strategy hinges on the ability to measure the beacon signal with sufficient timing accuracy.
+In the previous chapters, the overall performance of this strategy has been explored by using simulated waveforms.
+\\
+% ADC and filtering setup most important component.
+As mentioned in Chapter~\ref{sec:waveform}, the measured waveforms of a true detector will be influenced by characteristics of the antenna, the filter and the \gls{ADC}.
+Especially the filter and \gls{ADC} are important components to be characterised to compensate for possible systematic (relative) delays.
+This chapter starts an investigation into these systematic delays within \gls{GRAND}'s \gls{DU} V2.0\cite{GRAND:DU2}.
+\\
+
+%\section{GRAND DU}% <<<
+%\begin{figure}
+%	\begin{subfigure}{0.47\textwidth}
+%		\includegraphics[width=\textwidth]{grand/DU_board_encased}
+%	\end{subfigure}
+%	\hfill
+%	\begin{subfigure}{0.47\textwidth}
+%		\includegraphics[width=\textwidth]{grand/DU_board_nocase}
+%	\end{subfigure}
+%	\caption{
+%		\gls{GRAND}'s \acrlong{DU} V2.0 inside (\textit{left}) and outside (\textit{right}) its protective encasing.
+%	}
+%	\label{fig:grand_du}
+%\end{figure}
+
+% ADC
+At the base of every single antenna, a \gls{DU} is mounted.
+Its protective encasing has three inputs to which the different polarisations of the antenna are connected.
+These inputs are connected to their respective filterchains, leaving a fourth filterchain as spare.
+Each filterchain bandpasses the signal between $30\MHz$ and $200\MHz$.
+Finally, the signals are digitised by a four channel 14-bit \gls{ADC} sampling at $500\MHz$.
+%The input voltage ranges from $-900\mV$ to $+900\mV$.
+In our setup, the channels are read out after one of two internal ``monitoring'' triggers fire with the ten-second trigger (TD) linked to the 1~\acrlong{PPS} of the \gls{GNSS} chip and the other (MD) a variable randomising trigger.
+\\
+
+% timestamp = GPS + local oscillator
+%The \gls{DU} timestamps an event using a combination of the 1\gls{PPS} of a Trimble ICM 360 \gls{GNSS} chip and counting the local oscillator running at $500\MHz$.
+%At trigger time, the counter value is stored to obtain a timing accuracy of roughly $2\ns$.
+%The counter is also used to correct for fluctuating intervals of the 1\gls{PPS} by storing and resetting it at each incoming 1\gls{PPS}.
+
+
+\begin{figure}% <<<<
+	\centering
+		\includegraphics[width=0.5\textwidth]{grand/DU/1697110935017.jpeg}
+		\caption{
+			\gls{GRAND}'s \acrlong{DU} V2.0 inside its protective encasing.
+		}
+		\label{fig:grand_du}
+\end{figure}% >>>>
+
+% >>>
+%\section{Filterchain Relative Time Delays}% <<<
+Both the \gls{ADC} and the filterchains introduce systematic delays.
+Since each channel corresponds to a polarisation, it is important that relative systematic delays between the channels can be accounted for.
+\\
+
+\begin{figure}
+	\centering
+	\includegraphics[width=0.4\textwidth]{grand/setup/channel-delay-setup.pdf}
+	\caption{
+		Relative time delay experiment, a signal generator sends the same signal to two channels of the \gls{DU}.
+		The extra time delay incurred by the loop in the upper cable can be ignored by interchanging the cabling and doing a second measurement.
+	}
+	\label{fig:channel-delay-setup}
+\end{figure}
+Figure~\ref{fig:channel-delay-setup} illustrates a setup to measure the relative time delays of the filterchain and \gls{ADC}.
+Two \gls{DU}-channels receive the same signal from a signal generator where one of the channels takes an extra time delay $\Delta t_\mathrm{cable}$ due to extra cable length.
+In this ``forward'' setup, both channels are read out at the same time, and a time delay is derived from the channels' traces.
+Afterwards, the cables are interchanged and a second (``backward'') time delay is measured.
+\\
+The sum of the ``forward'' and ``backward'' time delays gives twice the relative time delay $\Delta t$ without needing to measure the time delays due to the cable lengths $t_\mathrm{cable}$ separately since
+\begin{equation}\label{eq:forward_backward_cabling}
+	\phantom{.}
+	\Delta t
+	= (t_\mathrm{forward} + t_\mathrm{backward})/2
+	= ([\Delta t + t_\mathrm{cable}] + [\Delta t - t_\mathrm{cable}])/2
+	.
+\end{equation}\\
+
+% setup: signal
+We used a signal generator to emit a single sine wave at frequencies from $50\MHz$ to $200\MHz$ at $200\;\mathrm{mVpp}$.
+Note that we measured the phases to determine the time delays for each channel.
+In Figure~\ref{fig:grand:signal} the time delay between the channels is clearly visible in the measured waveforms as well as in the phase spectrum.
+\\
+\begin{figure}% <<< fig:grand:signal
+	\begin{subfigure}{0.47\textwidth}
+		\includegraphics[width=\textwidth]{grand/split-cable/waveform_eid1_ch1ch2.pdf}
+		\label{fig:split-cable:waveform}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}{0.47\textwidth}
+		\includegraphics[width=\textwidth]{grand/split-cable/waveform_eid1_ch1ch2_spectrum.pdf}
+		\label{fig:split-cable:waveform:spectra}
+	\end{subfigure}
+	\caption{
+		Waveforms of the sine wave measured in the ``forward'' setup and their spectra around the testing frequency of $50\MHz$..
+		The sine wave was emitted at $200\;\mathrm{mVpp}$.
+	}
+	\label{fig:grand:signal}
+\end{figure}% >>>
+
+% Frequencies above 50mhz not true measurement
+In our setup, the cable length difference was $3.17-2.01 = 1.06\metre$, resulting in an estimated cable time delay of roughly $5\ns$.
+At a frequency of $50\MHz$, the difference between the forward and backward phase differences is thus expected to be approximately half a cycle.
+Figures~\ref{fig:grand:phaseshift:measurements} and~\ref{fig:grand:phaseshift} show this is in accordance with the measured delays.
+\\
+
+\begin{figure}% <<< fig:grand:phaseshift:measurements
+	\centering
+	\begin{subfigure}{0.47\textwidth}
+		\includegraphics[width=\textwidth]{grand/split-cable/sine-sweep/ch2ch4fig9-measurements.forward.50.pdf}
+		\caption{}
+		\label{fig:grand:phaseshift:measurements:forward}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}{0.47\textwidth}
+		\includegraphics[width=\textwidth]{grand/split-cable/sine-sweep/ch2ch4fig9-measurements.backward.50.pdf}
+		\caption{}
+		\label{fig:grand:phaseshift:measurements:backward}
+	\end{subfigure}
+	\caption{
+		The measured phase differences between channels 2 and 4 at $50\MHz$ converted to a time delay for the \subref{fig:grand:phaseshift:measurements:forward}~forward and \subref{fig:grand:phaseshift:measurements:backward}~backward setups.
+		The dashed vertical lines indicate the mean time delay, the errorbar at the bottom indicates the standard deviation of the samples.
+		Crosses are TD-triggered events, circles are MD-triggered.
+		The measurements are time-ordered within their trigger type.
+	}
+	\label{fig:grand:phaseshift:measurements}
+\end{figure}
+
+\begin{figure}% <<< fig:grand:phaseshift
+	\centering
+	\includegraphics[width=0.47\textwidth]{grand/split-cable/sine-sweep/ch2ch4fig8-histogram.50.pdf}
+	\caption{
+		Histogram of the measured phase differences in Figure~\ref{fig:grand:phaseshift:measurements}.
+		The relative signal chain time delay for the portrayed means is $0.2\ns$.
+	}
+	\label{fig:grand:phaseshift}
+\end{figure}% >>>
+
+\clearpage
+
+% Conclusion
+Figure~\ref{fig:channel-delays} shows the measured total time delays and the resulting signal chain time delays between both channels 1 and 2, and channels 2 and 4.
+Apart from two exceptional time delays upto $0.2\ns$, the signal chain time delays are in general below $0.05\ns$.
+\\
+Note that the reported signal chain time delays must be taken to be indications due to systematic behaviours (see below).
+\\
+Still, even when taking $0.2\ns$ as the upper limit of any relative signal chain time delay, the electric field at the antenna are reconstructable to a sufficient accuracy to use either the pulsed or sine beacon methods (see Figures~\ref{fig:pulse:snr_time_resolution} and~\ref{fig:sine:snr_time_resolution} for reference) to synchronise an array to enable radio interferometry.
+\\
+
+Note that at higher frequencies the phase differences are phase-wrapped due to contention of the used period and the cable time delay.
+Because it is symmetric for both setups, this should not affect the measurement of the signal chain time delay at the considered frequencies.
+Nevertheless, the result at these frequencies must be interpreted with some caution.
+\\
+
+% Discussion
+The time delays for both TD- and MD-triggered events in Figure~\ref{fig:grand:phaseshift:measurements} show a systematic behaviour of increasing total time delays for the forward setup.
+However, in the backward setup, this is not as noticable.
+\\
+This skewing of the channel time delays in one of the setups is also found at other frequencies (see Figures~\ref{fig:grand:phaseshift:ch1ch2} and~\ref{fig:grand:phaseshift:ch2ch4}), raising questions on the stability of the setup.
+Unfortunately, it is primarily visible in the larger datasets which correspond to measurements over larger timescales.
+As the number of these large datasets is limited, further investigation with the current datasets is prohibited.
+\\
+The skewing might also be an artifact of the short waveforms ($N\sim500\;\mathrm{samples}$) the data acquisition system was able to retrieve at the time of measurement.
+Since the data acquisition system is now able to retrieve the maximum size waveforms, this systematic behaviour can be investigated in a further experiment.
+\\
+
+\begin{figure}% <<<<
+	\centering
+	\begin{subfigure}{0.48\textwidth}
+		\includegraphics[width=\textwidth]{grand/split-cable/sine-sweep/ch1ch2fig2-combi-time-delays.pdf}
+		\caption{
+			Channels 1,2
+		}
+		\label{fig:channel-delays:1,2}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}{0.48\textwidth}
+		\includegraphics[width=\textwidth]{grand/split-cable/sine-sweep/ch2ch4fig2-combi-time-delays.pdf}
+		\caption{
+			Channels 2,4
+		}
+		\label{fig:channel-delays:2,4}
+	\end{subfigure}
+	\caption{
+		Total (\textit{upper}) and signal chain (\textit{lower}) time delays between \subref{fig:channel-delays:1,2} channels 1 and 2, and  \subref{fig:channel-delays:2,4} 2 and 4.
+		The dark grey vertical lines in the upper panes indicate the maximum measurable time delays at each frequency.
+		Due to systematic effects in the measurements and a low number of samples at certain frequencies, the signal chain time delays depicted here must be taken as indicative.
+		See text for discussion.
+	}
+	\label{fig:channel-delays}
+\end{figure}% >>>>
+% >>>
+\end{document}

documents/thesis/chapters/introduction.tex

@@ -40,7 +40,7 @@ It is therefore only at the very highest energies that the direction of an initi
 
 % CR: galaxy / extra-galactic
 The same argument (but in reverse) can be used to explain the steeper slope from the ``knee'' ($10^{6}\GeV$) onwards in Figure~\ref{fig:cr_flux}.
-The acceleration of cosmic rays equally requires strong and sizable magnetic fields.
+The acceleration of cosmic rays equally requires strong and sizeable magnetic fields.
 Size constraints on the Milky~Way lead to a maximum energy for which a cosmic ray can still be contained in our galaxy.
 It is thus at these energies that we can distinguish between galactic and extra-galactic origins.
 \\
@@ -67,7 +67,7 @@ Figure~\ref{fig:airshower:depth} shows the number of particles as a function of
 The atmospheric depth at which this number of particles reaches its maximum is called $\Xmax$.
 \pagebreak
 
-In Figure~\ref{fig:airshower:depth}, $\Xmax$ is different for the airshowers generated by a photon, a proton or an iron nucleus.
+In Figure~\ref{fig:airshower:depth}, $\Xmax$ is different for the air showers generated by a photon, a proton or an iron nucleus.
 Typically, heavy nuclei have their first interaction higher up in the atmosphere than protons, with photons penetrating the atmosphere even further.
 Therefore, accurate measurements of $\Xmax$ allow to statistically discriminate between photons, protons and iron nuclei.
 \\

documents/thesis/chapters/introduction.tex.bak

@@ -0,0 +1,188 @@
+% vim: fdm=marker fmr=<<<,>>>
+\documentclass[../thesis.tex]{subfiles}
+
+\graphicspath{
+	{.}
+	{../../figures/}
+	{../../../figures/}
+}
+
+\begin{document}
+\chapter{An Introduction to Cosmic Rays and Extensive Air Showers}
+\label{sec:introduction}
+%\section{Cosmic Particles}%<<<<<<
+%<<<
+% Energy and flux
+The Earth is bombarded with a variety of energetic, extra-terrestrial particles.
+The energies of these particles extend over many orders of magnitude (see Figure~\ref{fig:cr_flux}).
+The flux of these particles decreases exponentially with increasing energy.
+For very high energies, above $10^{6}\GeV$, the flux approaches one particle per~square~meter per~year, further decreasing to a single particle per~square~kilometer per~year for Ultra High Energies (UHE) at $10^{10}\GeV$.
+\\
+
+\begin{figure}%<<< fig:cr_flux
+	\centering
+	\includegraphics[width=0.9\textwidth]{astroparticle/The_CR_spectrum_2023.pdf}
+	\caption{
+		From \protect \cite{The_CR_spectrum}.
+		The diffuse cosmic ray spectrum (upper line) as measured by various experiments.
+		The intensity and fluxes can generally be described by rapidly decreasing power laws.
+		The grey shading indicates the order of magnitude of the particle flux, such that from the ankle on wards ($E>10^9\GeV$) the flux reaches $1$~particle per~square~kilometer per~year.
+	}
+	\label{fig:cr_flux}
+\end{figure}%>>>
+
+% CR: magnetic field
+At these high energies, the incoming particles are primarily cosmic rays\footnote{These are therefore known as \glspl{UHECR}.}, atomic nuclei typically ranging from protons ($Z=1$) up to iron ($Z=26$).
+Because these are charged, the various magnetic fields they pass through will deflect and randomise their trajectories.
+Of course, this effect is dependent on the strength and size of the magnetic field and the speed of the particle.
+It is therefore only at the very highest energies that the direction of an initial particle might be used to (conservatively) infer the direction of its origin.
+\\
+
+% CR: galaxy / extra-galactic
+The same argument (but in reverse) can be used to explain the steeper slope from the ``knee'' ($10^{6}\GeV$) onwards in Figure~\ref{fig:cr_flux}.
+The acceleration of cosmic rays equally requires strong and sizable magnetic fields.
+Size constraints on the Milky~Way lead to a maximum energy for which a cosmic ray can still be contained in our galaxy.
+It is thus at these energies that we can distinguish between galactic and extra-galactic origins.
+\\
+
+% Photons and Neutrinos
+Other particles at these energies include photons and neutrinos, which are not charged.
+Therefore, these particle types do not suffer from magnetic deflections and have the potential to reveal their source regions.
+Unfortunately, aside from both being much less frequent, photons can be absorbed and created by multiple mechanisms, while neutrinos are notoriously hard to detect due to their weak interaction.
+%\Todo{
+%	$\gamma + \nu$ production by CR,
+%	source / targets
+%}
+\\
+
+%>>>
+%\subsection{Air Showers}%<<<
+When a cosmic ray with an energy above $10^{3}\GeV$ comes into contact with the atmosphere, secondary particles are generated, forming an \gls{EAS}.
+This air shower consists of a cascade of interactions producing more particles that subsequently undergo further interactions.
+Thus, the number of particles rapidly increases further down the air shower.
+This happens until the mean energy per particle is sufficiently lowered from whereon these particles are absorbed in the atmosphere.
+\\
+
+Figure~\ref{fig:airshower:depth} shows the number of particles as a function of atmospheric depth where $0\;\mathrm{g/cm^2}$ corresponds with the top of the atmosphere.
+The atmospheric depth at which this number of particles reaches its maximum is called $\Xmax$.
+\pagebreak
+
+In Figure~\ref{fig:airshower:depth}, $\Xmax$ is different for the airshowers generated by a photon, a proton or an iron nucleus.
+Typically, heavy nuclei have their first interaction higher up in the atmosphere than protons, with photons penetrating the atmosphere even further.
+Therefore, accurate measurements of $\Xmax$ allow to statistically discriminate between photons, protons and iron nuclei.
+\\
+
+\begin{figure}%<<< airshower:depth
+	\centering
+	\vspace*{-10mm}
+	\includegraphics[width=0.5\textwidth]{airshower/shower_development_depth_iron_proton_photon.pdf}
+	\caption{
+		From H. Schoorlemmer.
+		Shower development as a function of atmospheric depth for an energy of $10^{19}\eV$.
+		Typically, iron- and proton-induced air showers have a difference in $\langle \Xmax \rangle$ of $100\;\mathrm{g/cm^2}$~\cite{Deligny:2023yms}.
+		For air showers from photons this is even further down the atmosphere.
+		They are, however, much more rare than cosmic rays.
+	}
+	\label{fig:airshower:depth}
+	\vspace*{-5mm}
+\end{figure}%>>>
+
+The initial particle type also influences the particle content of an air shower.
+Depending on the available interaction channels, we distinguish three components in air showers: the hadronic, electromagnetic and muonic components.
+Each component shows particular development and can be related to different observables of the air shower.
+\\
+For example, detecting a large hadronic component means the initial particle has access to hadronic interactions (creating hadrons such as pions, kaons, etc.) which is a typical sign of a cosmic ray.
+In contrast, for an initial photon, which cannot interact hadronicly, the energy will be dumped into the electromagnetic part of the air shower, mainly producing electrons, positrons and photons.
+\\
+
+Finally, any charged pions created in the air shower will decay into muons while still in the atmosphere, thus comprising the muonic component.
+The lifetime, and ease of penetration of relativistic muons allow them to propagate to the Earth's surface, even if other particles have decayed or have been absorbed in the atmosphere.
+These are therefore prime candidates for air shower detectors on the Earth's surface.
+\\
+
+% Radio measurements
+Processes in an air showers also generate radiation that can be picked up as coherent radio signals.
+%% Geo Synchro
+Due to the magnetic field of the Earth, the electrons in the air shower generate radiation.
+Termed geomagnetic emission in Figure~\ref{fig:airshower:polarisation}, this has a polarisation that is dependent on the magnetic field vector ($\vec{B}$) and the air shower velocity ($\vec{v}$).
+\\
+%% Askaryan / Charge excess
+An additional mechanism emitting radiation was theorised by Askaryan\cite{Askaryan:1961pfb}.
+Due to the large inertia of the positively charged ions with respect to their light, negatively charged electrons, a negative charge excess is created.
+In turn, this generates radiation that is polarised radially towards the shower axis (see Figure~\ref{fig:airshower:polarisation}).
+\\
+
+%% Cherenkov ring
+Due to charged particles moving relativistically through the refractive atmosphere, the produced radiation is concentrated on a cone-like structure.
+On the surface, this creates a ring called the Cherenkov-ring.
+On this ring, a peculiar inversion happens in the time-domain of the air shower signals.
+Outside the ring, radiation from the top of the air shower arrives earlier than radiation from the end of the air shower, whereas this is reversed inside the ring.
+Consequently, the radiation received at the Cherenkov-ring is maximally coherent, being concentrated in a small time-window.
+It is therefore crucial for radio detection to obtain measurements in this region.
+\\
+
+\begin{figure}%<<< airshower:polarisation
+	\centering
+	\begin{subfigure}{0.48\textwidth}
+		\includegraphics[width=\textwidth]{airshower/airshower_radio_polarisation_geomagnetic.png}%
+		\caption{
+			Geomagnetic emission
+		}
+		\label{fig:airshower:polarisation:geomagnetic}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}{0.48\textwidth}
+		\includegraphics[width=\textwidth]{airshower/airshower_radio_polarisation_askaryan.png}%
+		\caption{
+			Askaryan or charge-excess emission
+		}
+		\label{fig:airshower:polarisation:askaryan}
+	\end{subfigure}
+	\caption{
+		From \protect \cite{Schoorlemmer:2012xpa, Huege:2017bqv}
+		The Radio Emission mechanisms and the resulting polarisations of the radio signal: \subref{fig:airshower:polarisation:geomagnetic} geomagnetic and \subref{fig:airshower:polarisation:askaryan} charge-excess.
+		See text for explanation.
+	}
+	\label{fig:airshower:polarisation}
+	\vspace{-2mm}
+\end{figure}%>>>>>>
+%>>>>>>
+
+%\subsection{Experiments}%<<<
+As mentioned, the flux at the very highest energy is in the order of one particle per square kilometer per century (see Figure~\ref{fig:cr_flux}).
+Observatories therefore have to span huge areas to gather decent statistics at these highest energies on a practical timescale.
+In recent and upcoming experiments, such as the~\gls{Auger}\cite{Deligny:2023yms} and the~\gls{GRAND}\cite{GRAND:2018iaj}, the approach is typically to instrument a large area with a (sparse) grid of detectors to detect the generated air shower.
+With distances up to $1.5\;\mathrm{km}$ (\gls{Auger}), the detectors therefore have to operate in a self-sufficient manner with only wireless communication channels and timing provided by \gls{GNSS}.
+\\
+
+In the last two decades, with the advent of advanced electronics, the detection using radio antennas has received significant attention.
+Analysing air showers using radio interferometry requires a time synchronisation of the detectors to an accuracy in the order of $1\ns$\cite{Schoorlemmer:2020low} (see Chapter~\ref{sec:interferometry} for further details).
+Unfortunately, this timing accuracy is not continuously achieved by \glspl{GNSS}, if at all.
+For example, in the~\gls{AERA}, this was found to range up to multiple tens of nanoseconds over the course of a single day\cite{PierreAuger:2015aqe}.
+\\
+
+\pagebreak[1]
+
+% Structure summary
+This thesis investigates a relatively straightforward method (and its limits) to improve the timing accuracy of air shower radio detectors
+by using an additional radio signal called a beacon.
+It is organised as follows.
+\\
+
+First, an introduction to radio interferometry is given in Chapter~\ref{sec:interferometry}.
+This will be used later on and gives an insight into the timing accuracy requirements.
+\\
+Chapter~\ref{sec:waveform} reviews some typical techniques to analyse waveforms and to obtain timing information from them.
+\\
+
+In Chapter~\ref{sec:disciplining}, the concept of a beacon transmitter is introduced to synchronise an array of radio antennas.
+It demonstrates the achievable timing accuracy for a sine and pulse beacon using the techniques described in the preceding chapter.
+\\
+
+
+A degeneracy in the synchronisation is encountered when the timing accuracy of the \gls{GNSS} is in the order of the periodicity of a continuous beacon.
+Chapter~\ref{sec:single_sine_sync} establishes a method using a single sine wave beacon while using the radio interferometric approach to observe an air shower and correct for this effect.
+\\
+
+Finally, Chapter~\ref{sec:gnss_accuracy} investigates some possible limitations of the current hardware of \gls{GRAND} and its ability to record and reconstruct a beacon signal.
+\end{document}

documents/thesis/chapters/radio_measurement.tex

@@ -31,10 +31,10 @@ The $n$-th sample in this waveform is then associated with a time $t[n] = t[0] +
 
 % Filtering before ADC
 The sampling is limited by the \gls{ADC}'s Nyquist frequency at half its sampling rate.
-In addition, various frequency-dependent backgrounds can be reduced by applying a bandpass filter before digitisation.
+In addition, various frequency-dependent backgrounds can be reduced by applying a band-pass filter before digitisation.
 For example, in \gls{AERA} and in AugerPrime's radio detector \cite{Huege:2023pfb}, the filter attenuates all of the signal except for the frequency interval between $30 \text{--} 80\MHz$.
 \\
-In addition to a bandpass filter, more complex filter setups are used to remove unwanted components or introduce attenuation at specific frequencies.
+In addition to a band-pass filter, more complex filter setups are used to remove unwanted components or introduce attenuation at specific frequencies.
 For example, in \gls{GRAND} \cite{GRAND:2018iaj}, the total frequency band ranges from $20\MHz$ to $200\MHz$.
 such that the FM broadcasting band ($87.5\MHz \text{--} 108\MHz$) falls within this range.
 Therefore, notch filters have been introduced to suppress signals in this band.
@@ -49,7 +49,7 @@ Thus to reconstruct properties of the electric field signal from the waveform, b
 Different methods are available for the analysis of the waveform, and the antenna and filter responses.
 A key aspect is determining the frequency-dependent amplitudes (and phases) in the measurements to characterise the responses and, more importantly, select signals from background.
 With \glspl{FT}, these frequency spectra can be produced.
-This technique is especially important for the sinewave beacon of Section~\ref{sec:beacon:sine}, as it forms the basis of the phase measurement.
+This technique is especially important for the sine wave beacon of Section~\ref{sec:beacon:sine}, as it forms the basis of the phase measurement.
 \\
 The detection and identification of more complex time-domain signals can be achieved using the cross correlation,
 which is the basis for the pulsed beacon method of Section~\ref{sec:beacon:pulse}.
@@ -156,7 +156,7 @@ Since a complex plane wave can be linearly decomposed as
 		,
 	\end{aligned}
 \end{equation*}
-the above transforms can be decomposed into explicit real and imaginary parts aswell,
+the above transforms can be decomposed into explicit real and imaginary parts as well,
 i.e.,~\eqref{eq:fourier:dtft} becomes
 \begin{equation}
 	\phantom{.}
@@ -207,7 +207,7 @@ By missing the correct frequency bin for the sine wave, it estimates both a too
 
 
 % % Static sin/cos terms if f_s, f and N static ..
-When calculating the \gls{DTFT} for multiple inputs which share both an equal number of samples $N$ and equal sampling frequencies $f_s$, the $\sin$ and $\cos$ terms in \eqref{eq:fourier:dtft_decomposed} are the same for a single frequency $f$ upto an overall phase which is dependent on $t[0]$.
+When calculating the \gls{DTFT} for multiple inputs which share both an equal number of samples $N$ and equal sampling frequencies $f_s$, the $\sin$ and $\cos$ terms in \eqref{eq:fourier:dtft_decomposed} are the same for a single frequency $f$ up to an overall phase which is dependent on $t[0]$.
 Therefore, at the cost of an increased memory allocation, these terms can be precomputed, reducing the number of real multiplications to $2N+1$.
 
 % .. relevance to hardware if static frequency

documents/thesis/chapters/radio_measurement.tex.bak

@@ -0,0 +1,275 @@
+% vim: fdm=marker fmr=<<<,>>>
+\documentclass[../thesis.tex]{subfiles}
+
+\graphicspath{
+	{.}
+	{../../figures/}
+	{../../../figures/}
+}
+
+\begin{document}
+\chapter{Waveform Analysis Techniques}
+\label{sec:waveform}
+%Electric fields,
+%Antenna Polarizations,
+%Frequency Bandwidth,
+
+Radio antennas are sensitive to changes in their surrounding electric fields.
+The polarisation of the electric field that a single antenna can record is dependent on the geometry of this antenna.
+Therefore, in experiments such as \gls{Auger} or \gls{GRAND}, multiple antennas are incorporated into a single unit to obtain complementary polarisation recordings.
+Additionally, the shape and size of antennas affect how well the antenna responds to certain frequency ranges, resulting in different designs meeting different criteria.
+\\
+
+%Waveform time series data
+%Sampling,
+%Waveform + Time vector,
+In each radio detector, the antenna presents its signals to an \gls{ADC} as fluctuating voltages.
+In turn, the \gls{ADC} records the analog signals with a specified samplerate $f_s$ resulting in a sequence of digitised voltages or waveform.
+The $n$-th sample in this waveform is then associated with a time $t[n] = t[0] + n/f_s = t[0] + n\cdot\Delta t$ after the initial sample at $t[0]$.
+%In reality, the sampling rate will vary by very small amounts resulting in timestamp inaccuracies called jitter.
+\\
+
+% Filtering before ADC
+The sampling is limited by the \gls{ADC}'s Nyquist frequency at half its sampling rate.
+In addition, various frequency-dependent backgrounds can be reduced by applying a bandpass filter before digitisation.
+For example, in \gls{AERA} and in AugerPrime's radio detector \cite{Huege:2023pfb}, the filter attenuates all of the signal except for the frequency interval between $30 \text{--} 80\MHz$.
+\\
+In addition to a bandpass filter, more complex filter setups are used to remove unwanted components or introduce attenuation at specific frequencies.
+For example, in \gls{GRAND} \cite{GRAND:2018iaj}, the total frequency band ranges from $20\MHz$ to $200\MHz$.
+such that the FM broadcasting band ($87.5\MHz \text{--} 108\MHz$) falls within this range.
+Therefore, notch filters have been introduced to suppress signals in this band.
+\\
+
+% Filter and Antenna response
+From the above it is clear that the digitised waveform is a measurement of the electric field that is sequentially convolved with the antenna's and filter's response.
+Thus to reconstruct properties of the electric field signal from the waveform, both responses must be known.
+\\
+
+% Analysis, properties, frequencies, pulse detection, shape matching,
+Different methods are available for the analysis of the waveform, and the antenna and filter responses.
+A key aspect is determining the frequency-dependent amplitudes (and phases) in the measurements to characterise the responses and, more importantly, select signals from background.
+With \glspl{FT}, these frequency spectra can be produced.
+This technique is especially important for the sinewave beacon of Section~\ref{sec:beacon:sine}, as it forms the basis of the phase measurement.
+\\
+The detection and identification of more complex time-domain signals can be achieved using the cross correlation,
+which is the basis for the pulsed beacon method of Section~\ref{sec:beacon:pulse}.
+\\
+
+\section{Fourier Transforms}% <<<<
+\label{sec:fourier}
+\glspl{FT} allow for a frequency-domain representation of a time-domain signal.
+In the case of radio antennas, it converts a time-ordered sequence of voltages into a set of complex amplitudes that depend on frequency.
+By evaluating the \gls{FT} at appropriate frequencies, the frequency spectrum of a waveform is calculated.
+This method then allows to modify a signal by operating on its frequency components, i.e.~removing a narrow frequency band contamination within the signal.
+\\
+
+% DTFT from CTFT
+The continuous \acrlong{FT} takes the form
+\begin{equation}
+	\label{eq:fourier}
+	\phantom{.}
+	X(f) = \int_\infty^\infty \dif{t}\, x(t)\, e^{-i 2 \pi f t}
+	.
+\end{equation}
+It decomposes the signal $x(t) \in \mathcal{R}$ into plane waves with complex-valued amplitude $X(f)$ at frequency $f$.
+\\
+From the complex amplitude $X(f)$, the phase $\pTrue(f)$ and amplitude $A(f)$ are calculated as
+\begin{equation*}
+	\begin{aligned}
+		\phantom{.}
+		\pTrue(f) = \arg\left( X(f) \right), && \text{and} && A(f) = 2 \left| X(f) \right|
+		.
+	\end{aligned}
+\end{equation*}
+Note the factor $2$ in this definition of the amplitude.
+It is introduced to compensate for expecting a real valued input signal $x(t) \in \mathcal{R}$ and mapping negative frequencies to their positive equivalents.
+\\
+
+When $x(t)$ is sampled at discrete times, the integral of \eqref{eq:fourier} is discretized in time to result in the \gls{DTFT}:
+\begin{equation}
+	%\tag{DTFT}
+	\label{eq:fourier:dtft}
+	X(f) = \sum_{n=0}^{N-1}   x(t[n])\, e^{ -i 2 \pi f t[n]}
+\end{equation}
+where $x(t)$ is sampled a finite number of times $N$ at times $t[n]$.
+Note that the amplitude $A(f)$ will now scale with the number of samples~$N$, and thus should be normalised to $A(f) = 2 \left| X(f) \right| / N$.
+\\
+
+Considering a finite sampling size $N$ and periodicity of the signal, the bounds of the integral in \eqref{eq:fourier} have collapsed to $t[0]$ up to $t_{N-1}$.
+It follows that the lowest resolvable frequency is $f_\mathrm{lower} = 1/T = 1/(t_{N-1} - t[0])$.
+\\
+Additionally, when the sampling of $x(t)$ is equally spaced, the $t[n]$ terms can be written as a sequence, $t[n] - t[0] = n \Delta t = n/f_s$, with $f_s$ the sampling frequency.
+Here the highest resolvable frequency is limited by the Nyquist~frequency.
+\\
+
+% DFT sampling of DTFT / efficient multifrequency FFT
+Implementing the above decomposition of $t[n]$, \eqref{eq:fourier:dtft} can be rewritten in terms of multiples of the sampling frequency $f = k f_s/N$, becoming the \gls{DFT}
+\begin{equation*}
+	\label{eq:fourier:dft}
+	\phantom{.}
+	X(k)
+	% = \sum_{n=0}^{N-1}   x(t[n])\, e^{ -i 2 \pi f (t[0] + n/f_s)}
+	% = \sum_{n=0}^{N-1}   x(t[n])\, e^{ -i 2 \pi f t[0]}\, e^{ -i 2 \pi f n/f_s}
+	% = e^{ -i 2 \pi f t[0]}\, \sum_{n=0}^{N-1}   x(t[n])\, e^{ -i 2 \pi k n}
+	= e^{ -i 2 \pi f t[0]} \sum_{n=0}^{N-1}   x(t[n])\, \cdot e^{ -i 2 \pi \frac{k n}{N} }
+	.
+\end{equation*}
+
+The direct computation of this transform takes $2N$ complex multiplications and $2(N-1)$ complex additions for a single frequency $k$.
+When computing this transform for all integer $0 \leq k < N$, this amounts to $\mathcal{O}(N^2)$ complex computations.
+\acrlong{FFT}s (\acrshort{FFT}s) are efficient algorithms that derive all $X( 0 \leq k < N)$ in $\mathcal{O}( N \log N)$ calculations.
+
+\begin{figure}
+	\begin{subfigure}{0.49\textwidth}
+		\includegraphics[width=\textwidth]{methods/fourier/waveforms.pdf}%
+		%\caption{}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}{0.49\textwidth}
+		\includegraphics[width=\textwidth]{methods/fourier/spectrum.pdf}%
+		\label{fig:fourier:dtft_dft}
+		%\caption{}
+	\end{subfigure}
+	\caption{
+		\textit{Left:} A waveform sampling a sine wave with white noise.
+		\textit{Right:}
+			The frequency spectrum of the waveform.
+			Comparison of the \gls{DTFT} and \gls{DFT} of the same waveform.
+			The \gls{DFT} can be interpreted as sampling the \gls{DTFT} at integer multiple of the waveform's sampling rate $f_s$.
+	}
+	\label{fig:fourier}
+\end{figure}
+
+
+% Linearity fourier for real/imag
+In the previous equations, the resultant quantity $X(f)$ is a complex amplitude.
+Since a complex plane wave can be linearly decomposed as
+\begin{equation*}
+	\phantom{,}
+	\label{eq:complex_wave_decomposition}
+	\begin{aligned}
+		e^{-i x}
+		&
+		= \cos(x) + i\sin(-x)
+		%\\ &
+		= \Re\left(e^{-i x}\right) + i \Im\left( e^{-i x} \right)
+		,
+	\end{aligned}
+\end{equation*}
+the above transforms can be decomposed into explicit real and imaginary parts aswell,
+i.e.,~\eqref{eq:fourier:dtft} becomes
+\begin{equation}
+	\phantom{.}
+	\label{eq:fourier:dtft_decomposed}
+	\begin{aligned}
+		X(f)
+		&
+		= X_R(f) + i X_I(f)
+		%\\ &
+		\equiv \Re(X(f)) + i \Im(X(f))
+		\\ &
+		=     \sum_{n=0}^{N-1} \, x(t[n]) \, \cos( 2\pi f t[n] )
+		  - i \sum_{n=0}^{N-1} \, x(t[n]) \, \sin( 2\pi f t[n] )
+		.
+	\end{aligned}
+\end{equation}
+
+% FT term to phase and magnitude
+The normalised amplitude at a given frequency $A(f)$ is calculated from \eqref{eq:fourier:dtft} as
+\begin{equation}
+	\label{eq:complex_magnitude}
+	\phantom{.}
+	A(f)
+		\equiv \frac{2 \left| X(f) \right| }{N}
+		= \frac{ 2 \sqrt{ X_R(f)^2 + X_I(f)^2 } }{N}
+	.
+\end{equation}
+Likewise, the complex phase at a given frequency $\pTrue(f)$ is obtained by
+\begin{equation}
+	\label{eq:complex_phase}
+	\phantom{.}
+	\pTrue(f)
+		\equiv \arg( X(f) )
+		= \arctantwo\left( X_I(f), X_R(f) \right)
+	.
+\end{equation}
+
+% Recover A\cos(2\pi t[n] f + \phi) using above definitions
+Applying \eqref{eq:fourier:dtft_decomposed} to a signal $x(t) = A\cos(2\pi t f + \pTrue)$ with the above definitions obtains
+an amplitude $A$ and phase $\pTrue$ at frequency $f$.
+When the minus sign in the exponent of \eqref{eq:fourier} is not taken into account, the calculated phase in \eqref{eq:complex_phase} will have an extra minus sign.
+\\
+
+Figure~\ref{fig:fourier} shows the frequency spectrum of a simulated waveform that is obtained using either a \gls{DFT} or a \gls{DTFT}.
+It shows that the \gls{DFT} evaluates the \gls{DTFT} only at certain frequencies.
+By missing the correct frequency bin for the sine wave, it estimates both a too low amplitude and the wrong phase for the input function.
+\\
+
+
+% % Static sin/cos terms if f_s, f and N static ..
+When calculating the \gls{DTFT} for multiple inputs which share both an equal number of samples $N$ and equal sampling frequencies $f_s$, the $\sin$ and $\cos$ terms in \eqref{eq:fourier:dtft_decomposed} are the same for a single frequency $f$ upto an overall phase which is dependent on $t[0]$.
+Therefore, at the cost of an increased memory allocation, these terms can be precomputed, reducing the number of real multiplications to $2N+1$.
+
+% .. relevance to hardware if static frequency
+Thus, for static frequencies in a continuous beacon, the coefficients for evaluating the \gls{DTFT} can be put into the hardware of the detectors,
+opening the way to efficiently measuring the amplitude and phase in realtime.
+
+
+% >>>>
+
+\section{Cross-Correlation}% <<<<
+\label{sec:correlation}
+The cross-correlation is a measure of how similar two waveforms $u(t)$ and $v(t)$ are.
+By introducing a time delay $\tau$ in one of the waveforms it turns into a function of this time delay,
+\begin{equation}
+	\label{eq:correlation_cont}
+	\phantom{,}
+	\Corr(\tau; u, v) = \int_{-\infty}^{\infty} \dif t \, u(t)\, v^*(t-\tau)
+	,
+\end{equation}
+where the integral reduces to a sum for a finite amount of samples in either $u(t)$ or $v(t)$.
+Still, $\tau$ remains a continuous variable.
+% Figure example of correlation and argmax
+Figure~\ref{fig:correlation} illustrates how the best time delay $\tau$ between two waveforms can thus be found by finding the maximum cross-correlation.
+\\
+
+% Discrete \tau because of sampling
+In reality, both waveforms have a finite size, also reducing the time delay $\tau$ resolution to the highest sampling rate of the two waveforms.
+When the sampling rates are equal, the time delay variable is effectively shifting one waveform by a number of samples.
+\\
+% Upsampling? No
+Techniques such as upsampling or interpolation can be used to effectively change the sampling rate of a waveform up to a certain degree.
+\\
+
+% Approaching analog \tau; or zero-stuffing
+Since zero-valued samples do not contribute to the integral of \eqref{eq:correlation_cont}, they can be freely added (or ignored) to a waveform when performing the calculations.
+This means two waveforms of different sampling rates can be correlated when the sampling rates are integer multiples of each other, simply by zero-stuffing the slowly sampled waveform.
+This allows to approximate an analog time delay between two waveforms when one waveform is sampled at a very high rate as compared to the other.
+
+\begin{figure}
+	\centering
+	\begin{subfigure}{0.48\textwidth}
+		\includegraphics[width=\textwidth]{methods/correlation/waveforms.pdf}
+		%\caption{
+		%	Two waveforms.
+		%}%
+		\label{subfig:correlation:waveforms}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}{0.48\textwidth}
+		\includegraphics[width=\textwidth]{methods/correlation/correlation.pdf}
+		%\caption{
+		%	The correlation of two Waveforms as a function of time.
+		%}%
+		\label{subfig:correlation}
+	\end{subfigure}%
+	\caption{
+		\textit{Left:} Two waveforms to be correlated with the second waveform delayed by $5$.
+		\textit{Right:} The correlation of both waveforms as a function of the time delay $\tau$.
+		Here the best time delay (red dashed line) is found at $5$, which would align the maximum amplitudes of both waveforms in the left pane.
+	}
+	\label{fig:correlation}
+\end{figure}
+
+% >>>
+\end{document}

documents/thesis/chapters/single_sine_interferometry.tex

@@ -122,13 +122,16 @@ This falls into the dynamic setup previously mentioned.
 The best period defects must thus be recovered from a single event.
 \\
 When doing the interferometric analysis for a sine beacon synchronised array, waveforms can only be delayed by an integer amount of periods, thereby giving discrete solutions to maximising the interferometric signal.
-\Todo{add size of shower at plane vs period defects in meters}
 
 \clearpage
 \section{Air Shower simulation}
 % simulation of proton E15 on 10x10 antenna
 To test the idea of combining a single sine beacon with an air shower, we simulated a set of recordings of a single air shower that also contains a beacon signal.
-\footnote{\url{https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction/-/tree/main/airshower_beacon_simulation}}
+\footnote{
+\url{https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction/-/tree/main/airshower_beacon_simulation}
+or
+\url{https://etdeboone.nl/masters-thesis/airshower_beacon_simulation}
+}
 \\
 The air shower signal was simulated by \acrlong{ZHAireS}\cite{Alvarez-Muniz:2010hbb} on a grid of 10x10 antennas with a spacing of $50\,\mathrm{meters}$.
 Each antenna recorded a waveform of $500$ samples with a samplerate of $1\GHz$ for each of the X,Y and Z polarisations.
@@ -137,7 +140,7 @@ The air shower itself was generated by a $10^{16}\eV$ proton coming in under an
 \\
 Figure~\ref{fig:single:proton_grid} shows the maximum electric field measured at each of the antennas.
 The ring of antennas with maximum electric fields in the order of $25\uVperm$ at the center of the array is the Cherenkov-ring.
-The Cherenkov-ring forms due to the forward beaming of the radio emissions of the airshower.
+The Cherenkov-ring forms due to the forward beaming of the radio emissions of the air shower.
 Outside this ring, the maximum electric field quickly falls with increasing distance to the array core.
 As expected for a vertical shower, the projection of the Cherenkov-ring on the ground is roughly circular.
 \\
@@ -177,7 +180,7 @@ Of course, a gaussian white noise component is introduced to the waveform as a s
 		%\textit{Right:}
 		%Example of the earliest and latest recorded air shower waveforms in the array as simulated by ZHAireS
 		Excerpt of a fully simulated waveform ($N=10240\,\mathrm{samples}$) (blue) containing the air shower (a $10^{16}\eV$~proton), the beacon (green, $\fbeacon = 51.53\MHz$) and noise.
-		The part of the waveform between the vertical dashed lines is considered airshower signal and masked before measuring the beacon parameters.
+		The part of the waveform between the vertical dashed lines is considered air shower signal and masked before measuring the beacon parameters.
 		\textit{Right:}
 		Fourier spectra of the waveforms.
 		The green dashed lines indicate the measured beacon parameters.
@@ -196,12 +199,12 @@ Moreover, it falls in the order of magnitude of clock defects that were found in
 
 % separate air shower from beacon
 To correctly recover the beacon from the waveform, it must be separated from the air shower.
-Typically, a trigger sets the location of the airshower signal in the waveform.
-In our case, the airshower signal is located at $t=500\ns$ (see Figure~\ref{fig:single:proton}).
+Typically, a trigger sets the location of the air shower signal in the waveform.
+In our case, the air shower signal is located at $t=500\ns$ (see Figure~\ref{fig:single:proton}).
 Since the beacon can be recorded for much longer than the air shower signal, we mask a window of $500$ samples around the maximum of the trace as the air shower's signal.
 % measure beacon phase, remove distance phase
 The remaining waveform is fed into a \gls{DTFT} \eqref{eq:fourier:dtft} to measure the beacon's phase $\pMeas$ and amplitude.
-Note that due to explicitly including a time axis in a \gls{DTFT}, a number of samples can be omitted without introducing artifacts.
+Note that due to explicitly including a time axis in a \gls{DTFT}, a number of samples can be omitted without introducing artefacts.
 \\
 With the obtained beacon parameters, the air shower signal is in turn reconstructed by subtracting the beacon from the full waveform in the time domain.
 \\
@@ -212,6 +215,61 @@ The small clock defect $\tClockPhase$ is then finally calculated from the beacon
 From the above, we now have a set of reconstructed air shower waveforms with corresponding clock defects smaller than one beacon period $T$.
 Shifting the waveforms to remove these small clocks defects, we are left with resolving the correct number of periods $k$ per waveform (see Figure~\ref{fig:grid_power:repair_phases}).
 \\
+% >>>>
+\section{\textit{k}-finding} % <<<
+
+% unknown origin of air shower signal
+Up until now, the shower axis and thus the origin of the air shower signal have not been resolved.
+This means that the unknown propagation time delays for the air shower ($\tProp$) affect the alignment of the signals in Figure~\ref{fig:beacon_sync:period_alignment} in addition to the unknown clock period defects ($k T$).
+As such, both this origin and the clock defects have to be determined simultaneously.
+\\
+% radio interferometry
+If the antennas had been fully synchronised, radio interferometry as introduced in Chapter~\ref{sec:interferometry} can be applied to find the origin of the air shower signal, thus resolving the shower axis.
+Still, a (rough) first estimate of the shower axis might be made using this technique or by employing other detection techniques such as those using surface or fluorescence detectors.
+% (see Figure~\ref{fig:dynamic-resolve}).
+\\
+
+On the true shower axis, the waveforms would sum most coherently when the correct $k$'s are used.
+Therefore, around the estimated shower axis, we define a grid search to both optimise this sum and the location of the maximum power.
+In this process each waveform of the array is allowed to shift by a restricted amount of periods with respect to a reference waveform (taken to be the waveform with the highest maximum).
+\\
+
+Note that these grids are defined here in shower plane coordinates with $\vec{v}$ the true shower axis and $\vec{B}$ the local magnetic field.
+Searching a grid that is slightly misaligned with the true shower axis is expected to give comparable results.
+\\
+
+The below $k$-finding algorithm is an iterative process where the grid around the shower axis is redefined on each iteration.
+Discussion is found in the next Chapter.
+\begin{enumerate}[start=1, label={Step \arabic*.}, ref=\arabic*, topsep=6pt, widest={Step 1.}]
+	\label{algo:kfinding}
+	\item \label{algo:kfinding:grid}
+		Define a grid around the estimated shower axis, zooming in on each iteration.
+
+	\item  \label{algo:kfinding:optimisation}
+		$k$-optimisation: per grid point, optimise the $k$'s to maximise the sum of the waveforms (see Figure~\ref{fig:single:k-correlation}).
+
+	\item  \label{algo:kfinding:kfinding}
+		$k$-finding: find the grid point with the maximum overall sum (see Figure~\ref{fig:findks:maxima}) and select its set of $k$'s.
+	
+	\item \label{algo:kfinding:break}
+		Stop when the set of $k$'s is equal to the set of the previous iteration, otherwise continue.
+
+	\item \label{algo:kfinding:powermapping}
+		Finally, make a power mapping with the obtained $k$'s to re-estimate the shower axis (location with maximum power) (see Figure~\ref{fig:findks:reconstruction}), and return to Step~\ref{algo:kfinding:grid} for another iteration.
+\end{enumerate}
+\vspace*{2pt}
+
+Here, Step~\ref{algo:kfinding:optimisation} has been implemented by summing each waveform to the reference waveform (see above) with different time delays $kT$ and selecting the $k$ that maximises the amplitude of a waveform combination.\footnote{%<<<
+	Note that one could use a correlation method instead of a maximum to select the best time delay.
+	However, for simplicity and ease of computation, this has not been implemented.
+	Other improvements are discussed in the next Section.
+} %>>>
+As shown in Figure~\ref{fig:single:k-correlation}, the maximum possible period shift has been limited to $\pm 3\,\mathrm{periods}$.
+This corresponds to the maximum expected time delay between two antennas with a clock randomisation up to $30\ns$ for the considered beacon frequency.%
+\footnote{
+	Figure~\ref{fig:simu:error:periods} shows this is not completely true.
+	However, overall, it still applies.
+}
 
 \begin{figure}%<<<
 	\centering
@@ -223,51 +281,12 @@ Shifting the waveforms to remove these small clocks defects, we are left with re
 	\label{fig:single:k-correlation}
 \end{figure}%>>>
 
-% >>>>
-\section{\textit{k}-finding} % <<<
-
-% unknown origin of air shower signal
-Up until now, the shower axis and thus the origin of the air shower signal have not been resolved.
-This means that the unknown propagation time delays for the air shower ($\tProp$) affect the alignment of the signals in Figure~\ref{fig:beacon_sync:period_alignment} in addition to the unknown clock period defects ($k_j T$).
-As such, both this origin and the clock defects have to be determined simultaneously.
-\\
-% radio interferometry
-If the antennas had been fully synchronised, radio interferometry as introduced in Section~\ref{sec:interferometry} can be applied to find the origin of the air shower signal, thus resolving the shower axis.
-Still, a rough first estimate of the shower axis might be made using this or other techniques.% (see Figure~\ref{fig:dynamic-resolve}).
-\\
-
-Starting with an initial grid around this estimated axis, a two-step process zooms in on the shower axis while optimising the interferometric signal.
-In this process each waveform of the array is allowed to shift by a restricted amount of periods with respect to a reference waveform (taken to be the waveform with the highest maximum).
-\\
-Note that these grids are defined in shower plane coordinates with $\vec{v}$ the true shower axis and $\vec{B}$ the local magnetic field.
-Searching a grid that is slightly misaligned with the true shower axis is expected to give comparable results.
-\\
-
-The first step consists of finding the set of period shifts $k_j$ that maximises the overall maximum amplitude on the grid.
-At each location, after removing propagation delays, each waveform and the reference waveform are summed with different time delays $kT$ to find the maximum attainable amplitude of the combined trace.%
-\footnote{%<<<
-	Note that one could use a correlation method instead of a maximum to select the best time delay.
-	However, for simplicity and ease of computation, this has not been implemented.
-%As shown in Figure~\ref{fig:single:annotated_full_waveform}, the air shower signal has a length in the order of a few nanoseconds.
-%Since it is this peak that is of interest, it would have been possible to cut the waveforms such to only correlate the peaks.
-} %>>>
-\\
-As shown in Figure~\ref{fig:single:k-correlation}, the maximum possible period shift has been limited to $\pm 3\,\mathrm{periods}$.
-This corresponds to the maximum expected time delay between two antennas with a clock randomisation up to $30\ns$ for the considered beacon frequency.
-\\
-%
-This amplitude optimisation is iterated over the grid (see Figure~\ref{fig:findks:maxima}) resulting in a grid measurement of the maximum attainable amplitude and the corresponding set of period defects $k$.
-\\
-
-The second step applies the obtained period defects and measures the interferometric power on the same grid (see Figure~\ref{fig:findks:reconstruction}).
-Afterwards, a new grid zooms in on the power maximum and the process is repeated (Figures~\ref{fig:findks:maxima:zoomed},~\ref{fig:findks:reconstruction:zoomed}) until the set of period defects is constant between zooming grids.
-\\
 
 \begin{figure}%<<< findks
 	\begin{subfigure}[t]{0.45\textwidth}
 		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.maxima.run0.pdf}
 		\caption{
-			Maximum amplitudes obtainable by shifting the waveforms.
+			$k$-finding: optimise the $k$'s by shifting waveforms to find the maximum amplitude obtainable at each point.
 		}
 		\label{fig:findks:maxima}
 	\end{subfigure}
@@ -275,7 +294,7 @@ Afterwards, a new grid zooms in on the power maximum and the process is repeated
 	\begin{subfigure}[t]{0.45\textwidth}
 		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.reconstruction.run0.power.pdf}
 		\caption{
-			Power measurement with the $k$'s belonging to the overall maximum of the amplitude maxima.
+			Power measurement with the $k$'s belonging to the overall maximum of the tested amplitudes.
 		}
 		\label{fig:findks:reconstruction}
 	\end{subfigure}
@@ -283,7 +302,8 @@ Afterwards, a new grid zooms in on the power maximum and the process is repeated
 	\begin{subfigure}[t]{0.45\textwidth}
 		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.maxima.run1.pdf}
 		\caption{
-			Maximum amplitudes, zoomed to the location in \ref{fig:findks:reconstruction} with the highest amplitude.
+			$2^\mathrm{nd}$ $k$-finding iteration:
+				Zoom in on the location in \subref{fig:findks:reconstruction} with the highest amplitude and repeat algorithm.
 		}
 		\label{fig:findks:maxima:zoomed}
 	\end{subfigure}
@@ -299,7 +319,8 @@ Afterwards, a new grid zooms in on the power maximum and the process is repeated
 	\begin{subfigure}[t]{0.45\textwidth}
 		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.maxima.run2.pdf}
 		\caption{
-			Final test grid obtaining the same $k$'s as Figure~\ref{fig:findks:maxima:zoomed}.
+			$3^\mathrm{rd}$ $k$-finding iteration:
+			The same set of $k$'s has been found and we stop the algorithm.
 		}
 		\label{fig:findks:maxima:zoomed2}
 	\end{subfigure}
@@ -312,7 +333,7 @@ Afterwards, a new grid zooms in on the power maximum and the process is repeated
 		\label{fig:findks:reconstruction:zoomed2}
 	\end{subfigure}
 	\caption{
-		Iterative $k$-finding algorithm:
+			Iterative $k$-finding algorithm (see page~\pageref{algo:kfinding} for explanation):
 		First \subref{fig:findks:maxima}, find the set of period shifts $k$ per point on a grid that returns the highest maximum amplitude (blue cross).
 		The grid starts as a $8^\circ$ wide shower plane slice at $X=400\,\mathrm{g/cm^2}$, centred at the true shower axis (red cross).
 		Second \subref{fig:findks:reconstruction}, perform the interferometric reconstruction with this set of period shifts.
@@ -332,22 +353,20 @@ Phase synchronising the antennas gives a small increase in observed power, while
 
 The initial grid plays an important role here in finding the correct axis.
 Due to selecting the highest maximum amplitude, and the process above zooming in aggressively, wrong candidate axes are selected when there is no grid-location sufficiently close to the true axis.
-Figure~\ref{fig:findks:reconstruction} shows such a potential point near $(-1, 0.5)$ with a maximum that is comparable to the selected maximum near the true axis.\Todo{no longer shown}
-\\
 % premature optimisation / degeneracy
 Such locations are subject to differences in propagation delays that are in the order of a few beacon periods.
 The restriction of the possible delays is therefore important to limit the number of potential axis locations.
 \\
-
 % fall in local extremum, maximum
-In this analysis, the initial grid is defined as $8^\circ$ wide around the true axis.
+In this analysis, the initial grid is defined as a very wide $8^\circ$ around the true axis.
 As the number of computations scales linearly with the number of grid points ($N = N_x N_y$), it is favourable to minimise the number of grid locations.
-Unfortunately, the above process has been observed to fall into local maxima when a too coarse initial grid ($N_x < 13$ at $X=400\,\mathrm{g/cm^2}$) was used while restricting the time delays to $\left| k \right| \leq 3$.
+Unfortunately, the above process has been observed to fall into local maxima when a too coarse and wide initial grid ($N_x < 13$ at $X=400\,\mathrm{g/cm^2}$) was used while restricting the time delays to $\left| k \right| \leq 3$.
 \\
 
 % Missing power / wrong k
 As visible in the right side of Figure~\ref{fig:grid_power:repair_full}, not all waveforms are in sync after the optimisation.
-In this case, the period defects have been resolved incorrectly for two waveforms, lagging 1 and 3 periods respectively (see Figure~\ref{fig:simu:error:periods}).
+In this case, the period defects have been resolved incorrectly for two waveforms (see Figure~\ref{fig:simu:error:periods}) due to too stringent limits on the allowable $k$'s.
+Looking at Figure~\ref{fig:grid_power:repair_phases}, this was to be predicted since there are two waveforms peaking at $k=4$ from the reference waveform's peak (dashed line).
 As a result, the obtained power for the resolved clock defects is slightly less than the obtained power for the true clocks.
 \\
 
@@ -357,17 +376,6 @@ Figure~\ref{fig:grid_power:axis} shows the power mapping at four different atmos
 Except for the low power case at $X=800\,\mathrm{g/cm^2}$, the shower axis is found to be $<0.1^\circ$ of the true shower axis.
 \\
 
-
-% Future: at multiple depths
-Of course, this algorithm must be evaluated at relevant atmospheric depths where the interferometric technique can resolve the air shower.
-In this case, after manual inspection, the air shower was found to have \Xmax\ at roughly $400\,\mathrm{g/cm^2}$.
-The algorithm is expected to perform as long as a region of strong coherent power is resolved.
-This means that with the power in both Figure~\ref{fig:grid_power:axis:X200} and Figure~\ref{fig:grid_power:axis:X600}, the clock defects and air shower should be identified to the same degree.
-\\
-Additionally, since the true period shifts are static per event, evaluating the $k$-finding algorithm at multiple atmospheric depths allows to compare the obtained sets thereof to further minimise any incorrectly resolved period defect.
-\\
-
-
 \begin{figure}% fig:simu:error
 	\centering
 	%\includegraphics[width=\textwidth]{ZH_simulation/cb_report_measured_antenna_offsets.py.time-amplitudes-missing-k.residuals.pdf}
@@ -438,6 +446,26 @@ Additionally, since the true period shifts are static per event, evaluating the
 	\label{fig:grid_power_time_fixes}
 \end{figure}%>>>
 
+\pagebreak
+% Future: at multiple depths
+Of course, this algorithm must be evaluated at relevant atmospheric depths where the interferometric technique can resolve the air shower.
+In this case, after manual inspection, the air shower was found to have \Xmax\ at roughly $400\,\mathrm{g/cm^2}$.
+The algorithm is expected to perform as long as a region of strong coherent power is resolved.
+This means that with the power in both Figure~\ref{fig:grid_power:axis:X200} and Figure~\ref{fig:grid_power:axis:X600}, the clock defects and air shower should be identified to the same degree.
+\\
+
+Additionally, since the true period shifts are static per event, evaluating the $k$-finding algorithm at multiple atmospheric depths allows to compare the obtained sets thereof to further minimise any incorrectly resolved period defect.
+\\
+
+Further improvements to the algorithm are foreseen in both the definition of the initial grid (Step~\ref{algo:kfinding:grid}) and the optimisation of the $k$'s (Step~\ref{algo:kfinding:optimisation}).
+For example, the $k$-optimisation step currently sums the full waveform for each $k$ to find the maximum amplitude for each sum.
+Instead, the timestamp of the amplitude maxima of each waveform can be compared, directly allowing to compute $k$ from the difference.
+\\
+
+Finally, from the overlapping traces in Figure~\ref{fig:grid_power:repair_full}, it is easily recognisable that some period defects have been determined incorrectly.
+Inspecting Figure~\ref{fig:grid_power:repair_phases}, this was to be expected as there are two waveforms with the peak at $\left|k\right| = 4$ from the reference waveform.
+Therefore, either the $k$-optimisation should have been run with a higher limit on the allowable $k$'s, or, preferably, these waveforms must be optimised after the algorithm is finished with a higher maximum $k$.
+
 \begin{figure}%<<< grid_power:axis:X600
 	\vspace*{-5mm}
 	\centering

documents/thesis/chapters/single_sine_interferometry.tex.bak

@@ -0,0 +1,509 @@
+% vim: fdm=marker fmr=<<<,>>>
+\documentclass[../thesis.tex]{subfiles}
+
+\graphicspath{
+	{.}
+	{../../figures/}
+	{../../../figures/}
+}
+
+\begin{document}
+\chapter[Single Sine Synchronisation]{Single Sine Beacon Synchronisation and Radio Interferometry}
+\label{sec:single_sine_sync}
+
+% <<<
+As shown in Chapter~\ref{sec:disciplining}, both impulsive and sine beacon signals can synchronise air shower radio detectors to enable the interferometric reconstruction of extensive air showers.
+This chapter will focus on using a single sine beacon to synchronise an array due to the simple setup and analysis required for such a beacon.
+Additionally, at \gls{Auger}, a public TV-transmitter is broadcasting at $67.25\MHz$.
+This poses an opportunity to use a ``free'' beacon to synchronise the radio antennas of \gls{AERA} and \gls{AugerPrime}.
+\\
+
+Due to the periodicity of sine beacons, the ability to synchronise an array is limited up to the beacon period $T$.
+As previously mentioned, the correct periods can be ascertained by choosing a beacon period much longer than the estimated accuracy of another timing mechanism.\footnote{For reference, \gls{GNSS} timing is expected to be below $30\ns$}
+Likewise, this can be achieved using the beating of multiple frequencies such as the four frequency setup in \gls{AERA}, amounting to a total period of $>1\us$.
+\\
+
+In this chapter, a different method of resolving these period mismatches is investigated by recording an impulsive signal in combination with the sine beacon.
+Figure~\ref{fig:beacon_sync:sine} shows the steps of synchronisation using this combination.
+The extra signal declares a shared time $\tTrueEmit$ that is common to the stations, after which the periods can be counted.
+Note that the period mismatch term $\Delta k_{ij}$ in \eqref{eq:synchro_mismatch_clocks_periodic} will be referenced throughout this Chapter as $k$ since we can take station $i$ as reference ($k_i =0$).
+\\
+
+\begin{figure}%<<<
+	\centering
+	\begin{subfigure}[t]{0.47\textwidth}
+		\includegraphics[width=\textwidth]{beacon/beacon_sync.pdf}
+		\caption{
+			Phase alignment
+		}
+		\label{fig:beacon_sync:syntonised}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}[t]{0.47\textwidth}
+		\includegraphics[width=\textwidth]{beacon/beacon_sync_period.pdf}
+		\caption{
+			Period alignment
+		}
+		\label{fig:beacon_sync:period_alignment}
+	\end{subfigure}
+	\caption{
+		Synchronisation scheme for two antennas using a single sine wave beacon (orange) and an impulsive signal (blue).
+		Vertical dashed lines indicate periods of the sine wave.
+		\subref{fig:beacon_sync:syntonised} A small time delay $t_\varphi$ is derived from the phase difference of the beacon as measured by the antennas.
+		\subref{fig:beacon_sync:period_alignment} The period mismatch $k$ is determined from the overlap between the impulsive signals.
+		%Expecting the impulsive signals to come from the same source, the overlap between the impulsive signals is used to determine the period mismatch $k$.
+		Note that the impulsive signals do not coincide perfectly due to different propagation delays from the source to the antennas.
+	}
+	\label{fig:beacon_sync:sine}
+%	\begin{subfigure}{\textwidth}
+%		\centering
+%		\includegraphics[width=\textwidth]{beacon/08_beacon_sync_timing_outline.pdf}
+%		\caption{
+%			Measure two waveforms at different antennas at approximately the same local time (clocks are not synchronised).
+%		}
+%		\label{fig:beacon_sync:timing_outline}
+%	\end{subfigure}
+%	\begin{subfigure}{\textwidth}
+%		\centering
+%		\includegraphics[width=\textwidth]{beacon/08_beacon_sync_synchronised_outline.pdf}
+%		\caption{
+%			The beacon signal is used to remove time differences smaller than the beacon's period.
+%			The detector clocks are now an unknown amount of periods out of sync.
+%		}
+%		\label{fig:beacon_sync:syntonised}
+%	\end{subfigure}
+%	\begin{subfigure}{\textwidth}
+%		\centering
+%		\includegraphics[width=\textwidth]{beacon/08_beacon_sync_synchronised_period_alignment.pdf}
+%		\caption{
+%			Lifting period degeneracy ($k=n-m=7$ periods) using the optimal overlap between impulsive signals.
+%		}
+%		\label{fig:beacon_sync:period_alignment}
+%	\end{subfigure}
+%	\caption{
+%		Synchronisation scheme for two antennas using a continuous beacon and an impulsive signal, each emitted from a separate transmitter.
+%		Vertical dashed lines indicate periods of the beacon (orange),
+%		solid lines indicate the time of the impulsive signal (blue).
+%		\\
+%		\subref{fig:beacon_sync:syntonised}: The beacon allows to resolve a small timing delay ($\Delta \tClockPhase$).
+%		\\
+%		\subref{fig:beacon_sync:period_alignment}: Expecting the impulsive signals to come from the same source, the overlap between the two impulsive signals is used to lift the period degeneracy ($k=n-m$).
+%	}
+%	\label{fig:beacon_sync:sine}
+\end{figure}%>>>
+
+% Same transmitter / Static setup
+When the beacon transmitter is also used to emit the signal defining $\tTrueEmit$, the number of periods $k$ can be obtained directly from the signal.
+However, if this calibration signal is sent from a different location, its time delays differ from the beacon's time delays.
+\\
+% Dynamic setup
+For static setups, these time delays can be resolved by measuring the involved distances or by taking measurements of the time delays over time.
+In dynamic setups, such as for transient signals, the time delays change per event and the distances are not known a priori.
+The time delays must therefore be resolved from the information of a single event.
+\\
+
+% Beacon + Impulsive -> discrete
+As shown in Chapter~\ref{sec:beacon:array}, an impulsive signal allows to reconstruct the direction of origin in a single event depending on the timing resolution of the array.
+Synchronising the array with a sine beacon, any clock mismatch is discretized into a number of periods $k$.
+This allows to improve the reconstruction by iterating the discrete clock mismatches during reconstruction.
+\\
+Of course, a limit on the number of periods is required to prevent over-optimisation.
+In general, they can be constrained using estimates of the accuracy of other timing mechanisms (see below).
+\\
+With a restricted set of allowed period shifts, we can alternate optimising the calibration signal's origin and optimising the set of period time delays of the array.
+\\
+% >>>
+
+%\section{Lifting the Period Degeneracy with an Air Shower}% <<<
+% <<<<
+% Airshower gives t0
+In the case of radio detection of air showers, the very signal of the air shower itself can be used as the calibration signal.
+This falls into the dynamic setup previously mentioned.
+The best period defects must thus be recovered from a single event.
+\\
+When doing the interferometric analysis for a sine beacon synchronised array, waveforms can only be delayed by an integer amount of periods, thereby giving discrete solutions to maximising the interferometric signal.
+
+\clearpage
+\section{Air Shower simulation}
+% simulation of proton E15 on 10x10 antenna
+To test the idea of combining a single sine beacon with an air shower, we simulated a set of recordings of a single air shower that also contains a beacon signal.
+\footnote{\url{https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction/-/tree/main/airshower_beacon_simulation}}
+\\
+The air shower signal was simulated by \acrlong{ZHAireS}\cite{Alvarez-Muniz:2010hbb} on a grid of 10x10 antennas with a spacing of $50\,\mathrm{meters}$.
+Each antenna recorded a waveform of $500$ samples with a samplerate of $1\GHz$ for each of the X,Y and Z polarisations.
+The air shower itself was generated by a $10^{16}\eV$ proton coming in under an angle of $20^\circ$ from zenith.
+%Figure~\ref{fig:single:proton_waveform} shows the earliest and latest waveforms recorded by the array with their true time.
+\\
+Figure~\ref{fig:single:proton_grid} shows the maximum electric field measured at each of the antennas.
+The ring of antennas with maximum electric fields in the order of $25\uVperm$ at the center of the array is the Cherenkov-ring.
+The Cherenkov-ring forms due to the forward beaming of the radio emissions of the airshower.
+Outside this ring, the maximum electric field quickly falls with increasing distance to the array core.
+As expected for a vertical shower, the projection of the Cherenkov-ring on the ground is roughly circular.
+\\
+
+%% add beacon
+A sine beacon ($\fbeacon = 51.53\MHz$) was introduced at a distance of approximately $75\,\mathrm{km}$ northwest of the array, primarily received in the X polarisation.
+The distance between the antenna and the transmitter results in a phase offset with which the beacon is received at each antenna.%
+\footnote{%<<<
+	The beacon's amplitude is also dependent on the distance. Although simulated, this effect has not been incorporated in the analysis as it is negligible for the considered grid and distance to the transmitter.
+} %>>>
+The beacon signal was recorded over a longer time ($10240\,\mathrm{samples}$), to be able to distinguish the beacon and air shower later in the analysis.
+\\
+The final waveform of an antenna (see Figure~\ref{fig:single:proton}) was then constructed by adding its beacon and air shower waveforms and band-passing with relevant frequencies (here $30$ and $80\MHz$ are taken by default).
+Of course, a gaussian white noise component is introduced to the waveform as a simple noise model (see Figure~\ref{fig:sine:snr_time_resolution} for a treatise on the timing accuracy of a sine beacon).
+\\
+
+\begin{figure}% <<<
+	\centering
+	\includegraphics[width=0.5\textwidth]{ZH_simulation/array_geometry_shower_amplitude.pdf}
+	\caption{
+		The 10x10 antenna grid used for recording the air shower.
+		Colours indicate the maximum electric field recorded at the antenna.
+		The Cherenkov-ring is clearly visible as a circle of radius $100\metre$ centered at $(0,0)$.
+	}
+	\label{fig:single:proton_grid}
+\end{figure}% >>>
+\begin{figure}% <<<
+	\begin{subfigure}[t]{0.49\textwidth}
+		\includegraphics[width=\textwidth]{ZH_simulation/ba_measure_beacon_phase.py.A74.no_mask.zoomed.pdf}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}[t]{0.49\textwidth}
+		\includegraphics[width=\textwidth]{ZH_simulation/ba_measure_beacon_phase.py.A74.fourier.pdf}
+	\end{subfigure}
+	\caption{
+		\textit{Left:}
+		%\textit{Right:}
+		%Example of the earliest and latest recorded air shower waveforms in the array as simulated by ZHAireS
+		Excerpt of a fully simulated waveform ($N=10240\,\mathrm{samples}$) (blue) containing the air shower (a $10^{16}\eV$~proton), the beacon (green, $\fbeacon = 51.53\MHz$) and noise.
+		The part of the waveform between the vertical dashed lines is considered airshower signal and masked before measuring the beacon parameters.
+		\textit{Right:}
+		Fourier spectra of the waveforms.
+		The green dashed lines indicate the measured beacon parameters.
+		The amplitude spectrum clearly shows a strong component at roughly $50\MHz$.
+		The phase spectrum of the original waveform shows the typical behaviour for a short pulse.
+	}
+	\label{fig:single:proton}
+\end{figure}% >>>
+
+% randomise clocks
+After the creation of the antenna waveforms, the clocks are randomised by sampling a gaussian distribution with a standard deviation of $30\ns$. % (see Figure~\ref{fig:simu:sine:periods:repair_none}).
+At a beacon period of $\sim 20\ns$, this ensures that multiple antennas have clock defects of at least one beacon period.
+This in turn allows for synchronisation mismatches of more than one beacon period.
+Moreover, it falls in the order of magnitude of clock defects that were found in \gls{AERA}\cite{PierreAuger:2015aqe}.
+\\
+
+% separate air shower from beacon
+To correctly recover the beacon from the waveform, it must be separated from the air shower.
+Typically, a trigger sets the location of the airshower signal in the waveform.
+In our case, the airshower signal is located at $t=500\ns$ (see Figure~\ref{fig:single:proton}).
+Since the beacon can be recorded for much longer than the air shower signal, we mask a window of $500$ samples around the maximum of the trace as the air shower's signal.
+% measure beacon phase, remove distance phase
+The remaining waveform is fed into a \gls{DTFT} \eqref{eq:fourier:dtft} to measure the beacon's phase $\pMeas$ and amplitude.
+Note that due to explicitly including a time axis in a \gls{DTFT}, a number of samples can be omitted without introducing artifacts.
+\\
+With the obtained beacon parameters, the air shower signal is in turn reconstructed by subtracting the beacon from the full waveform in the time domain.
+\\
+The small clock defect $\tClockPhase$ is then finally calculated from the beacon's phase $\pMeas$ by subtracting the phase introduced by the propagation from the beacon transmitter.
+\\
+
+% introduce air shower
+From the above, we now have a set of reconstructed air shower waveforms with corresponding clock defects smaller than one beacon period $T$.
+Shifting the waveforms to remove these small clocks defects, we are left with resolving the correct number of periods $k$ per waveform (see Figure~\ref{fig:grid_power:repair_phases}).
+\\
+% >>>>
+\section{\textit{k}-finding} % <<<
+
+% unknown origin of air shower signal
+Up until now, the shower axis and thus the origin of the air shower signal have not been resolved.
+This means that the unknown propagation time delays for the air shower ($\tProp$) affect the alignment of the signals in Figure~\ref{fig:beacon_sync:period_alignment} in addition to the unknown clock period defects ($k T$).
+As such, both this origin and the clock defects have to be determined simultaneously.
+\\
+% radio interferometry
+If the antennas had been fully synchronised, radio interferometry as introduced in Chapter~\ref{sec:interferometry} can be applied to find the origin of the air shower signal, thus resolving the shower axis.
+Still, a (rough) first estimate of the shower axis might be made using this technique or by employing other detection techniques such as those using surface or fluorescence detectors.
+% (see Figure~\ref{fig:dynamic-resolve}).
+\\
+
+On the true shower axis, the waveforms would sum most coherently when the correct $k$'s are used.
+Therefore, around the estimated shower axis, we define a grid search to both optimise this sum and the location of the maximum power.
+In this process each waveform of the array is allowed to shift by a restricted amount of periods with respect to a reference waveform (taken to be the waveform with the highest maximum).
+\\
+
+Note that these grids are defined here in shower plane coordinates with $\vec{v}$ the true shower axis and $\vec{B}$ the local magnetic field.
+Searching a grid that is slightly misaligned with the true shower axis is expected to give comparable results.
+\\
+
+The below $k$-finding algorithm is an iterative process where the grid around the shower axis is redefined on each iteration.
+Discussion is found in the next Chapter.
+\begin{enumerate}[start=1, label={Step \arabic*.}, ref=\arabic*, topsep=6pt, widest={Step 1.}]
+	\label{algo:kfinding}
+	\item \label{algo:kfinding:grid}
+		Define a grid around the estimated shower axis, zooming in on each iteration.
+
+	\item  \label{algo:kfinding:optimisation}
+		$k$-optimisation: per grid point, optimise the $k$'s to maximise the sum of the waveforms (see Figure~\ref{fig:single:k-correlation}).
+
+	\item  \label{algo:kfinding:kfinding}
+		$k$-finding: find the grid point with the maximum overall sum (see Figure~\ref{fig:findks:maxima}) and select its set of $k$'s.
+	
+	\item \label{algo:kfinding:break}
+		Stop when the set of $k$'s is equal to the set of the previous iteration, otherwise continue.
+
+	\item \label{algo:kfinding:powermapping}
+		Finally, make a power mapping with the obtained $k$'s to re-estimate the shower axis (location with maximum power) (see Figure~\ref{fig:findks:reconstruction}), and return to Step~\ref{algo:kfinding:grid} for another iteration.
+\end{enumerate}
+\vspace*{2pt}
+
+Here, Step~\ref{algo:kfinding:optimisation} has been implemented by summing each waveform to the reference waveform (see above) with different time delays $kT$ and selecting the $k$ that maximises the amplitude of a waveform combination.\footnote{%<<<
+	Note that one could use a correlation method instead of a maximum to select the best time delay.
+	However, for simplicity and ease of computation, this has not been implemented.
+	Other improvements are discussed in the next Section.
+} %>>>
+As shown in Figure~\ref{fig:single:k-correlation}, the maximum possible period shift has been limited to $\pm 3\,\mathrm{periods}$.
+This corresponds to the maximum expected time delay between two antennas with a clock randomisation up to $30\ns$ for the considered beacon frequency.%
+\footnote{
+	Figure~\ref{fig:simu:error:periods} shows this is not completely true.
+	However, overall, it still applies.
+}
+
+\begin{figure}%<<<
+	\centering
+	\includegraphics[width=0.8\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.run0.i1.kfind.pdf}
+	\caption{
+		Finding the maximum correlation for integer period shifts (best $k=-1$) between two waveforms recording the same (simulated) air shower.
+		Randomising the antenna clocks up to $30\ns$ and $\fbeacon = 51.53\MHz$ corresponds to at most $3$ periods of time difference between two waveforms.
+	}
+	\label{fig:single:k-correlation}
+\end{figure}%>>>
+
+
+\begin{figure}%<<< findks
+	\begin{subfigure}[t]{0.45\textwidth}
+		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.maxima.run0.pdf}
+		\caption{
+			$k$-finding: optimise the $k$'s by shifting waveforms to find the maximum amplitude obtainable at each point.
+		}
+		\label{fig:findks:maxima}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}[t]{0.45\textwidth}
+		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.reconstruction.run0.power.pdf}
+		\caption{
+			Power measurement with the $k$'s belonging to the overall maximum of the tested amplitudes.
+		}
+		\label{fig:findks:reconstruction}
+	\end{subfigure}
+	\\
+	\begin{subfigure}[t]{0.45\textwidth}
+		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.maxima.run1.pdf}
+		\caption{
+			$2^\mathrm{nd}$ $k$-finding iteration:
+				Zoom in on the location in \subref{fig:findks:reconstruction} with the highest amplitude and repeat algorithm.
+		}
+		\label{fig:findks:maxima:zoomed}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}[t]{0.45\textwidth}
+		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.reconstruction.run1.power.pdf}
+		\caption{
+			Power measurement of the new grid.
+		}
+		\label{fig:findks:reconstruction:zoomed}
+	\end{subfigure}
+	\\
+	\begin{subfigure}[t]{0.45\textwidth}
+		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.maxima.run2.pdf}
+		\caption{
+			$3^\mathrm{rd}$ $k$-finding iteration:
+			The same set of $k$'s has been found and we stop the algorithm.
+		}
+		\label{fig:findks:maxima:zoomed2}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}[t]{0.45\textwidth}
+		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.reconstruction.run2.power.pdf}
+		\caption{
+			Final power measurement.
+		}
+		\label{fig:findks:reconstruction:zoomed2}
+	\end{subfigure}
+	\caption{
+			Iterative $k$-finding algorithm (see page~\pageref{algo:kfinding} for explanation):
+		First \subref{fig:findks:maxima}, find the set of period shifts $k$ per point on a grid that returns the highest maximum amplitude (blue cross).
+		The grid starts as a $8^\circ$ wide shower plane slice at $X=400\,\mathrm{g/cm^2}$, centred at the true shower axis (red cross).
+		Second \subref{fig:findks:reconstruction}, perform the interferometric reconstruction with this set of period shifts.
+		Zooming on the maximum power \subref{fig:findks:maxima:zoomed},\subref{fig:findks:reconstruction:zoomed} repeat the steps until the $k$'s are equal between the zoomed grids \subref{fig:findks:maxima:zoomed2},\subref{fig:findks:reconstruction:zoomed2}.
+	}
+	\label{fig:findks}
+\end{figure}%>>>
+% >>>
+
+\clearpage
+%\phantomsection
+\section{Strategy / Result} %<<<
+
+Figure~\ref{fig:grid_power_time_fixes} shows the effect of the various synchronisation stages on both the alignment of the air shower waveforms, and the interferometric power measurement near the true shower axis.
+Phase synchronising the antennas gives a small increase in observed power, while further aligning the periods after the optimisation process significantly enhances this power.
+\\
+
+The initial grid plays an important role here in finding the correct axis.
+Due to selecting the highest maximum amplitude, and the process above zooming in aggressively, wrong candidate axes are selected when there is no grid-location sufficiently close to the true axis.
+% premature optimisation / degeneracy
+Such locations are subject to differences in propagation delays that are in the order of a few beacon periods.
+The restriction of the possible delays is therefore important to limit the number of potential axis locations.
+\\
+% fall in local extremum, maximum
+In this analysis, the initial grid is defined as a very wide $8^\circ$ around the true axis.
+As the number of computations scales linearly with the number of grid points ($N = N_x N_y$), it is favourable to minimise the number of grid locations.
+Unfortunately, the above process has been observed to fall into local maxima when a too coarse and wide initial grid ($N_x < 13$ at $X=400\,\mathrm{g/cm^2}$) was used while restricting the time delays to $\left| k \right| \leq 3$.
+\\
+
+% Missing power / wrong k
+As visible in the right side of Figure~\ref{fig:grid_power:repair_full}, not all waveforms are in sync after the optimisation.
+In this case, the period defects have been resolved incorrectly for two waveforms (see Figure~\ref{fig:simu:error:periods}) due to too stringent limits on the allowable $k$'s.
+Looking at Figure~\ref{fig:grid_power:repair_phases}, this was to be predicted since there are two waveforms peaking at $k=4$ from the reference waveform's peak (dashed line).
+As a result, the obtained power for the resolved clock defects is slightly less than the obtained power for the true clocks.
+\\
+
+% directional reconstruction
+This does not impede resolving the shower axis.
+Figure~\ref{fig:grid_power:axis} shows the power mapping at four different atmospheric depths for the resolved clock defects.
+Except for the low power case at $X=800\,\mathrm{g/cm^2}$, the shower axis is found to be $<0.1^\circ$ of the true shower axis.
+\\
+
+\begin{figure}% fig:simu:error
+	\centering
+	%\includegraphics[width=\textwidth]{ZH_simulation/cb_report_measured_antenna_offsets.py.time-amplitudes-missing-k.residuals.pdf}
+	\includegraphics[width=0.5\textwidth]{ZH_simulation/cb_report_measured_antenna_time_offsets.py.time-periods.residuals.pdf}
+	\caption{
+		Errors in the resolved period defects with respect to the true period defects.
+	}
+	\label{fig:simu:error:periods}
+\end{figure}
+
+\begin{figure}%<<< grid power time fixes
+	%\vspace{-2cm}
+	\vspace*{-5mm}
+	\centering
+	\begin{subfigure}[t]{1\textwidth}
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.repair_none.scale4d.pdf}
+		\hfill
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/trace_overlap/on-axis/dc_grid_power_time_fixes.py.repair_none.axis.X400.trace_overlap.zoomed.repair_none.pdf}
+		\vspace*{-7mm}
+		\caption{
+			Randomised clocks
+		}
+		\label{fig:grid_power:repair_none}
+	\end{subfigure}
+	%\hfill
+	\\
+	\begin{subfigure}[t]{1\textwidth}
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.repair_phases.scale4d.pdf}
+		\hfill
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/trace_overlap/on-axis/dc_grid_power_time_fixes.py.repair_phases.axis.X400.trace_overlap.zoomed.repair_phases.pdf}
+		\vspace*{-7mm}
+		\caption{
+			Phase synchronisation
+		}
+		\label{fig:grid_power:repair_phases}
+	\end{subfigure}
+	\\
+	\begin{subfigure}[t]{1\textwidth}
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.repair_full.scale4d.pdf}
+		\hfill
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/trace_overlap/on-axis/dc_grid_power_time_fixes.py.repair_full.axis.X400.trace_overlap.zoomed.repair_full.pdf}
+		\vspace*{-7mm}
+		\caption{
+			Resolved clocks
+		}
+		\label{fig:grid_power:repair_full}
+	\end{subfigure}
+	\\
+	\begin{subfigure}[t]{1\textwidth}
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.no_offset.scale4d.pdf}
+		\hfill
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/trace_overlap/on-axis/dc_grid_power_time_fixes.py.no_offset.axis.X400.trace_overlap.zoomed.no_offset.pdf}
+		\vspace*{-7mm}
+		\caption{
+			True clocks
+		}
+		\label{fig:grid_power:no_offset}
+	\end{subfigure}
+	\vspace*{-7mm}
+	\caption{
+		Different stages of array synchronisation (unsynchronised, beacon synchronised, $k$-resolved and true clocks)
+		and
+		their effect on (\textit{right}) the alignment of the waveforms at the true axis
+			and (\textit{left}) the interferometric power near the simulation axis (red plus).
+		The maximum power is indicated by the blue cross.
+		In the right panes the vertical dashed line indicates the maximum of the reference waveform.
+	}
+	\label{fig:grid_power_time_fixes}
+\end{figure}%>>>
+
+\pagebreak
+% Future: at multiple depths
+Of course, this algorithm must be evaluated at relevant atmospheric depths where the interferometric technique can resolve the air shower.
+In this case, after manual inspection, the air shower was found to have \Xmax\ at roughly $400\,\mathrm{g/cm^2}$.
+The algorithm is expected to perform as long as a region of strong coherent power is resolved.
+This means that with the power in both Figure~\ref{fig:grid_power:axis:X200} and Figure~\ref{fig:grid_power:axis:X600}, the clock defects and air shower should be identified to the same degree.
+\\
+
+Additionally, since the true period shifts are static per event, evaluating the $k$-finding algorithm at multiple atmospheric depths allows to compare the obtained sets thereof to further minimise any incorrectly resolved period defect.
+\\
+
+Further improvements to the algorithm are foreseen in both the definition of the initial grid (Step~\ref{algo:kfinding:grid}) and the optimisation of the $k$'s (Step~\ref{algo:kfinding:optimisation}).
+For example, the $k$-optimisation step currently sums the full waveform for each $k$ to find the maximum amplitude for each sum.
+Instead, the timestamp of the amplitude maxima of each waveform can be compared, directly allowing to compute $k$ from the difference.
+\\
+
+Finally, from the overlapping traces in Figure~\ref{fig:grid_power:repair_full}, it is easily recognisable that some period defects have been determined incorrectly.
+Inspecting Figure~\ref{fig:grid_power:repair_phases}, this was to be expected as there are two waveforms with the peak at $\left|k\right| = 4$ from the reference waveform.
+Therefore, either the $k$-optimisation should have been run with a higher limit on the allowable $k$'s, or, preferably, these waveforms must be optimised after the algorithm is finished with a higher maximum $k$.
+
+\begin{figure}%<<< grid_power:axis:X600
+	\vspace*{-5mm}
+	\centering
+	\begin{subfigure}[t]{1\textwidth}
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X200.repair_full.scale2d.pdf}
+		\hfill
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X200.repair_full.scale02d.pdf}
+		\vspace*{-7mm}
+		\caption{$X=200\,\mathrm{g/cm^2}$}
+		\label{fig:grid_power:axis:X200}
+	\end{subfigure}
+	\begin{subfigure}[t]{1\textwidth}
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.repair_full.scale2d.pdf}
+		\hfill
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.repair_full.scale02d.pdf}
+		\vspace*{-7mm}
+		\caption{$X=400\,\mathrm{g/cm^2}$}
+		\label{fig:grid_power:axis:X400}
+	\end{subfigure}
+	\begin{subfigure}[t]{1\textwidth}
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X600.repair_full.scale2d.pdf}
+		\hfill
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X600.repair_full.scale02d.pdf}
+		\vspace*{-7mm}
+		\caption{$X=600\,\mathrm{g/cm^2}$}
+		\label{fig:grid_power:axis:X600}
+	\end{subfigure}
+	\begin{subfigure}[t]{1\textwidth}
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X800.repair_full.scale2d.pdf}
+		\hfill
+		\includegraphics[width=0.47\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X800.repair_full.scale02d.pdf}
+		\vspace*{-7mm}
+		\caption{$X=800\,\mathrm{g/cm^2}$}
+		\label{fig:grid_power:axis:X800}
+	\end{subfigure}
+	\vspace*{-6mm}
+	\caption{
+		Interferometric power for the resolved clocks (from Figure~\ref{fig:grid_power:repair_full})  at four atmospheric depths for an opening angle of $2^\circ$(\textit{left}) and $0.2^\circ$(\textit{right}).
+		The simulation axis is indicated by the red plus, the maximum power is indicated by the blue cross.
+		Except for \subref{fig:grid_power:axis:X800} where there is no power, the shower axis is resolved within $0.1^\circ$ of the true shower axis.
+	}
+	\label{fig:grid_power:axis}
+\end{figure}
+% >>>
+\end{document}

documents/thesis/thesis.tex

@@ -79,6 +79,9 @@
 \subfile{chapters/appendix-random-phasor.tex}
 
 %%% Print Bibliography
+
+\titlespacing{\chapter}{0pt}{*1}{*2}
+
 \backmatter
 \phantomsection
 \addcontentsline{toc}{chapter}{Bibliography}