Aide-memoire: 5月 2011

2011年5月28日土曜日

SciPyで数値積分 scipy.integrate.quad

#coding:sjis

'''
Created on 2011/03/18

@author: uchiyama
'''

"""
importの使い方・考え方は以下が分かりやすい。
http://1-byte.jp/2010/09/29/learning_python5/
そもそもこのシリーズは図解が多くて良さげなので後で読む。
"""
from math import cos,sin,pi, atan, exp, sqrt
from scipy import integrate

"""
普遍であるべき観測量
ここにグローバルな変数として書くのは正しいのだろうか。
"""
#postion of Sgr A* in the Galactic coordinate (degree)
L_GC=-0.056
B_GC=-0.046
# 1 pc in the unit of cm
#CM_PER_PC=3.09e18
#Position of the Earth from the plane (pc)
Z0=16.0
#Distance to Sgr A* (pc)
R0=8.5e3
#The rotation angle 
BD1=15.0
BD2=1.0

"""
座標変換に必要な関数。
本当はscipyの座標変換使えばもっと賢くできる気がする。
"""
def sin_d(degree):
    return sin(2*pi*degree/360.0);

def cos_d(degree):
    return cos(2*pi*degree/360.0);

def dl(l):
    return l-L_GC

def db_bulge(b):
    return b-B_GC-atan(Z0/R0)

def db(b):
    return b-B_GC

def r_f(s,l,b):
    l0=dl(l)
    b0=db(b)
    (R0**2+(s*cos_d(b0))**2-2*R0*s*cos_d(l0)*cos_d(b0))**0.5

def z_f(s, l, b):
    b0=db(b)
    return s*sin_d(b0)

def x_f_bulge(s, l, b):
    l0=dl(l)
    b0=db(b)
    return R0-s*cos_d(b0)*cos_d(l0)

def y_f_bulge(s, l, b):
    l0=dl(l)
    b0=db_bulge(b)
    s*cos_d(b0)*sin_d(l0)

def z_f_bulge(s, l, b):
    b0=db_bulge(b)
    return s*sin_d(b0)+Z0

def X_f(x, y, z):
    return x*cos_d(BD1)+y*sin_d(BD1)

def Y_f(x, y, z):
    return -1*x*sin_d(BD1)*cos_d(BD2)+y*cos_d(BD1)*cos_d(BD2)+z*sin_d(BD2);

def Z_f(x, y, z):
    return -1*x*sin_d(BD1)*sin_d(BD2)+y*cos_d(BD1)*sin_d(BD2)+z*cos_d(BD2);


"""
以下星質量モデル。
基本的にMuno+06に基づく…が数式に誤りたくさん+説明不足。
http://adsabs.harvard.edu/abs/2006ApJS..165..173M
Kent+91, Launhardt+02を元に情報を補足、整理。
http://adsabs.harvard.edu/abs/1991ApJ...378..131K
http://adsabs.harvard.edu/abs/2002A%26A...384..112L
更に詳しくは
http://www-utheal.phys.s.u-tokyo.ac.jp/~uchiyama/wordpress/wp-content/uploads/2011/05/SMM.pdf
"""

def NSC(x,l,b):
    norm1=3.3e6 # M_solar pc^{-3}
    norm2=norm1*8.98e7/3.3e6
    n1=2.0
    n2=3.0
    R_c=0.22 # pc
    R_lim=6.0
    R_cutoff=200.0

    r=r_f(x,l,b)
    z=z_f(x,l,b)
    """
    powや**は（defaultで）使えるのにsqrtはmathをimportしないと使えない。
    意味がよく分からないが、そもそもsqrt不要と言われるとそれが正しい気もする。
    """
    R=sqrt(pow(r,2)+pow(z,2))
    
    """
    元々は3項演算子を使っていたがpythonの3項演算子に慣れずネストの仕方が分からないのでif文で書く。
    """
    if R < R_lim : 
        return norm1/(1+pow((R/R_c),n1))
    elif R < R_cutoff : 
        return norm2/(1+pow((R/R_c),n2)) 
    else: 
        return 0   


def NSD(x,l,b):

    n1=-0.1
    n2=-3.5
    n3=-10.0

    r_1=120.0
    r_2=220.0
    r_3=2000.0
    r_d=1.0
    z_d=45.0
 
    norm1=300.0
    norm2=norm1*(r_1/r_d)**(-n2+n1)
    norm3=norm2*(r_2/r_d)**(-n3+n2)
    
    r=r_f(x,l,b)+1e-30
    z=z_f(x,l,b)
 
    
    if r <=0 :
        return 0.0
    elif r < r_1 :
        return norm1*(abs(r/r_d)**n1)*exp(-1.0*z/z_d)
    elif r < r_2 : 
        return norm2*(abs(r/r_d)**n2)*exp(-1.0*z/z_d)
    elif r < r_3 :
        return norm3*(abs(r/r_d)**n3)*exp(-1.0*z/z_d) 
    else:
        return 0.0
    


def GB(x,l,b):
        """
        Unit of value returned by this function
        M_solar pc^{-3}
        """

        #norm Unit: M_solar pc^{-3}
        norm=8.0

        #a_x,a_y,a_z : pc
        a_x=1100.0
        a_y=360.0
        a_z=220.0

        c_para=3.2
        c_perp=1.6      
        
        x_b=x_f_bulge(x, l, b);
        y=y_f_bulge(x, l, b);
        z=z_f_bulge(x, l, b);

        """
        Pythonはべき乘を表す演算子**もあるけどpowもある。
        そして小数乘もできるし、負の数を小数乗して複素数を
        返すこともできる。
        """
        R_perp=((abs(X_f(x_b,y,z))/a_x)**c_perp+(abs(Y_f(x_b,y,z))/a_y)**c_perp)**(1/c_perp) 
        Rs=(pow(abs(Z_f(x_b,y,z))/a_z,c_para)+pow(R_perp,c_para))**(1/c_para);

        return  norm*exp(-1*Rs);
    
    
def GD(x,l,b):
    """
    Unit of this function
    M_solar pc^{-3}
    """
    #norm Unit: M_solar pc^{-3}
    norm=5.0

    #r_d, z_d : pc
    r_d=2700.0;
    z_d=200.0;

    r=r_f(x, l, b)
    z=z_f(x, l, b)

    return  norm*exp(-1.0*r/r_d)*exp(-1.0*(z/z_d))


def Stellar_density_f(x, l, b):
        """
        Unit of this function
        M_solar pc^{-3}
        """
        return NSC(x,l,b)+NSD(x,l,b)+GB(x,l,b)+GD(x,l,b)
    

def SMM(l, b, s=16e3 , region='SMM'):

    """
    Pythonで関数はオブジェクトなので下記の様なことができてしまう。
    正直キモイ。ポインタを渡したくなる…。以下を参照。
    http://python.g.hatena.ne.jp/muscovyduck/20080707/p2
    """
    if region=='NSC':
        func=NSC
    elif region=='NSD':
        func=NSD
    elif region=='GD':
        func=GD
    elif region=='GB':
        func=GB
    elif region=='SMM':
        func=SMM
    else : 
        exit("3rd argument must be one of NSC/NSD/GD/GB/SMM.")

    """
    SciPy(さいぱい。すしぱいと読んでいた…)を使った積分。
    http://docs.scipy.org/doc/scipy/scipy-ref.pdf
    の1.4.1を見る。
    lambda式の意図は
    http://atkonn.blogspot.com/2008/02/python-python27-lambda.html
    をみるとなんとなく分かる。
    """
    
    res=integrate.quad(lambda x: func(x,l,b), 0.0, s)[0]
    
    """
    The unit of returned value is M_solar/pc^2/str
    """
    return res/(4.0*pi);

2011年5月25日水曜日

Tabを4スペースに変換

アホみたいだが

sed 's/(TABキー出力)/(4スペース)/g' $1

と書いたスクリプトを使う。

追記: Eclipseなら「ソース→タブをスペースに変換」でできる。

2011年5月24日火曜日

BloggerでPythonソースのシンタックスハイライト

クリボウの Blogger Tips: コードをハイライトする「Code Prettify」ウィジェット
上記の方のウィジェットでできます。

<pre class="prettyprint" id="python">
#coding:utf-8
import pylab
print "Hello, world!"
x=pylab.arange(0.0,1.0,0.01)
y=pylab.sin(2*pylab.pi*x)
figplot=pylab.subplot(111)
figplot.plot(x,y)
pylab.show()
</pre>

んで

#coding:utf-8
import pylab
print "Hello, world!"
x=pylab.arange(0.0,1.0,0.01)
y=pylab.sin(2*pylab.pi*x)
figplot=pylab.subplot(111)
figplot.plot(x,y)
pylab.show()

Guplot.pyでデータの図をプロット

Gnuplot.pyをプロッタとして使用。
使い慣れている分、確かに使い易い。

#coding:utf-8

"""
Gnuplot.pyが何故かeclipseで補完されない。
ライブラリに手動で設定したのだが。
何故か gnuplot_Suites が出てくる。
"""
import Gnuplot
import numpy as np

in_file='Data_all.txt'
in_txt_name='Data_GalacticPlane.txt'
out_file='Gnuplot_test.ps'
out_file2='Gnuplot_test2.ps'

"""
Pythonらしい?方法でファイル読み込み
"""
input_file=open(in_file,"r")
data_list=[]
Fe=[]
l=[]
line_data_cols=[4,7,10]
line_ids=[0,3,6]

for id in range(9):
    Fe.append([])

for input_line in input_file:
    if input_line.strip().split(',')[14]==' l':
        l.append(float(input_line.strip().split(',')[2]))
        for line_col_num in line_data_cols:
            """
            line_id:
            0-2 Fe I
            3-5 Fe XXV
            6-8 Fe XXVI 
            """
            line_id=line_col_num-4
      
            value=float(input_line.strip().split(',')[line_col_num])
            lower=float(input_line.strip().split(',')[line_col_num+1])
            upper=float(input_line.strip().split(',')[line_col_num+2])
            Fe[line_id].append(value*1e-7)
            Fe[line_id+1].append(lower*1e-7)
            Fe[line_id+2].append(upper*1e-7)

"""
numpyのloadtxtを使ってもう少し賢くやる
http://d.hatena.ne.jp/y_n_c/20091117/1258423212
を参考にする。
便利だが今回のようにデータを各行毎にチェックした上でplotしなくては
ならない場合は向かないのかもしれない。
l=-0.05のデータのみ抜き出したファイルを準備しておく。
"""
l_txt=np.loadtxt(in_txt_name, delimiter=',', usecols=[3],unpack=True)

"""
(np.)arrayとlistのどちらを使うべきかでよく混乱するが「自分の印象としては」
arrayは静的な数字そのものでlistは動的な変数。
・arrayにはappendはないけど、sin(array)みたいなことをやって、更にarrayを作れる。
・array*2=array[2*a,2*b]だが[a,b]*2=[a,b,a,b]
・そもそもarrayには数字以外は入らない。
"""

Fe_txt=[]
Fe_max_txt=[]
Fe_min_txt=[]

for id in range(3):
    Fe_txt.append([])
    Fe_max_txt.append([])
    Fe_min_txt.append([])

l_txt=np.loadtxt(in_txt_name, delimiter=',', usecols=[2],unpack=True)
for num in range(3):
    Fe_txt[num],Fe_min_txt[num],Fe_max_txt[num]=np.loadtxt(in_txt_name, delimiter=',', usecols=[3*num+4,3*num+5,3*num+6],unpack=True)


"""
Gnuplot.pyでプロット。
・Gnuplot.pyの正式なマニュアル
http://gnuplot-py.sourceforge.net/doc/
http://gnuplot-py.sourceforge.net/doc/Gnuplot/_Gnuplot/Gnuplot.html

http://www-acc.kek.jp/EPICS_Gr/products/gnuplot/Gnuplot_py.pdf
が日本語で詳しくて有り難い。
・withはpythonでは予約語なのでwith_を使う。
・Gnuplot.pyを使ってfittingも出来なくは無さそうだが、
今一つ頭が悪い感じ。少なくともpython内で定義した関数で
fitは出来そうにない。
http://osdir.com/ml/python.ipython.user/2003-07/msg00009.html
"""

"""
Python I/Oで読んだデータをプロット。
"""

gp1=Gnuplot.Gnuplot()

lw='3'
ps='1.5'

data1=[]
for num in range(3):
    data1.append([])    

for num in range(3):
    if num==0:
        line_name='Fe I K{/Symbol a}'
        pt='7'
        lt='0'
    elif num==1:
        line_name='Fe XXV K{/Symbol a}'
        pt='6'
        lt='1'
    elif num==2:
        line_name='Fe XXVI K{/Symbol a}'
        pt='9'
        lt='3'        
        
     
    """
    何故か titleだけ with_でなくてtitleでGnuplot.Dataに渡す。
    """
    with_option='e lw '+lw+' pt '+pt
    data1[num]=Gnuplot.Data(l,Fe[3*num],Fe[3*num+1],Fe[3*num+2],with_=with_option,title=line_name)

gp1('set key')
gp1('set log y')
gp1('set te unknown')
gp1.xlabel('Galactic longitude (degree)')
gp1.ylabel('Intensity (ph^{-1} cm^{-2} arcmin^{-2} s^{-1})')
gp1.title('Fe K{/Symbol a} profile along the Galactic plane')

gp1.plot(data1[0], xrange='[2.2:-3.2]', yrange='[1e-8:1e-5]')

for num in [1,2]:
    gp1.replot(data1[num])

gp1('set format y "10^{%L}"')
gp1('set te po co solid enhanced "" 24')
gp1('set ou '+'"'+out_file+'"')
gp1.replot()
gp1('set ou')

"""
loadtxtで読んだデータをプロット。
"""

data2=[]
for num in range(3):
    data2.append([])    

for num in range(3):
    if num==0:
        line_name='Fe I K{/Symbol a}'
        pt='7'
        lt='0'
    elif num==1:
        line_name='Fe XXV K{/Symbol a}'
        pt='6'
        lt='1'
    elif num==2:
        line_name='Fe XXVI K{/Symbol a}'
        pt='9'
        lt='3'   
        
    with_option='e lw '+lw+' pt '+pt
    data2[num]=Gnuplot.Data(l_txt,Fe_txt[num],Fe_min_txt[num],Fe_max_txt[num],with_=with_option,title=line_name)

gp1('set te unknown')
gp1.plot(data2[0], xrange='[2.2:-3.2]', yrange='[1e-1:1e2]')
for num in [1,2]:
    gp1.replot(data2[num])

gp1.ylabel('Intensity (10^{-7} ph^{-1} cm^{-2} arcmin^{-2} s^{-1})')
gp1('set te po mono dashed enhanced "" 24')
gp1('set ou '+'"'+out_file2+'"')
gp1.replot()
gp1('set ou')

出来た図は下。割と綺麗。

impotしたモジュールの元ファイル位置（インストールディレクトリ）を表示

__file__を使う。

In [515]: Gnuplot.__file__
Out[515]: '/Library/Frameworks/EPD64.framework/Versions/7.0/lib/python2.7/site-packages/Gnuplot/__init__.pyc'

のような感じ。

Pythonでモジュール一覧

自分世界やってみるのが第一歩::Pythonでインストール済みモジュール一覧を表示

help('modules')

2011年5月22日日曜日

Pylabでデータの図をプロット

要はGnuplot+bashスクリプトからの卒業を目指して。

#coding:utf-8

import pylab

input_file=open("../../Data_all.txt","r")
data_list=[]
Fe=[]
l=[]
line_data_cols=[4,7,10]
line_ids=[0,3,6]


for id in range(9):
    Fe.append([])

for input_line in input_file:
    if input_line.strip().split(',')[14]==' l':
        l.append(float(input_line.strip().split(',')[2]))
        for line_col_num in line_data_cols:
            """
            line_id:
            0-2 Fe I
            3-5 Fe XXV
            6-8 Fe XXVI 
            """
            line_id=line_col_num-4
      
            value=float(input_line.strip().split(',')[line_col_num])
            lower=float(input_line.strip().split(',')[line_col_num+1])
            upper=float(input_line.strip().split(',')[line_col_num+2])
            Fe[line_id].append(value*1e-7)
            Fe[line_id+1].append((value-lower)*1e-7)
            Fe[line_id+2].append((upper-value)*1e-7)




"""
Fontの設定の仕方は
http://www.scipy.org/Cookbook/Matplotlib/LaTeX_Examples
(rcParams.updateの使い方)
具体的にパラメータとして何が使えるかは
http://matplotlib.sourceforge.net/users/customizing.html#a-sample-matplotlibrc-file
とかを見る。
"""
fparams = {'backend': 'ps',
          'axes.labelsize': 20,
          'suptitle.fontsize': 20,
          'legend.fontsize': 20,
          'xtick.labelsize': 18,
          'ytick.labelsize': 18,
          'text.usetex': True,
          'font.family': 'serif',
          'font.serif': 'Times New Roman',
          'text.usetex' : True}

pylab.rcParams.update(fparams)

l_dist=pylab.subplot(111)


"""
ログスケールの図の作り方は下記。
http://matplotlib.sourceforge.net/examples/pylab_examples/log_demo.html#pylab-examples-example-code-log-demo-py
"""
l_dist.set_yscale("log", nonposy='clip')
l_dist.set_ylim(ymin=1e-8,ymax=1e-5)
l_dist.set_xlim(xmin=2.2,xmax=-3.3)



"""
TeXの書式を使う方法は下記。
http://www.scipy.org/Cookbook/Matplotlib/UsingTex
"""
for line_id in line_ids:
    if line_id==0:
        title=r'Fe I K$\alpha$'
    elif line_id==3:
        title=r'Fe XXV K$\alpha$'
    elif line_id==6:
        title=r'Fe XXVI K$\alpha$'



    """
    エラーのついた図の作り方は下記。
    http://matplotlib.sourceforge.net/api/pyplot_api.html#matplotlib.pyplot.errorbar
    """
    l_dist.errorbar(l,Fe[line_id], yerr=[Fe[line_id+1],Fe[line_id+2]], fmt='.', label=title, elinewidth=3 )

"""
凡例(Gnuplotで言うところのKey)の作り方は下記。
http://matplotlib.sourceforge.net/api/artist_api.html#matplotlib.legend.Legend
handletextpadは凡例中のマークと文章の間のスペース
"""
l_dist.legend(numpoints=1,frameon=False,markerscale=1.5,handletextpad=-0.5)

"""
suptitleのみ fontsizeが必要な理由が自分には意味不明。
"""
pylab.suptitle(r'Profile of the Fe K$\alpha$ lines along the Galactic longitude ($b=-0^{\circ}.05$)',size=20)
pylab.xlabel(r'Galactic longitude $l$ (degree)')
pylab.ylabel(r'Intensity (ph s$^{-1}$ cm$^{-2}$ arcmin$^{-2}$ )')

"""
表示した図を消さないためのpause
"""
raw_input()

pylab.savefig('GCXE.eps', format='eps')

読んだデータファイルはこれ。
出来た図は下。