# Creates fake photon evemnt lists (photenlist.dat) from the spectrum in the fak file
# The output is a FITS file with the event list prepended to the original FITS file
from astropy.io import fits
import numpy as np
from astropy.time import Time, TimeDelta
import datetime

def generate_event_file(input_spectrum,output_file,exptime=None):
    with fits.open(input_spectrum) as hdul:

        data = hdul[1].data
        chan=data['CHANNEL']
        cnt=data["RATE"]
        nchan=len(chan)

        cdf=np.cumsum(cnt)
        rate=cdf[-1]    # total photons per second
        if exptime==None:
            exptime=hdul[1].header['EXPOSURE']  # exposure time in seconds

        nphot=int(rate*exptime)
        s = np.random.uniform(0,rate,nphot)
        channellist=np.floor(np.interp(s,cdf,chan))
        t=  np.sort(np.random.uniform(0,exptime,nphot))

        
        
        newtable=fits.BinTableHDU.from_columns(fits.ColDefs([
            fits.Column(name='TIME',format='1D',unit='s',array=t),
            fits.Column(name='PI',format='1I',unit='chan',array=channellist),
            fits.Column(name='RAWX',format='1B',unit='pixel',array=0*t),
            fits.Column(name='RAWY',format='1B',unit='pixel',array=0*t),
            #fits.Column(name='PHA',format='1I',unit='chan',array=channellist),
            #fits.Column(name='PHA_FAST',format='1I',unit='chan',array=channellist),
            #fits.Column(name='DET_ID',format='1B',unit='chan',array=0*channellist+4),
            ]),name='EVENTS')
        
        hdul.insert(1,newtable)
        
        newtable=fits.BinTableHDU.from_columns(fits.ColDefs([
            fits.Column(name='START',format='1D',unit='s',array=[0]),
            fits.Column(name='STOP',format='1D',unit='s',array=[exptime]),
            #fits.Column(name='RAWX',format='1B',unit='pixel',array=0*t),
            #fits.Column(name='RAWY',format='1B',unit='pixel',array=0*t),
            #fits.Column(name='PHA',format='1I',unit='chan',array=channellist),
            #fits.Column(name='PHA_FAST',format='1I',unit='chan',array=channellist),
            #fits.Column(name='DET_ID',format='1B',unit='chan',array=0*channellist+4),
            ]),name='GTI')

        hdul.insert(2,newtable)

        t=Time(datetime.datetime.now())

        for i in range(3):
            hdul[i].header['INSTRUME']= ('TES','Instrument name')
            hdul[i].header['TELESCOP']= ('COLIBRI','Telescope (mission) name')
            hdul[i].header['OBJECT'] = ('SIMULATION')
            hdul[i].header['DATE-OBS']= (t.fits,'Start date of observations')
            hdul[i].header['DATE-END']= ((t+TimeDelta(exptime,format='sec')).fits,'End date of observations')
            hdul[i].header["TSTART"]=(0E0,'[s] time start')
            hdul[i].header["TSTOP"]=(exptime,'[s] time stop')
            hdul[i].header['MJDREFI'] = (int(t.mjd),'[d] MJD reference day %s' % t.fits)
            hdul[i].header['MJDREFF'] = (t.mjd-int(t.mjd),'[d] fractional part of reference day')
            hdul[i].header['RA_NOM'] =  (0E0,'[deg] R.A. of nominal aspect point [J2000]')
            hdul[i].header['DEC_NOM ']= (0E0,'[deg] Dec. of nominal aspect point [J2000]')
            hdul[i].header['RA_OBJ']=(0E0,'[deg] R.A. of target [J2000]')
            hdul[i].header['DEC_OBJ']= (0E0,'[deg] Dec. of target [J2000]')

        hdul[1].header['TCRPX3']= (0,'reference pixel')
        hdul[1].header['TCRVL3']= (0.0,'reference value')
        hdul[1].header['TCDLT3']= (1.0,'pixel step value')
        hdul[1].header['TCRPX4']= (0,'reference pixel')
        hdul[1].header['TCRVL4']= (0.0,'reference value')
        hdul[1].header['TCDLT4']= (1.0,'pixel step value')
        hdul[1].header['TLMIN']=(0,'minimum legal value in the column')
        hdul[1].header['TLMAX']=(73684,'maximum legal value in the column')
        hdul[1].header['TIMEPIXR'] = (0.0000000000000E+00, 'Bin time beginning=0 middle=0.5 end=1')
        hdul[1].header['TIMEDEL'] = (40.0E-9)
    
        for i in range(1,3):
            hdul[i].header['DATE']=(t.fits,'file creation date')
            hdul[i].header['EXPOSURE']=(exptime,'[s] Exposure time')
            hdul[i].header['ONTIME']=(exptime,'[s] Sum of stop-start times')
            hdul[i].header['TIMESYS'] = ('TT')
            hdul[i].header['DATAMODE'] = ('PHOTON')
            hdul[i].header["TIMEUNIT"]= ('s','unit for time keywords')
            hdul[i].header['TIERRELA']=(1.0E-8,'[s/s] relative errors expressed as rate')
            hdul[i].header['TIERABSO']=( 1.0, '[s] timing  precision  in seconds')
            hdul[i].header['EQUINOX']=(2000.0,'[yr] Equinox of celestial coord system')
            hdul[i].header['RADECSYS']= ('FK5','celestial coord system')            
            hdul[i].header['TIMEZERO']=(0E0,'Zero-point offset for TIME column')

        hdul.writeto(output_file,output_verify='warn')

    return

def main():
    import sys
    if len(sys.argv)<3:
        print("""Format:
python fakphot.py _Input_Fake_Spectum_ _Output_FITS_Event_File_ [
sure_Time]
""") 
    else:
        if len(sys.argv)>3:
            generate_event_file(sys.argv[1],sys.argv[2],exptime=float(sys.argv[3]))
        else:
            generate_event_file(sys.argv[1],sys.argv[2])
            
    
# nbins=1000.0
# spec,b=np.histogram(plist, bins=nbins)
# plt.plot(b[1:],spec,'k-')
# plt.plot(chan,cnt*mf*nchan/nbins,'b-')
# plt.show()

if __name__ == "__main__":
    # execute only if run as a script
    main()