2011/08/28

海外学振結果

不採用Aでした。

項目別評価は、
①研究業績 3.67
②研究計画 3.50
③外国の機関で研究することの意義 3.67
④総合評価 3.34
総合評価のTスコアは3.156。

特別な弱点も長所もなく、全体的に点が足りないということですね。
個人的には、乏しい研究業績(*)よりも研究計画の評価が低かったのは意外だったので、見直して今後に活かさなければなりません。といっても今年度の主だった書類はすでに申請ずみなのですが。

(*)=主著3報(PASJ×2, NIM A×1)+共著2報(ApJ 2nd×1, PASJ 3rd×1)+国際会議口頭×1+その他

2011/08/22

simxのinstallのtrouble shooting

1. fortranコンパイラをみつけられない問題

=configureで以下のwarningが出るとき。
configure: WARNING: cfitsio: == No acceptable f77 found in $PATH
configure: WARNING: cfitsio: == Cfitsio will be built without Fortran wrapper su


simx-1.1.2/configureの中で実行されるsimx-1.1.2/libsrc/cfitsio/configureはg77を探す。
おそらく多くのMac userはfortranコンパイラを持っていないか、HEASoftのコンパイルのために持っていたとしてもg77ではなくgfortran。てっとり早い対処はgfortranをg77として使わせること。
参考: http://monteboo.blogspot.com/2007/12/cfitsio-fortran-wrap.html
1-0. gfortranを入手。
1-0-1. finkを入れる。

    • Must be installed from source code at this time (2010/04)
    • After downloading a tarball,
      •    $ cd ~/Download
           $ tar -xvf fink-0.29.10.tar
           $ cd fink-0.29.10/
           $ ./bootstrap   
      • Choose [1: Default (mostly 32bit)], [1: Use sudo], [/sw] in the bootstrap.
      • Choose defaults for the others
    • After installation, add source /sw/bin/init.csh to .cshrc
    • and do $ sudo fink selfupdate once

          (なぜか外から閲覧できなくなってる?wikiから転載: http://byakko.scphys.kyoto-u.ac.jp:8830/MSworklog/iMacSetUpLog#General)

1-0-2. sudo fink install gcc4.4でgfortran4.4.1を入れる。
1-1. $PATHに/sw/lib/gcc4.4/lib追加
1-2. sudo ln -s /sw/bin/gfortran /sw/bin/g77
1-3. sudo ln -s /Users/sawada/Downloads/simx-1.1.2/libsrc/cfitsio/* /sw/lib/gcc4.4/lib
(1-3.が必須かどうか未確認。)

2. src/libsimx.a のtableができない問題

=makeで以下のerrorが出るとき。
ld: in ../src/libsimx.a, archive has no table of contents
collect2: ld returned 1 exit status
make[1]: *** [simx] Error 1
make: *** [all] Error 2


2-a. arのoptionを変更する
/Users/sawada/Downloads/simx-1.1.2/src/MakefileのL33を
$(AR) rv $@ $(OBJS)
から
$(AR) srv $@ $(OBJS)
に修正
参考: http://hidemon-memo.blogspot.com/2009/01/mac-os-ar-commandmode.html
2-b. (未試行) Randollの示した対処法
"cd src ; ranlib libsimx.a ; cd .."
参考: http://hea-www.harvard.edu/simx/install.html の末尾

2011/08/09

SPEX cie codeで計算に含まれる元素

SPEXのプラズマモデルはZ=1-30までの元素について組成比をパラメータとして持っていますが、実際のところCr, Mnなどは組成をどれだけ増やしてもHe-like Kα輝線とかが出てくることはありません。なぜならパラメータがあるだけで計算してないから!組成や電離度をちゃんと計算している元素がどれかは、ascdumpコマンドを叩けばわかります。

SPEX> com cie
 You have defined    1 component.
SPEX> calc
SPEX> ascd term 1 1 abund
 number of layer lines :        4076           0
 Atomic   Element  Abundance      Abundance     Average charge
 number   name     (solar units)  (absolute)
 -------------------------------------------------------------
  1        H       1.000          1.000             1.000
  2        He      1.000         9.7051E-02         2.000
  6        C       1.000         2.7733E-04         5.999
  7        N       1.000         8.1658E-05         6.996
  8        O       1.000         6.0534E-04         7.988
 10        Ne      1.000         1.2677E-04         9.923
 11        Na      1.000         2.2233E-06         10.82
 12        Mg      1.000         3.9719E-05         11.60
 13        Al      1.000         3.2584E-06         12.25
 14        Si      1.000         3.8548E-05         12.81
 16        S       1.000         1.6218E-05         14.21
 18        Ar      1.000         3.5727E-06         15.99
 20        Ca      1.000         2.3281E-06         17.87
 26        Fe      1.000         3.2659E-05         20.60
 28        Ni      1.000         1.8880E-06         22.38
てことで、N, Na, Alをのぞくodd Zの元素 (含 Mn)、Be, Ti, Crは計算してくれませんのでお気をつけを。

ちなみにlast columnのAverage chargeは楽しいですね。rtとか変えればもちろん値が変わっていきます。

SPEXのelectron densityについて

SPEXのいいところは (manualでも自画自賛の文章があるけど) データを再現する温度や元素組成比といった「表の」パラメータだけじゃなく、イオン存在比や各放射プロセスごとの寄与といった「裏の」プラズマの物理状態を知ることができる点だ。また単純な電離平衡プラズマでもAPECとは比べ物にならないほど高い自由度を持ってる。その分よくわかってないことに出くわす頻度も高く、遊んでるとなんだこりゃって思う部分も多々ある。

今日でくわしたのは密度に関する話。

cieを初めとするプラズマモデルでは電子密度(n_e)を陽に指定することができる。放射強度は電子密度とイオン密度の積(n_e * n_X)に比例するが、これとemissivityに効く電子温度(T_e)を一定に保っている限り、電子密度自体はCCDクオリティの分光ではスペクトルの形状をほんの少ししか変えない。(マイクロカロリメータ-クラスの分解能があればtripletのR ratioが変わるのが見える。)

でもよくよく考えると、n_e * n_X, T_e不変でn_eをいじるってことはn_Xいじってるんじゃね?ってことになって、そうだとすると電離の進行度とかが影響されちゃう。SPEXのascdumpコマンドではn_e / n_Hを教えてくれるので、n_eを変えてこの比がどう変わるか見てみた。すると、

SPEX> com cie
 You have defined    1 component.
SPEX> pa 1 1 ed s t
SPEX> pa 1 1 ed v 1e-14
SPEX> calc
SPEX> pa sh fr
--------------------------------------------------------------------------------------------------
sect comp mod  acro parameter with unit     value      status    minimum   maximum lsec lcom lpar

   1    1 cie  norm ne nX V (1E64/m**3)   1.000000     thawn     0.0      1.00E+20
   1    1 cie  t    Temperature (keV)     1.000000     thawn    5.00E-04  1.00E+03
   1    1 cie  ed   El  dens (1E20/m**3) 9.9999998E-15 thawn    1.00E-22  1.00E+10



--------------------------------------------------------------------------------
 Fluxes and restframe luminosities between   2.0000     and    10.000     keV

 sect comp mod   photon flux   energy flux nr of photons    luminosity
              (phot/m**2/s)      (W/m**2)   (photons/s)           (W)
    1    1 cie   3.732456E-03  1.568204E-18  4.690343E+42  1.970663E+27

SPEX> ascd term 1 1 plas
 number of layer lines :        4076           0
 Plasma parameters
 -----------------
 Electron temperature     :   1.0000     keV
 Electron density/1E20    :  1.00000E-14 /m**3
 Hydrogen density/1E20    :  8.30225E-15 /m**3
 Electron/Hydrogen density:   1.2045
SPEX> pa 1 1 ed v 1e-13
SPEX> calc
SPEX> pa sh fr
--------------------------------------------------------------------------------------------------
sect comp mod  acro parameter with unit     value      status    minimum   maximum lsec lcom lpar

   1    1 cie  norm ne nX V (1E64/m**3)   1.000000     thawn     0.0      1.00E+20
   1    1 cie  t    Temperature (keV)     1.000000     thawn    5.00E-04  1.00E+03
   1    1 cie  ed   El  dens (1E20/m**3) 9.9999998E-14 thawn    1.00E-22  1.00E+10



--------------------------------------------------------------------------------
 Fluxes and restframe luminosities between   2.0000     and    10.000     keV

 sect comp mod   photon flux   energy flux nr of photons    luminosity
              (phot/m**2/s)      (W/m**2)   (photons/s)           (W)
    1    1 cie   3.732456E-03  1.568204E-18  4.690343E+42  1.970663E+27

SPEX> ascd term 1 1 plas
 number of layer lines :        4076           0
 Plasma parameters
 -----------------
 Electron temperature     :   1.0000     keV
 Electron density/1E20    :  1.00000E-13 /m**3
 Hydrogen density/1E20    :  8.30225E-14 /m**3
 Electron/Hydrogen density:   1.2045

はい、変わってませんね。n_e/n_H不変。よくよく考えると放射強度はn_e * n_X * Vに比例なので、n_e増やしたらn_Xも同じだけ増やして、Vはその逆2乗で減らす、ということになてるんですね。イオンと電子の密度の比は不変に保たれるので、電離の進行度が変にn_eに左右されることはないと安心できます。

ところで、refパラメータをいじってH以外の元素(ここではHe)に対する組成比の計算に切り替えた上でHの組成比を0にしてやると、、
SPEX> pa 1 1 ed v 1e-14
SPEX> pa 1 1 ref v 02
SPEX> pa 1 1 01 v 0
SPEX> calc
SPEX> pa sh fr
--------------------------------------------------------------------------------------------------
sect comp mod  acro parameter with unit     value      status    minimum   maximum lsec lcom lpar

   1    1 cie  norm ne nX V (1E64/m**3)   1.000000     thawn     0.0      1.00E+20
   1    1 cie  t    Temperature (keV)     1.000000     thawn    5.00E-04  1.00E+03
   1    1 cie  ed   El  dens (1E20/m**3) 9.9999998E-15 thawn    1.00E-22  1.00E+10



--------------------------------------------------------------------------------
 Fluxes and restframe luminosities between   2.0000     and    10.000     keV

 sect comp mod   photon flux   energy flux nr of photons    luminosity
              (phot/m**2/s)      (W/m**2)   (photons/s)           (W)
    1    1 cie   2.904043E-03  1.210261E-18  3.649328E+42  1.520860E+27

SPEX> ascd term 1 1 plas
 number of layer lines :        4071           0
 Plasma parameters
 -----------------
 Electron temperature     :   1.0000     keV
 Electron density/1E20    :  1.00000E-14 /m**3
 Hydrogen density/1E20    :  4.89013E-14 /m**3
 Electron/Hydrogen density:  0.20449
ってなるんですが。。つまりascdumpで出てくるn_e/n_Hってほんとはn_e/(Σ{Z}{n_elem(Z)})だよね多分って突っ込み。実際計算してるのがどういう物理量なのかは現在問い合わせ中。。

2011/08/05

Python for shellscripters

解析フリークのみなさま、

High energy astro系の解析はNASAのHEASoftが主流だと思います。HEASoftに含まれるXSpecや各種ftoolsを使う上で、多くの人がshell scriptを書いていると思うのですが、個人的にshell scriptはクソだと思うので、即刻これを捨てて便利なpythonに移行しましょー!、というのが本ポストの趣旨です。

Pythonはちょいちょい使うものの、いざshell scriptを完全に置き換えたい、というとき具体的にどうすればいいんだー、と悩んだ個人的経験から、解析で頻出しそうなshell script的用法をまとめたいと思います。


 1. Shell commandや外部スクリプトの実行


os.systemを使います。
import os
os.system('cp aho1.txt aho2.txt')

()内はただのstringなのでpythonで使ってる変数も使えます。
import os
infile = 'aho1.txt'
outfile = 'aho2.txt'
os.system('cat '+infile+' > '+outfile)

これより複雑なことをしたければsubprocess.Popenを使います。

2. コマンドを実行して標準出力を返り値として得る


.communicate()を使います。
import subprocess
PIPE=subprocess.PIPE
eventfile = subprocess.Popen(['ls', 'temp.evt'], stdin=None, stdout=PIPE, shell=False).communicate()[0]
Popenの引数1つ目は実行するコマンドとその引数からなるリスト。shell=Trueにすれば、
import subprocess
PIPE=subprocess.PIPE
eventfile = subprocess.Popen('ls temp.evt', stdin=None, stdout=PIPE, shell=True).communicate()[0]
とベタ書きしてもok。また、shell=Trueではwild cardが使えるので、
eventfiles = subprocess.Popen('ls *_xis?.evt', stdin=None, stdout=PIPE, shell=True).communicate()[0]
なおこのままだと使い勝手がわるいし、最後に改行('\n')まで付いてきちゃうので、
eventfiles = subprocess.Popen('ls *_xis?.evt', stdin=None, stdout=PIPE, shell=True).communicate()[0].split()
としてリストにしておくと便利。こう来ると大抵
eventfiles = subprocess.Popen('ls *_xis?.evt', stdin=None, stdout=PIPE, shell=True).communicate()[0].split()
for eventfile in eventfiles:
とつづく。。

3. コマンド、プログラムを実行し、その標準入力とやり取りする

stdin.writeを使います。いわゆるhere documentです。XSpecやxselectを使うとき便利。
入力する各行末尾の改行を忘れずに。
import subprocess
PIPE=subprocess.PIPE
xsp = subprocess.Popen('xpec',stdin=PIPE,stdout=None,shell=False)
xsp.stdin.write(
    #ここにコマンド
    )
xsp.wait()
途中でfor回したりもできて便利。
import subprocess
PIPE=subprocess.PIPE
xse = subprocess.Popen('xselect',stdin=PIPE,stdout=None,shell=False)
xse.stdin.write(
    '\n'
    'no\n'
    'read event eventlist.txt\n'
    '.\n'
    'set image det\n'
    'set xybin 1\n'
    )

for seg in segs:

    subdir = 'xis'+xis+'_seg'+seg

    for grade in grades:
        qdpfile = str(tobs)+'_'+subdir+'_g'+grade+'.qdp'
        pifile = str(tobs)+'_'+subdir+'_g'+grade+'.pi'

        xse.stdin.write(
            'filter column "grade='+graderanges[grade]+' segment='+seg+':'+seg+' status=65536:131071 196608:262143 327680:393215 458752:524287"\n'
            'extract spec\n'
            'plot spec\n'
            'we '+subdir+'/'+qdpfile+'\n'
            'quit\n'
            'save spec '+subdir+'/'+pifile+'\n'
            'no\n'
            'clear column\n'
            )
xse.stdin.write(
    'exit\n'
    'no\n'
    '\n'
    )
xse.wait()

あとは得られたデータをnumpy, scipyで料理してmatplotlibでお絵描き、って流れです。
そこはまたおいおい。とりあえずここまでー。

さわだ

2011/08/03

Heasoft 6.11でのftools仕様変更

解析フリークのみなさま、

Heasoftを6.11にversion upしたところいくつか瑣末な仕様変更で面倒な目に遭遇したので情報共有します。

Heasoft 6.11 Release Note

1. Xspec->iplot->wenvironでのqdp出力

変更:
.pcoにいままでなかったViewから始まる行が追加される。

問題:
こいつはLocationと独立にwindowの図の表示範囲を指定するので、scriptでLocationのみでwinの調整している場合、予想と違う配置にされてしまう。

対処:
.pcoのViewの行を消すか、引数をかえて
View 0.1 0.1 0.9 0.9
とするとこれまでと同じ見た目になる。

2. mathphaのerrmeth, properrのデフォルト値

変更:
errmeth: POISS-1 -> POISS-0
properr: yes -> no

問題:
properrがnoだと誤差伝播しない上、演算した結果
負値を取った場合勝手に誤差を0に、bin quality
を5 (bad)にしてしまう。

対処:
mathphaを使用する際オプション
properr='yes'
をつける。
errmethに関しても新デフォルトではまずい場合は露に指定すべし。

仕様変更に気づかずに解析してしまっている恐れがあれば、fdumpなどで値を確認。
ちなみにfhelp mathphaで見られるヘルプファイルの記述にはこの変更が反映されてない。。

さわだ