Solve the delay of interferometer observation

Introduction

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.

Least squares

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.

Labor saving in calculation

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.

Summary code

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.

Return to the lottery

Recommended Posts

Solve the delay of interferometer observation
Solving the amplitude of interferometer observations
Solving the phase of interferometer observations
Solving the complex gain of interferometer observations
[At Coder] Solve the problem of binary search
The beginning of cif2cell
The meaning of self
Try to solve the problems / problems of "Matrix Programmer" (Chapter 1)
the zen of Python
The story of sys.path.append ()
Revenge of the Types: Revenge of types
Solve the verbal arithmetic
Try to solve the problems / problems of "Matrix Programmer" (Chapter 0 Functions)
Align the version of chromedriver_binary
Scraping the result of "Schedule-kun"
10. Counting the number of lines
Solve Copy-v0 of OpenAI Gym
The story of building Zabbix 4.4
Towards the retirement of Python2
Solve the Monty Hall problem
[Apache] The story of prefork
Compare the fonts of jupyter-themes
About the ease of Python
Get the number of digits
Explain the code of Tensorflow_in_ROS
Reuse the results of clustering
GoPiGo3 of the old man
Calculate the number of changes
Change the theme of Jupyter
The popularity of programming languages
Change the style of matplotlib
Visualize the orbit of Hayabusa2
About the components of Luigi
Connected components of the graph
Filter the output of tracemalloc
About the features of Python
Simulation of the contents of the wallet
The Power of Pandas: Python
Solve the problem of missing libcudart in Ubuntu 16.04 + CUDA 8.0 + Tensorflow environment
Try to solve the N Queens problem with SA of PyQUBO
Solve the initial value problem of ordinary differential equations with JModelica