Solving the phase of interferometer observations

Solve the phase

Phase $ \ phi $ is similar to Delay

\hat{\phi}_k = \phi_j - \phi_i \tag{4.1}

However, there are two points to note.

  1. The phase has $ 2n \ pi $ indefiniteness. It is not obvious whether the baseline-based phase is, for example, $ 2 \ pi / 3 $ or $ -4 \ pi / 3 $. So, $ \ hat {\ phi} _k = \ phi_j-\ phi_i \ pm 2n \ pi $, which is problematic to solve as well as delay.
  2. This relational expression can only be used when the closure phase is 0. If the celestial body is a point source, the closure phase will be 0, but this relational expression does not hold for a wide celestial body with a complicated structure.

Regarding 2., we will deal with it by selecting a point source as the phase calibration object. Regarding 1., in the case where the phase of each baseline is within $ \ pm \ pi / 2 $, it can be treated in the same way as the delay solution. Otherwise, it is generally better to solve with a complex solution.

Solve as a complex number

The phase term of the antenna gain is

e^{i\hat{\phi}_k} = e^{i\phi_j} e^{-i\phi_i} \tag{4.2}

It is expressed as. It satisfies equation (4.1) and has no indefiniteness of $ 2n \ pi $. However, since it is a non-linear equation, iteration by successive approximation is required. In the following, $ E_j = e ^ {i \ phi_j} $ is used to simplify the notation. Suppose the initial value $ E ^ 0_i $ is close enough to the optimal solution. The baseline-based residual is $ r_k = E_k --E ^ 0_j E_i ^ {0 *} $. The correction vector $ \ boldsymbol {c} $ given to $ \ phi ^ 0_j $ to minimize the residual norm $ \ sum_k r_k r ^ * _k $ is

\left(
\begin{array}{c}
r_0  \\
r_1  \\
r_2  \\
r_3  \\
r_4  \\
r_5  \\
\vdots
\end{array}
\right) = \left( \begin{array}{cccc}
 \frac{\partial {E_1 E^*_0}}{\partial \phi_1} &  \frac{\partial {E_1 E^*_0}}{\partial \phi_2} &  \frac{\partial {E_1 E^*_0}}{\partial \phi_3} & \cdots \\
 \frac{\partial {E_2 E^*_0}}{\partial \phi_1} &  \frac{\partial {E_2 E^*_0}}{\partial \phi_2} &  \frac{\partial {E_2 E^*_0}}{\partial \phi_3} & \cdots \\
 \frac{\partial {E_2 E^*_1}}{\partial \phi_1} &  \frac{\partial {E_2 E^*_1}}{\partial \phi_2} &  \frac{\partial {E_2 E^*_1}}{\partial \phi_3} & \cdots \\
 \frac{\partial {E_3 E^*_0}}{\partial \phi_1} &  \frac{\partial {E_3 E^*_0}}{\partial \phi_2} &  \frac{\partial {E_3 E^*_0}}{\partial \phi_3} & \cdots \\
 \frac{\partial {E_3 E^*_1}}{\partial \phi_1} &  \frac{\partial {E_3 E^*_1}}{\partial \phi_2} &  \frac{\partial {E_3 E^*_1}}{\partial \phi_3} & \cdots \\
 \frac{\partial {E_3 E^*_2}}{\partial \phi_1} &  \frac{\partial {E_3 E^*_2}}{\partial \phi_2} &  \frac{\partial {E_3 E^*_2}}{\partial \phi_3} & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{array} \right)
\left(
\begin{array}{c}
c_1  \\
c_2  \\
c_3  \\
\vdots
\end{array}
\right) \tag{4.3}

It is shown by the equation. Note that the phase of refant is set to 0 and generality is not lost, so $ c_0 $ is always 0, not unknown. $ \ frac {\ partial} {\ partial \ phi_0} $ is also 0. When I write out the elements of the matrix $ P $ in the middle,

P = \left( \begin{array}{cccc}
i{E_1 E^*_0} &  0 & 0 & \cdots \\
0 &  i{E_2 E^*_0} &  0 & \cdots \\
-i{E_2 E^*_0} &  i{E_2 E^*_0} &  0\\
0 &  0 & i{E_3 E^*_0}  & \cdots \\
-i{E_3 E^*_1} &  0 & i{E_3 E^*_1}  & \cdots \\
0 &  -i{E_3 E^*_2} & i{E_3 E^*_2}  & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{array} \right) \tag{4.4}

It looks like. The gain of refant is $ E_0 = E ^ * _0 = 1 $. Equation (4.3) is

\boldsymbol{r} = P \boldsymbol{c}

It is given by the simultaneous equations, so if you solve this for $ \ boldsymbol {c} $

P^{*T}\boldsymbol{r} = P^{*T} P \boldsymbol{c} \\
\boldsymbol{c} = (P^{*T} P)^{-1} P^{*T}\boldsymbol{r}

It will be. Add the $ \ boldsymbol {c} $ obtained in this way to the initial phase value $ \ boldsymbol {\ phi} ^ 0 $, and add $ \ boldsymbol {\ phi} ^ 1 \ leftarrow \ boldsymbol {\ phi} ^ Improve the phase as 0 + \ boldsymbol {c} $. If you repeat this, it will eventually converge. The convergence test can be judged by the norm of the correction vector $ \ boldsymbol {c} $.

Labor saving in calculation

Calculate $ P ^ {* T} P $ by hand in advance.

P^{*T} P = \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)

And the same matrix as for the delay. Therefore, the inverse matrix

(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)

Therefore, you can use Same matrix as when getting the solution of delay. The calculation of $ P ^ {* T} \ boldsymbol {r} $ is a bit tricky, but you can get it with the following Python code. It is unavoidable to end up in a loop with the number of antennas, but a loop with the number of baselines can be avoided.

PTY = np.zeros(antNum - 1, dtype=complex)
for ant_index in range(1,antNum):
  index0, index1 = np.where(ant0 == ant_index)[0].tolist(), np.where(ant1 == ant_index)[0].tolist()
  Y = np.zeros(blNum, dtype=complex)
  Y[index0] = -1.0j* antGain[ant_index].conjugate()* antGain[ant1[index0]]
  Y[index1] =  1.0j* antGain[ant_index]* antGain[ant0[index1]].conjugate()
  PTY[ant_index - 1] += Y.dot(resid)
#

Summary code

Summarizing the above, we show the function clphase_solve () that solves the antenna-based phase from the baseline-based visibility. The argument Vis is visibility, which is a vector of complex baseline numbers in Canonial order. IterNum is the number of iterations, which is set to 2 by default (if the S / N ratio is good, it will converge sufficiently). The return value is the antenna-based phase.

def clphase_solve(Vis, iterNum = 2):
    # Initial setup
    Vis = Vis / abs(Vis)    # Normalization (discard amplitude)
    blNum = len(Vis); antNum = Bl2Ant(blNum)[0]
    ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
    # Partial matrix
    PTP_inv = (np.diag(np.ones(antNum - 1)) + 1.0) / antNum
    antGain = np.append(1.0 + 0.0j, Vis[KERNEL_BL[0:antNum-1]])
    #---- Loop for iteration
    for iter_index in range(iterNum):
        resid = Vis - (antGain[ant0] * antGain[ant1].conjugate())
        PTY = np.zeros(antNum - 1, dtype=complex)
        #---- Calculate P*T r
        for ant_index in range(1,antNum):
            index0 = np.where(ant0 == ant_index)[0].tolist()
            index1 = np.where(ant1 == ant_index)[0].tolist()
            Y = np.zeros(blNum, dtype=complex)
            Y[index0] = -1.0j* antGain[ant_index].conjugate()* antGain[ant1[index0]]
            Y[index1] =  1.0j* antGain[ant_index]* antGain[ant0[index1]].conjugate()
            PTY[ant_index - 1] += Y.dot(resid)
        #
        #---- Correction for the antenna-based phase
        antPhs  = np.append(0, np.angle(antGain[1:antNum]) + PTP_inv.dot(PTY.real))
        antGain = np.cos(antPhs) + 1.0j* np.sin(antPhs)
    #
    return antPhs
#

that's all. Return to the lottery

Recommended Posts

Solving the phase of interferometer observations
Solving the amplitude of interferometer observations
Solving the complex gain of interferometer observations
Solve the delay of interferometer observation
Solving the equation of motion in Python (odeint)
The beginning of cif2cell
The meaning of self
○○ Solving problems in the Department of Mathematics by optimization
The story of sys.path.append ()
python chrome driver ver. Solving the problem of difference
Revenge of the Types: Revenge of types
Align the version of chromedriver_binary
Scraping the result of "Schedule-kun"
10. Counting the number of lines
The story of building Zabbix 4.4
Towards the retirement of Python2
[Apache] The story of prefork
Compare the fonts of jupyter-themes
About the ease of Python
Get the number of digits
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
Solving the Maze with Python-Supplement to Chapter 6 of the Algorithm Quick Reference-
Solving Knapsack Problems with Google's OR-Tools-Practicing the Basics of Combinatorial Optimization Problems