Aide-memoire: 2011

2011年8月31日水曜日

立体角

立体角をなかなか理解してもらえなかったので見つけてきたサイトとか。

----

立体角の比較的分かりやすい定義の説明です。

http://www.earliness-si-unit.com/hozyo_02.html

以下、電磁気の講義ノートetc.に現れる立体角の説明です。
いずれも分り易いと思うのでよく読んで下さい。
立体角はガウスの定理を習った時（マクスウェル方程式からクーロンの法則を導出した時）
に「必ず」習っている筈の基本概念なのでよく理解しておいて下さい。

http://www.phys.u-ryukyu.ac.jp/~maeno/cgi-bin/pukiwiki/index.php?%C5%C5%BC%A7%B5%A4%B3%D82007%C7%AF%C5%D9%C2%E84%B2%F3#q9a03723

http://www.erp.sophia.ac.jp/Projects/ocw/faculty/st/stCl/common_st/10SSSCT107E1/10012SCT107E1.pdf

http://hooktail.sub.jp/vectoranalysis/GaussSolidAngle/

2011年7月27日水曜日

2次元 SMM fit

#!/Library/Frameworks/EPD64.framework/Versions/Current/bin/python

SCRIPT_TYPE="InstrumentalityOfMankind" 
SCRIPT_SHELL="python" 
SCRIPT_NAME="Fit_infile_.py"
SCRIPT_HELP_ARG=""
SCRIPT_HELP_SENTENCE=""
SCRIPT_NUM_ARG=0
SCRIPT_VERSION=1.0

###IMPORT###
import SMM
from scipy.optimize import leastsq 
from numpy import *
import sys

###HELP###
if len(sys.argv[1:])!=SCRIPT_NUM_ARG:
    print 'Name: '+SCRIPT_NAME
    print 'Arguments: '+SCRIPT_HELP_ARG
    print 'Explanation: '+SCRIPT_HELP_SENTENCE
    sys.exit()


###MAIN###
def SMM_norm_free(norm_nsc,norm_nsd,norm_gb,norm_gd):
    return lambda l,b: norm_nsc * SMM.SMM(l,b,'NSC') + norm_nsd * SMM.SMM(l,b,'NSD') +  norm_gb * SMM.SMM(l,b,'GB') + norm_gd * SMM.SMM(l,b,'GD')

def SMM_norm_free_array(param,l_array,b_array):
    res=[]
    num=0
    for l in l_array :
        b=b_array[num]
        res.append(SMM_norm_free(*param)(l,b))
        num+=1                
    return res

def fit_SMM_norm_free(param,l_array,b_array,sb_array,sb_err_array):
    errorfunction = lambda p: (array(SMM_norm_free_array(p,l_array,b_array))-array(sb_array))/sb_err_array
    p=leastsq(errorfunction,param,full_output=True)
    return p

def demo(para=[1.0,2.0,3.0,4.0],para0=[1.1,1.2,1.3,1.4]) :
    l_in=linspace(0.0,5.0,10)
    b_in=linspace(0.0,5.0,10)
    data = SMM_norm_free_array(para,l_in, b_in)
    data_err=0.001*array(data)
    param=fit_SMM_norm_free(para0,l_in,b_in,data,data_err)
    
    val=param[0]
    covar=param[1]
    val_err=[]
    for num in range(0,4):
        val_err.append((covar[num][num])**0.5)
    
    print val,val_err
    

def fit_SMM_indata(infile):
    para0=(1.0,1.0,1.0,1.0)
    data=loadtxt(infile)
    l_in=data[:,0]
    b_in=data[:,1]
    sb_in=data[:,2]
    sb_in_err=data[:,3]
    param=fit_SMM_norm_free(para0,l_in,b_in,sb_in,sb_in_err)
    
    
    val=param[0]
    covar=param[1]
    val_err=[]
    for num in range(0,4):
        val_err.append((covar[num][num])**0.5)
    
    print val,val_err



fit_SMM_indata(str(sys.argv[1]))

2011年7月25日月曜日

XISとHXDのNormが食い違う理由についての考察?のようなもの。
後輩から質問を受けたのでまとめた。
間違えているかもしれないので注意。指摘していただけると助かります。

-------

以下、簡単に各XIS (0,1,3) とHXDのNormが食い違う理由になりうる較正上の
不定性について簡単に書いておきます。

1. XISのnon-X-ray backgroundの不定性
宇宙線由来のバックグラウンドは観測（衛星軌道）毎・検出器毎、更に検出器上の
場所毎に異なります。XIS 0,1,3間のnormalizationが観測毎に異なるのは主にこのNXBの
揺らぎによる筈です。暗い天体のハードな帯域 (5-10 keVあたり) のfluxはこの不
定性の影響を大きく受けます。 BGを視野内の暗い所から引くという解析では、
この不定性の影響が大きくなります。 xisnxbgenというftoolを使い、「(1) SRCとBG
の領域毎にNXBのスペクトルを作り、(2) SRCとBGの各々NXBを差し引いた上で更に
SRCからBGを引く」といったことを行えば、NXBの不定性は小さくなりますが少し
統計誤差が大きくなります。

2. XISのcontaminationの較正の不定性
XISには可視光遮断フィルターがついていますが、この表面にX線を吸収してしまう汚染
物質（ゴム等からアウトガス）が付着しつつあることが知られています。この汚染物質は
観測時期が進むにつれ増加しており、検出器毎、更に検出器上の場所毎に異なります。
この汚染物質量の較正の不定性はXIS 0,1,3間のnormalizationが異なる原因になります。
ソフトな帯域 (0.5-2 keVあたり) のfluxはこの不定性の影響を大きく受けます。

3. XRTの望遠鏡の較正の不定性
XRTの有効面積レスポンス (arf) を計算するモンテカルロシミュレーション
(xissimarfgen) 内では設計図どおりの理想的な望遠鏡の形状が仮定されてい
ます。しかし、実際には打ち上げ (+経年劣化?) によりXRTの形状がいくらか変化し
ていることが分かっており、XIS 0,1,3間、およびPINとのnormalizationの
食い違いの原因になります。この不定性はXIS 0,1,3で異なり、更に光軸（視野）
中心から離れるにつれて大きくなります。

4. 衛星の姿勢の揺れの不定性
観測中に「すざく」衛星の姿勢が揺れ、それが補正できないために天体が
光軸中心からずれてしまい、有効面積の見積もりを間違えてしまうことが
わかっています。XIS 0,1,3間、およびPINとのnormalizationの食い違いの
原因になります。この不定性は光軸（視野）中心から離れるにつれて大き
くなります。これまで割と良く補正されていたのですが、最近になってまた
姿勢の揺れが問題になりつつあるようです。aeattcor2というftoolを使うと少し
改善されるかもしれません。

2011年7月13日水曜日

FinderからFitsイメージをSnow Leopard用ds9 v6.2 X11で左クリックで開けるようにする

v6.1のaqua版を使え、というツッコミは正論すぎるので却下。
v6.1 aquaは何か不具合が生じたのでv6.2を使う。

1. Automaterで「アプリケーション」を作成。
2.「シェルスクリプトを実行」をドラック＆ドロップ
3. 以下みたいな感じ。-fileを付けないでうまく行かず詰まっていた。

4. 保存してApplicationフォルダに突っ込む。
（4+. できたAppを右クリックで「情報」を見て、一番上のアイコンのところにAqua版の”Appそのもの”をドラック＆ドロップするとアイコンを変えることができる）
5. Fitsイメージを右クリックして関連付けを変更

2011年7月2日土曜日

Declination angleとEuler angle2の関係

いつもいつも忘れるが
EA2=90°-Dec

2011年6月15日水曜日

Astro-phに投稿した時にEPSファイルがうまく表示されない時

割とアホだが

1. Previewで開いてjpgで保存
2. jpgを再びPreviewで開いてPDF->PSで保存
3. ps2epsでEPSに変換

でうまくいった。

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 でログのオン、オフ切り替え。

2011年4月12日火曜日

LEMOコネクタ? QLAコネクタ?

LEMOコネクタは商品名で所謂「LEMOコネクタ」以外のLEMOコネクタも存在するので正式名称を調べてみたのだが分からず。QLAコネクタも正式名称というよりは商品名らしいし…うーむ…。

2011年4月11日月曜日

GUIでネットリストを作成 on Mac OS X

GEDAというツールセットのgschemというコマンドを使うと
GUIでネットリストを作成できるようです。

http://ryusai-hp.web.infoseek.co.jp/gaf.htm
下の方に「次のページ　回路図エディタgschem」とかがあるので
順番に進んでみて下さい。

finkでgedaをinstallするには

fink install geda-gaf

で簡単にできます。

2011年1月28日金曜日

XISRMFPARAM

サワDに対するメールより

ae_xi?_rmfparam_*.fits
の中を見るとParam_s*_p* という32行のパラメータ列
が並んでいると思います。これらがレスポンスの
詳細情報（メインピークの幅とメインピークに対
する他コンポーネントの幅・強度）を示しています。
各行の意味については
src/function/xisRespUtil/xisRespUtil.c
を読んで下さい。（これ以外ドキュメントがなさ
そうです・・・）

現在はメインピークの幅のみを変更しています。
Param_s*_p* の
2行目=s12
3行目=s13
4行目=s14
と定義してE(keV)に対して

FWHM(eV)=(S12+S13*E+S14*E**2)**0.5

です。(Param_s*_p* の1行目は1/(3.65*2.35)でPIとSigmaに変換するための係数
です）Eの単位はkeV、計算結果はeVなので注意して下さい。

なぜ、Param_s*_p*が複数あるかというと
s*は各セグメント毎にp*は同一セグメント内の場所(4カ所)毎に
異なったレスポンスを設定できるからです。
現在(ae_xi?_rmfparam_20080901.fits)は同一センサー
内では16列に全て同じ値を詰めています。

2011年1月12日水曜日

Xspecで異なるレスポンスかけてfit

diffuseや暗い点源の解析をやっているとCXBのモデルには異なるレスポンス(arf)をかけたくなることがあるが、Xspecでもそれは出来るようだ。
実際にやってみないと分からないが取りあえずマニュアルの対応箇所をペーストしておく。

変なマークアップが残ったままだけどまあいいか・・・

2011年1月6日木曜日

MITによるSpaced-row Charge Injection (SCI) の本論文

みんな Uchiyama et al. (2007) or (2009) を引いてくれるのを嬉しいのだが、あれはsaw-tooth correctionの論文であって、「すざく」の機上SCI論文としては下記を引くべき。

CCD Charge Injection Structure at Very Small Signal Levels

Prigozhin, Gregory; Burke, Barry; Bautz, Marshall; Kissel, Steve; Lamarr, Beverly
IEEE Transactions on Electron Devices, vol. 55, issue 8, pp. 2111-2120 (2008)

\bibitem[Prigozhin et al.(2008)]{2008ITED...55.2111P} Prigozhin, G., Burke, B., Bautz, M., Kissel, S., Lamarr, B.\ 2008, IEEE Transactions on Electron Devices, 55, 2111

MITの方が主著者で、査読付きのIEEE論文になっている、「筋としては」SCIの機上運用の本論文として引くべき一番の論文。…CIそのものについてが主でSCIの記述は少ないけど。

Mitigating CCD radiation damage with charge injection: first flight results from Suzaku

Bautz, M. W.; LaMarr, B. J.; Miller, E. D.; Kissel, S. E.; Prigozhin, G. Y.; Burke, B. E.; Gregory, J. A.; Uchiyama, H.; Hyodo, Y.; Yamaguchi, H.; Mori, H.; Tsuru, T.; Matsumoto, H.; Koyama, K.; Torii, K.; Katsuda, S.; Hasuike, K.; Nakajima, H.; Hayashida, K.; Tsunemi, H.; Murakami, H.; Ozaki, M.; Dotani, T.

UV, X-Ray, and Gamma-Ray Space Instrumentation for Astronomy XV. Edited by Siegmund, Oswald H. Proceedings of the SPIE, Volume 6686, pp. 66860Q-66860Q-11 (2007)

\bibitem[Bautz et al.(2007)]{2007SPIE.6686E..23B} Bautz, M.~W., LaMarr,
B.~J., Miller, E.~D., et al.\ 2007, \procspie, 6686

「すざく」でのSCI機上運用についての情報、見せたい（見せるべき）図などがまとまったBautz先生のSPIE論文。SCIの機上運用の論文としては内容的に最適。学生さんにも読んでもらいたいところ。…SPIE論文な上、Astro-phにも上がっていないので入手しにくいのが玉に瑕…。

でもUchiyama et al.も引いてくれると嬉しい。
\bibitem[Uchiyama et al.(2009)]{2009PASJ...61S...9U} Uchiyama, H., et al. 2009, \pasj, 61, 9