フラックス測定観測の方法と解析

パワーメーターを用いて天体のフラックスを求める観測の方法と解析手順を示す。スクリプト一式は自分のPCに落として使えばよいが、一部のスクリプトにはasptropyライブラリがインストールされたpython環境が必要。以下ではcorr3を使った実行例を示す。

観測スケジュール作成

作業ディレクトリを作成しそこに移動する。

[oper@corr3 ~]$ cd ~/flux_measure/
[oper@corr3 ~]$ mkdir F19001A ← 適当な観測コードを設定してディレクトリ作成
[oper@corr3 ~]$ cd F19001A

アンテナ駆動スケジュールファイルをPythonスクリプト make_skd_flux_measure.py で作成する。スケジュール作成のたびに、スクリプト自体を編集する。

[oper@corr3 F19001A]$ vi ~/bin/make_skd_flux_measure.py

EXPER_prefix      = "F" # F of F19001A
EXPER_suffix      = "A" # A of F19001A
rec_start_UTC     = datetime.datetime(2019, 12, 10, 18, 30, 00)
num_seq           = 1   # Number of observation sequences. One sequence consists of multiple pointings.
pointing_gap_s    = 10.0 # Gap time (in second) between the pointings
on_source_s       = 10.0 # On-source time (in second)
az_default_offset = 0    # in arc-minutes
el_default_offset = 0

target_info = '3C48    01 37 41.2995845985    +33 09 35.079126038'

azel_arcmin = ( # Cross-scanning
    (-10, 0), # OFF-POINTING REQUIRED
    (-2, 0), (-1, 0), (0, 0), (+1, 0), (+2, 0),
    (+10, 0),
    (0, -10), 
    (0, -2), (0, -1), (0, 0), (0, +1), (0, +2),  
    (0, +10) # OFF-POINTING REQUIRED
)

各変数は以下の通り。

rec_start_UTC: 観測開始時刻 (UTC)。例では 2019-12-10 18:30:00 UTC に指定している。
num_seq: 観測シーケンスの数。ここでいうシーケンスとは、例えば5点法観測を1回することであり、1シーケンスで、その仰角におけるフラックス測定は完了する。例えばシーケンス数を10回とすると、N点観測によるフラックス測定を10回連続的に実施する。
pointing_gap_s: あるポインティングから次のポインティングに移るときの猶予時間 (秒単位)。望遠鏡の移動が完了して駆動が安定するのに十分な時間でなければならない。
on_source_s: 天体を観ている時間 (秒単位)。望遠鏡の安定性、人工電波の影響、天体の仰角変化などを考慮した値でなくてはならない。この値が大きすぎると、フラックス測定の間に天体の仰角が変わりすぎてまともな測定にならない。
target_info: 天体の名前、赤経 (hour:min:sec)、赤緯 (deg:arcmin:arcsec)。
azel_arcmin: ポインティングオフセットの座標 (AZ, EL) を分角単位で指定。例では、望遠鏡が中心だと思っている指向位置から (AZ, EL) を (-10', 0') だけずらした位置から観測を開始し、次に (-2', 0'), ..., (0', +10') という位置を観測する。この時、最初と最後の指向位置は、必ずオフ点 (sky) でなければならない。

スクリプトを編集したら実行する。するとカレントディレクトリにskdファイルが作成される。

[oper@corr3 F19001A]$ make_skd_flux_measure.py
Saved: /???/F19001A/F19001Aa.skd

観測する

作成したスケジュールファイルを用いてアンテナを駆動し、パワーメーターを用いて電力測定する。パワーのデータファイル名は任意 (e.g., F19001A-power.txt)。 skdファイルを32m制御PCに転送するコマンドは次の通り。

[oper@corr3 F19001A]$ scp  F19001Aa.skd  ymg32:~/F19001Aa.skd

データファイルの回収

解析に必要な次の2つのファイルを回収する:

アンテナ駆動スケジュールファイル (e.g., F19001Aa.skd)
パワーデータファイル (e.g., F19001A-power.txt)

上記作業により、カレントディレクトリには既にskdがある。ここにパワーメーターデータファイルを下記の通り転送してくる。

[oper@corr3 F19001A]$ scp  powermonit:/Users/astro/VEE_Programs/F19001A-power.txt  ./

解析する

スケジュールファイルとパワーデータファイルから、ビームパターンを求める。後述するスクリプト群を一括して使うインターフェーススクリプト flux_measure.sh を使って、とりあえず次のように実行するとビームパターンが求まり、フラックスを算出できる。

[oper@corr3 F19001A]$ flux_measure.sh
使い方が表示される。

このシェルスクリプトで行っている内容は下記の通り。

パワーデータファイルのフォーマットを変換する

スクリプト power_data_formatter.py を用いて、茨城大学で使用しているパワーデータのフォーマットに変換する。

$ power_data_formatter.py
使い方が表示される。

各スキャン毎の電力を求める

スクリプト flux_measure_1_power_vs_scan.py を用いて、スケジュールファイル *.skd とパワーデータファイル *.txt を照合し、スキャンごとの電力値を求める。

$ flux_measure_1_power_vs_scan.py
使い方が表示される。

スキャンごとの電力値として、スキャン中の電力の平均値と中央値を求めるが、以降の解析では中央値を使うと都合が良いことが多い。

ビームパターンを求める

スクリプト flux_measure_2_fit.py を用いて、スキャンごとの電力値からそのデータにフィットするビームパターンを求める。モデル関数として、ガウス関数とエアリー関数の2つを出力できる。実行にはAstropyライブラリが必要。

$ flux_measure_2_fit.py
使い方が表示される。

解析にかけるとまず、各スキャンごとの電力データを時系列で表示し、また最初と最後のデータ点を用いたベースライン除去を行う。

グラフが、観測時のスキャン方法に合った形状であることを確認する。 AZ方向とEL方向に十字スキャンした場合は、図のような2つの山をもった時系列データになる。 エラーバーが極端に大きい場合は、人工電波の干渉など何らかの原因で観測失敗しているので観測し直す。

表示されたグラフウィンドウを×ボタンで消すと、次に関数フィットした3つのグラフが表示される。 2つの立体グラフのうち、Cyanで色付けされた関数は2次元ガウス関数、Magentaで色付けされた関数はエアリーディスク関数であり、系列名欄にそれぞれの関数パラメータが示されている。これらのグラフはPDFファイルとしても出力されており、そのファイル名は標準出力に書かれている。 1つの平面グラフは、フィットしたモデル関数との残差プロットであり、系列名欄に、その残差の統計量として平均値、偏差、歪度を示している。

より詳細な関数パラメータは標準出力に出力される。

Fitted Gaussian2D Parameters:
    amplitude (= Ta / Tsys)     = 2.259089e-02
    x_mean (AZ OFFSET) / arcmin = -0.086471
    y_mean (EL OFFSET) / arcmin = 0.365159
    x_stddev / arcmin           = 2.094082
    y_stddev / arcmin           = 2.020382
    rotation angle / rad        = -0.054848
    x_fwhm (AZ FWHM) / arcmin   = 4.931186
    y_fwhm (EL FWHM) / arcmin   = 4.757636
    Goodness of fit: (chi-square, degree of freedom, upper probability) = (9.80, 20, 97.17%)
    Residual / Tsys: (mean, stdev, skew) = (-1.059103e-04, 8.251155e-04, 9.160260e-01)

Saved: /.../3C345_20191001T100000Z_gauss.pdf

Fitted AiryDisk2D Parameters:
    amplitude (= Ta / Tsys)       = 2.217251e-02
    x_0 (AZ OFFSET) / arcmin      = -0.077773
    y_0 (EL OFFSET) / arcmin      = 0.360073
    radius of first null / arcmin = 5.965962
    FWHM = 0.84 radius / arcmin   = 5.027056
    Goodness of fit: (chi-square, degree of freedom, upper probability) = (4.63, 22, 100.00%)
    Residual / Tsys: (mean, stdev, skew) = (2.517859e-05, 7.532696e-04, 1.151687e+00)

Saved: /.../3C345_20191001T100000Z_airydisk.pdf

出力文に書かれている通り、アンテナ温度 $T_\text{A}$ は \[ T_\text{A} = \textit{amplitude} \times T_\text{sys} \] によって得られる。その誤差は、とりあえずresidualの偏差 (stdev値) で与えても良いかもしれない: \[ T_\text{A, err} = \textit{residual-stdev} \times T_\text{sys}. \] フラックス $S$ は、観測対象が点源とみなせるならば、 \[ S = \frac{2 k_\text{B} T_\text{A}}{\eta_\text{a} A_\text{p}} \] として得られる。ここで $k_\text{B}$ はボツルマン定数 $1.38 \times 10^{-23}\: \text{J}/\text{K}$, $A_\text{p}$ は物理開口面積 (e.g., $3.14 \times 16^2 \: \text{m}^2$), $\eta_\text{a}$ は開口能率であり、その仰角依存性は2004年時の8 GHz帯のデータを使用すれば \[ \eta_\text{a} \times 100 = 57.543 + 0.815 \times \left(\frac{\textit{EL}}{\text{deg}}\right) - 0.014\: 042 \times \left(\frac{\textit{EL}}{\text{deg}}\right)^2 +0.000\: 064\: 673 \times \left(\frac{\textit{EL}}{\text{deg}}\right)^3 \] で与えられる。

より正確を期すなら、光学的厚み $\tau$ の大気による放射を加味しつつ大気による吸収を補正した (つまり大気の影響を加味しつつ宇宙空間にアンテナを置いたときの) システム雑音温度 $T_\text{sys}^* = e^\tau T_\text{sys}$ を用いて、宇宙空間においたアンテナで観測されるアンテナ温度 $T'_\text{A}$ を \[ T'_\text{A} = \textit{amplitude} \times T_\text{sys}^* \] で与える。さらにこれに各種効率を考慮すると、天体の放射輝度温度に対する most useful な推定値として \[ T_\text{R}^* = \frac{T'_\text{A}}{\eta_\text{r} \eta_\text{rss} \eta_\text{fss}} \] を得る。これはKutner & Ulich, 1981 の式 (10) であり、詳細は

で述べられ、それらを端的にまとめた資料がNRAOのページに上がっている。山口のR-Sky法ではそもそも大気込みシステム雑音温度 $T_\text{sys}^*$ を測っている。