Canonical Orderingによる干渉計の基線番号付け

# はじめに
干渉計の基線に番号を振る方法のうち、よく使われるcanonical ordering（調和順列）について解説します。Canonical orderingでない並びには、例えばlexical orderingがあり、ALMAデータのASDMはlexical orderingを採用しているのですが、較正処理などのデータ解析のコードが複雑になるので canonical orderingに変換してから較正処理します。

# 記号

* $N_{\rm a}$ : アンテナ数
* $N_{\rm b} = N_{\rm a} (N_{\rm a} - 1) / 2$ : 基線数（自己相関を含まない）
* $i < j$ : アンテナ番号
* $i = \{0, 1, 2, \dots, N_{\rm a}-2\}$
* $j = \{i+1, i+2, \dots, N_{\rm a}-1 \}$
* $k$ : 基線番号, $k = \{0, 1, 2, \dots,  N_{\rm b}-1 \}$
* refant : $i=0$が指すアンテナ
* kernel baseline : refantを含む基線

# Canonical ordering とは

|      | $i=0$ | $i=1$ | $i=2$ | $i$ | $i=N_{\rm a}-2$ |
|------|------:|------:|------:|--------:|----------------:|
|$j=1$ |     0 | 
|$j=2$ |     1 |      2|
|$j=3$ |     3 |      4|      5|
|$j$|$j(j-1)/2$|$j(j-1)/2+1$|$j(j-1)/2+2$|$j(j-1)/2+i$|
|$j=N_{\rm a}-1$ |$(N_{\rm a}-1)(N_{\rm a}-2)/2$ |$(N_{\rm a}-1)(N_{\rm a}-2)/2 + 1$| $(N_{\rm a}-1)(N_{\rm a}-2)/2 + 2$| $(N_{\rm a}-1)(N_{\rm a}-2)/2 + i$ |$N_{\rm b}-1$|

基線番号を$k = j(j-1)/2 + i$ というように順番付けます。Canonical orderingの特長は、順番付の式にアンテナ総数$N_{\rm a}$が含まれないので、アンテナ数によらず$i, j, k$の関係は不変であることです。だから、任意のアンテナ数のアレイにおけるアンテナ番号と基線番号の対応は、よりアンテナ数の大きいアレイにおける対応の部分集合になります。

# Python Codes
## Canonical orderingでのアンテナ番号と基線番号とを対応させるPythonコード

```Python
import numpy
from numpy import *
KERNEL_BL = arange(64)*arange(1,65)/2
#---- Antenna -> baseline index (without autocorr)
def Ant2Bl(ant0, ant1): 
  antenna0, antenna1 = max(ant0, ant1), min(ant0, ant1)
  return KERNEL_BL[antenna0 - 1] + antenna1
#
#---- Antenna -> baseline index and direction (True if inverted)
def Ant2BlD(ant0, ant1):
  antenna0 = max(ant0, ant1); antenna1 = min(ant0, ant1)
  return KERNEL_BL[antenna0 - 1], (ant0 < ant1)
#
#---- Baseline -> antenna index (without autocorr)
def Bl2Ant(bl_index):
  ant0 = max(np.where(KERNEL_BL<= bl_index)[0]) + 1
  return ant0, bl_index - KERNEL_BL[ant0 - 1]
#
```

上記のコードでKERNEL_BLはrefantを含む基線番号のnumpy arrayで、便利なのでグローバルに使います。64はALMAの最大アンテナ数で、想定する干渉計に応じて適当に調整するとよいでしょう。
関数Ant2Bl(ant0, ant1)は、２つのアンテナ番号入力すると対応する基線番号を返します。ant0 $>$ ant1が呼び出し側で保証されているなら、max, minの部分は省略してもよいでしょう。
関数Ant2BlD(ant0, ant1)はAnt2Bl(ant0, ant1)に加えて、基線反転フラグも返します。ant0 $<$ ant1が入力されると返り値の2番目はTrue (基線が反転している）、ant0 $>$ ant1ならFalseです。基線反転フラグがTrueの場合はVisibilityが複素共役になり、さらに偏波観測の場合には$\left< XY \right>$と$\left< YX \right>$が入れ替わります。
関数Bl2Ant(bl_index)は、基線番号を入力すると2つのアンテナ番号をtupleで返します。アンテナ番号は降順になります。

関数のままにしておくと高頻度に呼び出されてその度にオーバーヘッドが発生しますから、あらかじめlistを作っておくと便利です。

```Python
ANT0 = []; ANT1 = []     # List the BL -> antenna indexing
for bl_index in range(2016):    # Maximum number of baseline
  ants = Bl2Ant(bl_index)
  ANT0.append(ants[0])        # bl -> ant0 (baseline-end antenna) mapping
  ANT1.append(ants[1])        # bl -> ant1 (baseline-begin antenna) mapping
#
```

ここで2016はALMAの最大基線数です。想定する干渉計によって適当に調整するといいでしょう。このように作ったANT0およびANT1は、基線番号からアンテナ番号のmappingの役割をするので、例えば基線番号のリストを Bl = [0,4,6,8] と与えたときに ANT0[Bl] とか ANT1[Bl] とかすればその基線に使われるアンテナのリストがBlの順に得られます。

## Lexical ordering から Canonical orderingへ変換するPythonコード
ALMAで出力されるASDMでは、Lexical orderingでビジビリティデータが格納されています。Lexical orderingとは

|      | $i=0$ | $i=1$ | $i=2$ | $i$ | $i=N_{\rm a}-2$ |
|------|------:|------:|------:|--------:|----------------:|
|$j=1$ |     0 | 
|$j=2$ |     1 |     $N_{\rm a}-1$ |
|$j=3$ |     2 |    $N_{\rm a}$| $2N_{\rm a}-3$|
|$j$|$j-1$|$N_{\rm a}+j-3$|$2N_{\rm a} + j-6$|$N_{\rm a}i + j - i(i+3)/2-1$|
|$j=N_{\rm a}-1$ |$N_{\rm a}-２$ |$2N_{\rm a}-4$| $3N_{\rm a} - 7$| $N_{\rm a}(i+1)-i(i+3)/2-2$ |$N_{\rm b}-1$|

のように、$k = N_{\rm a}i + j - i(i+3)/2-1$ で表される並び順です。refantを含むkernel baselineが最初に並ぶ事が利点ですが、アンテナ番号と基線番号とを対応させるのにアンテナ数$N_{\rm a}$を与える必要があります。
なお、ALMAのASDMには自己相関も一緒に並んでいます。これをCanonical orderingとして処理するには並び替えが必要です。

``` Python
## Need to get the total number of antennas, such as : antNum = len(GetAntName(msfile)) where msfile is the file name of MS
#---- Convert lexical ordering (with autocorr) into canonical ordering
def Ant2Bla_RevLex(ant0, ant1, antNum):    # Reverse Lexical, with autcorr
  antenna0 = min(ant0, ant1); antenna1 = max(ant0, ant1)
  kernel = antNum* antenna0 - antenna0* (antenna0 - 1)/2
  return kernel + antenna1 - antenna0
#
blNum = antNum* (antNum - 1)/2
xcorr_index = range(blNum)
for bl_index in range(blNum): xcorr_index[bl_index] = Ant2Bla_RevLex(ANT0[bl_index], ANT1[bl_index], antNum)
#
```

このようにして作ったxcorr_indexというリストは、lexical orderingからcanonical orderingへ変換するマッピングに使えます。