Aide-memoire: Python

ラベル Python の投稿を表示しています。すべての投稿を表示

2017年12月7日木曜日

Pyorbitalで方位仰角、Rangeの計算

Pyorbitalで方位仰角、Rangeの計算。こんなのがさっくりできるのは、とてもありがたい...。どのくらい信用できるのかはわからないけど。

from pyorbital import orbital　＃これを追加すると動く（外すと動かない）のは何故
from pyorbital import orbital
import pyorbital
import numpy
import datetime

#Sat info
satname='ITUPSAT 1'
f=437.325
tlef='TLE.txt'

#Position Info
#静岡大学
lon=138.431704
lat=34.965720
alt=90

#Obs Info
start_jst=datetime.datetime(2018,02,16,15,25,00)
end_jst=datetime.datetime(2018,02,16,15,45,00)

#Conf
c=299792.458 # km/s
dt=datetime.timedelta(0,30)

def jst2utc(time) :
    return time-datetime.timedelta(0,9.*60*60,0)

def obsf_ICR6(obsf) :
    f0=obsf*1000.
    f1=int(obsf*100.)*10
    df=f0-f1
    if df <2.5 :
        return f1/1000.
    elif df <7.5 : 
        return (f1+5.)/1000. 
    else :
        return (f1+10.)/1000. 
    
tle=pyorbital.tlefile.read(satname,tlef)
orb=pyorbital.orbital.Orbital(satname,tlef)

print "TLE Epoch:",",",str(datetime.datetime(2000+int(tle.epoch_year),1,1)+datetime.timedelta(days=tle.epoch_day)).split(' ')[0]
print "Date:",",",str(start_jst).split(' ')[0]
print "Lon:",",",lon
print "Lat:",",",lat
print "Alt:",",",alt
print ""

jst=start_jst
while jst < end_jst :

    now=jst2utc(jst)
    obs_pos=numpy.array(pyorbital.astronomy.observer_position(now,lon,lat,alt)[0])
    sat_pos,sat_vel=orb.get_position(now,False)

    range_dis=numpy.linalg.norm(sat_pos-obs_pos)
    range_vec=(sat_pos-obs_pos)/range_dis
    rang_rate=numpy.dot(sat_vel,range_vec)

    Az,El=orb.get_observer_look(now,lon,lat,alt)
    obsf=f*(1-rang_rate/c)
    
    print str(jst).split(' ')[1],",",Az,",",El,",",obsf_ICR6(obsf),",",obsf
    jst+=dt

2016年3月31日木曜日

Arf_BG_corrector.py

お察し下さい (pyfits練習) 第2弾。前よりは多少ファイルI/Oの書き方が賢くなっている...かも。あくまで個人的に作ったスクリプト。

SCRIPT_SHELL="python" 
SCRIPT_NAME="Arf_BG_corrector.py"
SCRIPT_HELP_ARG="1) Input Rate Bg PI file made by mathpha 2) Src Arf 3) Bg Arf 4) Output Bg PI file"
SCRIPT_HELP_SENTENCE="Modiefy Bg PI Channel with Src / Bg  arf ratio."
SCRIPT_NUM_ARG=4

###IMPORT###
import sys
import pyfits
import os
###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###
ev_per_ch=3.65

bg_pif=str(sys.argv[1])
src_arf=str(sys.argv[2])
bg_arf=str(sys.argv[3])
out_pif=str(sys.argv[4])

def main() :
    
    
    if os.path.isfile(out_pif) :
        os.remove(out_pif)
    
    data_src_arf = (pyfits.open(src_arf))['SPECRESP'].data
    data_bg_arf = (pyfits.open(bg_arf))['SPECRESP'].data
   
    ratio=[] 
    for n in range(len(data_src_arf)) :
        ratio.append([data_src_arf['ENERG_LO'][n],data_src_arf['ENERG_HI'][n],data_src_arf['SPECRESP'][n] / data_bg_arf['SPECRESP'][n]])
    
    hdulist=pyfits.open(bg_pif)
    
    sp=hdulist['SPECTRUM'].data
    ch_num=len(hdulist['SPECTRUM'].data)
    
    hdulist['SPECTRUM'].header['COMMENT']='Modefied with Arf_BG_corrector.py'
    hdulist['SPECTRUM'].header['COMMENT']='Src arf : '+src_arf
    hdulist['SPECTRUM'].header['COMMENT']='Bg arf : '+bg_arf
    hdulist['SPECTRUM'].header['COMMENT']='Bg pi : '+bg_pif

        
    for n in xrange(ch_num) :
        enegy=ev_per_ch*n/1e3
        for erange in ratio :
            if erange[0] <= enegy and enegy < erange[1] :
                sp['RATE'][n]=sp['RATE'][n]*erange[2]
                sp['STAT_ERR'][n]=sp['STAT_ERR'][n]*erange[2]
                break

    hdulist.writeto(out_pif)

main()

2016年3月27日日曜日

openpyxlを使ってPythonでExcelファイルを読み書きする練習

openpyxlを使ってPythonでExcelファイルを読み書きする練習。

###IMPORT###

import sys
import openpyxl
from openpyxl.cell import get_column_letter

###MAIN###

def main() :
    inbook_file="Book1.xlsx"
    csv_file="CSV.csv"
    out_file="Book2.xlsx"
        
    #Read XLSX            
    wb1=openpyxl.load_workbook(inbook_file)
    ws1=wb1.active
    AB={}  
    
    for u in ws1.rows :
        AB[str(u[2].value)]=[u[0].value,u[1].value]
    
    #Read CSV
    AL=[]
    for line in open(csv_file, 'r'):
        AL.append((line.split(","))[0])
    
    #Open WB for Output
    wb2=openpyxl.workbook.Workbook()
    ws2=wb2.active     
    for i,ad in enumerate(AL) :
        ws2['A'+str(i+1)].value=ad
        if AB.has_key(ad) :
            ws2['B'+str(i+1)]=AB[ad][0]
            ws2['C'+str(i+1)]=AB[ad][1]    
        else :
            ws2['B'+str(i+1)]=" "
            ws2['C'+str(i+1)]=" "   
  
    wb2.save(filename=out_file)

main()

2016年2月24日水曜日

PI Shifter for Suzaku XIS

お察し下さい。個人的に作ったスクリプト。pyfitsの練習というか。

SCRIPT_SHELL="python" 
SCRIPT_NAME="PiShift.py"
SCRIPT_HELP_ARG="1) input Count PI file made by mathpha 2)output PI file 3)c0 4)c1 5)c2"
SCRIPT_HELP_SENTENCE="Modify PI Channel with f(x)=c2*x**2+c1*x+c0. Only for Suzaku XIS."
SCRIPT_NUM_ARG=5

###IMPORT###
import sys
import pyfits
import os
###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###
inf=str(sys.argv[1])
outf=str(sys.argv[2])
c0=float(sys.argv[3])
c1=float(sys.argv[4])
c2=float(sys.argv[5])


ch_num=4096

hdulist = pyfits.open(inf)

if os.path.isfile(outf) :
    os.remove(outf)

sp=hdulist[1].data

hdulist[1].header['COMMENT']='Modefied with PiSift.py'
hdulist[1].header['COMMENT']='c0='+str(c0)
hdulist[1].header['COMMENT']='c1='+str(c1)
hdulist[1].header['COMMENT']='c2='+str(c2)

def f(x) :
    x=float(x)
    return c2*x**2+c1*x+c0


def main() :
    new_counts=[]
    for n in range(ch_num) :
        new_counts.append(0.0)  
        
    old_counts=sp['COUNTS']
    
    for n in range(ch_num) :
        new_ch0=f(n)
        new_ch1=f(n+1)
        if new_ch1 < new_ch0 :
            sys.exit('Correction function is not monotone.')
        elif new_ch0 > 0 and new_ch1 < ch_num :
            for m in range(int(new_ch0+1),int(new_ch1)) :
                new_counts[m]+=old_counts[n]/(new_ch1-new_ch0)

            new_counts[int(new_ch0)]+=old_counts[n]/(new_ch1-new_ch0)*(int(new_ch0+1)-new_ch0)
            new_counts[int(new_ch1)]+=old_counts[n]/(new_ch1-new_ch0)*(new_ch1-int(new_ch1))
         
  

    for n in range(ch_num) :
        sp['COUNTS'][n]=int(round(new_counts[n]))
        sp['STAT_ERR'][n]=float(new_counts[n]**0.5)

    hdulist.writeto(outf)

main()

2015年10月18日日曜日

Pythonで文字列と数値を含んだテキストファイルの読み込み

こんなファイルを読みたい。本当はnumpyのloadtxtでdtypeを指定する (or 読みたい列だけ読む) のが賢い気は果てしなくする。

data=map(lambda str:str.split(), open("FeDist_Screened.plt").readlines())

とはいえ、これはメンドくさすぎると感じるのは自分だけ?dtypeどうにかならんのか。 dtype=('f8','f8','S16')とかダメなのか...。

data=np.loadtxt('FeDist_Screened.plt', delimiter=' ', usecols=[0,2,11],skiprows=1,dtype={'names':('l','6.4','id'), 'formats':('f8','f8','S16')})

一応やりたかったことをやってみた結果。awkだと瞬間なんだけど...。

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

###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###

infile="FeDist_Screened.plt"

data=map(lambda str:str.split(), open("FeDist_Screened.plt","r").readlines())

f=1/(3.14*18**2)

enes=['6.4','6.7']
for ene in enes :
 if ene=="6.4" :
     line_num=2
 elif ene=="6.7" :
     line_num=5
 
 outname='Fe'+ene+'Dist_GRXE.txt'
 line_emin_num=line_num+1
 line_emax_num=line_num+2
 
 outfile=open(outname,"w")
 
 for data_row in data[1:] :
     if float(data_row[0])**2>10**2:
            l,b=data_row[0:2]
            ids=data_row[11]+' '+data_row[12]
                val, val_min, val_max=map(float,data_row[line_num:line_num+3]) 
                val=f*val
                val_err=f*(val_max-val_min)*0.5/1.64
                out_str=l+' '+b+' '+str(val)+' '+str(val_err)+' '+ids+'\n'
                outfile.write(out_str)

 outfile.close()

書き込む前にstrにするのとかが果てしなく面倒くさい。結論 : awkで良いじゃん (ヲイ) or 適材適所。

2015年7月7日火曜日

N人のクラスで同じ誕生日の人がM人以上いる確率

学生さんのレポートでの質問に答えるために作った。40人クラスで3人以上同じ誕生日の確率は6.7%くらい。

###IMPORT###

import sys
import random

###MAIN###

def dice():
    return int(random.uniform(0,365))

def main():
    trial=1e5
    count=[]
    same=0
    N=40
    M=3
    
    for day in range(0,365) :
        count.append(0)
    
    for num1 in range(0,int(trial)):
        for day in range(0,365) :
            count[day]=0
    
        for num2 in range(0,N) :
            count[dice()]+=1
    
        for day in range(0,365) :
            if count[day] > (M-1) :
                same+=1
                break
        
    print same/trial

main()

2015年7月1日水曜日

篠本先生のじゃんけんマシン in Python

篠本先生の教科書で昔勉強したパーセプトロンのじゃんけんマシン (ソースはここを参照しました) をPythonで書き直してみた。もし情報処理の授業を持つ事になったら、Pythonで機械学習とか面白いかと思ったのだが、ハードルが高い (&篠本先生の教科書は既に入手困難) か...。 GUIをつけてみたいところ。

###MAIN###
N=5

def perceptron(m,x,w,v):
    prec=[]
    
    if m<=0:
        return -1
                
    for k in range(0,3) :
        prec.append(-1)
        
    prec[m-1]=1
  
    for k in range(0,3):
            if prec[k]*v[k]<=0 :
                for j in range(0,3*N+1) :
                    w[(3*N+1)*k+j]+=prec[k]*x[j]
    
    for i in range(0,3*N-3) :
        x[3*N-1-i]=x[3*N-4-i]
        
    for i in range(0,3):
        x[i]=prec[i]

    for k in range(0,3):
        v[k]=0

    for k in range(0,3):
        for j in range(0,3*N+1):
            v[k]+=w[(3*N+1)*k+j]*x[j]
            
    vmax=-1000000
    for k in range(0,3) :
        if v[k] >=vmax :
            vmax=v[k]
            kmax=k
            
    return kmax+1
    
def main() :
    
    total=0
    v=[]
    fw=[]
    x=[]
    w=[]
    pred=0
    
    for i in range(0,3):
        v.append(0)
        fw.append(0)
        
    for i in range(0,3*N):
        x.append(0)
    x.append(-1)
    
    for i in range(0, 9*N+3):
        w.append(0)
    
    print "x:",x
    print "w:",w
    print "v:",v       
    m=1
    
    while m>0:
        pred=perceptron(m,x,w,v)

    
        while True:
            try:
                m = int(raw_input("1 Gu 2 Choki 3 Pa :"))
                break
            except ValueError:
                print "Try again..."    
        
        if m>3 :
            m=3
        
        print "You",m,"vs Machine",(pred+1)%3+1 
        if pred==m :
            print "Machine win."
            fw[0]+=1
        elif pred%3==m-1:
            print "You win"
            fw[1]+=1
        else:
            print "Draw"
            fw[2]+=1
            
        total+=1        
        print "You :",fw[1], "Machine:", fw[0], "Draw :",fw[2],"Total:", total
        #print "x:",x
        #print "w:",w
        #print "v:",v
     
    return 0

main()

2015年1月8日木曜日

ついカッとして作った

今は反省しています...。...辞書の使い方の練習、とか(汗)。

import random

def main():
    vals=['ST','CON','LK','DEX','IQ','CHR','SPD','MON']
    ch={}
    
    for val in vals:
        ch[val]=int(three_dices())
    
    bonus=cal_bonus(ch)
    print ch,bonus

def cal_bonus(ch) :
    vals=['ST','LK','DEX']
    bonus=0
    for val in vals :
        if ch[val]>12 :
            bonus=bonus+ch[val]-12
        elif ch[val]<9 :
            bonus=bonus+ch[val]-9

    return bonus

def dice():
    random.seed()
    return random.randint(1, 6)

def three_dices():
    s=0
    for num in [0,1,2]:
        s=s+dice()
    return s

main()

2014年12月24日水曜日

ヘルスライスでバーサーク

「『カザンの闘技場』でもらえるヘルスライス(42D)でバーサークするとどの位ヒットが出るのか?」という子供の頃の疑問に答えるスクリプトを誘惑に負けて作成。関数を再帰的に使うことの練習ということで(汗)。

100万回トライして
平均 : 357
標準偏差 : 23
最大 : 501
最小 : 268

非バーサーク時(期待値 147)の2.4倍。

分布でなんか変なヒゲが生えているのは何なのだろう...?

ハイパーバーサークだと1万回トライして
平均 : 499
標準偏差 : 51
最大 : 765
最小 : 348

...仕事します。

import random
import numpy
import matplotlib.pyplot as plt

def normal_dice(dice_num) :
    vals = []
    total_val = 0
    
    random.seed()
    for num in range(0,dice_num):
        vals.append(random.randint(1, 6))
    
    #print vals
    
    total_val = sum(vals)
    
    return total_val


def berserk_dice(dice_num) :
    vals = []
    total_val = 0
    
    random.seed()
    for num in range(0,dice_num):
        vals.append(random.randint(1, 6))
    
    #print vals
    
    total_val = total_val + sum(vals)
   
    for num in range(1,6):
        num_of_same = vals.count(num)
        if num_of_same > 1 :
            #print "Bersrk by %d dices of %d marks" % (num_of_same,num) 
            total_val=total_val+berserk_dice(num_of_same)
    
    return total_val

def hyper_berserk_dice(dice_num) :
    vals = []
    total_val = 0
    
    random.seed()
    for num in range(0,dice_num):
        vals.append(random.randint(1, 6))
    
    #print vals
    
    total_val = total_val + sum(vals)
   
    for num in range(1,6):
        num_of_same = vals.count(num)
        if num_of_same > 1 :
           # print "Bersrk by %d dices of %d marks" % (num_of_same,num) 
            total_val=total_val+hyper_berserk_dice(num_of_same+1)
    
    return total_val
    
def main():
    dice_num = input("Input dice_num :")
    trial = input("Input Num. of Trials :")
    vals=[]

    for num in range(0, trial) :
        vals.append(berserk_dice(dice_num)) 
        if num%10000==0:
            print num

    vals=numpy.array(vals)
    print "Mean : ", numpy.mean(vals)
    print "SD : ", numpy.std(vals)
    print "Max : ", vals.max()
    print "Min : ", vals.min()
   
    plt.rc('font', **{'family': 'serif'})
    fig = plt.figure()
    
    ax = fig.add_subplot(111)
    ax.hist(vals, bins=100)
    ax.set_title('Berserk with %d dices! (%d trials)' % (dice_num, trial))
    ax.set_xlabel('Hits!')
    ax.set_ylabel('Frequency')

    plt.show()    
    return None 
    
main()

2014年12月23日火曜日

中心極限定理をテストするスクリプト

学生さんに宿題で出したので自分でも書いてみた。
Python (matplotlib, numpy) 便利。
Enthought Canopyも便利。

# -*- coding: utf-8 -*-
import random
import numpy
import matplotlib.pyplot as plt

'''
あまり格好は良くないけど、中心極限定理をテストするスクリプト。
歪んだベータ分布から正規分布が出てくる。
以下を参考にした。
http://akiyoko.hatenablog.jp/entry/2013/06/07/213448
http://docs.python.jp/2/library/random.html
'''

def main():
    random_val=[]
    
    print("This is a test of central limit theorem by beta functions.")
    alpha = input("Input alpha :")
     
    beta = input("Input beta :")
    
    N = input("Input N :")
    trial = input("Input Num. of Trials :")

    random.seed()

    for num in range(0, N):
        random_val.append(random.betavariate(alpha,beta))
        
    random_val=numpy.array(random_val)

        
    mean=numpy.mean(random_val)    
    print "mean=%.3f" % mean

    sigma=numpy.std(random_val)
    print "sigma=%.3f" % sigma
    
    vals=[]
    
    for num0 in range(0,trial):
        val=0.0
        for num1 in range(0, N):
            val=val+random.betavariate(alpha,beta)
        vals.append(val) 
                 
    vals=numpy.array(vals)

    plt.rc('font', **{'family': 'serif'})
    fig = plt.figure()
    
    ax0 = fig.add_subplot(311)
    ax0.hist(random_val, bins=100, range=(0, 1), normed=False, facecolor='r', alpha=0.8)
    ax0.set_xlim(0, 1)

    y,x=numpy.histogram(random_val, bins=100, range=(0, 1), normed=False)

    ax0.set_title('Beta distribution')
    ax0.set_xlabel('Random variable')
    ax0.set_ylabel('Frequency')
    ax0.text(0.8, y.max()*0.8,r'alpha=%.2f' % alpha)
    ax0.text(0.8, y.max()*0.7,r'beta=%.2f' % beta)
    ax0.text(0.8, y.max()*0.6,r'N=%d' % N)

    ax1 = fig.add_subplot(313)
    ax1.hist(vals, bins=100, normed=False, facecolor='g', alpha=0.8)
    ax1.set_title('Distribution of sums Trial=%d' % trial)
    ax1.set_xlabel('Sum of random variables')
    ax1.set_ylabel('Frequency')
  
    plt.show()
    
    return None

main()

2013年7月1日月曜日

Twitterで勝手に掛け合いをするBot的な何か

ついカッとなって、週末作った。
今は反省している。 (これをcrontabで定期的に走らせたりすると・・・ゴホゴホ)
勢いでやったからよく覚えていないがこれを使ったっぽい。

#!/Library/Frameworks/Python.framework/Versions/7.2/bin/python
# -*- coding: utf-8 -*-

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

###IMPORT###

import sys
import tweepy
import random
import time

###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###

CONSUMER_KEY = ''
CONSUMER_SECRET = ''
logfile="temp.log"




###ふぁんくしょん###

def Tweets(ack, acs, text) : 
    CONSUMER_KEY = ''
    CONSUMER_SECRET = ''
 
    auth = tweepy.OAuthHandler(CONSUMER_KEY, CONSUMER_SECRET)
    auth.set_access_token(ack, acs)
    api = tweepy.API(auth, api_root='/1.1')
    api.update_status(text)

     for p in tweepy.Cursor(api.user_timeline).items(1):
          lastid=p.id
    print text
    del api  
    return lastid


def Tweets_as_A(text) :
    ACCESS_KEY1 = ''
    ACCESS_SECRET1 = ''
    return Tweets(ACCESS_KEY1, ACCESS_SECRET1, text)


def Tweets_as_B(text) :
    ACCESS_KEY2 = ''
    ACCESS_SECRET2 = '' 
    return Tweets(ACCESS_KEY2, ACCESS_SECRET2, text)


def Reply(ack, acs, text, lastid, to_usr) : 
    CONSUMER_KEY = ''
    CONSUMER_SECRET = ''
 
    auth = tweepy.OAuthHandler(CONSUMER_KEY, CONSUMER_SECRET)
    auth.set_access_token(ack, acs)
    api = tweepy.API(auth, api_root='/1.1')

    name_to_reply = to_usr
    status_id_to_reply = lastid
    text2 = u'@%s %s' % (name_to_reply, text)
    api.update_status(text2, in_reply_to_status_id = status_id_to_reply)
    print text2
    for p in tweepy.Cursor(api.user_timeline).items(1):
          a=p.id
    del api
    return a


def Reply_as_A(text,lastid) :
    ACCESS_KEY1 = ''
    ACCESS_SECRET1 = ''
    return Reply(ACCESS_KEY1, ACCESS_SECRET1, text, lastid, "B")


def Reply_as_B(text, lastid) :
    ACCESS_KEY2 = ''
    ACCESS_SECRET2 = '' 
    return Reply(ACCESS_KEY2, ACCESS_SECRET2, text, lastid, "A")

 

####おだいけってい###

num_wadai=2

f1=open(logfile,"r")
old_nums=[]
for line in f1:
   old_nums.append(line)
f1.close()

prev_num=int(old_nums[-1])

num_now=random.randint(1, num_wadai)
while prev_num == num_now  : 
    num_now=random.randint(1, num_wadai)
    print num_now, prev_num

f1=open(logfile,"w")
old_nums.append(str(num_now)+"\n")
f1.writelines(old_nums)
f1.close()


####おしゃべり#####


if num_now==1 :
    text=u"はーちみつ！"
    sid=Tweets_as_A(text)
 
    text=u"レモン!"
    sid=Reply_as_B(text, sid)
 
 
elif num_now==2 :
    text=u"レモンはー！"
    sid=Tweets_as_B(text)
 
    text=u"ハチミツ！"
    sid=Reply_as_A(text, sid)
 
    text=u"なんでやねん！"
    sid=Reply_as_B(text, sid)

2012年10月10日水曜日

importしてあるモジュール名から、そのモジュールがどのファイルからロードされているのかを調べる方法

http://blog.taikomatsu.com/2008/02/28/python-%E3%83%A2%E3%82%B8%E3%83%A5%E3%83%BC%E3%83%AB%E3%81%AE%E5%A0%B4%E6%89%80/

module名.__file__

に書いてあるそうです

2012年4月28日土曜日

Shellスクリプトからの卒業?を目指してPythonでXSPECの解析用スクリプトを書いてみる。
サワDに教えてもらった方法を使用。

やりたいことは「たくさんのRegionのスペクトルをビミョーに条件を変えながらモデルFit、エラーを計算」。読み込んでいるModel_Basic.xcmはこんな感じ。

やってみた感じとしては、（便利なところも多いが）今回の様に「たくさんのスペクトルをモデルFit、エラー計算」するだけなら、Bashスクリプトで十分（かつ便利）な気がする。
numpyとかscipyとか連携させたい時は良いけど、xspecを操作する部分までpythonでやる必要があるかは謎。

xsp.stdin.write( Command_Fit )とかは、もはや「お約束」なので、関数&モジュール化して再利用しやすくしておくと楽になるのかもしれない。

PyXspecを使えば、別次元の扉が開ける気もするが、import するとFatal Python error・・・何故だ・・・。

他の人がどうやっているのか & sherpaを使った場合を知りたいところ。

#coding:utf-8

SCRIPT_TYPE = "InstrumentalityOfMankind"
SCRIPT_SHELL = "python"
SCRIPT_NAME = "tempFit.py"
SCRIPT_HELP_ARG = ""
SCRIPT_HELP_SENTENCE = "This is a sample script for"
SCRIPT_NUM_ARG = 0
SCRIPT_VERSION = 1.0

###IMPORT###

import sys
import os
import subprocess

###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###
range_min = 2.3
range_max = 10.0
Basic_Model = 'Model_Basic.xcm'
Command_error = 'error max 1000 1-29 \n'

#デバック用
#Command_error = '\n'


for Reg in 'abcd':
    ###Input file###
    filename = 'sp.Position_' + Reg + '.FI.rebin.pi'
    Command_file = \
    'data 1 ' + filename + '\n ' + \
    '@' + Basic_Model + '\n'


    ###Temperature###
    for Temp in ['1Temp', '2Temp']:
        if Temp == '1Temp':
            Command_Temperature = """
            newpar 6 1 -1 \n 
            newpar 8 0 -1 \n
            newpar 9 0 -1 \n
            newpar 10 0 -1 \n
            """
            Abunds = ['1Solar', 'Free']
            Abss = ['Same']
        elif Temp == '2Temp':
            Command_Temperature = """
            thaw 6  \n 
            newpar 1 \n
            thaw 8  \n
            thaw 9  \n
            thaw 10  \n
            """
            Abunds = ['1Solar', 'Free', 'Free_Respectively']
            Abss = ['Same', 'Free_Respectively']

        ###Abund###
        Command_Abund1 = 'newpar 7=2 \n'
        for Abund in Abunds:
            if Abund == '1Solar':
                Command_Abund = Command_Abund1 + ' newpar 2 1 -1 \n'
            elif Abund == 'Free':
                Command_Abund = Command_Abund1 + ' thaw 2 \n'
            elif Abund == 'Free_Respectively':
                Command_Abund = Command_Abund1 + """
                untie 2,7 \n
                thaw 2 \n
                thaw 7 \n
                """
            ###Absorption###
            for Abs in Abss:
                if Abs == 'Same':
                    Command_Abs = 'newpar 10 = 5 \n'
                elif Abs == 'Free_Respectively':
                    Command_Abs = """
                    thaw 5 \n
                    thaw 10 \n
                    """

                ###Output files###
                output_id = 'Region_' + Reg + '_' + Temp + '_Abund_' + Abund + '_Abs_' + Abs
                logfile = output_id + '.log'
                xcmfile = output_id + '.xcm'
                if os.path.exists( logfile ):
                    os.remove( logfile )
                if os.path.exists( xcmfile ):
                    os.remove( xcmfile )

                ###Fit###
                Command_Fit = \
                'ignore **-' + str( range_min ) + '\n ' + \
                'ignore ' + str( range_max ) + '-** \n' + \
                """
                
                query yes \n 
                renorm \n 
                fit \n
                """ + \
                Command_error + \
                'log ' + logfile + '\n' + \
                'show \n' + \
                Command_error + \
                'log none \n' + \
                'save all ' + xcmfile + ' \n'

                xsp = subprocess.Popen( 'xspec', stdin = subprocess.PIPE )
                xsp.stdin.write( Command_file )
                xsp.stdin.write( Command_Temperature )
                xsp.stdin.write( Command_Abund )
                xsp.stdin.write( Command_Abs )
                xsp.stdin.write( Command_Fit )
                xsp.stdin.write( 'quit \n \n' )
                xsp.wait()

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年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月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に到達する前に挫折。やれやれ。

登録: 投稿 (Atom)