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年12月26日水曜日

LabViewで拡大

は、できないそうだがMacだとcontrol+マウスホイールで幸せになれる。
投稿者 時刻:
ラベル:

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日土曜日

PythonでXspec

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()
投稿者 時刻:
ラベル: ,

2012年4月27日金曜日

WebExが"Connecting..."で停止する

TV会議で使っているCiscoのWebExが"Connecting..."で停止してしまって困っていたが ~/Library/Application Support/WebEx Folderを丸ごと削除すると問題解決。
投稿者 時刻:
ラベル:

2012年4月25日水曜日

set table

新しいgnuplotで"set terminal table"が使えなくて困っていたが"set table"で同じことができるようだ。

http://www.gnuplot.info/docs_4.2/gnuplot.html#x1-27000043.61
投稿者 時刻:
ラベル:

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 (0,1,3) とHXDのNormが食い違う理由

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論文としては下記を引くべき。


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の記述は少ないけど。
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

投稿者 時刻:
ラベル:

2010年12月26日日曜日

gnuplotで太陽を表す記号(\odot)の表示方法

http://objectmix.com/graphics/140030-astronomical-symbols-gnuplot.html
set term post eps enhanced color "Helvetica" 20 fontfile "cmsy10.pfb"
set label 1 "1 M_{/CMSY10 \014}" とか。
cmsy10.pfbをcurrent dirにおく。finkなら/sw/share/texmf-dist/fonts/type1/bluesky/cm/の下にある。
投稿者 時刻:
ラベル:

2010年12月18日土曜日

caioのインストール時、configure, ciao.bashでコケる。

何故かPATHを追加するだけでなく書き換えてしまってうまく行かない。
仕方ないので一旦finkと/usr/localへのPATH以外はすべて消去してconfigure後、.bash_profileの一番先頭でciao.bashを走らせてやる（その後PATHを設定）とうまく行く。
Sherpaでfittingが出来るようになった。幸せ。
投稿者 時刻:
ラベル: ,

2010年12月7日火曜日

Scientific linuxでキーアサインを日本語に変更

Preferences -> Keyboard -> Layouts -> add でJapanを加えたあと、U.S. Englishを削除。


こんな感じ。
投稿者 時刻:
ラベル:

2010年11月20日土曜日

Root チュートリアル

Rootが使えずに（そしてRootのマニュアルを読んでみたがあまりに詳細すぎて）往生していたが、良いチュートリアルを発見。
  • Root tutorial
http://www.nevis.columbia.edu/~seligman/root-class/

By B. Seligman (コロンビア大学 Nevis Lab)
マクロの作り方とか実際に使う上で欲しい情報がまとまっていて良いです。
投稿者 時刻:
ラベル:

2010年10月15日金曜日

眠気覚ましの手ツボ

中衝
http://www.genki39.com/archives/64675132.html
がとても良く効きました。
投稿者 時刻:
ラベル:

行動ログノート

Kokuyoのノー242Cが便利。
投稿者 時刻:
ラベル:

2010年10月14日木曜日

PASJでの集録などの参考文献

PASJの原稿作成の手引きには何故か会議集録を参考文献にする時の方法が明示されていない。
仕方ないので天文月報の手引きと書籍を引用する際の例を参考に書いてみる。

例えば、これなら

\bibitem[Yusef-Zadeh et al.(1996)]{Yu96}
Yusef-Zadeh, F., Robinson, B.~T., Roberts, D.~A., Goss, W.~M., Frail, D.~A.,
\& Green, A.\ 1996, in The Galactic Center: 4th Eso/Ctio Workshop, ed.
R. Gredel (San Francisco Astronomical: Society of the Pacific), 151

とか？
出版社（とその所在）が分からない時はISBNから調べて見るとか。
投稿者 時刻:
ラベル:

ThunderbirdでAPOP

京都のメールサーバーはAPOPを使わないとメールを取りにいけないが、thunderbird v3.1ではAPOPを選択できない。あきらめてv3.0にダウングレードすると成功。うーん・・・。
投稿者 時刻:
ラベル: ,

Macのショートカット

  • 「Command」+「Option」を押しながらドラッグ＆ドロップするとエイリアスを作成
  • スクリーンショット
    • 指定範囲 コマンド＋ Shift ＋ 4 キーを押してから、ピクチャに取り込みたい領域をドラッグして選択
    • 全画面  を撮影するときは、コマンド＋ Shift ＋ 3 キー
投稿者 時刻:
ラベル:

2010年9月29日水曜日

Gmailで未読メールのみ検索

is:unread label:Inbox
投稿者 時刻:
ラベル:

2010年9月24日金曜日

すざくのデータプロセッシングの流れ

USCにて「RPTって何の略だっけ？」と思った際に見つけた尾崎さんの解説。
勉強になります。Raw Packet Telemetry。 
投稿者 時刻:
ラベル:

2010年8月12日木曜日

Gnuplotで作ったPSファイルの色を弄る時

http://t16web.lanl.gov/Kawano/gnuplot/postscript.html#8.4


を参考にRGB辺りを弄れば良いワケだが

http://www.scollabo.com/banban/lectur/websafe.html

を255で割って入れてやれば望みの出力得られる。

51/255=0.2
投稿者 時刻:
ラベル:

2010年4月29日木曜日

古いAide-memoire

    QDPで上付き・下付き文字

    • \uで１レベル上へ、\dで１レベル下へ。以降、変更しない限り、そのレベルを維持する。
    • なので、Counts s^{-1} keV^{-1} とやりたい時は Counts s\u-1\d keV\u-1\dとする。

    QDPでギリシア文字

    • \gを付ければOK。\gaはα。\ggはγ。\gGでΓ。

    SCI-on時STATUS ColumnによるCut

    • CI行の2nd tailing row (good pixel) がB18でキャルマスクのB16より上位なので割と込み入ってます。以下オマジナイ。
    • calsourecを含んだgood pixel : status=0:524287
    • calsourceを抜いたgood pixel : status=0:65535 131072:196607 262144:327679 393216:458751
    • calsource内ののgood pixel : status=65536:131071 196608:262143 327680:393215 458752:524287

    Fedora Core Linux (少なくとも4) の管理コマンド

    • system-config-(色々) networkとかprinterとか。

    Rsyncでバックアップ

    • rsync -avz （バックアップ元マシン）:（バックアップ元DIR) (バックアップ先マシン）：（バックアップ先DIR）

    Gnuplotで特殊記号

    • とりあえずpostsrcript enhancedで文字コード直接入力。困ったらps_guide.psを見る。以下使いそうなもの。
    • × : {/Symbol \264}, ° : {/Symbol \260}, ・: {/Symbol \267}

    Gnuplotで〜を表示

    • 上記の方法で{/Symbol \176}とすると何故かうまくいかない。"\\〜(半 角)"と入力するとうまくいく。

    Tarで圧縮

    • tar zcvf backup.tar.gz backup_dir

    日本語環境(LANG=JP)でgeditを使用すると異常に文字が小さくなる

    • 慌てず騒がずgnome-font-properties

    日本語環境(LANG=JP)でds9を使用すると異常に不安定になる

    • alias ds9 'env LANG=POSIX ds9' とか.cshrcに書いておく。

    「すざく」XISデータのCALDBの掛け直し方

    • いろんな人に割とよく質問される。
    • CALDBを掛け直したunfilteredイベントの作り方はこ れを読んで下さい。(宇宙研の森さん謹製。「すざく」公式ページにも同じものがある筈だが見つからない・・・）
    • unfilteredイベントからcleanedイベントを作る方法はこ こ。
    • xispiを掛け直す時、便利かもしれないス クリプトとか。ただし使用は自己責任で。

    SCI時makepiファイルのログ (20080825まで)

    • こ れ。しかし、このファイルも公開用に準備したのにどこかに消えてしまった・・・。

    GNU Scientific Library (GSL)

    • 今更だが、とてもとても便利。
    • こちら(日本語での解説、インス トールガイド）を見ると使える様になります(外部ページ）。

    中性子イメージング教科書

    • 最近、うちの研究室でX・ガンマともに流行って(？)いる中性子イメージングについての日本アイソトープ協会による解説。
    • 中性子イメージング技術 の基礎と応用
    • コンバータを使うのはその方が圧倒的に断面積がお得だから。

英語・論文書式

もういい加減、英語で怒られるのはなんとかしたいです。
投稿者 時刻:
ラベル:
登録: 投稿 (Atom)