シミュレーションでの星形成の扱い
基本
ある銀河の中で生まれる星の量をシミュレーションするには，星の元になる\(H_2\)分子の量( \(f_{H_2}\))と，そのgasが星に変換される割合を計算する必要がある．
\(H_2\)はdustを触媒として生成され，UVの輻射場によって破壊されるが，self-shieldingの効果などを考慮しなくてはならない．
zoom-in シミュレーションでgas densityやmetallicity, UV backgroundなどの関数として\(f_{H_2}\)を計算し，他のシミュレーションに組み込んでいる．
星形成の基本はKennicutt-Schmidt law（星形成面密度とgas面密度の関係式）から来ている．
numerical simulationの場合
基本的なレシピは\(\dot{\rho_{\star}} = \frac{\epsilon_{\star} \rho_{gas}}{t_{ff}} \propto \rho^{1.5}_{gas}\)という式で星形成密度がgas密度の1.5乗に比例するというものである．これはfree fall timeが\(\rho^{-0.5}\)に比例するため．
密度がある閾値を超えると星が形成されるとする．閾値はfree parameterである．
しかし，disk銀河のシミュレーションでこの単純なモデルを使うと簡単にgasがcollopseしてしまい，clumpyな銀河になってしまうという問題がある（おそらくgridの荒さの問題）．人工的に圧力を高めることでこの問題は解決されているようである．
セミアナの場合
上の式を簡単にした\(\dot{m_{\star}} = \epsilon_{\star} \)\( \frac{m_{cold}}{t_{\star}}\)が基本的な式．\(\dot{m_{\star}} \)はSFR，\(m_{cold}\)はcold gas mass である．
最近のセミアナでは分子，原子を区別しているものもあり，\(H_2\)が星形成に寄与する．
参考文献
Somerville and Dave 2015 3.1章
Sextractorをインストールする
ずっとmatplotlibについて書いてきたのでたまには天文らしい記事にします笑。
天体検出ソフトのSextractorのインストールの仕方についてです。
Sextractorはfits画像から天体を検出して測光、簡単なサイズ測定なんかをやってくれる便利なソフトです。
HSTの画像に走らせるとこんな感じで検出してくれます。
インストール手順
Sextractorのページ(http://www.astromatic.net/software/sextractor)に行ってDownloads the latest stable versionのsourceから最新版をダウンロードできます。現在v2.19.5が最新ですが、自分の環境では後でmakeした際にエラーを吐かれました汗。なので現在v2.5.0を入れてます。うまくインストールできなかったら古いバージョンを試してみてください。
さてダウンロードしたファイルを展開するとsextractor-2.5.0という名前のディレクトリになるはずなので、これを適当なディレクトに置きます。どこでもいいです。
後はターミナルでsextractor-2.5.0のディレクトリに入って、
./configure
make
sudo make install
とコマンドを打ちます。
エラーがでなければこれでインストール完了です。
うまくインストールできているか確認するにはターミナルでsexと打ってみて下のような画面がでるか確認してみてください。ちゃんと表示されればおしまいです。お疲れ様でした。Sextractorの使い方はそのうち書きます。
あれ？表示されてるバージョン違うような・・・。
matplotlibでヒストグラムの作成
matplotlibでヒストグラムを作ってみたのでメモ。
matplotlibを使えばヒストグラムが非常に簡単に作れます。gnuplotよりもおすすめです。
plt.histに配列を与えればそれだけでヒストグラムにしてくれます。
まず乱数で平均１、分散１のガウシアンの分布を発生させてみます。最近知ったんですがpythonのrandomはメルセンヌツイスタで発生させてくれているそうです。さすが気が利いてますね(Cのrandとは・・・）。
最も簡単にやろうと思えばplt.hist(配列,binの数)を指定するだけで書けますが一応オプションの説明をします。
bins : ビンの数
range : ビンの範囲（デフォルトは配列の最小値と最大値)
alpha : 透過率
weight : 重み(第一引数と同じサイズの配列を与えて、重みを付けます。ヒストグラムなので通常１の重みを持つわけですが、配列に１以外を入れることで重み付けすることが可能です。詳しくは下。）
normed : Trueにすると合計1に規格化してくれる
import numpy as np import matplotlib.pyplot as plt from pylab import * from math import * from scipy.stats import norm nor = np.random.normal(0,1,500) #ガウシアンに従う乱数500個発生 x = np.arange(-5,5,0.01) plt.xlabel("x") plt.ylabel("N") plt.title("平均0,分散1の正規分布に従って生成した乱数のヒストグラム") plt.hist(nor,bins=20,range=(-4,4)) plt.plot(x,norm.pdf(x)*200,lw=3,color="r") plt.show()
確率密度関数に従う乱数発生法についてはそのうち書きます。
これで終わりなんですがオプションについてちょっと解説。
normedについて
これは便利です。自動で規格化してくれます。自分で底面積がいくつで・・・とか考えなくてよいです。
import numpy as np import matplotlib.pyplot as plt from pylab import * from math import * from scipy.stats import norm nor = np.random.normal(0,1,500) #ガウシアンに従う乱数500個発生 x = np.arange(-5,5,0.01) plt.xlabel("x") plt.ylabel("N") plt.title("平均0,分散1の正規分布に従って生成した乱数のヒストグラム") plt.hist(nor,bins=20,range=(-4,4),normed=True) plt.plot(x,norm.pdf(x),lw=3,color="r") plt.show()
みごと規格化されました。確率密度関数とぴったりです。
weightについて
[0.5,1.5,2.5]の配列を1刻みでヒストグラムにすると、普通はすべてのbinで高さ1になりますが、weightの配列を[1,2,3]とか与えてやると、高さが1,2,3になります。
import numpy as np import matplotlib.pyplot as plt from pylab import * from math import * from scipy.stats import norm nor = np.arange(1.0,10,0.1) #ヒストグラムにする配列 wht = np.arange(1.0,10,0.1) #weightの配列 plt.xlabel("x") plt.ylabel("N") plt.hist(nor,bins=9,weights=wht,label="weight") #weightあり plt.hist(nor,bins=9,weights=None,label="weight None") #weightなし plt.legend(shadow=True,numpoints=1,prop={"size":15},loc=2) plt.show() [f:id:cranethree:20150729192449p:plain]
matplotlibで領域の図示
matplotlibで領域に色を付けてみたのでメモ。
fill_between関数fill関数を使うと大体できます。
sin(x)とcos(x)の間の領域に色を付ける。
まずはsin(x)とcos(x)の間に色を付けてみます。
fill_betweenは3つの変数を取ります。
x: x軸の範囲指定
y1,y2: この2つの間を塗りつぶす
あとはオプションでfacecolorは色、alphaは透過度なんかを指定できます。
import numpy as np import matplotlib.pyplot as plt from math import * x = np.arange(0,10,0.01) y1 = np.sin(x) y2 = np.cos(x) plt.fill_between(x,y1,y2,facecolor='y',alpha=0.5) plt.show()
条件を指定して図示
図示の条件がついた場合は少し面倒になります。
1条件の場合(y1>y2)などでいい場合はオプションwhere=y1>y2で追加すればいいですが、2条件(y1>y2 かつ xの範囲指定)になるとwhereでは書けません。
下のコードではグラフを書く部分と色を塗る部分を完全に分けています。
import numpy as np import matplotlib.pyplot as plt from math import * x = np.arange(0,10,0.01) y1 = np.sin(x) y2 = np.cos(x) plt.plot(x,y1,color="k") plt.plot(x,y2,color="k") x1 = np.arange(0,5,0.01) z1 = np.sin(x1) z2 = np.cos(x1) plt.fill_between(x1,z1,z2,where=z1>z2,facecolor='y',alpha=0.5) plt.show()
さらに複雑な場合
3つ以上の関数で領域に色を塗る場合にはminimum関数やmaximum関数を使います。(minimum,maximumの使用は(from pylab import *)が必要です。
\[y > 5x,\ \ y < x,\ \ y < -x^2 + 5\]
この条件を満たす領域に色を塗ってみます。
import numpy as np import matplotlib.pyplot as plt from math import * from pylab import * x1 = np.arange(0, 2, 0.001) y2 = - x1 ** 2 + 5 y3 = x1 y4 = 5 * x1 y5 = minimum(y2,y4) plt.fill_between(x1,y3,y5,where=y5>y3,facecolor='y',alpha=0.5) plt.plot(x1,y2,color="k") plt.plot(x1,y3,color="k") plt.plot(x1,y4,color="k") plt.show()
ちょっとごちゃごちゃしてますがこんな感じでしょうか。
おまけ
座標が分かっていてその内部を塗りつぶす場合はfill関数が楽です。長方形を図示したいような場合はこっちがいいです。
例(1,0),(2,0.3),(2,1),(1.4,0.7),(1,1)の内部の塗りつぶし。
import matplotlib.pyplot as plt testx = [1,2,2,1.4,1] testy = [0,0.3,1,0.7,1] plt.fill(testx,testy,color="y",alpha=0.5) plt.show()
はてなブログでフォントの変更
前々からフォントが読みにくい（特に数字）と思っていたので、はてなブログでフォントの変更をしてみました。
トミー (id:tosssaurus)さんの記事を参考にさせていただきました。
http://itkeijyoshilog.hatenablog.com/entry/2015/02/16/201000
トミーさんありがとうございます。
フォントの変更は意外と簡単です。
デザイン→カスタマイズ→デザインCSSといってコードを貼るだけです。
初めてのmatplotlib(対数グラフ)
グラフを対数にするにはplt.xscale("log"),plt.yscale("log")というコマンドを使用します。
例としてx軸が対数で直線のグラフを描いてみましょう。
\[ y = log(x)\]をプロットしてみましょう。
ちなみにnumpyで作ったリストに対して数学関数をつかうにはnp.~とする必要があります。ここではnumpyで作ったxの配列に対してlogを取りたいのでnp.log(x)としています。
# -*- coding: utf-8 -*- #x軸がlog-scaleでy=log(x)のグラフをプロットする import numpy as np import matplotlib.pyplot as plt plt.xscale("log") x = np.arange(0.1,10,0.01) y = np.log(x) plt.plot(x,y) plt.show()
y軸もlogにしたい場合は、plt.yscale("log")を追加して下さい。
また、グリッド線を入れたい場合は、plt.grid(which="both")とすると下の図のようにできます。
