#stazeni souboru
import sys,os
import pyfits
import numpy as np
resp2=pyfits.open("spec.rmf")
arf2=pyfits.open("spec.arf")
#from matplotlib import pyplot as pl


#prevedeme ridkou matici na hustou
llen=800
ematff=np.zeros(llen)
gmat2ff=[]
for i in range(0,2000):
    q=resp2[1].data[i]
    if q['N_GRP']==0: continue
    if q['N_GRP']>1:
        smat=np.concatenate([np.arange(a,a+b) for a,b in np.r_[[q['F_CHAN'],q['N_CHAN']]].T])
    else:
        smat=np.arange(q['F_CHAN'][0],q['F_CHAN'][0]+q['N_CHAN'][0])
    ematff[:]=0
    ematff[smat]=np.log(q['MATRIX'])
    gmat2ff.append(ematff.copy())
gmat2ff=np.array(gmat2ff)

x1b=(resp2[1].data['ENERG_LO'][:2000]+resp2[1].data['ENERG_HI'][:2000])/2.
smatff=np.exp(gmat2ff)
smatff[gmat2ff==0]=0

#zprumerujeme matici v jednom smeru - stredni odezva na jeden kanal
x1a=resp2[1].data['ENERG_HI'][:2000:]
widb=(x1a-x1b)*2 #sirka binu  (x1b jsou stredy binu)
cnorm=smatff.sum(0) #normalizace na instr. kanal
cxhan=x1b.dot(smatff)/cnorm # stredni energie pro dany kanal
wnorm=(smatff*widb.reshape(2000,1)).sum(0) #stredni sirka

#pl.semilogy(wnorm)
#pl.ylim(0.01,0.1)
#pl.grid()

myspec="agn_grpspec.pha"
sp1=pyfits.open(myspec)[1].data['COUNTS']

# data pro fitovani:
#  * energie -- cxhan
#  * spektrum -- sp1
#
#for i in range(len(cxhan)):
#    print("{0} {1}".format(cxhan[i],sp1[i]))
#pl.semilogxy(cxhan,sp1)