原始惑星系円盤・原始星

原始惑星系円盤「輝線キューブイメージからガス円盤の回転を知る」

ここでは作成されたイメージデータをもとにサイエンスを行うための例として、幾つかの典型的な解析を紹介します。


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.1 (Python 3) で動作確認しています。 使用している CASA 関係のコマンドは CASA 5 (Python 2) でも問題ないはずですが、一部の Python コマンドが Python 2 で問題が起こるかもしれません。

イメージデータ

今回使用するデータは、ALMA で観測された原始惑星系円盤を持つ「HL Tau」の H2CO 分子輝線キューブイメージデータです。
ALMA Science Archive より、以下を検索してください。
  • Project code: 2018.1.01037.S
  • Source Name: HL_Tau
JVO のサイトからも イメージ FITS ALMA01424659_00_00_00.fits を直接ダウンロードできます。

FITS イメージ表示

タスク viewer を使用します。 
# In CASA
viewer
または、 
# In command line
% casaviewer &
で、viewer が起動します。

"Viewer Display Panel" の他に、"Data Manager" というウインドウが開きます。 ここで開きたいデータを選択します。 CASA 形式のイメージデータの他に、イメージ FITS ファイルも直接開くことが出来ます。

"Viewer Display Panel" に、選択したイメージデータの画像が表示されます。

ここで扱うイメージデータはキューブデータですので、周波数 (速度) 方向の情報を持ちます。 周波数方向のチャンネルを変える場合は以下の様に操作してください。
  • チャンネル移動 (ムービー)
  • チャンネル移動 (1ch のみ)
  • ムービーの停止
  • 最小・最大チャンネルへの移動
スライダー をマウスでドラッグしてもチャンネルを変えることが出来ます。

"Viewer Display Panel" 上部のパネルの各アイコンをクリックすると各種操作ができます。
  • 拡大・縮小
  • 移動 (画面上でマウスをドラッグ)
  • カラートーン変更 (画面上でマウスをドラッグ)


右図は例として、249 チャンネル (6.11555 km/s) を拡大して表示させた HL Tau の H2CO 輝線イメージです。

viewer の扱い方については、CASA docs の "Viewing Images" のページも参考にして下さい。

イメージ FITS 読み込み

# In CASA
inimfits= 'ALMA01424659_00_00_00.fits'
imname  = 'HL_Tau.H2CO.cube.pbcor'
importfits(fitsimage=inimfits, 
           imagename=imname, overwrite=True)
タスク importfits でイメージ FITS を CASA 形式の画像ファイルに変換しています。

ヘッダーから合成ビームサイズを確認

# In CASA
imname = 'HL_Tau.H2CO.cube.pbcor'
imhd_A_res = imhead(imname)
bmaj  =imhd_A_res['restoringbeam']['major']['value']
bmajut=imhd_A_res['restoringbeam']['major']['unit']
bmin  =imhd_A_res['restoringbeam']['minor']['value']
bminut=imhd_A_res['restoringbeam']['minor']['unit']
bmpa  =imhd_A_res['restoringbeam']['positionangle']['value']
bmpaut=imhd_A_res['restoringbeam']['positionangle']['unit']
print('Synthesized Beam:')
print(f'   Major Axis = {bmaj:.7f} {bmajut}')
print(f'   Minor Axis = {bmin:.7f} {bminut}')
print(f'   Pos. Angle = {bmpa:.7f} {bmpaut}')

Synthesized Beam:
   Major Axis = 0.3242316 arcsec
   Minor Axis = 0.2745839 arcsec
   Pos. Angle = -2.3669729 deg
イメージデータのヘッダから合成ビームサイズが確認できます。 タスク imhead を使用し、Python の辞書オブジェクトとして様々な情報が読み込まれます。 これにより、ビームサイズの他に、イメージデータの各軸の座標情報・周波数情報なども確認できます。

平均スペクトルを表示

タスク viewer を使って、キューブデータから平均スペクトルを作成することが出来ます。
viewer でキューブデータを開いた後、 アイコンのどれかを使って視野中心方向に領域 (Region) を指定し、上部パネルの アイコンをクリックして "Spectral Profile Tool" を開いて下さい。 縦軸は Mean, Median, Sum, Flux Density、横軸は channel, radio velocity [km/s], frequency [GHz] などから選んで変更できます。

視野中心より半径 9" の円形内の平均スペクトルです。ダブルピークの輝線が検出されていることが分かります。

ラインのある範囲とない範囲、ラインの中心速度 (システム速度 (原始星系の速度) と同等と予測される) がどのあたりになるかに注意してください。

イメージ RMS (及び輝線ピーク) を確認

# In CASA
imname = 'HL_Tau.H2CO.cube.pbcor'
ctx=400 ; cty=400 # center pixs
stats1=imstat(imagename=imname, 
          region='circle[['+str(ctx)+'pix,'+str(cty)+'pix],9.0arcsec],range=[30pix,180pix]')
stats2=imstat(imagename=imname, 
          region='circle[['+str(ctx)+'pix,'+str(cty)+'pix],9.0arcsec],range=[300pix,450pix]')
imgrms=(stats1['rms'][0]+stats2['rms'][0])/2
print(f'RMS [mJy/beam] of {imname} : {imgrms*1000:.7f}')
stats3=imstat(imagename=imname, 
          region='circle[['+str(ctx)+'pix,'+str(cty)+'pix],4.5arcsec],range=[200pix,280pix]')
imgmax=stats3['max'][0]
print(f'Max [mJy/beam] of {imname} : {imgmax*1000:.7f}')

RMS [mJy/beam] of HL_Tau.H2CO.cube.pbcor : 1.8930310
Max [mJy/beam] of HL_Tau.H2CO.cube.pbcor : 53.1570651
タスク imstat で、パラメータ region で指定された範囲内のイメージ RMS やピーク値などの統計量を計算しています。
領域 "region" として、
  • "circle" は中心位置と半径を設定して円形領域を指定します。
  • "range" は周波数チャンネルの範囲 (始点, 終点) を指定します。
上記のスクリプトの場合、"circle" で視野中心 ((x,y) = (400pix,400pix)) から半径 9" (RMS) または半径 4.5" (Max) の円形領域を、"range" で速度 (周波数) 方向の範囲をピクセルで指定しています。

region の指定方法について、詳しくは CASA docs の "Region File Format" をご覧ください。

イメージ RMS (ノイズの大きさ) を測るときは、天体からの放射のない領域で測定する様にします。 上記では、"range" で輝線の検出されていない周波数チャンネルを指定して RMS を求めています。

求められた統計量は Python の辞書オブジェクト形式として読み込まれ、RMS の他にフラックス密度の最大値 & 最小値 (及びその位置)・合計・平均なども知る事が出来ます。
全ての統計量を表示したい場合は、先のコマンド実行後に 
# In CASA
print(stats1)
としてください。

これらの統計量は viewer 上でも、region の指定に制限がありますが、簡単に確認することが出来ます。

まず、調べたい周波数 (速度) チャンネルを選んで表示します。 イメージ RMS を調べるのなら、輝線のないチャンネルで測定します。 楕円形領域選択ボタン ( の中央) をクリックし、視野のなるべく広い領域を画像上でマウスをドラッグさせる事によって region を指定します。 詳細な領域設定は右側の "Regions" から "Properties" タブを開いて指定することもできます (中心位置・X 方向 Y 方向の直径を指定)。

region を指定した後、"Statistics" タブを開くと、選択した領域内の RMS を含む統計量を確認することが出来ます。

天体からの輝線放射を取り囲む region を選択すれば、ピークのフラックス密度などが分かります。

だたし、viewer で総計量を求める場合、周波数 (速度) チャンネルは表示されている 1 channel のみに対応し、領域設定は四角形・楕円形 (ただしいずれも回転はなし) 領域と、ポリゴンでの領域指定のみが可能です。

モーメントマップ作成

スペクトル線データから作成されるイメージは、2 つの空間軸と 1 つの周波数 (又は速度) 軸の 3 次元のイメージキューブになり、これらの 3 つのイメージ平面を操作してモーメントマップが作成できます。 
# In CASA
imname = 'HL_Tau.H2CO.cube.pbcor'
imname_m0 = imname+'.m0'
imname_m1 = imname+'.m1'
imname_m2 = imname+'.m2'
imname_m8 = imname+'.m8'
imgrms=1.8930310/10**3      # Jy/beam  RMS
ctx=402 ; cty=391 ; mwd=180 # center pixs & image size
os.system('rm -rf '+imname+'.m*')
# moment0 & moment8
immoments(imagename = imname, outfile = imname_m0,
       axis = 'spectral', moments = 0,
       chans = '200~280',
       region = 'centerbox[['+str(ctx)+'pix,'+str(cty)+'pix],['+str(mwd)+'pix,'+str(mwd)+'pix]]')
immoments(imagename = imname, outfile = imname_m8,
       axis = 'spectral', moments = 8,
       chans = '200~280',
       region = 'centerbox[['+str(ctx)+'pix,'+str(cty)+'pix],['+str(mwd)+'pix,'+str(mwd)+'pix]]')
# moment1 & moment2
immoments(imagename = imname, outfile = imname_m1,
       axis = 'spectral', moments = 1,
       chans = '200~280', includepix=[imgrms*5.0, 9999.],
       region = 'centerbox[['+str(ctx)+'pix,'+str(cty)+'pix],['+str(mwd)+'pix,'+str(mwd)+'pix]]')
immoments(imagename = imname, outfile = imname_m2,
       axis = 'spectral', moments = 2,
       chans = '200~280', includepix=[imgrms*5.0, 9999.],
       region = 'centerbox[['+str(ctx)+'pix,'+str(cty)+'pix],['+str(mwd)+'pix,'+str(mwd)+'pix]]')
ここで作成したモーメントマップは、moment 0 は積分強度マップ、moment 1 は速度場マップ、moment 2 は速度分散マップ、moment 8 はピーク強度マップになります。
なお、ここではモーメントマップを作成する際、円盤の中心付近 (中心位置を (x,y) = (402pix,391pix) とする) の領域のみを切り出しています。 ctx, cty, mwd はそれぞれ円盤中心位置と切り出すイメージサイズをピクセルで指定しています。 また、モーメント 1 及びモーメント 2 マップについてはノイズの影響を避けるため 5σ 以上のみのデータを使って作成しています。 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
モーメントマップについて、詳しくは以下を参考にして下さい。

連続波イメージとの比較

viewer を使って、モーメント 0 マップ (積分強度図) を連続波マップと重ねて表示することも出来ます。 それぞれのイメージは原則、同じ cell size・image size・イメージ中心 (phasecenter) で作られている事を想定しています。 連続波マップは「連続波イメージよりダスト円盤の物理量を導出する」で作った CASA 形式のイメージデータ (HL_Tau.cont.pbcor) を使います。

モーメント 0 マップのカラー画像に、連続波マップをコントアイメージで重ねたものです

この図より、ダスト円盤とガス円盤の位置がだいたい一致していることが分かります。 また、ここでは若干ガス円盤の方が広がっている様にも見えます。

ガス円盤の見た目の位置角の推定

大雑把な見積もりで良ければ、viewer の "物差しツール" ( アイコン) を使って、先に作ったモーメント 0 マップの長軸の X 方向の長さ・Y 方向の長さが得られますので、それらの値から長軸の傾きを推定できます。

また、赤方偏移成分と青方偏移成分について別々にモーメント 0 マップ (積分強度図) を作成し、それらの位置のずれから長軸の傾きを推定することも出来ます。
赤方偏移成分・青方偏移成分のモーメント 0 マップは、以下の様に作成します。 
# In CASA
imname = 'HL_Tau.H2CO.cube.pbcor'
imname_m0a = imname+'.m0_red'
imname_m0b = imname+'.m0_blue'
ctx=402 ; cty=391 ; mwd=180 # center pixs & image size
sysch=243 # channel of system velocity
offch=13  # offset channel for red/blue shift components
os.system('rm -rf '+imname+'.m0_red '+imname+'.m0_blue')
immoments(imagename = imname, outfile = imname_m0a,
      axis = 'spectral', moments = 0, chans = '200~'+str(sysch-offch),
      region = 'centerbox[['+str(ctx)+'pix,'+str(cty)+'pix],['+str(mwd)+'pix,'+str(mwd)+'pix]]')
immoments(imagename = imname, outfile = imname_m0b,
      axis = 'spectral', moments = 0, chans = str(sysch+offch)+'~280',
      region = 'centerbox[['+str(ctx)+'pix,'+str(cty)+'pix],['+str(mwd)+'pix,'+str(mwd)+'pix]]')
上記では、システム速度のピクセルを sysch、システム速度の前後のオフセットを offch で指定し、そのプラス側・マイナス側を赤方・青方偏移成分としてモーメントマップを作成しています。 ここでは sysch=243, offch=13 としていますが、赤方偏移成分として 200~230 channel を、青方偏移成分として 256~280 channel を足し合わせてモーメント 0 マップ (積分強度図) を作っている事になります。 これにより HL_Tau.H2CO.cube.pbcor.m0_red, HL_Tau.H2CO.cube.pbcor.m0_blue の 2 つのマップが作成されます。

赤方偏移・青方偏移の 2 つのモーメント 0 マップの中心位置を、viewer によるガウシアンフィットを使って以下の手順で測定して下さい。
  1. イメージを読み込む
  2. 天体が分かる様に拡大表示・位置調整
  3. 対象の天体を取り囲む様にregionを指定
  4. 右下の "Regions" から "Fit" タブを開き "gaussfit" をクリック
  5. 中心位置の pixel 値 (Xcen, Ycen) を確認する
それぞれのガウシアンフィットの結果からピーク位置を測定すると、赤方偏移成分のピークは (x,y) = [97.5716, 97.4801] pixels、青方偏移成分のピークは (x,y) = [81.3155, 81.0002] pixels と分かります。 したがって、それぞれの中心位置のずれから球面三角関数 () より、長軸の位置角 PA (Position Angle) が以下の様に推定できます。 
# In CASA
import math
incrx = 0.054 # arcsec : X 方向の pixel 間隔 (イメージング時の cell size) 
incry = incrx # arcsec : Y 方向の pixel 間隔
cnt_rd = [97.5716, 97.4801] # pixel : 赤方偏移成分のピーク位置
cnt_bl = [81.3155, 81.0002] # pixel : 青方偏移成分のピーク位置
dxx = cnt_rd[0] - cnt_bl[0] # pixel
dyy = cnt_rd[1] - cnt_bl[1] # pixel
#
gdkpa =  90. + math.degrees(math.atan2( math.sin(math.radians(dxx*incrx)), math.tan(math.radians(dyy*incry)) ))
print(f'Gas disk PA (position angle) = {gdkpa:.5f} deg')

Gas disk PA (position angle) = 134.60488 deg
また、viewer を用いて、前述の全速度範囲のモーメント 0 マップをカラーで、今回の赤方/青方偏移成分のモーメント 0 マップをコントアで重ねるとこのようになります。

長軸方向に円盤の両側で、右上が赤方偏移成分 (システム速度からプラス速度側)、左下が青方偏移成分 (システム速度からマイナス速度側) が分布していることが分かります。 この事から、円盤が短軸方向に近い軸を中心に回転運動をしていることが推測できます。

PV マップ作成

イメージキューブデータは、空間方向 XY の 2 成分と周波数 (速度) 成分の 3 次元の情報を含みます。 通常、viewer 等を使ってキューブデータを表示させると空間方向 2 成分が表示されます。 一方、任意の空間方向と周波数 (速度) 方向の強度分布を、位置速度図 (PVマップ; Position-Velocity map) と言います。 円盤の中心を基準に長軸・短軸の方向にカットしたPVマップから、ガス円盤のダイナミクスが分かります。

PV マップは、viewer を使って以下の様に簡単に作成できます。
キューブデータを表示した後、"P/V ツール" ( アイコン) をクリックし、PVマップを作成したい方向に始点から終点までマウスをドラッグして下さい。 右下の "Region" に "pV" タブが開きますので、"Generate P/V" をクリックすると、新しいウインドウが開き PV マップが表示されます。 また、この "pV" タブで、PV マップを作るためのパラメータを任意に指定することも出来ます。 averaging width は、スライスする方向と垂直な方向へ平均する幅をピクセル単位で表します。 この値を大きくすると、PV マップの SN が良くなるはずです。

まず、長軸方向にスライスした PV マップを作成します。 右の図では、背景はカラーのキューブイメージと、コントアでモーメント 0 マップを表示しています。


新しいウインドウが開き、作成された PV マップが表示されます。 なお、周波数 (速度) 方向にはキューブの全範囲が指定されて作成されます。

拡大するとこのように見えます。 横軸が位置 (オフセット)、縦軸が速度になります。

短軸方向にスライスした PV マップを作成すると、

新しいウインドウが開き、短軸方向に作成された PV マップが表示されます。

拡大するとこのように見えます。

また、PV マップをタスク impv を使って作成すると以下の様になります。
# In CASA
imname = 'HL_Tau.H2CO.cube.pbcor'
pvname1 = 'HL_Tau.H2CO.cube.pv01'
pvname2 = 'HL_Tau.H2CO.cube.pv02'
ctx=402 ; cty=391 # center pixs
pvpa  = 134.60 # PA (deg)
pleng = 7.0 # arcsec
vlrag = '160~310' # vel.ch.range [ch]
# pv maps
os.system('rm -rf '+pvname1+'*')
impv(imagename = imname, outfile = pvname1, 
     mode = 'length', center = [ctx, cty], 
     length = str(pleng)+'arcsec', pa = str(pvpa)+'deg',
     width = 5, chans = vlrag)
os.system('rm -rf '+pvname2+'*')
impv(imagename = imname, outfile = pvname2, 
     mode = 'length', center = [ctx, cty], 
     length = str(pleng)+'arcsec', pa = str(pvpa+90.)+'deg',
     width = 5, chans = vlrag)
# fits
exportfits(imagename=pvname1, fitsimage=pvname1+'.fits', 
     velocity=True, dropstokes=True, overwrite=True)
上記の impv では、中心位置 ctx, cty(pixels) を基準に、円盤の長軸の位置角 (PA) pvpa(deg) の方向とそれに垂直な方向 (短軸方向) へスライスして、長さ pleng (arcsec) の PV マップを作成します。 また、速度方向へは vlrag(channels) の範囲で作成するように制限しています。 これで、HL_Tau.H2CO.cube.pv01, HL_Tau.H2CO.cube.pv02 の 2 つの PV マップが作成されます。

ケプラー回転を仮定すると、予想されるガスの回転曲線は以下の様に計算されて PV 図上 (長軸の PV マップ) に表示できます。 ケプラー回転する薄く軸対称なガス円盤を考えると、中心から半径 r の回転速度は と表せます。 また、ガス円盤の傾き (face on より) を incl とすると、回転速度 Vrot の視線方向の成分は になります。 この視線方向の速度成分が PV 図上での速度に対応します。 なお、ここでは中心星の質量を 1.3 Mo、星系のシステム速度を Vsys = 7.2 km/s、ガス円盤の傾きを incl = 42.5 deg と仮定しています。

PV 図上のフラックス分布は、ここで予想したケプラー回転の位置-速度プロファイルとよく一致して見えます。 なお図中の黒点線はシステム速度 Vsys と、中心からのオフセットが 0 の位置を示しています。

このイメージを作成するスクリプトとは以下の通りになります。 以下のスクリプトでは、計算された予想される回転曲線を PV マップに重ねるために、Python の matplotlib (図作成) ・astropy (FITS 読込み) のモジュールを使用しています。 前のステップで PV マップのイメージ FITS (HL_Tau.H2CO.cube.pv02.fits) を作成した上で、matplotlib・astropy をインストールした Python 環境で実行して下さい。 
# In Python
import math
from matplotlib import pyplot as plt
from astropy.io import fits
au  = 1.49597870700*10**(11) # m/AU
pc  = 3.085677581*10**(16)   # m/pc
Msol= 1.98847*10**(33)       # g        Solar mass
Grv = 6.67430*10**(-14)      # m3/g s2  Newtonian Constant of Gravitation

# for HL Tau
dist = 140.  # pc
stmass = 1.3 # Mo
Vsys = 7.2   # km/s
incl= 42.5   # deg

pvname1 = 'HL_Tau.H2CO.cube.pv01'
ftdat = fits.open(pvname1+'.fits')
ft_img = ftdat[0].data
ft_hdr = ftdat[0].header
ftdat.close()

NN = int(ft_hdr['NAXIS1']/2.)
dt_rx1 = [[] for _ in range(NN)] ; dt_rx2 = [[] for _ in range(NN)]
dt_vt1 = [[] for _ in range(NN)] ; dt_vt2 = [[] for _ in range(NN)]
for rrr in list(range(NN)):
    rrr += 1
    radu = dist*pc*math.tan(math.radians(rrr*ft_hdr['CDELT1']/3600))/2 # m
    dt_rx1[rrr-1] = -rrr*ft_hdr['CDELT1'] # arcsec
    dt_rx2[rrr-1] =  rrr*ft_hdr['CDELT1']
    dt_vt1[rrr-1] = Vsys - math.sqrt(Grv*(stmass*Msol)/radu)/1000 * math.cos(math.radians(90.-incl)) # km/s
    dt_vt2[rrr-1] = Vsys + math.sqrt(Grv*(stmass*Msol)/radu)/1000 * math.cos(math.radians(90.-incl))

ofs_st=(ft_hdr['CRVAL1']+(              1 -ft_hdr['CRPIX1'])*ft_hdr['CDELT1'])
ofs_ed=(ft_hdr['CRVAL1']+(ft_hdr['NAXIS1']-ft_hdr['CRPIX1'])*ft_hdr['CDELT1'])
vel_st=(ft_hdr['CRVAL2']+(              1 -ft_hdr['CRPIX2'])*ft_hdr['CDELT2'])/1000
vel_ed=(ft_hdr['CRVAL2']+(ft_hdr['NAXIS2']-ft_hdr['CRPIX2'])*ft_hdr['CDELT2'])/1000

ff, aa = plt.subplots(1,1, figsize=(4,5))
img = aa.imshow(ft_img, origin='lower', extent=(ofs_st,ofs_ed,vel_st,vel_ed), 
                vmin=-0.015, vmax=0.05, cmap='jet')
aa.set_xlabel("offset[arcsec]") ; aa.set_xlim(ofs_st,ofs_ed)
aa.set_ylabel("velocity[km/s]") ; aa.set_ylim(vel_st,vel_ed)
aa.set_aspect(0.5)
cbar=plt.colorbar(img, ax=aa)
cbar.set_label("flux density ["+ft_hdr['BUNIT']+"]",size=9)
aa.set_title('pv map')
plt.vlines(   [0], vel_st, vel_ed, "black", linewidth = 0.7, linestyles='dashed')
plt.hlines([Vsys], ofs_st, ofs_ed, "black", linewidth = 0.7, linestyles='dashed')
imr = aa.plot(dt_rx1, dt_vt1, "-", color='white', linewidth = 2.0, label="rotation curve")
imr = aa.plot(dt_rx2, dt_vt2, "-", color='white', linewidth = 2.0)
aa.legend(loc="upper left")
plt.savefig('line_pvmap_model.png')

ページ先頭に戻る

Last Update: 2024.02.15