輻射輸送計算入門¶
ここでは、RADMC-3Dという輻射輸送シミュレーションを行うためのプログラムの使い方を示しています。このページは2020年2月に実施された天文学研究体験2020というプログラムをベースに作られています。なお、このページはjupyter notebookという環境で作成されています。使い方に困ったら jupyter notebookとググってください。
(下記計算コードの一部はCornelis Dullemondや土井聖明さんが書いたものを改変しています。また、観測データは塚越崇さんと橋本淳さんの提供していただいています。下記ではクレジットが明記されておらず申し訳ありませんが、この場を借りて感謝いたします。）
そもそも -- 使う計算機の仕組み(?)の大まかな説明¶
天文学の研究ではUnixというシステム上で動くプログラムを使います。Unixは、例えば「端末」というアプリケーションから操作できます。コマンドと呼ばれる文字列を打ち込んで、操作します。例えば ls と打ち込んでエンターを押すと、今いるディレクトリ（フォルダ）のファイルリストが出てきます。別のフォルダの中身を見たい場合は、（普通のパソコンならダブルクリックすれば済む話ですが）Unixでは 「cd (開きたいフォルダ名)」というコマンドを打ち込みます。こんな感じで進めて行きます。これらのコマンドも、困ったら適宜ググってみてください。
必要なソフトのインストール1 -- jupyter notebook¶
これは一般的に使われるソフトウェアの一つです。webブラウザを用いてプログラムを動かしたり描画をすることができます。コマンドラインに「jupyter notebook」と打ち込んで表示されればあなたのパソコンには既にインストールされています(anaconda等をインストールしていると自動的に入っている)。そうでない場合は、グーグル検索等を駆使して自力でインストールしておいてください。(例えば「anaconda install mac」や「anaconda install windows」のように検索すると良いと思います)
必要なソフトのインストール2 -- RADMC-3D¶
こちらは天文学において輻射輸送計算を用いるときに使うソフトなので、普通インストールはされていません。下記ページに行って
Download→current releaseと進んでダウンロードしてください。おそらくダウンロードフォルダにダウンロードされると思います。つづいて、端末（ターミナル）からこのダウンロードフォルダに移動します。このダウンロードフォルダがどこにあるかはパソコンによるのですが、殆どの場合ホームディレクトリ直下にダウンロードフォルダがあります。「cd」を打つと自動的にホームディレクトリに移りますので、「ls」を押してそこにあるディレクトリリストを表示しましょう。もしそこに「Downloads」のようなディレクトリがあれば、「cd Downloads」と打って、ダウンロードフォルダに言ってください。そこに先程ダウンロードしたファイルが有るはずです。
ダウンロードしたファイルを、インストールしたい場所に移動します。例えばホームディレクトリに移動する場合は、「mv [ダウンロードしたファイル] ~」とすればホームディレクトリに移動できます。
次に圧縮されているファイルを展開します。zipファイルであればunzipというコマンドで展開できます。
展開したディレクトリの中にsrcというディレクトリがあるはずです。これはソースコードが全て入っているディレクトリです。ここに行って(cdコマンドを利用してこのディレクトリに行って)、make とタイプして、その後にmake install をします。こうすればインストールが完了します。
より詳しい方法は、manualの中のradmc-3d_v0.41.pdfというファイルを見てください。
インストールされたかどうか確認するため、端末上で「radmc3d」と打ち込んでみてください。command not foundと出たらインストールできていません。下記のような表示が出れば成功です。
================================================================
WELCOME TO RADMC-3D: A 3-D CONTINUUM AND LINE RT SOLVER プログラムを動かす準備をする。¶
いよいよプログラムを動かす準備をしましょう。まず、作業したいディレクトリを作成します。例えば、ホームディレクトリに行って（「cd」と打てば行ける）、「mkdir test1」などと打てばtest1というディレクトリができます。「cd test1」と打てばそのディレクトリに移動できます。
必要なファイルを入手し展開します。
コマンド「jupyter notebook」と打つと、ブラウザにjupyter notebookが立ち上がりますので、planetexperience2020_ver1.0.ipynbというファイルをクリックしてみましょう。
さて、見た目の変わらないページが出てきたと思います。が、このページはプログラムを動かせるページです。ここから直接プログラムを実行できます。仕組みをもっと理解したい人は、jupyter notebookとググってください。
0. 始める前に¶
はじめに断っておきますが、今回の研究体験の数日間で研究のすべてを体験することはできません。ここでは、実際の研究で使われる輻射輸送計算を用いて、実際の観測データを再現するためにはどんな物理条件が必要かを探っていきます。その一方で、得られた条件が一体惑星形成についてどういうふうに面白いかを考えていくところはプログラム化されていません。みなさんは、得られた条件がどう惑星形成に結びつくか、他の参加者や講師・TAと議論しながらどんどん理解を深めていってください。
おさらい
- 観測データについて
- 観測天体の名前はなんですか？
- 観測に使われた望遠鏡はなんですか？
- 観測データは、各ピクセルに明るさの値が入っています。この値の単位は何でしょうか？
- 観測データの動径プロファイルの横軸、縦軸の単位はなんでしょうか？
- 観測波長はいくつですか？（単位も含めて）
- 輻射輸送シミュレーションについて
- いまからやることの目的は何でしょうか。
- 最終的にどんな図を作ればよいでしょうか。
1. Problem setup -- 密度・温度ファイル等を作って輻射輸送計算の準備をする¶
まず、自分のいるディレクトリにどんなファイルがあるかを確認します。lsというコマンドで確認しましょう（このjupyter notebookから呼び出すときは最初に「!」をつけます）
In [1]:
!ls
「planetexperience2020_ver1.0.ipynb」 というファイルが今編集しているファイルです。次のBoxで、これから計算する天体の構造などの情報を作成し、ファイルとして書き出します。書き出すファイルは主に「.inp」という拡張子です。次のボックスにカーソルを合わせてクリックし、Ctrl+Enterを押してみてください。あるいは、上の方にあるRunというボタンを押してみてください。
In [2]:
# problem setup
# このボックス自体は、今からシミュレーションしたい密度温度分布を計算し、
# それらをファイルとして書き出すものです。
# 出力されたファイルをradmc-3dが読み込んで輻射輸送計算を行います。
import numpy as np
import math
# ----------------------------------
# Physical constants in cgs unit
# ----------------------------------
# ここでは、後で使うときに便利な定数を先に定義しています。
# なお、天文学ではcgs単位系を使い続けています。
# 長さの単位は[cm]、重さの単位は[g]です。ご容赦ください。
au = 1.49598e13 #地球と太陽の距離[cm]
M_sun = 2.0e33 # 中心星質量[g]
G = 6.67e-8 # 万有引力定数
mu = 2.34 # 平均分子量
kb = 1.38e-16 # boltzmann定数
NA = 6.02e23 # アボガドロ定数
mH = 1.0/NA # 水素の質量
pc = 3.08572e18 # Parsec [cm]
ms = 1.98892e33 # Solar mass [g]
ts = 5.78e3 # Solar temperature [K]
ls = 3.8525e33 # Solar luminosity [erg/s]
rs = 6.96e10 # Solar radius [cm]
# --------------------------------------------------------------------
# Make the grid
# 3次元のグリッドを作り、amr_grid.inpというファイルに書き込む。
# --------------------------------------------------------------------
#
# grid parameters
#
nr = 32 # 動径方向のgrid数。初期設定は32。初期設定は荒い設定。
ntheta = 32 # 仰角方向のgrid数。初期設定は32。初期設定は荒い設定。
nphi = 1 # 方位角方向のgrid数。初期設定は1。すなわち軸対称。
rin = 10*au # 円盤の内縁。これより内側は密度を0とする [cm]
rout = 100*au # 円盤の外縁。これより外側は密度を0とする [cm]
r0 = 1.0*au # べき分布を考える際の基準点。[cm]
thetaup = np.pi*0.5-0.7# 仰角のシミュレーションの範囲(rad)。 天頂から引き算部分の角度(rad)分だけ計算する。他は計算しない。
#
# Make the coordinates
#
#ri = np.logspace(np.log10(rin),np.log10(rout),nr+1)
ri = np.linspace(rin,rout,nr+1)
thetai = np.linspace(thetaup,0.5e0*np.pi,ntheta+1)
phii = np.linspace(0.e0,np.pi*2.e0,nphi+1)
rc = 0.5 * ( ri[0:nr] + ri[1:nr+1] )
thetac = 0.5 * ( thetai[0:ntheta] + thetai[1:ntheta+1] )
phic = 0.5 * ( phii[0:nphi] + phii[1:nphi+1] )
#
#
# Make the grid
#
qq = np.meshgrid(rc,thetac,phic,indexing='ij')
rr = qq[0]
tt = qq[1]
zr = np.pi/2.e0 - qq[1]
#
# Write the grid file
#
with open('amr_grid.inp','w+') as f:
f.write('1\n') # iformat
f.write('0\n') # AMR grid style (0=regular grid, no AMR)
f.write('100\n') # Coordinate system: spherical
f.write('0\n') # gridinfo
f.write('1 1 0\n') # Include r,theta coordinates
f.write('%d %d %d\n'%(nr,ntheta,1)) # Size of grid
for value in ri:
f.write('%21.14e\n'%(value)) # X coordinates (cell walls)
for value in thetai:
f.write('%21.14e\n'%(value)) # Y coordinates (cell walls)
for value in phii:
f.write('%21.14e\n'%(value)) # Z coordinates (cell walls)
# ----------------------------------
# Star parameters
# 中心星のパラメータを決めて stars.inp ファイルを作る
# ----------------------------------
mstar = ms
rstar = rs
tstar = ts
pstar = np.array([0.,0.,0.])
#本来は星のスペクトルを知らないと星の放射は再現できない。
#ここでは、簡単のため、星の放射は表面温度tstarの黒体放射を仮定する。
# --------------------------------------------------------------------
# Temperature profile
#
# 温度分布はべき乗分布を仮定する。下記のような分布を考える。
# T = T_0 * ( r /1 au) ^ q
# すなわち、温度は軌道長半径にのみ依存し鉛直方向は一様温度とする。
# --------------------------------------------------------------------
t0 = 100.0 #r0における温度[K]。*今回はこのように指定するが、本当は中心星の明るさから自動的に決まるはず。
pltemp = -0.5 #べき乗分布のべき。無次元。
tdust = t0 * (rr/r0)**pltemp #温度分布の配列[K]
#
# Write the temperature file
#
with open('dust_temperature.dat','w+') as f:
f.write('1\n') # Format number
f.write('%d\n'%(nr*ntheta*nphi)) # Nr of cells
f.write('1\n') # Nr of dust species
data = tdust.ravel(order='F') # Create a 1-D view, fortran-style indexing
data.tofile(f, sep='\n', format="%21.14e")
f.write('\n')
# ----------------------------------
#
# Setting the scale height
#
# ダストの鉛直方向の典型的な高さであるスケールハイトを計算しておく。後の密度分布計算に必要。
# この値は、温度分布で決まっている。
#
# ----------------------------------
f_set = 1.0 # Dust settling factor *ダストがガスに対してどれくらい沈殿しているかを表すパラメータ
hh = np.sqrt(kb*tdust/(mH*2.34)) / np.sqrt(G*mstar/rr**3)/f_set #[cm]
hhr = hh / rr #スケールハイトをその場所の軌道長半径で割ったもの。アスペクト比と呼ぶ。[無次元]
# ----------------------------------
#
# Make the dust density model
# ダストの空間密度分布を作り、dust_density.inpというファイルに書き込む。
#
# ----------------------------------
rhod = rr*0 #rhodは空間密度\rho_{dust}のこと。単位は[g/cm^3]。各グリッドにこの値を書き込む必要がある。
# power-law profile
#
# 円盤のダスト面密度[g/cm^2]をSigma_{dust}、軌道長半径をr、pを定数だとして、下記のような分布を考える。
# Sigma_{dust} = Sigma_{dust,0} * ( r /1 au) ^ p
# このような分布をべき乗分布(power-law profile)と呼ぶ
# ここでSigma_{dust,0}=0とすれば、もちろん全部０になる。
sigmad0 = 100 # Sigma_{dust,0}のこと。軌道長半径がr0の位置での面密度。[g/cm^2]
plsig = -1.0 # pのこと。面密度の軌道長半径依存性のべき。
sigmad = sigmad0 * (rr/r0)**plsig
rhod +=( sigmad / (np.sqrt(2.e0*np.pi)*hh) ) * np.exp(-(zr**2/hhr**2)/2.e0)
#rhodは空間密度[g/cm^3]、sigmadは[g/cm^2]。最後にrhodに入れるときに、鉛直方向の構造を仮定した。
#今回は高さ方向に、スケールハイトhhを持つガウシアン分布を仮定した。
# ring-like profile
# 動径方向にGaussianの形をしたリングを想定する。
sigmadc = 0.6 #Sigma_{dust, center}リング中心のダスト面密度[g/cm^2]
ringcau = 50.0 #R_{ring, center}リング中心の位置をau単位で表したもの。
ringwau = 3.0 #R_{ring, width} リングの幅をau単位で表したもの。
ringc = ringcau*au #上記au単位で指定したものを[cm]単位に変換
ringw = ringwau*au #上記au単位で指定したものを[cm]単位に変換
sigmadring = sigmadc * np.exp(-((rr-ringc)**2/ringw**2)/2.0) #動径方向にガウシアンを仮定した面密度分布[g/cm^2]
rhod += ( sigmadring / (np.sqrt(2.e0*np.pi)*hh) ) * np.exp(-(zr**2/hhr**2)/2.e0)
#
# Write the density file
#
with open('dust_density.inp','w+') as f:
f.write('1\n') # Format number
f.write('%d\n'%(nr*ntheta*nphi)) # Nr of cells
f.write('1\n') # Nr of dust species
data = rhod.ravel(order='F') # Create a 1-D view, fortran-style indexing
data.tofile(f, sep='\n', format="%21.14e")
f.write('\n')
# ----------------------------------
#
# Write the wavelength_micron.inp file
# 波長の情報を持ったファイルを作る。
#
# ----------------------------------
lam1 = 0.1e0
lam2 = 7.0e0
lam3 = 25.e0
lam4 = 1.0e4
n12 = 20
n23 = 100
n34 = 30
lam12 = np.logspace(np.log10(lam1),np.log10(lam2),n12,endpoint=False)
lam23 = np.logspace(np.log10(lam2),np.log10(lam3),n23,endpoint=False)
lam34 = np.logspace(np.log10(lam3),np.log10(lam4),n34,endpoint=True)
lam = np.concatenate([lam12,lam23,lam34])
nlam = lam.size
#
# Write the wavelength file
#
with open('wavelength_micron.inp','w+') as f:
f.write('%d\n'%(nlam))
for value in lam:
f.write('%21.14e\n'%(value))
#
#
# Write the stars.inp file
#
with open('stars.inp','w+') as f:
f.write('2\n')
f.write('1 %d\n\n'%(nlam))
f.write('%21.14e %21.14e %21.14e %21.14e %21.14e\n\n'%(rstar,mstar,pstar[0],pstar[1],pstar[2]))
for value in lam:
f.write('%21.14e\n'%(value))
f.write('\n%21.14e\n'%(-tstar))
#これ以外にも必要な輻射輸送計算のコントロールファイルを作る。
nphot = 10 # これはthermal monte carlo のときに使用するものなので今回は不要
#
# Dust opacity control file
#
with open('dustopac.inp','w+') as f:
f.write('2 Format number of this file\n')
f.write('1 Nr of dust species\n')
f.write('============================================================================\n')
f.write('1 Way in which this dust species is read\n')
f.write('0 0=Thermal grain\n')
f.write('silicate Extension of name of dustkappa_***.inp file\n')
f.write('----------------------------------------------------------------------------\n')
#
# Write the radmc3d.inp control file
#
with open('radmc3d.inp','w+') as f:
f.write('nphot = %d\n'%(nphot))
f.write('scattering_mode_max = 0\n') #
f.write('iranfreqmode = 1\n')
では、lsでどんなファイルができたか確認しましょう。
In [3]:
!ls
例えば、dust_density.inp というファイルには、これから輻射輸送計算を走らせる上でどんな密度構造をおいたかが記録されています。dust_temperature.datというファイルには温度構造が記録されています。次の輻射輸送計算では、これらのファイルを読み込んで計算します。
2. 輻射輸送計算本番 -- RADMC-3Dを動かす¶
ここまでで輻射輸送計算の準備ができました。次に輻射輸送計算を行います。次のボックスをrunしてください(!マークを入れると、ここでシェルコマンドを動かせます。)
(補足) シェルとは？：「端末」や「ターミナル」を立ち上げて「cd」とか「ls」とコマンドを打ち込んでいるそこを、シェルと呼んでいます。シェルで打ち込むコマンドのことをシェルコマンドと言います。「cd」や「ls」はシェルコマンドの一種です。シェルにはいくつか種類があります（bashやzsh）。計算機室のターミナルでは、デフォルトのbashを使っています。¶
ここで使用しているパラメータがなにかについては、RADMC-3Dのマニュアルを見てください。「9.1 Basics of image making with RADMC-3D」です。 例えば、 lambda 1300 は波長1300micronであることを表します。 npixは 一辺のピクセルの数です。
In [4]:
!radmc3d image lambda 1300 npix 200 sizeau 400 incl 0 posang 30
image.outというファイルができていれば成功です。試しにlsでファイルができているかみてみましょう。(ちなみに、上の出力の中のDust mass totalという欄は見ておくと良いかもしれません。例えば円盤のダスト質量が太陽質量より大きかったら、非物理的な感じがしますよね）
In [5]:
!ls
出来上がったimage.outというファイルには、二次元画像データがテキストとして保存されています。これが、輻射輸送計算の結果です。しかし、このデータは数字の羅列であり、まだ人間の読める形にはなりきれていません。そこで次に、結果を二次元カラープロットして確認してみます。
3. 解析¶
3-1. 図を書く --- テキストデータとして出力されたimage.outを読み込んで絵にする¶
'image.out'というファイルを読み込んで、それを2次元カラープロットする計算です。とりあえず下のボックスをRunしてみましょう。
In [6]:
import numpy as np
# ------------------------------------------------
#
# reading the image.out file
#
# ファイルを読み込みます。
# まず、読み込むファイルの名前を指定します。直前の「radmc3d image ...」というコマンドで作成された image.outというファイルを使います。
# pixsizeは一つのピクセルの物理的な長さです。これは「radmc3d image ...」のコマンドで、
# npix (一辺を何ピクセルに分割するか)とsizeau（一辺の物理長さを何auにするか）の２つを指定していますので、
# そこで決まっています。（いまはその値をファイルから読み込むようになっている。）
# image.outのファイル構造を知りたい人は、radmc3dのマニュアルA.16章を御覧ください。
#
# ------------------------------------------------
fnameread = 'image.out'
imrangeau = [-200,200,-200,200]
with open(fnameread,'r') as f:
s = f.readline()
im_n = f.readline()
im_nx = int(im_n.split()[0])
im_ny = int(im_n.split()[1])
nlam = int(f.readline())
pixsize = f.readline()
pixsize_x = float(pixsize.split()[0])
pixsize_y = float(pixsize.split()[1])
wavelength_mic = np.zeros(nlam)
im = np.zeros((im_nx,im_ny,nlam))
rr_au = np.zeros((im_nx,im_ny,nlam))
phi = np.zeros((im_nx,im_ny,nlam))
for ilam in range(nlam): wavelength_mic[ilam] = f.readline()
for ilam in range(nlam):
dummy = f.readline() #blank line
for iy in range(im_ny):
for ix in range(im_nx):
image_temp = f.readline()
im[ix,iy,ilam] = image_temp
rr_au[ix,iy,ilam] = np.sqrt(((ix-im_nx/2)*pixsize_x)**2+((iy-im_ny/2)*pixsize_y)**2)/au
phi[ix,iy,ilam] = np.arctan2((iy-im_ny/2)*pixsize_y,(ix-im_nx/2)*pixsize_x)
sizex_au = im_nx*pixsize_x/au
sizey_au = im_ny*pixsize_y/au
# ------------------------------------------------
#
# Beam convolution
#
# ガウシアンビームでなまします。
# ------------------------------------------------
'''
from scipy.ndimage import gaussian_filter
im_smoothed = gaussian_filter(im,10)
'''
# ------------------------------------------------
#
# Conversion of the intensity unit
#
# image.outから呼んだ画像情報は、いまimという配列に格納されています。
# imという配列には、各ピクセルに天体の明るさ（輝度=Intensity）がcgs 単位系で入っています。
# (書き下すと、「erg/s/cm/cm/Hz/ster」です)
# これを、観測画像と比べるときには観測された輝度の単位と合わせる必要があります。
# 観測された輝度の単位は望遠鏡にもよるので、ここではいくつかの選択肢を例示します。
# ------------------------------------------------
# ------------------------------------------------
# parameters:
#
# distance to the object
# 地球に届く天体の明るさを計算するときは、天体までの距離の情報が必要です。
# (天体から地球までの距離が遠ければ遠いほど暗いので)
d_pc = 100.0
# beam size
Bmaj_as = 0.0727 # Major axis of the beam in arcsec
Bmin_as = 0.0478 # Minor axis of the beam in arcsec
wavelength = 0.13 # wavelength in unit of cm
# ------------------------------------------------
# im_Jyppix : intensity with a unit of [Jy/pixel]
# Conversion from erg/s/cm/cm/Hz/ster to Jy/pixel
#
# まず、Jy(ジャンスキー)という電波望遠鏡でよく使う単位に変換します。
# また、とりあえず「1ピクセルあたり何Jyか」という単位である「Jy/pixel」を使いましょう。
'''
conv = pixsize_x * pixsize_y / (d_pc*pc)**2. * 1e23
im_Jyppix = np.zeros((im_nx,im_ny,nlam))
im_Jyppix = im.copy()*conv
im_forplot = im_Jyppix
pltunit = '[Jy/pix]'
'''
# im_Jypas2 : intensity with a unit of [Jy/arcsec^2]
# 1 arcsecは、距離d_pc[pc]において d_pc [au]。
'''
conv = pixsize_x * pixsize_y / (d_pc*pc)**2. * 1e23 * (pixsize_x/au*d_pc)*(pixsize_y/au*d_pc)
print('plotting with the unit of [Jy/arcsec^2]')
print('ref: arcsec^2 in the unit of au^2 = ', (pixsize_x/au*d_pc)*(pixsize_y/au*d_pc))
im_Jypas2 = np.zeros((im_nx,im_ny,nlam))
im_Jypas2 = im.copy()*conv
im_forplot = im_Jypas2
pltunit = '[Jy/arcsec^2]'
'''
# im_Jypbeam : intensity with a unit of [Jy/beam]
# 電波望遠鏡ではビームサイズという実効的空間分解能を表す単位を使います。
# FWHM(full width at half maximum)という単位を使うので、
# ガウシアンの標準偏差とは少し異なります。
'''
Bmajinpix = Bmaj_as*pixsize_x/au*d_pc
Bmininpix = Bmin_as*pixsize_x/au*d_pc
pixelinbeam = (np.pi/4.0/np.log(2.0))*Bmajinpix*Bmininpix
print('plotting with the unit of [Jy/beam]')
print('Bmaj, Bmin = ', Bmaj_as, ',', Bmin_as, 'in arcseconds')
conv = pixsize_x * pixsize_y / (d_pc*pc)**2. * 1e23 * pixelinbeam
print('ref: how many pixels in one beam?:', pixelinbeam)
im_Jypbeam = np.zeros((im_nx,im_ny,nlam))
im_Jypbeam = im.copy()*conv
im_forplot = im_Jypbeam
pltunit = '[Jy/beam]'
'''
# im_BT_RJ : intensity with a unit of brightness temperature in RJ limit [K]
# 天体の輝度を、温度何Kの黒体輻射に対応するかで表すのが輝度温度です。
# ここではレイリージーンズ測を仮定します。
# 単位変換は例えばここを参考にしてください。
# https://science.nrao.edu/facilities/vla/proposing/TBconv
#'''
conv = wavelength*wavelength/kb/2.0
im_TB_RJ = np.zeros((im_nx,im_ny,nlam))
im_TB_RJ = im.copy()*conv
im_forplot = im_TB_RJ
pltunit = '[K]'
#'''
# ------------------------------------------------
#
# Plotting the image...
#
# 以下、図をプロットするための情報です。
# intmax,intminは表示画像の最大値・最小値を表します。
# もし画像に何も見えない場合はここを変更してみてください。
# ------------------------------------------------
fontsize = 12
fontfamily = 'serif'
interpolation = 'bicubic'
cmap = 'Oranges'
intmax = np.max(im_forplot)
intmin = np.max(im_forplot)*0.01
tickcolor = 'white'
xlabel = 'x [AU]'
ylabel = 'y [AU]'
norm = None
drawpolvect = False
title = 'default image'
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
%matplotlib inline
matplotlib.rcParams.update({'font.size': fontsize,'font.family': fontfamily})
fig=plt.figure(num=None, figsize=(8,6), dpi=400, facecolor='w', edgecolor='k')
ax1=fig.add_subplot(111)
norm=LogNorm()
#norm=None
plt.imshow(im_forplot[:,:,0].T,interpolation=interpolation, cmap=cmap, vmin=intmin, vmax=intmax,norm=norm,origin='norm',extent= [-sizex_au/2.0,sizex_au/2.0,-sizey_au/2.0,sizey_au/2.0])
plt.xlabel(xlabel)
plt.ylabel(ylabel)
plt.xlim(imrangeau[0],imrangeau[1])
plt.ylim(imrangeau[2],imrangeau[3])
plt.title(title)
plt.tick_params(which='major', color=tickcolor,width=2.0,length=6.0)
#
# Setting the colorbar
#
cbar=plt.colorbar()
cbar.set_label(pltunit)
plt.savefig('2dplot.pdf')
#次のボックスでの計算用に以下の配列を準備しておきます。
im_for1dplot=im_forplot.copy()
この画像を見ながら、一つ前の 'radmc3d image' を走らせたときのパラメータを確認しましょう。例えば、incl 70 は円盤を70度傾いた場所から見ることに対応しています。
3-1. 単位を変換する¶
RADMC-3Dは、天体の輝度分布をcgs単位系で出力しています。観測画像と比較するため、この単位を観測の単位と揃える必要があります。これに必要な単位の変換が、上記ボックス内でコメントアウトしてありますので、コメントアウトを外したりしながら別単位で表示をしてみてください。
3-2. ビームでなます¶
望遠鏡は解像度（あるいは空間分解能）という概念があります。視力と同じようなもので、望遠鏡の口径（あるいは基線長）がながければ長いほど、より高解像度の画像を得ることができます。すなわち、実際の望遠鏡観測データは、望遠鏡の大きさによって少し「ぼけて」見えています。これを再現するために、輻射輸送計算の結果を、望遠鏡の大きさに対応してぼかす必要があります。
Beam convolution という場所がコメントアウトしてあります。コメントアウトを外して再度runしてみてください。
3-4. 輻射輸送計算結果の1次元輝度分布を書く --- 観測と直接比較のために¶
さて、輻射輸送計算の結果を眺めることに成功しました。しかし、これが実際に観測とどれくらい整合的なのか、このままでは判断が難しいです。そこで、観測結果と直接比較するための一つの手法として、ここでは動径プロファイルを用います。中心からとある扇形の方向の輝度を、中心からの距離の関数として書くことで、横軸中心からの距離、縦軸輝度の図を書くことができます。こうすれば、観測の図を重ねて、直接比較することができます。
いま、im_for1dplotという配列に、各点の輝度の情報が入っています。この中のとある領域の輝度を1次元プロットすることを考えてみます。以下のボックスをRunしてみてください。
In [7]:
#注意点：このプログラムは、通常のx,y軸における極座標を使っています。
#すなわち、単位円上で(x,y)=(1,0)の場所を角度0度として、(x,y)=(0,1)をπ/2, (x,y)=(0,-1)を-π/2としています。
#値の範囲は-π<θ<πとなるようにしていますが、詳しくはnumpy.arctan2を調べてください。
#一方で、天文学における場所の角度はPosition Angleと呼ばれ北を0度としたものです。
#通常北を上にして図を作るので、極座標では90度だと思うところが0度になります。
#観測の図を重ねるときはとても注意してください。
#輝度・中心からの距離・その点の角度を1次元配列に入れる。
im_1d = np.ravel(im_for1dplot[:,:,0])
rr_au_1d = np.ravel(rr_au[:,:,0])
phi_au_1d = np.ravel(phi[:,:,0])
#ある角度領域のヒストグラムを書く。
phi0 = np.pi/2.0 # ここでは角度90度のあたりのprofileを書く
mask_condition = np.abs(phi_au_1d-phi0)<10.0*np.pi/180.0 # phi0の周り10度くらいの扇形のプロファイルを書く
im_1d_masked = im_1d[mask_condition]
rr_au_masked = rr_au_1d[mask_condition]
from scipy import stats
bin_means, bin_edges, binnumber=stats.binned_statistic(rr_au_masked,im_1d_masked,statistic='mean',bins=np.arange(1,100,5))
#結果をプロットする
fig=plt.figure(num=None, figsize=(8,6), dpi=400, facecolor='w', edgecolor='k')
plt.subplot(211)
plt.scatter(rr_au_masked,im_1d_masked,label='scatter plot before average')
plt.hlines(bin_means, bin_edges[:-1], bin_edges[1:], colors='red', lw=1,label='binned statistic of data')
plt.legend()
plt.title('radial profile in '+pltunit+' along theta = 90 degrees')
plt.xlim([0,100])
plt.ylim([0,np.max(im_forplot)*1.2])
#別の角度領域のヒストグラムを書く。
phi0 = 0 # ここでは角度0度のあたりのprofileを書く
mask_condition = np.abs(phi_au_1d-phi0)<10.0*np.pi/180.0 # phi0の周り10度くらいの扇形のプロファイルを書く
im_1d_masked = im_1d[mask_condition]
rr_au_masked = rr_au_1d[mask_condition]
bin_means, bin_edges, binnumber=stats.binned_statistic(rr_au_masked,im_1d_masked,statistic='mean',bins=np.arange(1,100,5))
plt.subplot(212)
norm=LogNorm()
#norm=None
plt.scatter(rr_au_masked,im_1d_masked,label='scatter plot before average')
plt.hlines(bin_means, bin_edges[:-1], bin_edges[1:], colors='red', lw=1,label='binned statistic of data')
plt.title('radial profile in '+pltunit+' along theta = 0 degrees')
plt.xlim([0,100])
plt.ylim([0,np.max(im_forplot)*1.2])
plt.tight_layout()
plt.savefig('1dplot.pdf')
3-5. 観測から得られたプロファイルを重ねて書く¶
では、観測プロファイルを重ねて書いてみましょう。'twhya.major.dat'というファイルにテスト用の観測結果が入っています。これを読み込んで重ねて書いてみます。なお、テストデータは Tsukagoshi et al. 2016, ApJL, 829, L35 の結果です。ビームサイズは72.7*47.8 milliarcsecondsです。
In [8]:
#結果をプロットする
fig=plt.figure(num=None, figsize=(8,6), dpi=400, facecolor='w', edgecolor='k')
plt.subplot(111)
plt.scatter(rr_au_masked,im_1d_masked,label='RADMC-3D results '+pltunit)
plt.hlines(bin_means, bin_edges[:-1], bin_edges[1:], colors='red', lw=1,label='binned statistic of data')
plt.title('radial profile')
plt.legend()
#観測結果を読み込んでプロットする
data_twhya_major=np.loadtxt('twhya.major.dat')
#plt.plot(data_twhya_major[:,0],data_twhya_major[:,1],label='TWHya profile [Jy/beam]')
conv_TB_RJ = 1.36e3*wavelength*wavelength/Bmaj_as/Bmin_as
plt.plot(data_twhya_major[:,0],data_twhya_major[:,1]*conv_TB_RJ,label='TWHya profile [K]')
plt.legend()
plt.xlim([0,100])
plt.ylim([1e-2,1e2])
plt.yscale('log')
plt.tight_layout()
plt.savefig('1dplot_withobservation.pdf')
これでようやく、観測結果と輻射輸送計算を直接比較することができるようになりました。例えば「40auの位置では、明るさが3桁足りないなどうしよう。全体の面密度を上げればいいだろうか？それとも円盤の面密度の動径依存性を帰るべきだろうか？」などです。これらの直感をもとに、パラメータを変えて、観測を再現できるモデルを探していきます。なお、描画の単位が違う場合は3-1まで戻って単位の変換をやり直してみてください。
4. パラメータを変える --- そもそも何をしていたかを一つずつ理解する。¶
ここまで、言われるがままにコマンドを打ち込んできたことでしょう。何やってるかわからなくても、一通りの流れを理解できたと思います。ここから、上記で何をやっていたかを少しずつ理解しながら、観測量を説明するシミュレーションパラメータを探してもらいます。なお、この際radmc-3dのディレクトリのどこかにあるmanualも参考にしてください。下記には、少し自分でいじるとき用に変えられるパラメータを例示しています。これ以外にもいくらでも方向はあると思います。
- 円盤を観測する角度を変えたらどうなるでしょう？ 例：radmc3dを動かしたときのパラメータを incl 40から、incl 80にしてみる。
- リングの幅を太くするにはどの数字を変えればいいでしょう？
- リングの位置を変えるにはどの数字を変えればいいでしょう？
- リングの幅を1auにしましょう。ぼやけて見えますか？もっと解像度を上げたい場合、どのパラメータを変えるべきでしょうか。
- リングを2つおいてみましょうか。
- リングではなく、スムーズな円盤を再現するにはどうするべきでしょうか？
- いまの観測波長は何マイクロメートルでしょうか？それは赤外線？電波？波長を変えるにはどうする？波長を変えたら結果はどう変わる？
4-1. パラメータを変える例: 波長を変える --- 特に赤外線の結果¶
先程は波長1300マイクロメートル、すなわち電波領域でのシミュレーションを行いました。これは、例えばアルマ望遠鏡の観測できる波長に対応します。ここではパラメータを変えて、近赤外線でのシミュレーションをする例を示します。近赤外線は例えばすばる望遠鏡の観測できる波長に対応します。
なお、この先のボックスを動かすと、*.inpファイル等がすべて上書きされてしまいます。上書きを避けたい場合は、ディレクトリを分けることをおすすめします。例えば下記のようにしてみると良いと思います。
In [9]:
!mkdir test2
In [10]:
!cp *.ipynb test2
In [11]:
!cp dustkapscatmat_pyrmg70.inp test2
そして、ブラウザのHomeタブから、test2を選択、その中の.ipynbファイルをクリックして、再開しましょう。
ところで、赤外線では偏光データを利用しています。それに伴い、下記ではこっそり以下のような変更をしています。変更した箇所は、[変更点]と書いてありますので確認してみてください。
- thetaのグリッド配置を上下対称を仮定しないものにした。
- 散乱ありの計算にするため、opacityファイルに散乱の情報を入れた。
- それに伴い scattering_mode_max = 5にした。（フル散乱）
- 描画も、polarized intensityを書き足した。
In [12]:
# problem setup
# このボックス自体は、今からシミュレーションしたい密度温度分布を計算し、
# それらをファイルとして書き出すものです。
# 出力されたファイルをradmc-3dが読み込んで輻射輸送計算を行います。
import numpy as np
import math
# ----------------------------------
# Physical constants in cgs unit
# ----------------------------------
# ここでは、後で使うときに便利な定数を先に定義しています。
# なお、天文学ではcgs単位系を使い続けています。
# 長さの単位は[cm]、重さの単位は[g]です。ご容赦ください。
au = 1.49598e13 #地球と太陽の距離[cm]
M_sun = 2.0e33 # 中心星質量[g]
G = 6.67e-8 # 万有引力定数
mu = 2.34 # 平均分子量
kb = 1.38e-16 # boltzmann定数
NA = 6.02e23 # アボガドロ定数
mH = 1.0/NA # 水素の質量
pc = 3.08572e18 # Parsec [cm]
ms = 1.98892e33 # Solar mass [g]
ts = 5.78e3 # Solar temperature [K]
ls = 3.8525e33 # Solar luminosity [erg/s]
rs = 6.96e10 # Solar radius [cm]
# --------------------------------------------------------------------
# Make the grid
# 3次元のグリッドを作り、amr_grid.inpというファイルに書き込む。
# --------------------------------------------------------------------
#
# grid parameters
#
nr = 32 # 動径方向のgrid数。初期設定は32。初期設定は荒い設定。
ntheta = 32 # 仰角方向のgrid数。初期設定は32。初期設定は荒い設定。
nphi = 1 # 方位角方向のgrid数。初期設定は1。すなわち軸対称。
rin = 10*au # 円盤の内縁。これより内側は密度を0とする [cm]
rout = 100*au # 円盤の外縁。これより外側は密度を0とする [cm]
r0 = 1.0*au # べき分布を考える際の基準点。[cm]
thetaup = np.pi*0.5-0.7# 仰角のシミュレーションの範囲(rad)。 天頂から引き算部分の角度(rad)分だけ計算する。他は計算しない。
#
# Make the coordinates
#
#ri = np.logspace(np.log10(rin),np.log10(rout),nr+1)
ri = np.linspace(rin,rout,nr+1)
#thetai = np.linspace(thetaup,0.5e0*np.pi,ntheta+1)
thetai = np.linspace(thetaup,np.pi-thetaup,ntheta+1)#[変更点]鏡対称の条件を外した。
phii = np.linspace(0.e0,np.pi*2.e0,nphi+1)
rc = 0.5 * ( ri[0:nr] + ri[1:nr+1] )
thetac = 0.5 * ( thetai[0:ntheta] + thetai[1:ntheta+1] )
phic = 0.5 * ( phii[0:nphi] + phii[1:nphi+1] )
#
#
# Make the grid
#
qq = np.meshgrid(rc,thetac,phic,indexing='ij')
rr = qq[0]
tt = qq[1]
zr = np.pi/2.e0 - qq[1]
#
# Write the grid file
#
with open('amr_grid.inp','w+') as f:
f.write('1\n') # iformat
f.write('0\n') # AMR grid style (0=regular grid, no AMR)
f.write('100\n') # Coordinate system: spherical
f.write('0\n') # gridinfo
f.write('1 1 0\n') # Include r,theta coordinates
f.write('%d %d %d\n'%(nr,ntheta,1)) # Size of grid
for value in ri:
f.write('%21.14e\n'%(value)) # X coordinates (cell walls)
for value in thetai:
f.write('%21.14e\n'%(value)) # Y coordinates (cell walls)
for value in phii:
f.write('%21.14e\n'%(value)) # Z coordinates (cell walls)
# ----------------------------------
# Star parameters
# 中心星のパラメータを決めて stars.inp ファイルを作る
# ----------------------------------
mstar = ms
rstar = rs
tstar = ts
pstar = np.array([0.,0.,0.])
#本来は星のスペクトルを知らないと星の放射は再現できない。
#ここでは、簡単のため、星の放射は表面温度tstarの黒体放射を仮定する。
# --------------------------------------------------------------------
# Temperature profile
#
# 温度分布はべき乗分布を仮定する。下記のような分布を考える。
# T = T_0 * ( r /1 au) ^ q
# すなわち、温度は軌道長半径にのみ依存し鉛直方向は一様温度とする。
# --------------------------------------------------------------------
t0 = 100.0 #r0における温度[K]。*今回はこのように指定するが、本当は中心星の明るさから自動的に決まるはず。
pltemp = -0.5 #べき乗分布のべき。無次元。
tdust = t0 * (rr/r0)**pltemp #温度分布の配列[K]
#
# Write the temperature file
#
with open('dust_temperature.dat','w+') as f:
f.write('1\n') # Format number
f.write('%d\n'%(nr*ntheta*nphi)) # Nr of cells
f.write('1\n') # Nr of dust species
data = tdust.ravel(order='F') # Create a 1-D view, fortran-style indexing
data.tofile(f, sep='\n', format="%21.14e")
f.write('\n')
# ----------------------------------
#
# Setting the scale height
#
# ダストの鉛直方向の典型的な高さであるスケールハイトを計算しておく。後の密度分布計算に必要。
# この値は、温度分布で決まっている。
#
# ----------------------------------
f_set = 1.0 # Dust settling factor *ダストがガスに対してどれくらい沈殿しているかを表すパラメータ
hh = np.sqrt(kb*tdust/(mH*2.34)) / np.sqrt(G*mstar/rr**3)/f_set #[cm]
hhr = hh / rr #スケールハイトをその場所の軌道長半径で割ったもの。アスペクト比と呼ぶ。[無次元]
# ----------------------------------
#
# Make the dust density model
# ダストの空間密度分布を作り、dust_density.inpというファイルに書き込む。
#
# ----------------------------------
rhod = rr*0 #rhodは空間密度\rho_{dust}のこと。単位は[g/cm^3]。各グリッドにこの値を書き込む必要がある。
# power-law profile
#
# 円盤のダスト面密度[g/cm^2]をSigma_{dust}、軌道長半径をr、pを定数だとして、下記のような分布を考える。
# Sigma_{dust} = Sigma_{dust,0} * ( r /1 au) ^ p
# このような分布をべき乗分布(power-law profile)と呼ぶ
# ここでSigma_{dust,0}=0とすれば、もちろん全部０になる。
sigmad0 = 0 # Sigma_{dust,0}のこと。軌道長半径がr0の位置での面密度。[g/cm^2]
plsig = -1.0 # pのこと。面密度の軌道長半径依存性のべき。
sigmad = sigmad0 * (rr/r0)**plsig
rhod +=( sigmad / (np.sqrt(2.e0*np.pi)*hh) ) * np.exp(-(zr**2/hhr**2)/2.e0)
#rhodは空間密度[g/cm^3]、sigmadは[g/cm^2]。最後にrhodに入れるときに、鉛直方向の構造を仮定した。
#今回は高さ方向に、スケールハイトhhを持つガウシアン分布を仮定した。
# ring-like profile
# 動径方向にGaussianの形をしたリングを想定する。
sigmadc = 0.6 #Sigma_{dust, center}リング中心のダスト面密度[g/cm^2]
ringcau = 50.0 #R_{ring, center}リング中心の位置をau単位で表したもの。
ringwau = 3.0 #R_{ring, width} リングの幅をau単位で表したもの。
ringc = ringcau*au #上記au単位で指定したものを[cm]単位に変換
ringw = ringwau*au #上記au単位で指定したものを[cm]単位に変換
sigmadring = sigmadc * np.exp(-((rr-ringc)**2/ringw**2)/2.0) #動径方向にガウシアンを仮定した面密度分布[g/cm^2]
rhod += ( sigmadring / (np.sqrt(2.e0*np.pi)*hh) ) * np.exp(-(zr**2/hhr**2)/2.e0)
#
# Write the density file
#
with open('dust_density.inp','w+') as f:
f.write('1\n') # Format number
f.write('%d\n'%(nr*ntheta*nphi)) # Nr of cells
f.write('1\n') # Nr of dust species
data = rhod.ravel(order='F') # Create a 1-D view, fortran-style indexing
data.tofile(f, sep='\n', format="%21.14e")
f.write('\n')
# ----------------------------------
#
# Write the wavelength_micron.inp file
# 波長の情報を持ったファイルを作る。
#
# ----------------------------------
lam1 = 0.1e0
lam2 = 7.0e0
lam3 = 25.e0
lam4 = 1.0e4
n12 = 20
n23 = 100
n34 = 30
lam12 = np.logspace(np.log10(lam1),np.log10(lam2),n12,endpoint=False)
lam23 = np.logspace(np.log10(lam2),np.log10(lam3),n23,endpoint=False)
lam34 = np.logspace(np.log10(lam3),np.log10(lam4),n34,endpoint=True)
lam = np.concatenate([lam12,lam23,lam34])
nlam = lam.size
#
# Write the wavelength file
#
with open('wavelength_micron.inp','w+') as f:
f.write('%d\n'%(nlam))
for value in lam:
f.write('%21.14e\n'%(value))
#
#
# Write the stars.inp file
#
with open('stars.inp','w+') as f:
f.write('2\n')
f.write('1 %d\n\n'%(nlam))
f.write('%21.14e %21.14e %21.14e %21.14e %21.14e\n\n'%(rstar,mstar,pstar[0],pstar[1],pstar[2]))
for value in lam:
f.write('%21.14e\n'%(value))
f.write('\n%21.14e\n'%(-tstar))
#これ以外にも必要な輻射輸送計算のコントロールファイルを作る。
nphot = 10 # これはthermal monte carlo のときに使用するものなので今回は不要
#
# Dust opacity control file
#
with open('dustopac.inp','w+') as f:
f.write('2 Format number of this file\n')
f.write('1 Nr of dust species\n')
f.write('============================================================================\n')
f.write('10 Way in which this dust species is read\n')#[変更点]散乱の情報の入ったopacityに変更
f.write('0 0=Thermal grain\n')
f.write('pyrmg70 Extension of name of dustkappa_***.inp file\n')#[変更点]散乱の情報の入ったopacityに変更
f.write('----------------------------------------------------------------------------\n')
#
# Write the radmc3d.inp control file
#
with open('radmc3d.inp','w+') as f:
f.write('nphot = %d\n'%(nphot))
f.write('scattering_mode_max = 5\n') #[変更点]散乱を解くためmode変更
f.write('iranfreqmode = 1\n')
In [13]:
!radmc3d image lambda 1.0 npix 200 sizeau 400 incl 40 posang 30 nphot_scat 100000 stokes setthreads 8
In [15]:
import numpy as np
# ------------------------------------------------
#
# reading the image.out file
#
# ファイルを読み込みます。
# まず、読み込むファイルの名前を指定します。直前の「radmc3d image ...」というコマンドで作成された image.outというファイルを使います。
# pixsizeは一つのピクセルの物理的な長さです。これは「radmc3d image ...」のコマンドで、
# npix (一辺を何ピクセルに分割するか)とsizeau（一辺の物理長さを何auにするか）の２つを指定していますので、
# そこで決まっています。（いまはその値をファイルから読み込むようになっている。）
# image.outのファイル構造を知りたい人は、radmc3dのマニュアルA.16章を御覧ください。
#
# ------------------------------------------------
fnameread = 'image.out'
imrangeau = [-200,200,-200,200]
with open(fnameread,'r') as f:
s = f.readline()
im_n = f.readline()
im_nx = int(im_n.split()[0])
im_ny = int(im_n.split()[1])
nlam = int(f.readline())
pixsize = f.readline()
pixsize_x = float(pixsize.split()[0])
pixsize_y = float(pixsize.split()[1])
wavelength_mic = np.zeros(nlam)
im = np.zeros((im_nx,im_ny,4,nlam)) #[変更点]Stokes I,Q,U,Vの４パラメータ格納用配列
rr_au = np.zeros((im_nx,im_ny,4,nlam)) #[変更点]Stokes I,Q,U,Vの４パラメータ格納用配列
phi = np.zeros((im_nx,im_ny,4,nlam)) #[変更点]Stokes I,Q,U,Vの４パラメータ格納用配列
for ilam in range(nlam): wavelength_mic[ilam] = f.readline()
for ilam in range(nlam):
dummy = f.readline() #blank line
for iy in range(im_ny):
for ix in range(im_nx):
image_temp = f.readline()
im[ix,iy,:,ilam] = [image_temp.split()[istokes] for istokes in [0,1,2,3]] #[変更点]Stokes I,Q,U,Vの４パラメータ格納用配列
rr_au[ix,iy,ilam] = np.sqrt(((ix-im_nx/2)*pixsize_x)**2+((iy-im_ny/2)*pixsize_y)**2)/au
phi[ix,iy,ilam] = np.arctan2((iy-im_ny/2)*pixsize_y,(ix-im_nx/2)*pixsize_x)
sizex_au = im_nx*pixsize_x/au
sizey_au = im_ny*pixsize_y/au
#変更輝度をim_PIに格納
im_PI = np.zeros((im_nx,im_ny,nlam,4))
im_PI = ((im[:,:,1,0])**2.0 + (im[:,:,2,0])**2.0)**0.5
# ------------------------------------------------
#
# Beam convolution
#
# ガウシアンビームでなまします。
# ------------------------------------------------
'''
from scipy.ndimage import gaussian_filter
im_smoothed = gaussian_filter(im,10)
'''
# ------------------------------------------------
#
# Conversion of the intensity unit
#
# image.outから呼んだ画像情報は、いまimという配列に格納されています。
# imという配列には、各ピクセルに天体の明るさ（輝度=Intensity）がcgs 単位系で入っています。
# (書き下すと、「erg/s/cm/cm/Hz/ster」です)
# これを、観測画像と比べるときには観測された輝度の単位と合わせる必要があります。
# 観測された輝度の単位は望遠鏡にもよるので、ここではいくつかの選択肢を例示します。
# ------------------------------------------------
# ------------------------------------------------
# parameters:
#
# distance to the object
# 地球に届く天体の明るさを計算するときは、天体までの距離の情報が必要です。
# (天体から地球までの距離が遠ければ遠いほど暗いので)
d_pc = 100.0
# ------------------------------------------------
# im_Jyppix : intensity with a unit of [Jy/pixel]
# Conversion from erg/s/cm/cm/Hz/ster to Jy/pixel
#
# まず、Jy(ジャンスキー)という電波望遠鏡でよく使う単位に変換します。
# また、とりあえず「1ピクセルあたり何Jyか」という単位である「Jy/pixel」を使いましょう。
conv = pixsize_x * pixsize_y / (d_pc*pc)**2. * 1e23
im_forplot = im[:,:,:,0].copy()*conv
im_PI_forplot = im_PI.copy()*conv
pltunit = '[Jy/pix]'
# ------------------------------------------------
#
# Plotting the image...
#
# 以下、図をプロットするための情報です。
# intmax,intminは表示画像の最大値・最小値を表します。
# もし画像に何も見えない場合はここを変更してみてください。
# ------------------------------------------------
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
%matplotlib inline
matplotlib.rcParams.update({'font.size': fontsize,'font.family': fontfamily})
# Plotting the Stokes I image...
fontsize = 12
fontfamily = 'serif'
interpolation = 'bicubic'
cmap = 'Blues'
intmax = np.max(im_forplot[np.where(im_forplot<np.max(im_forplot)*0.1)])
#解説：一番明るい点の0.1倍以下の明るさの天の中で最も明るい点を最大値とする
#一番明るい点は星のことです。円盤の輝度はそれよりずっと暗いので、円盤の最も明るいところを最大値とするように設定しています。
intmin = intmax*0.01
tickcolor = 'white'
xlabel = 'x [AU]'
ylabel = 'y [AU]'
norm = None
drawpolvect = False
title = 'Stokes I image'
fig=plt.figure(num=None, figsize=(8,6), dpi=400, facecolor='w', edgecolor='k')
ax1=fig.add_subplot(111)
norm=LogNorm()
#norm=None
plt.imshow(im_forplot[:,:,0].T,interpolation=interpolation, cmap=cmap, vmin=intmin, vmax=intmax,norm=norm,origin='norm',extent= [-sizex_au/2.0,sizex_au/2.0,-sizey_au/2.0,sizey_au/2.0])
plt.xlabel(xlabel)
plt.ylabel(ylabel)
plt.xlim(imrangeau[0],imrangeau[1])
plt.ylim(imrangeau[2],imrangeau[3])
plt.title(title)
plt.tick_params(which='major', color=tickcolor,width=2.0,length=6.0)
#
# Setting the colorbar
#
cbar=plt.colorbar()
cbar.set_label(pltunit)
# Plotting the Polarized Intensity image...
fontsize = 12
fontfamily = 'serif'
interpolation = 'bicubic'
cmap = plt.cm.magma
cmap.set_over('white', 1.0)
cmap.set_under('black', 1.0)
cmap.set_bad('black', 1.0)
intmax = np.max(im_PI_forplot)
intmin = np.max(im_PI_forplot)*0.01
tickcolor = 'white'
xlabel = 'x [AU]'
ylabel = 'y [AU]'
norm = None
drawpolvect = False
title = 'Polarized Intensity'
fig2=plt.figure(num=None, figsize=(8,6), dpi=400, facecolor='w', edgecolor='k')
ax2=fig2.add_subplot(111)
norm=LogNorm()
#norm=None
plt.imshow(im_PI_forplot[:,:].T,interpolation=interpolation, cmap=cmap, vmin=intmin, vmax=intmax,norm=norm,origin='norm',extent= [-sizex_au/2.0,sizex_au/2.0,-sizey_au/2.0,sizey_au/2.0])
plt.xlabel(xlabel)
plt.ylabel(ylabel)
plt.xlim(imrangeau[0],imrangeau[1])
plt.ylim(imrangeau[2],imrangeau[3])
plt.title(title)
plt.tick_params(which='major', color=tickcolor,width=2.0,length=6.0)
#
# Setting the colorbar
#
cbar=plt.colorbar()
cbar.set_label(pltunit)
plt.savefig('irplot.pdf')
簡単な解説：赤外線では、中心星の明るさが円盤に対してとても明るいです。もし現実の望遠鏡で観測すると、おそらく星の光が広がって見えてしまい(point spread funciton)、それによって円盤が隠されてしまうでしょう。それを避けるために、偏光強度の観測を行います。星の光は無偏光なのに対し円盤の光は偏光しているので、円盤の成分だけを取り出すことができます。
観測された円盤を再現する（以下は研究体験2020で使用したもの）¶
ではここから、自分の担当する天体を輻射輸送計算で再現し、その構造を議論してもらいます。下記に、必要にいなるかもしれないパラメータをリストしました。適宜参考にしてください。新しい計算を始めるときは、計算ごとにディレクトリを新しくつくり、そこに別のjupyter notebookを立ち上げ、適宜コピーアンドペーストして進めていってください。
TWHya¶
- distance = 60pc
- inclination = 7.0deg
- Mstar = 0.8Msun (Boekel+2016)
- Rstar = 1.2Rsun (Boekel+2016)
- Teff = 3800K (Boekel+2016)
- Lstar = 0.26Lsun (Boekel+2016)
- Tdisk = 26[K]*(R/1au)**(-0.4) (Tsukagoshi+2016)
PDS70¶
- distance = 113pc(Muller+2018)
- inclination = 52deg (Isella+2019)
- Mstar = 0.76Msun (Muller+2018)
- Rstar = 1.26 (Muller+2018)
- Teff = 3972K (Muller+2018)
- Lstar = 0.36Lsun (Muller+2018)
- Tdisk = 19[K](R/22au)(-0.5) (Keppler+2019)
J1604¶
- distance = 150pc
- inclination = (Mayama+2018)
- Mstar = 1Msun (Preibisch+2002)
- Rstar = 1.41 Rsun (Preibisch+2002)
- Teff = 4550K (Preibisch+2002)
- Lstar = 0.76 Lsun (Preibisch&Zinnercker1999)
- Tdisk = 15[K]*(R/100au)**(-0.64) (Zhang+2014, gas temperatureのモデルから1/3として算出)