Résolution de la phase des observations de l'interféromètre

Résolvez la phase

La phase $ \ phi $ est similaire à Delay

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

Cependant, il y a deux points à noter.

  1. La phase est indéfinie $ 2n \ pi $. Il n'est pas évident que la phase basée sur la ligne de base soit, par exemple, $ 2 \ pi / 3 $ ou $ -4 \ pi / 3 $. Donc, $ \ hat {\ phi} _k = \ phi_j- \ phi_i \ pm 2n \ pi $, ce qui est problématique à résoudre ainsi que le délai.
  2. Cette relation ne peut être utilisée que lorsque la phase de clôture est 0. Si le corps céleste est une source ponctuelle, la phase de fermeture sera 0, mais cette expression relationnelle ne vaut pas pour un corps céleste large avec une structure compliquée.

Concernant 2., nous allons y remédier en sélectionnant une source ponctuelle comme corps céleste d'étalonnage de phase. Concernant 1., dans le cas où la phase de chaque ligne de base est à l'intérieur de $ \ pm \ pi / 2 $, elle peut être traitée de la même manière que la solution de retard. Sinon, il est généralement préférable de résoudre avec une solution complexe.

Résoudre comme un nombre complexe

Le terme de phase du gain d'antenne est

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

Il est exprimé comme. Il satisfait l'équation (4.1) et n'a pas d'indétermination de $ 2n \ pi $. Cependant, puisqu'il s'agit d'une équation non linéaire, il est nécessaire de la répéter par approximation séquentielle. Dans ce qui suit, $ E_j = e ^ {i \ phi_j} $ est utilisé pour simplifier la notation. Supposons que la valeur initiale $ E ^ 0_i $ soit suffisamment proche de la solution optimale. Le résidu basé sur la ligne de base est $ r_k = E_k --E ^ 0_j E_i ^ {0 *} $. Le vecteur de correction $ \ boldsymbol {c} $ donné à $ \ phi ^ 0_j $ pour minimiser la norme résiduelle $ \ sum_k r_k r ^ * _k $ est

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

Cela est illustré par l'équation. Notez que la phase de refant est mise à 0 et que la généralité n'est pas perdue, donc $ c_0 $ est toujours 0, pas inconnu. $ \ frac {\ partial} {\ partial \ phi_0} $ vaut également 0. Quand j'écris les éléments de la matrice du milieu $ P $,

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}

On dirait. Le gain de refant est $ E_0 = E ^ * _0 = 1 $. L'équation (4.3) est

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

Il est donné par les équations simultanées, donc si vous résolvez ceci pour $ \ boldsymbol {c} $

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

Ce sera. Ajoutez le $ \ boldsymbol {c} $ obtenu de cette manière à la valeur de phase initiale $ \ boldsymbol {\ phi} ^ 0 $, et ajoutez $ \ boldsymbol {\ phi} ^ 1 \ leftarrow \ boldsymbol {\ phi} ^ Améliorez la phase avec 0 + \ boldsymbol {c} $. Si vous répétez ceci, il finira par converger. Le jugement de convergence peut être jugé par la norme du vecteur de correction $ \ boldsymbol {c} $.

Économie de travail dans le calcul

Calculez manuellement $ P ^ {* T} P $ à l'avance.

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)

Et la même matrice que pour le delay. Par conséquent, la matrice inverse

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

Vous pouvez donc utiliser la même matrice que lorsque vous obtenez la solution pour le retard. Le calcul de $ P ^ {* T} \ boldsymbol {r} $ est un peu délicat, mais vous pouvez l'obtenir avec le code Python suivant. Il est inévitable de se retrouver dans une boucle avec le nombre d'antennes, mais une boucle avec le nombre de lignes de base peut être évitée.

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

Code récapitulatif

Pour résumer ce qui précède, nous montrons la fonction clphase_solve () qui résout la phase basée sur l'antenne à partir de la visibilité basée sur la ligne de base. L'argument Vis est la visibilité, qui est un vecteur dans lequel des nombres complexes de lignes de base sont alignés dans l'ordre canonique. IterNum est le nombre d'itérations, qui est fixé à 2 par défaut (convergence suffisante si le rapport S / N est bon). La valeur de retour est la phase basée sur l'antenne.

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
#

c'est tout. Retour à Mokuji

Recommended Posts

Résolution de la phase des observations de l'interféromètre
Résolution de l'amplitude des observations d'interféromètre
Résolution du gain complexe des observations d'interféromètre
Résoudre le retard d'observation de l'interféromètre
Résolution d'équations de mouvement en Python (odeint)
Le début de cif2cell
Le sens de soi
○○ Résolution de problèmes dans le département de mathématiques avec optimisation
L'histoire de sys.path.append ()
La vengeance des types: la vengeance des types
Aligner la version de chromedriver_binary
Grattage du résultat de "Schedule-kun"
10. Compter le nombre de lignes
L'histoire de la construction de Zabbix 4.4
Vers la retraite de Python2
Comparez les polices de jupyter-themes
Obtenez le nombre de chiffres
GoPiGo3 du vieil homme
Calculez le nombre de changements
Changer le thème de Jupyter
La popularité des langages de programmation
Changer le style de matplotlib
Visualisez la trajectoire de Hayabusa 2
À propos des composants de Luigi
Composants liés du graphique
Filtrer la sortie de tracemalloc
À propos des fonctionnalités de Python
Simulation du contenu du portefeuille
Le pouvoir des pandas: Python
Résolution du labyrinthe avec Python-Supplément au chapitre 6 de la référence rapide de l'algorithme-
Résolution des problèmes de sac à dos avec les outils OR de Google - Pratiquer les bases des problèmes d'optimisation combinée