Merge branch 'main' into step-up-2

documents/thesis/chapters/beacon_discipline.tex

@@ -76,7 +76,7 @@ The setup of an additional in-band synchronisation mechanism using a transmitter
 \\
 
 % time delay
-The distances between the transmitter $T$ and the antennas $A_i$ incur a time delay $(\tProp)_i$ caused by the finite propagation speed of the radio signal over these distances.
+The distance between the transmitter $T$ and the antenna $A_i$ incur a time delay $(\tProp)_i$ caused by the finite propagation speed of the radio signal over these distances.
 Since the signal is an electromagnetic wave, its instantanuous velocity $v$ depends solely on the refractive index~$n$ of the medium as $v = \frac{c}{n}$.
 In general, the refractive index of air is dependent on factors such as the pressure and temperature of the air the signal is passing through and the frequencies of the signal.
 However, in many cases, the refractive index can be taken constant over the trajectory to simplify models.
@@ -192,7 +192,7 @@ The dead time in turn, allows to emit and receive strong signals such as a singl
 Schemes using such a ``ping'' can be employed between the antennas themselves.
 Appointing the transmitter role to differing antennas additionally opens the way to calibrating the antennas in the array.
 \\
-Note the following method works fully in the time-domain.
+Note the following method works fully within the time-domain.
 
 % conceptually simple + filterchain response
 The detection of a pulse is conceptually simple.
@@ -337,19 +337,58 @@ This changes the synchronisation mismatches in \eqref{eq:synchro_mismatch_clocks
 	.
 \end{equation}
 
+\begin{figure}
+	\begin{subfigure}{\textwidth}
+		\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}
+		\includegraphics[width=\textwidth]{beacon/08_beacon_sync_synchronised_outline.pdf}
+		\caption{
+			Phase alignment syntonising the antennas using the beacon.
+		}
+		\label{fig:beacon_sync:syntonised}
+	\end{subfigure}
+	\begin{subfigure}{\textwidth}
+		\includegraphics[width=\textwidth]{beacon/08_beacon_sync_synchronised_period_alignment.pdf}
+		\caption{
+			Lifting period degeneracy ($k=m-n=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.
+		Grey dashed lines indicate periods of the beacon (orange),
+		full lines indicate the time of the impulsive signal (blue).
+		\\
+		Middle panel: The beacon allows to resolve a small timing delay ($\Delta t_\phase$).
+		\\
+		Lower panel: 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=m-n$).
+	}
+	\label{fig:beacon_sync:sine}
+	\todo{
+		Redo figure without xticks and spines,
+		rename $\Delta t_\phase$,
+		also remove impuls time diff?
+	}
+\end{figure}
+
 % lifting period multiplicity -> long timescale
 Synchronisation is possible with the caveat of being off by an unknown integer amount of periods $\Delta k_{ij}$.
 In phase-locked systems this is called syntonisation.
 There are two ways to lift this period degeneracy.
 \\
-First, if the timescale of the beacon is much longer than the estimated accuracy of another timing mechanism (such as \gls{GNSS}),
+First, if the timescale of the beacon is much longer than the estimated accuracy of another timing mechanism (such as a \gls{GNSS}),
  one can be confident to have the correct period.
 In AERA \cite{PierreAuger:2015aqe} for example, the total beacon repeats only after $\sim 1 \us$ (see Figure~\ref{fig:beacon:pa}).
 With an estimated accuracy of the \gls{GNSS} below $50 \ns$ the correct beacon period can be determined, resulting in a unique $\tTrueEmit$ transmit time\todo{reword}.
 \\
 
 % lifing period multiplicity -> short timescale counting +
-Another scheme is using an additional discrete signal to declare a unique $\tTrueEmit$.
+Another scheme is using an additional discrete signal to declare a unique $\tTrueEmit$ (see Figure~\ref{fig:beacon_sync:sine}).
 This relies on the ability of counting how many beacon periods have passed since the discrete signal has been recorded.
 
 \begin{figure}
@@ -370,9 +409,11 @@ In the following section, the scenario of a (single) sine wave as a beacon is wo
 It involves the tuning of the signal strength to attain the required accuracy.
 Later, a mechanism to lift the period degeneracy using an airshower as discrete signal is presented.
 
+
 %%
 %% Phase measurement
 \subsection{Phase measurement} % <<<
+% <<<
 A continuous beacon can syntonise an array of antennas by correcting for the measured difference in beacon phases $(\Delta \pMeasArriv)_{ij}$.
 They are derived by applying a \gls{FT} to the traces of each antenna.
 
@@ -416,11 +457,12 @@ These aspects are examined in the following section.
 	\label{fig:beacon:ttl_sine_beacon}
 \end{figure}
 % >>>
+%
 % DTFT
 \subsubsection{Discrete Time Fourier Transform}% <<<
 % FFT common knowledge ..
-The typical \gls{FT} to obtain spectral information from periodic data is the \gls{FFT} (a fast implementation of the \gls{DFT} \eqref{eq:fourier:dft}).
+The typical method to obtain spectral information from periodic data is the \gls{FFT} (a fast implementation of the \gls{DFT} \eqref{eq:fourier:dft}).
-Such an algorithm efficiently finds the magnitudes and phases within a trace $x$ at specific frequencies $f = f_s \tfrac{k}{N}$ determined solely by the number of samples $N$ ($0 \leq k < N$) and the sampling frequency $f_s$.
+Such an algorithm efficiently finds the amplitudes and phases within a trace $x$ at specific frequencies $f = f_s \tfrac{k}{N}$ determined solely by the number of samples $N$ ($0 \leq k < N$) and the sampling frequency $f_s$.
 \\
 
 % .. but we require a DTFT
@@ -431,78 +473,208 @@ Especially when a single frequency is of interest, a shorter route can be taken
 \\
 
 % DTFT from CTFT
-Spectral information in data can be obtained using a \acrlong{FT}.
+The continuous formulation of the \acrlong{FT} takes the following form,
 \begin{equation}
 	\label{eq:fourier}
-	X(f) = \frac{1}{2\pi} \int_\infty^\infty \dif{t}\, x(t)\, e^{i 2 \pi f t}
+	\phantom{.}
+	X(f) = \int_\infty^\infty \dif{t}\, x(t)\, e^{-i 2 \pi f t}
+	.
 \end{equation}
+It decomposes the signal $x(t)$ into complex-valued plane waves $X(f)$ of frequency $f$.
 
-
+When $x(t)$ is sampled at discrete times, the integral of \eqref{eq:fourier} is discretized in time to result in the \acrlong{DTFT}:
-The general (continuous) \gls{FT} \eqref{eq:fourier} can be discretized in time to result in the \acrlong{DTFT}:
 \begin{equation}
+	\tag{DTFT}
 	\label{eq:fourier:dtft}
-	X(f) = \frac{1}{2\pi N} \sum_{n=0}^{N-1}   x(t[n])\, e^{i 2 \pi f t[n]}
+	X(f) = \sum_{n=0}^{N-1}   x(t[n])\, e^{ -i 2 \pi f t[n]}
 \end{equation}
-where $X(f)$ is the transform of $x(t)$ at frequency $f$, sampled at $t[n]$.
+where $x(t) \in \mathcal{R} $ is sampled at times $t[n]$.
+Considering a finite sampling size $N$ and periodicity of the signal, the bounds of the integral in \eqref{eq:fourier} collapse to $t[0]$ up to $t[N]$.
 \\
+From this it follows that the lowest resolvable frequency is $f_\mathrm{lower} = \tfrac{1}{T} = \tfrac{1}{t[N] - 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 = \tfrac{n}{f_s}$, with $f_s$ the sampling frequency.
+The highest resolvable frequency, known as the Nyqvist frequency, is limited by this sampling frequency as $f_\mathrm{nyqvist} = \tfrac{f_s}{2}$.
+\\
+
+% 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 $k$ of the sampling frequency, becoming the \acrlong{DFT}
+\begin{equation*}
+	\label{eq:fourier:dft}
+	\phantom{,}
+	X(k) = \sum_{n=0}^{N-1}   x[n]\, \cdot e^{ -i 2 \pi {\frac{k n}N} }
+	,
+\end{equation*}
+with $k = \tfrac{f}{f_s}$.
+For integer $0 \leq k < N $, efficient algorithms exist that derive all $X( 0 \leq k < N )$ in $\mathcal{O}( N \log N )$ calculations, a~\acrlong{FFT}, sampling a subset of the frequencies.\Todo{citation?}
 
 \bigskip
 
-% DFT sampling of DTFT / efficient multifrequency FFT
+% Linearity fourier for real/imag
-When the sampling of $x(t)$ is equally spaced, the $t[n]$ terms can be decomposed as a sequence, $t[n] = \tfrac{n}{f_s}$ such that \eqref{eq:fourier:dtft} becomes the \acrlong{DFT}:
+In the previous equations, the resultant quantity $X(f)$ is a complex value.
+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}
-	\label{eq:fourier:dft}
 	\phantom{.}
-	X(k) = \frac{1}{N} \sum_{n=0}^{N-1}   x[n]\, \cdot e^{ i 2 \pi {\frac{k n}N} }
+	\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[n] \, \cos( 2\pi f t[n] )
+		  - i \sum_{n=0}^{N-1} \, x[n] \, \sin( 2\pi f t[n] )
 		.
+	\end{aligned}
 \end{equation}
 
 % FT term to phase and magnitude
-\bigskip
+The normalised amplitude at a given frequency $A(f)$ is calculated from \eqref{eq:fourier:dtft} as
-The magnitude of at frequency $f$ 
-
-
-\bigskip
-% Beacon frequency known -> single DTFT run
-When the beacon frequency is known, a single \gls{DTFT} needs to be evaluated.
-From this $X(f)$, the magnitude $A$ and phase $\pTrue$ are derived using
 \begin{equation}
-	\label{eq:magnitude_and_phase}
+	\label{eq:complex_magnitude}
-	\phantom.
+	\phantom{.}
-	A(f) = {\left|X(f)\right|}^2
+	A(f) \equiv \frac{ 2 \sqrt{ X_R(f)^2 + X_I(f)^2 } }{N}
-	\hfill
-	\pTrue(f) = \arctantwo\left(\Re(X(f)), \Im(X(f))\right)
 	.
 \end{equation}
-The decomposition of $X(f)$ into a real and imaginary part
+Likewise, the complex phase at a given frequency $\pTrue(f)$ is obtained by
+\begin{equation}
+	\label{eq:complex_phase}
+	\phantom{.}
+	\pTrue(f) \equiv \arctantwo\left( X_I(f), X_R(f) \right)
+	.
+\end{equation}
+\\
 
-With a constant beacon frequency, the coefficients for evaluating the \gls{DTFT} can be put into the hardware of the detectors.
+% 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[n] 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.
 
-% Beacon frequency unknown -> either zero-padding FFT or DTFT grid search
+\bigskip
+% 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$.
+Therefore, these can be precomputed ahead of time, reducing the number of calculations to $2N$ multiplications.
+
+% .. 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 phases in realtime.\Todo{figure?}
+
+
+
+% Beacon frequency known -> single DTFT run
+% Beacon frequency unknown -> either zero-padding FFT or, DTFT grid search
+%When the beacon frequency is known, a single \gls{DTFT} needs to be evaluated.
 
 
 % Removing the beacon from the signal trace
 
+% >>>
+%
 % >>>
 % Signal to noise
 \subsubsection{Signal to Noise}% <<<
 
 % Gaussian noise
-The traces will contain noise from various sources, both internal (e.g. LNA) and external (e.g. radio communications) to the detector.
+The traces will contain noise from various sources, both internal (e.g.~LNA~noise) and external (e.g.~radio~communications) to the detector.
-Adding gaussian noise to the traces in simulation gives a simple noise model, associated to many random noise sources.
+A simple noise model is given by gaussian noise in the time-domain, associated to many independent random noise sources.
 Especially important is that this simple noise model will affect the phase measurement depending on the strength of the beacon with respect to the noise level.
+
+In the following, this aspect is shortly described in terms of two frequency-domain phasors;
+the noise phasor written as $\vec{m} = a \, e^{i\pTrue}$ with phase $-\pi < \pTrue \leq \pi$ and amplitude $a \geq 0$,
+and the signal phasor written as $\vec{s} = s \, e^{i\pTrue_s}$, but rotated such that its phase $\pTrue_s = 0$.
+\Todo{reword; phasor vs plane wave}
+Further reading can be found in Ref.~\cite{goodman1985:2.9}.
 \\
 
+% Phasor concept
+\begin{figure}
+	\label{fig:phasor}
+	\caption{
+		Phasors picture
+	}
+\end{figure}
+
 \bigskip
 
-
+% Noise phasor description
-Phasor concept
+The noise phasor is fully described by the joint probability density function
-\cite{goodman1985:2.9}
-
-Known phasor $\vec{s}$ + random phasor $\vec{m} = a e^{i\pTrue}$ with $-\pi < \pTrue < \pi$ and $a \geq 0$.
-
 \begin{equation}
-	\label{eq:random_phasor_pdf}
+	\label{eq:noise:pdf:joint}
+	\phantom{,}
+	p_{A\PTrue}(a, \pTrue; \sigma)
+	=
+	\frac{a}{s\pi\sigma^2} e^{-\frac{a^2}{2\sigma^2}}
+	,
+\end{equation}
+for $-\pi < \pTrue \leq \pi$ and $a \geq 0$.
+\\
+
+Integrating \eqref{eq:noise:pdf:joint} over the amplitude $a$, it follows that the phase is uniformly distributed.
+
+Likewise, the amplitude follows a Rayleigh distribution
+\begin{equation}
+	\label{eq:noise:pdf:amplitude}
+	\label{eq:pdf:rayleigh}
+	\phantom{,}
+	p_A(a; \sigma)
+	%= p^{\mathrm{RICE}}_A(a; \nu = 0, \sigma)
+	= \frac{a}{\sigma^2} e^{-\frac{a^2}{2\sigma^2}}
+	,
+\end{equation}
+for which the mean is $\bar{a} = \sigma \sqrt{\frac{\pi}{2}}$ and the standard~deviation is given by $\sigma_{a} = \sigma \sqrt{ 2 - \tfrac{\pi}{2} }$.
+
+\begin{figure}
+	\begin{subfigure}{0.45\textwidth}
+		\includegraphics[width=\textwidth]{beacon/pdf_noise_phase.pdf}
+		\caption{
+			The phase of the noise is uniformly distributed.
+		}
+		\label{fig:noise:pdf:phase}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}{0.45\textwidth}
+		\includegraphics[width=\textwidth]{beacon/pdf_noise_amplitude.pdf}
+		\caption{
+			The amplitude of the noise is Rayleigh distribution \eqref{eq:noise:pdf:amplitude}.
+		}
+		\label{fig:noise:pdf:amplitude}
+	\end{subfigure}
+	\caption{
+		Marginal distribution functions of the noise phasor.
+		Rayleigh and Rice distributions.
+		\Todo{expand captions}
+	}
+	\label{fig:noise:pdf}
+\end{figure}
+
+\bigskip
+
+% Random phasor sum
+
+In this work, the addition of the signal phasor to the noise phasor will be named ``Random Phasor Sum''.
+The addition shifts the mean in \eqref{eq:noise:pdf:joint}
+from $\vec{a}^2 = a^2 {\left( \cos \pTrue + \sin \pTrue \right)}^2$
+to ${\left(\vec{a} - \vec{s}\right)}^2 = {\left( a \cos \pTrue -s \right)}^2 + {\left(\sin \pTrue \right)}^2$
+,
+resulting in a new joint distribution
+\begin{equation}
+	\label{eq:phasor_sum:pdf:joint}
+	\phantom{.}
 	p_{A\PTrue}(a, \pTrue; s, \sigma)
 	= \frac{a}{2\pi\sigma^2}
 		\exp[ -
@@ -513,45 +685,63 @@ Known phasor $\vec{s}$ + random phasor $\vec{m} = a e^{i\pTrue}$ with $-\pi < \p
 				2 \sigma^2
 			}
 		]
+	.
 \end{equation}
-requiring $ -\pi < 0 \leq pi $ and $a > 0$, otherwise $p_{A\PTrue} = 0$.
+\\
 
-\bigskip
+Integrating \eqref{eq:phasor_sum:pdf:joint} over $\pTrue$ one finds
-
+a Rice (or Rician) distribution for the amplitude,
-Noise only Amplitude:
-Rayleigh distribution
 \begin{equation}
-	\label{eq:amplitude_pdf:rayleigh}
+	\label{eq:phasor_sum:pdf:amplitude}
-	p_A(a; s=0, \sigma)
+	\label{eq:pdf:rice}
-	= p^{\mathrm{RICE}}_A(a; \nu = 0, \sigma)
+	\phantom{,}
-	= \frac{a}{\sigma^2} e^{-\frac{a^2}{2\sigma^2}}
+	p_A(a; s, \sigma)
-\end{equation}
-with $\sigma = \frac{\mu_1}{\sqrt{\frac{\pi}{2}}}$ and $\mu_2 = \frac{ 4 - \pi }{2}\sigma^2$.
-
-\bigskip
-Gaussian distribution
-\begin{equation}
-	\label{eq:amplitude_pdf:gauss}
-	p_A(a; \sigma) = \frac{1}{\sqrt{2\pi}} \exp[-\frac{{\left(a + s\right)}^2}{2\sigma^2}]
-\end{equation}
-
-
-Rician distribution ( 2D Gaussian at $\nu$ with $\sigma$ spread)
-\begin{equation}
-	\label{eq:amplitude_pdf:rice}
-	p^{\mathrm{RICE}}_A(a; s, \sigma)
 	= \frac{a}{\sigma^2}
 		\exp[-\frac{a^2 + s^2}{2\sigma^2}]
 		\;
 		I_0\left( \frac{a s}{\sigma^2} \right)
+	,
 \end{equation}
-with $I_0(z)$ the modified Bessel function of the first kind with order zero.\\
+where $I_0(z)$ is the modified Bessel function of the first kind with order zero.
-No signal $\mapsto$ Rayleigh ($s = 0$);\\
-Large signal $\mapsto$ Gaussian ($s \gg a$)
 
+For the Rician distribution, two extreme cases can be highlighted (as can be seen in Figure~\ref{fig:phasor_sum:pdf:amplitude}).
+In the case of a weak signal ($s \ll a$), \eqref{eq:phasor_sum:pdf:amplitude} behaves as a Rayleigh distribution~\eqref{eq:noise:pdf:amplitude}.
+Meanwhile, it approaches a gaussian distribution around $s$ when a strong signal ($s \gg a$) is presented.
+
+\begin{equation}
+	\label{eq:strong_phasor_sum:pdf:amplitude}
+	p_A(a; \sigma) = \frac{1}{\sqrt{2\pi}} \exp[-\frac{{\left(a - s\right)}^2}{2\sigma^2}]
+\end{equation}
+
+\begin{figure}
+	\begin{subfigure}{0.45\textwidth}
+		\includegraphics[width=\textwidth]{beacon/pdf_phasor_sum_phase.pdf}
+		\caption{
+			The Random Phasor Sum phase distribution \eqref{eq:phasor_sum:pdf:phase}.
+		}
+		\label{fig:phasor_sum:pdf:phase}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}{0.45\textwidth}
+		\includegraphics[width=\textwidth]{beacon/pdf_phasor_sum_amplitude.pdf}
+		\caption{
+			The Random Phasor Sum amplitude distribution \eqref{eq:phasor_sum:pdf:amplitude}.
+		}
+		\label{fig:phasor_sum:pdf:amplitude}
+	\end{subfigure}
+	\caption{
+			A signal phasor's amplitude in the presence of noise will follow a Rician distribution.
+			For strong signals, this approximates a gaussian distribution, while for weak signals, this approaches a Rayleigh distribution.
+		\Todo{expand captions}
+	}
+	\label{fig:phasor_sum:pdf}
+\end{figure}
 
 \bigskip
-Random Phasor Sum phase distribution: uniform (low $s$) + gaussian (high $s$)
+Like the amplitude distribution \eqref{eq:phasor_sum:pdf:amplitude}, the marginal phase distribution of \eqref{eq:phasor_sum:pdf:joint} results in two extremes cases;
+weak signals correspond to the uniform distribution for \eqref{eq:noise:pdf:joint}, while strong signals are well approximated by a gaussian distribution.
+
+The analytic form takes the following complex expression,
 \begin{equation}
 	\label{eq:phase_pdf:random_phasor_sum}
 	p_\PTrue(\pTrue; s, \sigma) =
@@ -565,28 +755,27 @@ Random Phasor Sum phase distribution: uniform (low $s$) + gaussian (high $s$)
 		\right)}{2}
 		\cos{\pTrue}
 \end{equation}
-with
+where
 \begin{equation}
 	\label{eq:erf}
+	\phantom{,}
 	\erf{\left(z\right)} = \frac{2}{\sqrt{\pi}} \int_0^z \dif{t} e^{-t^2}
+	,
 \end{equation}
-.
+is the error function.
-
-\bigskip
-Phase distribution: gaussian
-\begin{equation}
-	\label{eq:phase_pdf:gaussian}
-	p_\PTrue(\pTrue; s, \sigma) = \frac{1}{\sqrt{2} \sigma} \exp\left(- \frac{s^2}{2\sigma^2}  \right)
-\end{equation}
 
 \begin{figure}
 	\includegraphics[width=0.5\textwidth]{beacon/time_res_vs_snr.pdf}
-	\caption{Measured Time residuals vs Signal to Noise ration}
+	\caption{
+		Measured Time residuals vs Signal to Noise ratio
+	}
 	\label{fig:time_res_vs_snr}
 \end{figure}
 
 % Signal to Noise >>>
 
+% Phase measurement >>>
+%
 \subsection{Period degeneracy}% <<<
 % period multiplicity/degeneracy
 

figures/radio_interferometry/.gitignore

@@ -0,0 +1 @@
+rit_schematic_*.*

figures/radio_interferometry/Makefile

@@ -0,0 +1,42 @@
+SUBDIRS := $(subst Makefile,,$(wildcard */Makefile))
+
+.PHONY: all dist dist-clean $(SUBDIRS)
+
+all: dist $(SUBDIRS)
+
+dist: dist.png dist.pdf
+	#
+
+.PHONY: dist.png
+dist.png: \
+	rit_schematic_base.png \
+	rit_schematic_true.png \
+	rit_schematic_close.png \
+	rit_schematic_far.png \
+	#
+
+.PHONY: dist.pdf
+dist.pdf: \
+	rit_schematic_base.pdf \
+	rit_schematic_true.pdf \
+	rit_schematic_close.pdf \
+	rit_schematic_far.pdf \
+	#
+
+$(SUBDIRS):
+	@$(MAKE) -C $@
+
+dist-clean:
+	rm -v rit_schematic_*
+
+rit_schematic_base.%: src/rit_scheme.py
+	$< 'base' $@
+
+rit_schematic_true.%: src/rit_scheme.py
+	$< 'true' $@
+
+rit_schematic_close.%: src/rit_scheme.py
+	$< 'closeby' $@
+
+rit_schematic_far.%: src/rit_scheme.py
+	$< 'far-away' $@

figures/radio_interferometry/src/rit_scheme.py

@@ -0,0 +1,153 @@
+#!/usr/bin/env python3
+
+__doc__ = \
+"""
+Show geometry and time delay between radio antennas 
+and a source (true, and expected).
+"""
+
+import matplotlib.pyplot as plt
+from matplotlib.path import Path
+import matplotlib.patches as patches
+
+import numpy as np
+
+def antenna_path(loc, size=1, height=(.5)**(.5), width=1, stem_height=None):
+    stem_height = stem_height if stem_height is not None else height
+
+    vertices = [
+        (0, 0), # center, middle
+        ( width/2*size,  height/2*size), #right, top
+        ( width/2*size, -height/2*size), #right, bottom
+        (-width/2*size,  height/2*size), #left, top
+        (-width/2*size, -height/2*size), #left, bottom
+        (0, 0), # back to center
+        (0, -stem_height), # stem
+        (0, 0), # back to center
+    ]
+
+    codes = [
+            Path.MOVETO,
+            Path.LINETO,
+            Path.LINETO,
+            Path.LINETO,
+            Path.LINETO,
+            Path.LINETO,
+            Path.LINETO,
+            Path.CLOSEPOLY,
+    ]
+
+    # modify vertices
+    for i, v in enumerate(vertices):
+        vertices[i] = ( loc[0] + v[0], loc[1] + v[1] )
+
+    return Path(vertices, codes)
+
+
+def radio_interferometry_figure(emit_loc=(3,8), resolve_loc=None, N_antenna=4, ant_size=0.5, **fig_kwargs):
+    fig, ax = plt.subplots(**fig_kwargs)
+
+    annot_kwargs = dict(
+            color='red',
+            fontsize=15,
+    )
+
+    ray_kwargs = dict(
+            marker=None,
+            ls='solid',
+            color='red',
+            alpha=0.8
+    )
+    antenna_patch_kwargs = dict(
+            edgecolor='k',
+            facecolor='none',
+            lw=2
+    )
+
+    antenna_patches = []
+    dx_antenna = 1.5
+    stem_height = ant_size
+    for i in range(N_antenna):
+        path = antenna_path( (4+dx_antenna*i, stem_height), size=ant_size, stem_height=stem_height)
+        patch = patches.PathPatch(path, **antenna_patch_kwargs)
+        ax.add_patch(patch)
+        antenna_patches.append(patch)
+
+        if i == N_antenna - 1:
+            ant_loc = path.vertices[0]
+            ax.annotate("$\\vec{a_i}$", (ant_loc[0]+0.4, 0.8), va='top', ha='left', **{**annot_kwargs, **dict(color='k')})
+
+    # ground level
+    ax.axhline(0, color='k')
+
+    # indicate antenna signal
+    ant_loc = antenna_patches[int(2/4*N_antenna)].get_path().vertices[0]
+    ax.annotate("$S_i(t)$", (ant_loc[0], ant_loc[1]-stem_height-0.2), va='top', ha='center', **annot_kwargs)
+
+
+    if emit_loc or resolve_loc:
+        # resolve_loc
+        if resolve_loc is None:
+            resolve_loc = emit_loc
+
+        if resolve_loc:
+            for i, antenna in enumerate(antenna_patches):
+                    ant_loc = antenna.get_path().vertices[0]
+
+                    # rays from the antenna to resolve_loc
+                    ax.plot( (ant_loc[0], resolve_loc[0]), (ant_loc[1], resolve_loc[1]), **ray_kwargs)
+                    if i == N_antenna - 1:
+                        ax.annotate("$\Delta_i$", ( (ant_loc[0]+resolve_loc[0])/2, (ant_loc[1]+resolve_loc[1])/2 ), va='bottom', ha='left', **annot_kwargs)
+
+            ax.plot(*resolve_loc, 'ro')
+            ax.annotate("$S(\\vec{x}, t)$", resolve_loc, ha='left', va='bottom', **annot_kwargs)
+
+        # emit loc
+        if emit_loc:
+            ax.plot(*emit_loc, 'ko')
+
+            ax.annotate('$S_0$', emit_loc, ha='right', va='top', **{**annot_kwargs, **dict(color='k')})
+
+    ax.set_xlim(-1, 10)
+    ax.set_ylim(-1, 10)
+
+    ax.axis('off')
+
+    fig.tight_layout()
+
+    return fig
+
+
+if __name__ == "__main__":
+    emit_loc = (2.5, 8)
+    figsize = (6,6)
+
+    from argparse import ArgumentParser
+    import os.path as path
+
+    parser = ArgumentParser(description=__doc__)
+    parser.add_argument('scenario', choices=['base', 'true', 'closeby', 'far-away'], default='true')
+    parser.add_argument("fname", metavar="path/to/figure[/]", nargs="?", help="Location for generated figure, will append __file__ if a directory. If not supplied, figure is shown.")
+
+    args = parser.parse_args()
+
+    if args.fname is not None and path.isdir(args.fname):
+        args.fname = path.join(args.fname, path.splitext(path.basename(__file__))[0] + ".pdf")
+
+    ###
+    emit_loc = None
+    resolve_loc = None
+    if args.scenario != 'base':
+        emit_loc = (2.5, 8)
+
+        if args.scenario == 'closeby':
+            resolve_loc = (3, 8)
+        elif args.scenario == 'far-away':
+            resolve_loc = (8, 8)
+    
+    fig = radio_interferometry_figure(emit_loc=emit_loc, resolve_loc=resolve_loc, figsize=figsize)
+
+    if args.fname is not None:
+        plt.savefig(args.fname)
+    else:
+        plt.show()