calculate SXS sensitivity (sxs2.py)

#!/usr/bin/env python
import math

myname = 'sxs2.py'
sxs_area = 'sxs-area-gbgd.txt'
# paste sxs-area.qdp apec02kev1solarNorm150.txt | awk '{print $1, $2, $5}' >| sxs-area-gbgd.txt
# 0625-1m.xcm

# general para.
sigma = 3.0
arcmin2 = 8.46e-8 # sr
time = 3e5 # sec

# SXS
omega_arcmin2 = 9 # arcmin2
nxb = 0.7e-3 # /s/sxs/keV from QR at 6 keV
e_reso = 5.0/1000 # keV

omega = omega_arcmin2 * arcmin2 # sr
#print omega,

print '! created by %s¥n' % myname
print '! energy (keV) NXB-limit (LU) CXB-limit (LU) photon-limit (LU) total-limit (LU0 ! Area (cm2)¥n'

for line in open(sxs_area):            # use file iterators
#    line.rstrip()
#    print line
    lines = line.split(' ')

    if len(lines) != 3 : continue
    energy = float(lines[0]);
    area = float(lines[1]);
    cxb_galactic = float(lines[2]);
    # /s/cm2/keV/sr
    if energy <=0.0 or area <= 0.0:
        continue

    nxb_per_so = nxb / area / omega # cts/s/cm2/sr/keV
    cxb_per_so = cxb_galactic
    f = area * omega * time

    limit_bgd_nxb = sigma * math.sqrt (nxb_per_so * e_reso/f)
    limit_bgd_cxb = sigma * math.sqrt (cxb_per_so * e_reso/f)
    limit_photon = sigma ** 2 * (f) ** -1

    limit_sum = math.sqrt(limit_bgd_nxb **2 + limit_photon **2 + limit_bgd_cxb **2)

    print "%6.4f " % energy,
    print "%4.2e " % limit_bgd_nxb,
    print "%4.2e " % limit_bgd_cxb,
    print "%4.2e " % limit_photon,
    print "%4.2e " % limit_sum,
    print "!",
    print "%4.1e " % area,
    print

print '!end of sxs1.py'