The delay $ \ tau $ has a relationship of $ \ phi = \ phi_0 + 2 \ pi \ nu \ tau $ between the phase and the angular frequency $ 2 \ pi \ nu $. The delay $ \ hat {\ tau} $ obtained from the visibility per baseline is for the delay $ \ tau $ per antenna.
\hat{\tau}_k = \tau_j - \tau_i \tag{3.1}
It is represented by. $ i, j, k $ shall comply with canonical ordering. This relationship holds even if a common delay offset is given to all antennas, so setting the delay of the reference antenna (refant) to 0 does not lose generality.
The equation that solves the delay $ \ tau_i $ for each antenna from the delay $ \ hat {\ tau} _k $ measured for each baseline by the least squares method is
\left(
\begin{array}{c}
\hat{\tau}_0 \\
\hat{\tau}_1 \\
\hat{\tau}_2 \\
\hat{\tau}_3 \\
\hat{\tau}_4 \\
\hat{\tau}_5 \\
\vdots
\end{array}
\right) = \left( \begin{array}{cccc}
1 & 0 & 0 & \cdots \\
0 & 1 & 0 & \cdots \\
-1 & 1 & 0 & \cdots \\
0 & 0 & 1 & \cdots \\
-1 & 0 & 1 & \cdots \\
0 & -1 & 1 & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{array} \right)
\left(
\begin{array}{c}
\tau_1 \\
\tau_2 \\
\tau_3 \\
\vdots
\end{array}
\right) \tag{3.2}
It will be like that. Since we set $ \ tau_0 = 0 $, the number of unknowns is $ N_a -1 $. The middle $ (N_a --1) \ times N_b $ size matrix is $ P $, the vector of the observable $ \ hat {\ tau} _k $ is $ \ boldsymbol {Y} $, and the unknown number $ \ tau_i $ If the vector is $ \ boldsymbol {X} $, the simultaneous equations that obtain the delay for each antenna in the sense of least squares are
(P^TP) \boldsymbol{X} = P^T \boldsymbol{Y} \tag{3.3}
It will be. Therefore, $ \ boldsymbol {X} = (P ^ TP) ^ {-1} P ^ T \ boldsymbol {Y} $, so $ \ boldsymbol {X} $ can be obtained.
As the number of antennas and the number of baselines increases, the amount of calculation for $ P $ becomes enormous. ALMA has a maximum of 64 antennas and a maximum of 2016 baselines, so it is difficult to solve this. Fortunately, the calculation of $ P ^ TP $ is simplified as follows when canonical ordering is adopted.
P^TP = \left( \begin{array}{ccccc}
N_a - 1 & -1 & -1 & -1 & \cdots \\
-1 & N_a - 1 & -1 & -1 & \cdots \\
-1 & -1 & N_a-1 & -1 & \cdots \\
-1 & -1 & -1 & N_a-1 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{3.4}
In short, the diagonal component is $ N_a-1 $ and the other components are $ -1 $. Since this is a real symmetric matrix, it can be solved by using the Cholesky decomposition to obtain the lower triangular matrix $ L $ which can be expressed as $ LL ^ T = P $. Furthermore, the inverse matrix of $ (P ^ TP) $ can be easily calculated directly,
(P^TP)^{-1}= \frac{1}{N_a}\left( \begin{array}{ccccc}
2 & 1 & 1 & 1 & \cdots \\
1 & 2 & 1 & 1 & \cdots \\
1 & 1 & 2 & 1 & \cdots \\
1 & 1 & 1 & 2 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{3.5}
It will be. In short, a matrix of $ (N_a ―― 1) \ times (N_a ―― 1) $ size where the diagonal component is $ \ frac {2} {N_a} $ and the other components are $ \ frac {1} {N_a} $ is. If you write it in Python code, you can use antNum as the number of antennas in one line.
import numpy as np
PTP_inv = (np.diag(np.ones(antNum - 1)) + 1.0) / antNum
All that remains is the calculation of $ P ^ T \ boldsymbol {Y} $. Although we avoid calculating $ P ^ T $, which is a huge sparse matrix, we want to avoid loops with the number of baselines, so we perform vector operations within the loop with the number of antennas as shown below.
ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
PTY = np.zeros(antNum - 1)
for ant_index in range(1, antNum):
index0 = np.where(ant0 == ant_index)[0].tolist()
index1 = np.where(ant1 == ant_index)[0].tolist()
PTY[ant_index - 1] += np.sum(Y[index0]) - np.sum(Y[index1])
#
The lists ANT0 and ANT1 are prepared in advance by referring to this article Baseline numbering of interferometers by Canonical Ordering. index0 has "ant_index as the baseline end point" By taking the inner product of $ (P ^ TP) ^ {-1} $ and $ P ^ TY $ created in this way, an unknown vector $ X $ is obtained.
Delay_ant = [0.0] + (PTP_inv.dot(PTY)).tolist()
0.0 is inserted at the beginning of the list because it is a list including refant. By doing this, the delay for each baseline will be reduced.
Delay_bl = Delay_ant[ANT1] - Delay_ant[ANT0]
You can get it like this.
Here's a Python function that gets an antenna-based delay from a baseline-based delay. The argument bl_delay is a numpy array of baseline-based delays, which must be aligned in Canonical Ordering. The return value is an antenna-based delay, which is 0 because it starts with refant.
def cldelay_solve(bl_delay):
blNum = len(bl_delay); antNum = Bl2Ant(blNum)[0]
ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
PTP_inv = (np.diag(np.ones(antNum - 1)) + 1.0) / antNum
PTY = np.zeros(antNum - 1)
for ant_index in range(1, antNum):
index0 = np.where(ant0 == ant_index)[0].tolist()
index1 = np.where(ant1 == ant_index)[0].tolist()
PTY[ant_index - 1] += np.sum(bl_delay[index0]) - np.sum(bl_delay[index1])
#
return np.array( [0.0] + (PTP_inv.dot(PTY)).tolist() )
#
that's all.
Recommended Posts