- Maps to Science
- 原始惑星系円盤・原始星
- 星形成領域
- 銀河
- 近傍銀河: M100 CO 分子輝線データ
- 遠方銀河: BR1202 [CII] 輝線データ
- 活動銀河核
- (準備中)
銀河
遠方銀河「BR1202 [CII] 輝線データ」
ここでは作成されたイメージデータをもとにサイエンスを行うための例として、幾つかの典型的な解析を紹介します。Jupyter Notebook 版スクリプト
ここで紹介する解析について Jupyter Notebook 形式の Python スクリプトも用意しました。 Jupyter Notebook 形式のファイルは以下からダウンロードできます。ipynb ファイル (2024/02/09 改訂)
Jupyter Notebook 版スクリプトを実行するためには Jupyter Notebook がインストールされている必要があります。 また、Python への CASA モジュールのインストールも必要になります (tarball でインストールされた CASA には対応しません)。 Jupyter Notebook 版スクリプトは例えば Google Colaboratory (※) 上で実行する事も出来ます (Google Colab 環境かを識別して 自動的に CASA モジュールのインストールを行います)。
なお Jupyter Notebook 版スクリプトでは、以降で説明する解析とタスク
viewer を使用しない等の点で異なります。
※ Google Colaboratory (略称: Colab) は、Google の提供するブラウザから Python を記述、実行できるサービスです。 Jupyter Notebook 形式で書かれたスクリプトを読み込み、実行することが出来ます。
以降の解析は CASA 6.4 (Python 3) で動作確認しています。 使用している CASA 関係のコマンドは CASA 5 (Python 2) でも問題ないはずですが、一部の Python コマンドが Python 2 で問題が起こるかもしれません。
イメージデータ
ここで使用するデータは、ALMA Band 7 で観測された遠方銀河「BR1202」の [CII] 輝線キューブイメージデータです。Science Verification Data の一つとして取得されたデータ (「6. BR1202-0725 Band 7」を参照) から作成されたイメージを使用します。 また、初期解析については CASA ガイドも参考になります (CASA Guide：BR1202 SV Band7 Calibration notes)。
JVO のサイトからもイメージ FITS を直接ダウンロードできます。 輝線キューブイメージは
ALMA00000012_00_00_00.fitsから、
連続波イメージは
ALMA00000011_00_00_00.fitsからダウンロードして下さい。FITS イメージ表示
タスクviewer を使用します。
# In CASA
viewer
# In command line
% casaviewer &
viewer が起動します。
"Viewer Display Panel" の他に、"Data Manager" というウインドウが開きます。
ここで開きたいデータを選択して、"raster image" をクリックして下さい。
CASA 形式のイメージデータの他に、イメージ FITS ファイルも直接開くことが出来ます。
"Viewer Display Panel" に、選択したイメージデータの画像が表示されます。
右図は一例として、70 チャンネル目を表示させた BR1202 の 輝線イメージです。
(速度表示はイメージングの際に指定した周波数 (通常観測中心周波数) からの速度オフセットなので無視して下さい)
キューブイメージデータを扱う場合、周波数 (速度) 方向の情報を持ちます。
周波数方向のチャンネルを変える場合は以下の様に操作してください。
チャンネル移動 (ムービー)
チャンネル移動 (1ch のみ)
ムービーの停止
最小・最大チャンネルへの移動
をマウスでドラッグしてもチャンネルを変えることが出来ます。
"Viewer Display Panel" 上部のパネルの各アイコンをクリックすると各種操作ができます。
拡大・縮小
移動 (画面上でマウスをドラッグ)
カラートーン変更 (画面上でマウスをドラッグ)
こちらは連続波イメージになります。
なお、連続波イメージの場合、周波数方向には 1 チャンネルしかないのでチャンネル操作パネルは表示されません。
viewer の扱い方については、CASA docs の
"Viewing Images"
のページも参考にして下さい。
イメージ FITS 読み込み
# In CASA
# cube image for line
orgimfits='ALMA00000012_00_00_00.fits'
orgimgdat='BR1202_line_cube.image'
importfits(fitsimage=orgimfits,
imagename=orgimgdat, overwrite=True)
importfits でイメージ FITS を CASA 形式の画像ファイルに変換しています。
fitsimage には入力イメージ FITS 名を、imagename には出力される CASA 形式のファイル名 (実体はディレクトリ) を指定します。
# In CASA
# continuum image
orgimfitsct='ALMA00000011_00_00_00.fits'
orgimgdatct='BR1202_cont.image'
importfits(fitsimage=orgimfitsct,
imagename=orgimgdatct, overwrite=True)
ヘッダーから合成ビームサイズ等を確認
# In CASA
orgimgdat='BR1202_line_cube.image'
orgimgdatct='BR1202_cont.image'
for imname in [orgimgdat, orgimgdatct]:
bmaj =imhead(imname, mode='get', hdkey='beammajor')['value']
bmajut=imhead(imname, mode='get', hdkey='beammajor')['unit']
bmin =imhead(imname, mode='get', hdkey='beamminor')['value']
bminut=imhead(imname, mode='get', hdkey='beamminor')['unit']
bmpa =imhead(imname, mode='get', hdkey='beampa')['value']
bmpaut=imhead(imname, mode='get', hdkey='beampa')['unit']
if bmpaut=='rad': bmpa=bmpa*180/math.pi; bmpaut='deg'
print('Synthesized Beam of '+imname+':')
print(f' Major Axis = {bmaj:.7f} {bmajut}')
print(f' Minor Axis = {bmin:.7f} {bminut}')
print(f' Pos. Angle = {bmpa:.7f} {bmpaut}')
Synthesized Beam of BR1202_line_cube.image:
Major Axis = 1.7564977 arcsec
Minor Axis = 1.3982995 arcsec
Pos. Angle = 87.3220291 deg
Synthesized Beam of BR1202_cont.image:
Major Axis = 1.7226058 arcsec
Minor Axis = 1.3101563 arcsec
Pos. Angle = -86.6755066 deg
imhead を用いて、イメージデータのヘッダーから合成ビームサイズが確認できます。
# In CASA
orgimgdat ='BR1202_line_cube.image'
orgimgdatct='BR1202_cont.image'
imhd_A1 = imhead(orgimgdat)
print(imhd_A1)
imhd_B1 = imhead(orgimgdatct)
print(imhd_B1)
imhead を用いて Python の辞書オブジェクトとして読み込まれ、keyword を指定して呼び出すことが出来ます。
これによりイメージデータの各軸の座標情報・周波数情報なども確認できます。
# In CASA
import numpy as np
orgimgdat ='BR1202_line_cube.image'
imhd_A1 = imhead(orgimgdat)
print(imhd_A1['axisnames'])
print(imhd_A1['axisunits'])
print(imhd_A1['shape'])
pixinc=imhd_A1['incr'][0:2]
if imhd_A1['axisunits'][0]=='rad': pixinc=pixinc*180/np.pi
print('ref. freq. :', imhd_A1['refval'][2]/10**9, 'GHz')
print('pixel incr.:', pixinc*3600, 'arcsec/pixel')
print('ch width :', abs(imhd_A1['incr'][2])/10**6, 'MHz/pixel')
['Right Ascension' 'Declination' 'Frequency' 'Stokes']
['rad' 'rad' 'Hz' '']
[256 256 128 1]
ref. freq. : 334.9602625 GHz
pixel incr.: [-0.25 0.25] arcsec/pixel
ch width : 15.625 MHz/pixel
# In CASA
orgimgdatct='BR1202_cont.image'
imhd_B1 = imhead(orgimgdatct)
print(imhd_B1['axisnames'])
print(imhd_B1['axisunits'])
print(imhd_B1['shape'])
pixinc=imhd_B1['incr'][0:2]
if imhd_B1['axisunits'][0]=='rad': pixinc=pixinc*180/np.pi
print('ref. freq. :', imhd_B1['refval'][2]/10**9, 'GHz')
print('pixel incr.:', pixinc*3600, 'arcsec/pixel')
['Right Ascension' 'Declination' 'Frequency' 'Stokes']
['rad' 'rad' 'Hz' '']
[256 256 1 1]
ref. freq. : 341.9 GHz
pixel incr.: [-0.25 0.25] arcsec/pixel
これら情報より、このイメージデータは Right Ascension, Declination, Frequency, Stokes の 4 軸を持ち (ただし Stokes 軸は total intensity (Stokes I) しか持たない)、 空間方向はいずれも 256 × 256 [pixel]、キューブは周波数方向に 128 [pixel] あり、空間方向のピクセル間隔はいずれも 0.25 [arcsec]、キューブの周波数方向のピクセル間隔は 15.625 [MHz] である事が分かります。
イメージ RMS (及び輝線ピーク) を確認
タスクimstat で、パラメータ region で指定された空間・周波数範囲内のピクセル値から、イメージ RMS やピーク値などの統計量を計算できます。領域 "
region" の指定は、空間・周波数範囲を以下の様な設定が可能です。- 空間方向：
- "
centerbox" は中心位置と縦横方向の長さで四角形の領域を指定します。 - "
circle" は中心位置と半径で円形の領域を指定します。 - "
annulus" は中心位置と内縁・外縁の半径でリング状の領域を指定します。 - 周波数方向：
- "
range" は周波数 (速度) チャンネルの範囲 (始点, 終点) を指定します。
空間方向の領域指定方法：
# In CASA
# cube image
imname = 'BR1202_line_cube.image'
imhd_A = imhead(imname)
radipx = 100 # pixel
cetpx = [str(int((imhd_A['shape'][0]-1)/2)), str(int((imhd_A['shape'][1]-1)/2))]
ofpx = 10 ; wdpx = 10
nchan = imhd_A['shape'][2]
stats1=imstat(imagename=imname,
region='circle[['+cetpx[0]+'pix,'+cetpx[1]+'pix],'+str(radipx)+'pix],\
range=['+str(ofpx)+'pix,'+str(ofpx+wdpx-1)+'pix]')
stats2=imstat(imagename=imname,
region='circle[['+cetpx[0]+'pix,'+cetpx[1]+'pix],'+str(radipx)+'pix],\
range=['+str(nchan-ofpx-wdpx)+'pix,'+str(nchan-ofpx-1)+'pix]')
stats3=imstat(imagename=imname,
region='circle[['+cetpx[0]+'pix,'+cetpx[1]+'pix],'+str(radipx)+'pix],\
range=['+str(ofpx+wdpx)+'pix,'+str(nchan-ofpx-wdpx-1)+'pix]')
imgrms=(stats1['rms'][0]+stats2['rms'][0])/2.0
imgpek=stats3['max'][0]
print(f'RMS of {imname}: {imgrms*1000:.7f} [mJy/beam]')
print(f'Peak of {imname}: {imgpek*1000:.7f} [mJy/beam]')
RMS of BR1202_line_cube.image: 3.5801296 [mJy/beam]
Peak of BR1202_line_cube.image: 24.2306907 [mJy/beam]
region として、空間方向に "circle" でイメージ中心 (ヘッダーから参照) から半径 100 [pix] の円形領域内を指定し、
"range" で速度 (周波数) 方向の範囲をイメージ RMS では両バンド端から 10 [pix] 除いて 10 [pix] 幅 (10 [pix] ～ 19 [pix] 及び 108 [pix] ～ 117 [pix]) を、
ピーク値ではバンド端 20 [pix] を除いたその内側 (20 [pix] ～ 107 [pix]) を指定して算出しています。
# In CASA
# continuum image
imname = 'BR1202_cont.image'
imhd_B = imhead(imname)
inrdipx = 30 # pixel; inner radius
otrdipx = 100 # pixel; outer radius
cetpx = [str(int((imhd_B['shape'][0]-1)/2)), str(int((imhd_B['shape'][1]-1)/2))]
stats1=imstat(imagename=imname,
region='annulus[['+cetpx[0]+'pix,'+cetpx[1]+'pix],['+str(inrdipx)+'pix,'+str(otrdipx)+'pix]]')
stats2=imstat(imagename=imname,
region='circle[['+cetpx[0]+'pix,'+cetpx[1]+'pix],'+str(inrdipx)+'pix]')
imgrmsct=stats1['rms'][0]
imgpekct=stats2['max'][0]
print(f'RMS of {imname}: {imgrmsct*1000:.7f} [mJy/beam]')
print(f'Peak of {imname}: {imgpekct*1000:.7f} [mJy/beam]')
RMS of BR1202_cont.image: 0.2586689 [mJy/beam]
Peak of BR1202_cont.image: 10.6590521 [mJy/beam]
region として空間方向に、
イメージ RMS では "annulus" でイメージ中心 (ヘッダーから参照) から半径 30 [pix] ～ 半径 100 [pix] のリング状の空間領域を、
ピーク値では "circle" で半径 30 [pix] の円形領域を指定しています。
イメージ RMS (ノイズの大きさ) を測るときは、天体からの放射のない領域で測定する必要があります。 上記において、キューブイメージでは "
range" で輝線の検出されていない両端周辺の周波数チャンネルを指定し、連続波イメージでは中心付近の放射を避けて "annulus" でリング状の空間領域を指定して RMS を求めています。
パラメータ
region で指定される "領域" は、タスク imstat だけでなく様々なタスクで使用されます。
領域 "region" の指定方法について、詳しくは CASA docs の
"Region File Format"
をご覧ください。
求められた統計量は Python の辞書オブジェクト形式として読み込まれ、RMS の他に指定領域内のフラックス密度の最大値と最小値 (及びその位置)・合計・平均などの統計値も知る事が出来ます。
全ての統計量を表示したい場合は、先のコマンド実行後に (例えばキューブイメージのピーク値の場合)
# In CASA
print(stats3)
これらの統計量は
viewer 上でも、領域 "region" の指定に制限がありますが、簡単に確認することが出来ます。
まず、調べたい周波数 (速度) チャンネルを選んで表示します。
イメージ RMS を調べるのなら、輝線のないチャンネルで測定します (右図は 15 チャンネルの場合)。
円形領域選択ボタン (
の真ん中) をクリックし、画像上でマスクのドラッグによって region を指定できます。
詳細な領域設定は右側の "Regions" から "Properties" タブを開いて指定することもできます (中心位置・X 方向 Y 方向の幅を指定)。
領域を指定した後に、 "Regions" の "Statistics" タブを開くと、選択した領域内の RMS を含む統計量を確認することが出来ます。
天体の輝線が検出されている領域を取り囲む領域を指定すれば、ピークのフラックス密度などが分かります (右図は 70 チャンネルの場合)。
だたし、
viewer で統計量を求める場合、周波数 (速度) 方向は表示されている 1 チャンネルのみに対応し、空間方向の領域設定は四角形・楕円形 (いずれも回転なし) 領域と、ポリゴンでの領域指定 () のみが可能です。
サブイメージ作成
# In CASA
# cube image
orgimgdat='BR1202_line_cube.image'
myimages ='BR1202_line_cube.subimage'
imsubimage(imagename=orgimgdat, outfile=myimages, dropdeg=True,
region='centerbox[[127pix,127pix],[80pix,80pix]]')
# continuum image
orgimgdatct='BR1202_cont.image'
myimagesct ='BR1202_cont.subimage'
imsubimage(imagename=orgimgdatct, outfile=myimagesct, dropdeg=True,
region='centerbox[[127pix,127pix],[80pix,80pix]]')
imsubimage で天体成分のあるイメージ中心付近について切り出すことにします。
切り出す領域 "region" は、"centerbox" により四角形の中心位置及び縦横の幅をピクセル (pixel) で指定しています。
上記では、イメージ中心を中心とした 80 [pix] × 80 [pix] の範囲を切り出しています。
切り出されたイメージについて、70 チャンネル目を viewer で表示させた BR1202 の 輝線イメージです。
こちらは切り出された連続波イメージです。
スペクトルを表示
スペクトル線データから作成される 3 次元のキューブイメージにおいて、一部の領域内のピクセル値を空間方向に足し合わせる (又は平均する) 事で、周波数方向のスペクトルを作成する事が出来ます。ここではタスク
viewer を使って、キューブデータからスペクトルを表示します。viewer でキューブデータを開いた後、
アイコンのどれかを使って見たい領域 "region" を指定し、上部パネルの 
アイコンをクリックして "Spectral Profile Tool" を開いて下さい。
スペクトルの縦軸は Mean, Median, Sum, Flux Density、横軸は channel, radio velocity [km/s], frequency [GHz] などから選んで変更できます。
中心部の 35 [pix] × 35 [pix] の四角形を "region" として指定しています。この空間方向の範囲で合成されたスペクトルを表示します。
アイコンをクリックすると、この "Spectral Profile" ウインドウが開き、
上記で指定された "region" 内の平均スペクトル (Mean) が表示されます。
下部で横軸を frequency [GHz] に、縦軸を合成ビームサイズを考慮した Flux Density [mJy] を選択して表示したスペクトルです。
輝線が検出されているのは分かりますが、あまり顕著ではありません。 輝線が十分強くない事と、速度構造の異なる複数のコンパクトな成分の足し合わせのためスペクトルが平均化されている考えられます。 詳しくスペクトルを見るには、後述の様に空間成分ごとにスペクトルを作る事が必要になると考えられます。
モーメントマップ作成
スペクトル線データから作成されるイメージは、2 つの空間軸と 1 つの周波数 (又は速度) 軸の 3 次元のイメージキューブになり、これらの 3 つのイメージ平面を操作してモーメントマップが作成できます。# In CASA
imname = 'BR1202_line_cube.subimage'
imname_m0 = imname+'.m0'
imname_m1 = imname+'.m1'
imname_m2 = imname+'.m2'
imname_m8 = imname+'.m8'
imgrms = 3.5801296/10**3 # [Jy/beam] : RMS
cutrms = 3.0 # intensity cut off for moment1 & moment2
mochs = [25, 112] # channel ranges for moment maps
# moment0 & moment8
immoments(imagename = imname, outfile = imname_m0, moments = 0,
axis = 'spectral', chans=str(mochs[0])+'~'+str(mochs[1]))
immoments(imagename = imname, outfile = imname_m8, moments = 8,
axis = 'spectral', chans=str(mochs[0])+'~'+str(mochs[1]))
# moment1 & moment2
immoments(imagename = imname, outfile = imname_m1, moments = 1,
axis = 'spectral', chans=str(mochs[0])+'~'+str(mochs[1]),
includepix=[imgrms*cutrms, 9999.])
immoments(imagename = imname, outfile = imname_m2, moments = 2,
axis = 'spectral', chans=str(mochs[0])+'~'+str(mochs[1]),
includepix=[imgrms*cutrms, 9999.])
immoments で作成出来ます。
ここで作成したモーメントマップは、moment 0 は積分強度マップ、moment 1 は速度場マップ、moment 2 は速度分散マップ、moment 8 はピーク強度マップに相当します (immoments のパラメータ moments で種別指定)。
上記では "mochs"で指定した周波数チャンネル範囲 (25~112 チャンネル) でモーメントマップが作成されます (
immoments のパラメータ chans でチャンネル範囲指定)。
また、モーメント 1 及びモーメント 2 マップについてはノイズの影響を避けるため、ある閾値 "cutrms" σ (上記では 3.0 σ) 以上のデータのみを使って作成しています。 "imgrms" で先に測定したイメージ RMS を指定して下さい。
作成されたモーメントマップを viewer で表示するとこのようになります。
モーメントマップを複数開き でアニメーション表示することが出来ます。
モーメントマップには、他に以下の様な種類があります。 それぞれタスク
immoments でパラメータ moments に以下の値を指定すれば作ることが出来ます。
- moments=-1 - mean value of the spectrum
- moments=0 - integrated value of the spectrum
- moments=1 - intensity weighted coordinate; traditionally used to get 'velocity fields'
- moments=2 - intensity weighted dispersion of the coordinate; traditionally used to get 'velocity dispersion'
- moments=3 - median of I
- moments=4 - median coordinate
- moments=5 - standard deviation about the mean of the spectrum
- moments=6 - root mean square of the spectrum
- moments=7 - absolute mean deviation of the spectrum
- moments=8 - maximum value of the spectrum
- moments=9 - coordinate of the maximum value of the spectrum
- moments=10 - minimum value of the spectrum
- moments=11 - coordinate of the minimum value of the spectrum
- CASA Docs: "Spectral Analysis - Computing the Moments of an Image Cube"
- CASA Toolkit Reference Manual: "image.moments - Function"
連続波イメージとの比較
viewer を使って、モーメント 0 マップ (積分強度図) を連続波マップと重ねて表示することも出来ます。
それぞれのイメージは原則、同じ cell size・image size・イメージ中心 (phasecenter) で作られている事を想定しています
(今回のイメージデータでそれぞれ同じであることはヘッダー情報から確認)。
"Data Manager" ウインドウで、ラインのモーメント 0 マップを選んで "raster image" をクリック、連続波イメージを選んで "contour map" をクリックします。
この様に、モーメント 0 マップのカラー画像に、連続波マップをコントアイメージが重ねて表示されます。
なお、右のイメージではコントア間隔等は調整されています。
コントア表示の細かな調整は、上部パネルの 
アイコンをクリックして開かれる下図の様な "Data Display Options" で行う事が出来ます。
こちらが "Data Display Options" ウインドウです。
上部の対象イメージのタブに切り替えて、"basic settings" から"Relative Contour Levels", "Base Contour Level", "Unit Contour Level" を修正してコントアイメージを調節する事が出来ます。
例えば、以下の様に入力します。
- Relative Contour Levels : [1,3,5,7,9,11,13,15]
- Base Contour Level : 0
- Unit Contour Level : 0.0006467 (~ rms x 2.5)
ピーク位置と輝線フラックスの算出
モーメント 0 マップ (及び連続波イメージ) から、2 つの放射成分がある事が分かります。 そこで、下記の様にviewer で表示されたモーメント 0 マップ上の 2 つの赤円の領域に対して2次元ガウシアンフィットを行い、各領域における輝線フラックス・ピーク値・ピークの位置を算出出来ます。
モーメント 0 マップを表示して、それぞれの放射成分に対して領域を設定します。
試したい領域を選択した状態で、右側の "Regions" から "Fit" タブを開いて下部にある "gaussfit" ボタンをクリックします。
成分 A (左下) への 2 次元ガウシアンフィットの結果が、右側の "Fit" タブに表示されます。
中心位置の座標・長軸短軸サイズと傾き・ピーク値と積分したフラックス等が表示されています。
また、イメージ中には緑色で、フィティングされたガウシアンの中心位置と長軸方向の向きが表示されます。
こちらは成分 B (右上) の結果です。指定した領域にマウスを持ってくと領域が切り替わり、右側の結果の表示も切り替わります。
また、2 次元ガウシアンフィットを行うと CASA (または
viewer) を開いたターミナルに、フィッティングの結果の詳細が出力されます。
下記は、その結果の抜粋になります (左側の "...." には実行された時刻とメッセージの種類が表記されます)。こちらは成分 A (左下) の結果、
.... Region::getLayerCenter Centering results:
.... Region::getLayerCenter --- Ra_J2000: 12:05:23.127
.... Region::getLayerCenter --- Dec_J2000: - 7.42.33.04
.... Region::getLayerCenter --- W-Majorax: 2.02805 arcsec
.... Region::getLayerCenter --- W-Minorax: 1.61319 arcsec
.... Region::getLayerCenter --- W-Posang: 92.7321 deg
.... Region::getLayerCenter --- Xcen: 40.5293 pix
.... Region::getLayerCenter --- Ycen: 38.4568 pix
.... Region::getLayerCenter --- I-Majorax: 8.11218 pix
.... Region::getLayerCenter --- I-Minorax: 6.45276 pix
.... Region::getLayerCenter --- I-Posang: 2.73208 deg
.... Region::getLayerCenter --- Radeg: 181.346 deg
.... Region::getLayerCenter --- Decdeg: -7.70918 deg
.... Region::getLayerCenter --- IntegrFlux: 11.9647 Jy.km/s
.... Region::getLayerCenter --- PeakFlux: 8.98227 Jy/beam.km/s
.... Region::getLayerCenter --- Skylevel: -1.02083
.... Region::getLayerCenter --- Converged: YES
.... Region::getLayerCenter Centering results:
.... Region::getLayerCenter --- Ra_J2000: 12:05:22.996
.... Region::getLayerCenter --- Dec_J2000: - 7.42.29.84
.... Region::getLayerCenter --- W-Majorax: 2.0747 arcsec
.... Region::getLayerCenter --- W-Minorax: 1.61217 arcsec
.... Region::getLayerCenter --- W-Posang: 109.058 deg
.... Region::getLayerCenter --- Xcen: 48.3114 pix
.... Region::getLayerCenter --- Ycen: 51.2268 pix
.... Region::getLayerCenter --- I-Majorax: 8.29879 pix
.... Region::getLayerCenter --- I-Minorax: 6.44867 pix
.... Region::getLayerCenter --- I-Posang: 19.058 deg
.... Region::getLayerCenter --- Radeg: 181.346 deg
.... Region::getLayerCenter --- Decdeg: -7.70829 deg
.... Region::getLayerCenter --- IntegrFlux: 13.0052 Jy.km/s
.... Region::getLayerCenter --- PeakFlux: 9.54993 Jy/beam.km/s
.... Region::getLayerCenter --- Skylevel: -0.637841
.... Region::getLayerCenter --- Converged: YES| 中心位置 [RA, Dec] | 中心位置 [RA, Dec](pixel) | PeakFlux [Jy/beam.km/s] | IntegrFlux [Jy.km/s] | |
|---|---|---|---|---|
| 成分 A | 12:05:23.127, -07.42.33.04 | 40.5293, 38.4568 | 8.98227 | 11.9647 |
| 成分 B | 12:05:22.996, -07.42.29.84 | 48.3114, 51.2268 | 9.54993 | 13.0052 |
イメージへの 2 次元ガウシアンフィットは、タスク
imfit を使って行う事も出来ます
(Jupyter Notebook 版スクリプトをご覧ください)。
輝線の赤方偏移と速度構造
前述の様にモーメント 0 マップ上 (連続波イメージ上でも) で 2 つの放射成分あることが分かります。 イメージ上では明確に 2 つの成分を切り分ける事ができ、それぞれ独立した「銀河」に由来する事が期待されます。 そこで、2 つの銀河の成分を別々に輝線スペクトルを作成して、それぞれの周波数又は速度情報の算出を行います。
viewer でキューブイメージを開いた後、それぞれの放射成分に対して領域を設定します (2 次元ガウシアンフィットの時と同様)。
試したい領域を選択した状態で、上部パネルの アイコンをクリックして "Spectral Profile" ウインドウを開いて下さい。
成分 A (左下) のスペクトルです。
先に説明した通り、縦軸・横軸は変更可能です (ここでは、Frequency [GHz] vs FLux Density [mJy])。
こちらは、成分 B (右上) のスペクトルになります。
指定した領域にマウスを持ってくとスペクトルも切り替わります (少し切り替わるのに時間がかかる事があります)。
viewer を使わずに、タスク specflux を使用してスペクトルのデータを作成する事も出来ます。
specflux ではスペクトルデータを、パラメータ logfile で指定したファイルにテキストデータとして書き出します。
下記は specflux で作成したスペクトルデータを読み込んで、Python (matplotlib モジュール; CASA 上で使用可能) を使って表示するスクリプトです。
# In CASA
import numpy as np
from matplotlib import pyplot as plt
imname ='BR1202_line_cube.subimage'
radipx = 7 # pixels
# Comp.A (South)
dtext1 = "Comp.A"
splogf1 = imname+'.specA.log'
cetpx1 = [ 40.5293, 38.4568 ] # pixs; peak position
specflux(imagename=imname,
region = 'circle[['+str(cetpx1[0])+'pix,'+str(cetpx1[1])+'pix],'+str(radipx)+'pix]',
function='flux density', unit='km/s',
logfile = splogf1, overwrite=True)
dt1_ax,dt1_npx,dt1_frq,dt1_vl,dt1_fl = np.loadtxt(splogf1, comments='#', unpack=True)
# Comp.B (North)
dtext2 = "Comp.B"
splogf2 = imname+'.specB.log'
cetpx2 = [ 48.3114, 51.2268 ] # pixs; peak position
specflux(imagename=imname,
region = 'circle[['+str(cetpx2[0])+'pix,'+str(cetpx2[1])+'pix],'+str(radipx)+'pix]',
function='flux density', unit='km/s',
logfile = splogf2, overwrite=True)
dt2_ax,dt2_npx,dt2_frq,dt2_vl,dt2_fl = np.loadtxt(splogf2, comments='#', unpack=True)
dt_flcb = np.r_[dt1_fl, dt2_fl]
ymin = min(dt_flcb); ymax = max(dt_flcb)
xmin = min(dt1_frq); xmax = max(dt1_frq)
ff, aa = plt.subplots(1,1, figsize=(8,5))
imr1 = aa.step(dt1_frq, dt1_fl, color="r", label=dtext1, where="mid")
imr2 = aa.step(dt2_frq, dt2_fl, color="b", label=dtext2, where="mid")
aa.set_xlim(xmin, xmax) ; aa.set_xlabel("Frequency [MHz]")
aa.set_ylim(ymin*1.2, ymax*1.2) ; aa.set_ylabel("Flux density [Jy]")
plt.hlines([0], xmin, xmax, "blue", linestyles='dashed')
aa.legend(loc="upper left")
plt.savefig('BR1202_line_spectrum2.png', bbox_inches="tight", pad_inches=0.05)
plt.show()
region は、各成分毎に "cetpx1", "cetpx2" で指定される中心位置ピクセルを中心に、半径 "radipx" ピクセルで指定される円形で定義されています。
中心位置ピクセルについては、モーメント 0 マップに対する 2 次元ガウシアンフィットで求めた各成分の中心位置 (pixel) を参考にして下さい。
上記スクリプトで作成されたプロットです。
以上で作成された成分毎のスペクトルに対して 1 次元ガウシアンフィットを行い、各成分の周波数 (速度) 情報を算出する事が出来ます。
viewer から開かれるスペクトルを表示している "Spectral Profile" ウインドウにおいて、ウインドウ下部の "Spectral-Line Fitting" タブにある "fit" ボタンをクリックすると、スペクトルに対する 1 次元ガウシアンフィットが行われます。
成分 A (左下) スペクトルへのガウシアンフィットの結果です。
フィッティングの際、"Spectral-Line Fitting" タブにある "from:", "to: " にフィッティングを行う周波数範囲を指定して下さい。
輝線が一つの場合はデフォルト (全域指定) でも問題ないはずです。
右側で、フィッティング関数を指定します。単一ガウス関数の場合は、"gauss", "no poly" を選択してください。
"poly" は線形関数を意味します。連続波成分などが含まれると予想される場合は "poly 1" を選ぶ事も出来ます。
成分 B (右上) スペクトルへのガウシアンフィットの結果です。
また、スペクトルへのガウシアンフィットを行うと CASA (または
viewer) を開いたターミナルに、
下記の様なフィッティングの結果の詳細が出力されます (左側の "...." には実行された時刻とメッセージの種類が表記されます)。
こちらは成分 A (左下) スペクトルの結果、
.... SpecFitOptical::specLineFit Successful fit!
.... SpecFitOptical::specLineFit No. of iterations: 18
.... SpecFitOptical::specLineFit Chi-square: 2250.41
.... SpecFitOptical::specLineFit Gauss amplitude: 20.2042+-1.18062 [mJy]
.... SpecFitOptical::specLineFit Gauss centre: 333.775+-0.0117115 [GHz]
.... SpecFitOptical::specLineFit Gauss FWHM: 0.408726+-0.0275785 [GHz]
.... SpecFitOptical::specLineFit Gaussian area: 8.79033+-0.554814 [mJy*GHz]
.... SpecFitOptical::specLineFit Successful fit!
.... SpecFitOptical::specLineFit No. of iterations: 22
.... SpecFitOptical::specLineFit Chi-square: 3602.27
.... SpecFitOptical::specLineFit Gauss amplitude: 13.5433+-1.10258 [mJy]
.... SpecFitOptical::specLineFit Gauss centre: 333.987+-0.0299612 [GHz]
.... SpecFitOptical::specLineFit Gauss FWHM: 0.750683+-0.0706369 [GHz]
.... SpecFitOptical::specLineFit Gaussian area: 10.8221+-0.952163 [mJy*GHz]
| ピーク強度 [mJy] | 中心周波数 [GHz] | 周波数幅 [GHz] | 成分A | 20.2042+-1.18062 | 333.775+-0.0117115 | 0.408726+-0.0275785 |
|---|---|---|---|
| 成分B | 13.5433+-1.10258 | 333.987+-0.0299612 | 0.750683+-0.0706369 |
これより輝線の中心周波数の赤方偏移 (redshift; z) 及び輝線の速度幅 (半値全幅: FWHM) を以下の様に算出できます。
# In CASA
popt1v = [ 20.2042, 333.775, 0.408726 ] # [mJy , GHz, GHz]
popt2v = [ 13.5433, 333.987, 0.750683 ]
ccc = 299792458*10**-3 # [km/s] : speed of light
rfrq = 1900.539 # [GHz] : rest frequency for [CII] line
print("# Compornent A")
print(f'Redshift (z) : {rfrq/popt1v[1]-1:.7f}')
print(f'Line Width(FWHM): {ccc/popt1v[1]*popt1v[2]:.7f} km/s')
print("# Compornent B")
print(f'Redshift (z) : {rfrq/popt2v[1]-1:.7f}')
print(f'Line Width(FWHM): {ccc/popt2v[1]*popt2v[2]:.7f} km/s')
# Compornent A
Redshift (z) : 4.6940724
Line Width(FWHM): 367.1124925 km/s
# Compornent B
Redshift (z) : 4.6904580
Line Width(FWHM): 673.8259326 km/s
また、ガウシアンフィットで求められた輝線周波数幅 (半値全幅: FWHM) は、視線速度の速度幅として表せます。 周波数幅 (Δνobs) から電波観測における速度幅 (ΔVradio) は以下の様に計算しています。 νref にはガウシアンフィットから求められた輝線中心周波数を採用します。 (マイナスは周波数と速度で正負の向きが逆であることを示していて、幅がマイナスになるわけではありません)
上記の結果をまとめると、
| 中心周波数 [GHz] | 赤方偏移 (z) | 速度幅 [km/s] | |
|---|---|---|---|
| 成分 A | 333.775 | 4.69407 | 367.1125 |
| 成分 B | 333.987 | 4.69046 | 673.8259 |
それぞれの銀河の成分の輝線スペクトルとフィッティング結果 (速度構造) が分かりやすい様に、2 つの成分別々に輝線スペクトル及びフィッティングカーブを重ねて、中心周波数からのオフセットとして横軸を速度で表示します。
# In CASA
import numpy as np
from matplotlib import pyplot as plt
ccc = 299792458*10**-3 # [km/s] : speed of light
def gaussian_func(x, A, mu, sigma):
C = 0.0
return A * np.exp( - (x - mu)**2 / (2 * sigma**2)) + C
# Comp.A (South)
splogf1 = 'BR1202_line_cube.subimage.specA.log'
dt1_ax,dt1_npx,dt1_frq,dt1_vl,dt1_fl = np.loadtxt(splogf1, comments='#', unpack=True)
# Comp.B (North)
splogf2 = 'BR1202_line_cube.subimage.specB.log'
dt2_ax,dt2_npx,dt2_frq,dt2_vl,dt2_fl = np.loadtxt(splogf2, comments='#', unpack=True)
popt1v = [ 20.2042, 333.775, 0.408726 ] # [mJy , GHz, GHz]
popt2v = [ 13.5433, 333.987, 0.750683 ]
popt1 = [ popt1v[0]/(10**3), popt1v[1]*(10**3), popt1v[2]*(10**3)/(2*np.sqrt(2*np.log(2))) ] # [Jy , MHz, MHz]
popt2 = [ popt2v[0]/(10**3), popt2v[1]*(10**3), popt2v[2]*(10**3)/(2*np.sqrt(2*np.log(2))) ]
fig,aa = plt.subplots(nrows=1, ncols=2, figsize=(10,4.5), sharex=False)
xd1 = np.arange(dt1_frq.min(), dt1_frq.max(), 0.01)
xd2 = np.arange(dt2_frq.min(), dt2_frq.max(), 0.01)
estimated_curve1 = gaussian_func(xd1, popt1[0], popt1[1], popt1[2])
estimated_curve2 = gaussian_func(xd2, popt2[0], popt2[1], popt2[2])
vof1 = (popt1[1]-dt1_frq)/popt1[1]*ccc
vec1 = (popt1[1]-xd1)/popt1[1]*ccc
aa[0].step(vof1, dt1_fl, color="r")
aa[0].plot(vec1, estimated_curve1, color="black", linestyle='dashed')
vof2 = (popt2[1]-dt2_frq)/popt2[1]*ccc
vec2 = (popt2[1]-xd2)/popt2[1]*ccc
aa[1].step(vof2, dt2_fl, color="b")
aa[1].plot(vec2, estimated_curve2, color="black", linestyle='dashdot')
aa[0].set_title("Component A",size=12)
aa[1].set_title("Component B",size=12)
aa[0].set_xlabel("Velocity offset [km/s]",size=10)
aa[1].set_xlabel("Velocity offset [km/s]",size=10)
aa[0].set_ylabel("Flux density [Jy]",size=10)
aa[0].hlines([0], min(vof1), max(vof1), "blue", linestyles='dashed')
aa[1].hlines([0], min(vof1), max(vof1), "blue", linestyles='dashed')
plt.tight_layout()
plt.savefig('BR1202_line_spectrum4_b.png', bbox_inches="tight", pad_inches=0.05)
plt.show()
それぞれの輝線スペクトル及びフィッティングカーブです。
横軸は中心周波数からの速度オフセットになっています。
なお Jupyter Notebook 版では、
scipy モジュールで CASA とは独立にガウシアンフィットを行っています。
ページ先頭に戻る
Last Update: 2024.02.15