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')

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

しかし、pylabを使ってしまっているせいで結局どのオブジェクトについて扱っているのか、分かりにくくなってしまっている。
このサイトを参考に書きなおしてみる。
あと、ここをみると使いなれたGnuplotを使うのも悪くなさそう。
しかし、一番の目的であるFittingに到達する前に挫折。やれやれ。

Pylabの凡例 (legend)について

matplotlib artists — Matplotlib v1.0.1 documentation

defaultがnumpoints=2(1でない)なことにキレそうになった時にどうぞ。

2011年5月21日土曜日

Pylabでlog scaleの使い方

pylab_examples example code: log_demo.py — Matplotlib v1.0.1 documentation

2011年5月20日金曜日

ファイルからデータを読んでlistに突っ込む

for line in test_file:
temp_list = float(line.strip().split(',')[0])
test_list.append(temp_list)

Python で型を調べる

Python で型を調べる – types モジュール | すぐに忘れる脳みそのためのメモ

type(object)

でオブジェクトの型を返してくれる。

サンプルコード勉強中等は激しく便利。

良くも悪くも宣言しない時は何やっているのか訳が分からなくなる。

Python/インタラクティブシェルでの操作をファイルに保存する方法

Python/インタラクティブシェルでの操作をファイルに保存する方法 - TOBY SOFT wiki:
より

①前に入力した分も含めてまるごと保存

import readline
readline.write_history_file('test.py')

これはこれで便利

②IPythonのログ機能を使う方法
XSPCEライク?にデータを保存して行く方法

%logstartであらかじめログファイルを指定

%logstart 'test.py'
%logstopでログ記録終了

%logstop
%logon,%logoff でログのオン、オフ切り替え。

Aide-memoire

2011年5月28日土曜日

SciPyで数値積分 scipy.integrate.quad

2011年5月25日水曜日

Tabを4スペースに変換

2011年5月24日火曜日

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

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

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

Pythonでモジュール一覧

2011年5月22日日曜日

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

Pylabの凡例 (legend)について

2011年5月21日土曜日

Pylabでlog scaleの使い方

2011年5月20日金曜日

ファイルからデータを読んでlistに突っ込む

Python で型を調べる

Python/インタラクティブシェルでの操作をファイルに保存する方法

プロフィール

ブログアーカイブ

ラベル

このブログを検索

2011年5月28日土曜日

2011年5月25日水曜日

2011年5月24日火曜日

2011年5月22日日曜日

2011年5月21日土曜日

2011年5月20日金曜日

プロフィール

ブログ アーカイブ

ラベル

このブログを検索

ブログアーカイブ