Surface brightness of the Galactic diffuse X-ray emission

系内、特に銀河面付近のソースを扱うＸ線屋さんが避けて通れない背景拡散放射 (GDXE)。
その輝度の見積もりがむずかしくbackground estimationで苦労することもしばしば。

すざくの良質なスペクトルを元にH. Uchiyamaさんがその表面輝度をモデル化しています。
H. Uchiyama PhD (2010) PDFはこちら

Sgr A*を中心とする２つの2-D exponentialの和を用いたモデル (eq.4.2) は特に便利なので、

銀河座標を放り込むと表面輝度を表示するかんたんなスクリプトを書きました。

座標を２セットいれると２領域間での輝度の比を表示します。

無駄にnumpyつかってるのはご愛嬌。。

※初期版にひどいバグ:
l0, b0で絶対値とるの忘れてました！
12/5に参照した方は修正した下記をご覧下さい。

#!/usr/bin/env python

'''
The model expression and the parameters are respectively 
taken from eq.4.2 and table 4.11 of H. Uchiyama 2010, PhD Thesis
http://repository.tksc.jaxa.jp/pl/dr/IS8000028000/en
The parameters were measured with the K-alpha lines from FeXXV 
and FeXXVI and the 5-10 keV continuum bands.
The unit of surface brightness is 1e-7 ph/s/cm2/am2 for 
Fe Ka lines and 1e-13 erg/s/cm2/am2 for 5-10keV continuum.
'''

import numpy as np

gpos_SgrAs = np.array([-0.056, -0.046])

a1 = (26.3, 10.0, 3.4)
a2 = (1.6, 0.25, 0.20)

hl1 = (0.47, 0.51, 0.48)
hb1 = (0.22, 0.22, 0.19)

hl2 = (30, 47, 55)
hb2 = (1.7, 17, 2.5)

band2colnum = {
    'FeXXV':0,
    'FeXXVI':1,
    'hard':2
    }

def calc_sb(l, b, band):

    i = band2colnum[band]

    l0, b0 = np.abs(np.array([l, b]) - gpos_SgrAs)
    sb1 = a1[i] * np.exp(-l0/hl1[i]) * np.exp(-b0/hb1[i])
    sb2 = a2[i] * np.exp(-l0/hl2[i]) * np.exp(-b0/hb2[i])
    sb = sb1 + sb2

    return(sb)

def calc_ratio(l1,b1, l2,b2, band):

    sb1 = calc_sb(l1, b1, band)
    sb2 = calc_sb(l2, b2, band)
    ratio = sb1 / sb2

    return(ratio)

if __name__ == '__main__':

    import sys

    argvs = sys.argv

    try:
        l, b, band = argvs[1:4]
        cnum = band2colnum[band]
        sb = calc_sb(float(l), float(b), band)
        print(sb)
    except:
        try:
            l1, b1, l2, b2, band = argvs[1:6]
            cnum = band2colnum[band]
            rt = calc_ratio(float(l1), float(b1), float(l2), float(b2), band)
            print(rt)
        except:
            print('Usage1 calc sb : calcgdxe.py l b band')
            print('Usage2 calc sb ratio : calcgdxe.py l1 b1 l2 b2 band')
            print('band should be selected from FeXXV, FeXXVI, or hard')
            exit()

追記。

標準で入ってないないnumpyというライブラリを使っているので、
import numpyでこける場合は以下の３点を修正すれば動くと思います。
めんどくさいので動作みかくにん。

1)
import numpy as np
gpos_SgrAs = np.array([-0.056, -0.046])
↓
import math
l_SgrAs = -0.056
b_SgrAs = -0.046

2)
l0, b0 = np.abs(np.array([l, b]) - gpos_SgrAs)
↓
l0 = abs(l - l_SgrAs)
b0 = abs(b - b_SgrAs)

3)
np.exp()
↓
math.exp()