[PYTHON] Résoudre le retard d'observation de l'interféromètre

introduction

Le retard $ \ tau $ a une relation de $ \ phi = \ phi_0 + 2 \ pi \ nu \ tau $ entre la phase et la fréquence angulaire $ 2 \ pi \ nu $. Le retard $ \ hat {\ tau} $ obtenu à partir de la visibilité pour chaque ligne de base est le retard $ \ tau $ pour chaque antenne.

\hat{\tau}_k = \tau_j - \tau_i \tag{3.1}

Il est représenté par. $ i, j, k $ seront soumis à un ordre canonique. Cette relation est valable même si un décalage de retard commun est donné à toutes les antennes, de sorte que le réglage du retard de l'antenne de référence (refant) à 0 ne perd pas sa généralité.

Méthode du carré minimum

L'équation qui résout le retard $ \ tau_i $ pour chaque antenne à partir du retard $ \ hat {\ tau} _k $ mesuré pour chaque ligne de base par la méthode du carré minimum est

\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}

Ce sera comme ça. Puisque nous posons $ \ tau_0 = 0 $, le nombre d'inconnues est $ N_a -1 $. La matrice de taille moyenne $ (N_a --1) \ times N_b $ est $ P $, le vecteur de montant observé $ \ hat {\ tau} _k $ est $ \ boldsymbol {Y} $, nombre inconnu $ \ tau_i $ Si le vecteur est $ \ boldsymbol {X} $, les équations simultanées qui obtiennent le retard pour chaque antenne au sens du carré minimum sont

(P^TP) \boldsymbol{X} = P^T \boldsymbol{Y} \tag{3.3}

Ce sera. Par conséquent, $ \ boldsymbol {X} = (P ^ TP) ^ {-1} P ^ T \ boldsymbol {Y} $, donc $ \ boldsymbol {X} $ peut être obtenu.

Économie de travail dans le calcul

À mesure que le nombre d'antennes et de lignes de base augmente, la quantité de calcul pour $ P $ devient énorme. ALMA a un maximum de 64 antennes et un maximum de lignes de base 2016, il est donc difficile de résoudre ce problème. Heureusement, le calcul de $ P ^ TP $ est simplifié comme suit lors de l'utilisation de l'ordre canonique.

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}

En bref, la composante diagonale est $ N_a-1 $ et les autres composantes sont $ -1 $. Comme il s'agit d'une matrice symétrique réelle, elle peut être résolue en utilisant la décomposition de Cholesky pour obtenir la matrice triangulaire inférieure $ L $ qui peut être exprimée comme $ LL ^ T = P $. De plus, la matrice inverse de $ (P ^ TP) $ peut être facilement calculée directement,

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

Ce sera. En bref, une matrice de taille $ (N_a --1) \ times (N_a --1) $ dont les composantes diagonales sont $ \ frac {2} {N_a} $ et les autres composantes sont $ \ frac {1} {N_a} $. est. Si vous l'écrivez en code Python, vous pouvez utiliser antNum comme nombre d'antennes sur une ligne.

import numpy as np
PTP_inv = (np.diag(np.ones(antNum - 1)) + 1.0) / antNum

Il ne reste plus qu'à calculer $ P ^ T \ boldsymbol {Y} $. Même si nous évitons de calculer $ P ^ T $, qui est une énorme matrice creuse, nous voulons éviter les boucles avec le nombre de lignes de base, nous effectuons donc des opérations vectorielles dans la boucle avec le nombre d'antennes comme indiqué ci-dessous.

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])
#

Les listes ANT0 et ANT1 sont créées à l'avance en se référant à cet article Numérotation de la ligne de base des interféromètres par commande canonique. index0 a "ant_index comme point final de la ligne de base" En prenant le produit interne de $ (P ^ TP) ^ {-1} $ et $ P ^ TY $ créé de cette manière, un vecteur inconnu $ X $ peut être obtenu.

Delay_ant = [0.0] + (PTP_inv.dot(PTY)).tolist()

0.0 est inséré au début de la liste car il s'agit d'une liste comprenant le refant. En faisant cela, le délai pour chaque ligne de base

Delay_bl = Delay_ant[ANT1] - Delay_ant[ANT0]

Vous pouvez l'obtenir comme ça.

Code récapitulatif

Vous trouverez ci-dessous une fonction Python qui obtient un retard basé sur l'antenne à partir d'un retard basé sur la ligne de base. L'argument bl_delay est un tableau numpy de retards basés sur la ligne de base, qui doivent être alignés dans Canonical Ordering. La valeur de retour est un retard basé sur l'antenne, qui est 0 car il commence par 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() )
#

c'est tout.

Retour à Mokuji