XIS sensitivity xis-limit.py

#!/usr/bin/env python
import math

myname = 'xis-test1.py'

area_file = 'fi2-area.qdp' # for 2 FI
nxb = 'fi1-nxb-fov.qdp' # for XIS1, cts/s/keV/fov/1FI

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

# XIS
omega_fov_arcmin2 = 17.0 * 17.0 # arcmin2
omega_arcmin2 = 9 # arcmin2

omega = omega_arcmin2 * arcmin2 # sr
omega_fov = omega_fov_arcmin2 * arcmin2 # sr

# functions
def eff_area(energy):
#    area = 2.0 * energy
    for line in open(area_file):
            lines = line.split(' ')
            try:
                in_energy = float(lines[0]);
                in_area = float(lines[1]);
            except ValueError: continue
            if in_energy <=0.0 or in_area <= 0.0: continue

            if in_energy >= energy :
                area = in_area
                break
    return area
# for 2 FI area (cm2)  

print '! created by %s¥n' % myname
print '! for 2FI ¥n'
print '! energy (keV) NXB-limit (LU) photon-limit (LU) ! Area (cm2)¥n'

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


    if len(lines) != 4 : continue

    try:
        energy = float(lines[0]);
        nxb_cts = float(lines[2]);
    except ValueError: continue

    if energy <=0.0 or nxb_cts <= 0.0:
        continue

    area = eff_area(energy) # for 2 FI
    e_reso = 130.0/1000.0 * (energy/6.0) ** 0.5

    nxb_per_so = nxb_cts * 2.0 / area / omega_fov # cts/s/cm2/sr/keV for 2 FI

    cxb_norm = 10.0 # cts/cm2/sr/keV
    cxb_per_so = cxb_norm * (energy) ** -1.4

    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.3f " % energy,
    print "%4.1e " % limit_bgd_nxb,
    print "%4.1e " % limit_bgd_cxb,
    print "%4.1e " % limit_photon,
    print "%4.1e " % limit_sum,
    print "!",
    print "%4.1e " % area,
    print "%4.1e " % nxb_cts,
    print "%4.1e " % e_reso,
#    print "%4.1e " % nxb_per_so,
    print

print '!end '