#!/usr/bin/env python

version_num    = "1.1"		# Version number to be updated

#
# VERSION = '1.1';
# author  = 'Jose Vicente Perea';
# date    = '2017-04-04';
#
#
# ChangeLog
# =========
#
# Version 1.1 - 2017-04-04 (JVP)
# ------------
#
# + Include matplotlib.use('Agg') to avoid issue:
#   'This application failed to start because it could not find or load the Qt platform plugin "xcb"'
#
# Version 1.0 - 2017-03-28 (JVP)
# ------------
#
# + Initial version
#   Python script to get the Flaring Background filtering threshold cut
#   from the Source Time Series. That cut value will be used to produce
#   the Source Spectrum.
# + In addition a diagnostic plot is created. Default: No
#
# ------------------------------------------------------------------------------
#
#	usage: srcBkgGti.py [-h] [--version] --srctslist SourceTimeSeries.ds [--outfile OUTFILE] --outformat {png,pdf}
#
#	Get the Flaring Background filtering threshold cut from the Source Time Series.
#	
#	optional arguments:
#	  -h, --help     	show this help message and exit
#	  --version, -V  	show program's version number and exit
#
#	  --srctslist INFILE	Input Source Time Series FITS file. Default: SourceTimeSeries.ds
#	  --outfile OUTFILE	Output figure file. Default: figure
#	  --outformat {png,pdf}	Output figure format. Default: png
#
# ------------------------------------------------------------------------------
#
import sys,os
import warnings
import argparse



###################
# Inut Options
#
parser = argparse.ArgumentParser(description='Get the Flaring Background filtering threshold cut from the Source Time Series', prog='srcBkgGti.py')

#parser.add_argument('--version', '-v', action='version', version='%(prog)s (%(prog)s-1.0) [-]')
parser.add_argument('--version', '-v', action='version', version='%(prog)s (%(prog)s-{}) [-]'.format(version_num))

parser.add_argument('--srctslist', help='Input Source Time Series FITS file. Default: SourceTimeSeries.ds', default='SourceTimeSeries.ds')

parser.add_argument('--diagplot'  , help='Create a diagnostic plot. Default: No', default='No')
parser.add_argument('--outfile'  , help='Output figure file. Default: figure. If the output filename already has a valid extension, it will be kept', default='figure')
parser.add_argument('--outformat', help='Output figure format. Default: png', choices=['png', 'pdf'], default='png')

args = parser.parse_args()
#
###################

from astropy.io import fits
from astropy.table import Table
from astropy.coordinates import SkyCoord
from astropy import units as u

import numpy as np

import matplotlib
matplotlib.use('Agg')	# To avoid problems with DISPLAY variable
			# raise RuntimeError('Invalid DISPLAY variable')
			# RuntimeError: Invalid DISPLAY variable

import matplotlib.pyplot as plt


srcTSList = args.srctslist
ext       = args.outformat

diagPlot  = args.diagplot
diagPlotTest = 0
yes = ('yes', 'y')
no  = ('no', 'n')
if diagPlot.lower() in yes:
  diagPlotTest = 1
elif diagPlot.lower() in no:
  print('No plot produced.')
else:
  print('Wrong option for --diagplot. No plot produced.')


def isNotNaN(num):
    return num == num


# If input file name already has a valid extension (png or pdf): keep it 
if args.outfile.lower().endswith(('.png', '.pdf')):
  #print('\n ----> outfile name already has a valid format')
  figure = args.outfile
else:
  figure = args.outfile + '.' + ext
#print(' ----> figure = ', figure, '\n')


if os.path.isfile(srcTSList) and os.access(srcTSList, os.W_OK):
    #print("\nInput file " + srcTSList + " exists and has write persmissions.")
    pass
else:
    print("\nEither file " + srcTSList + " is missing or has not write persmissions.\n")
    sys.exit(1)
###################


#
# Read input FITS file
#
with warnings.catch_warnings():
  warnings.simplefilter("ignore", category=Warning)	# Ignore warnings about non-standard FITS units
  tsr = Table.read(srcTSList, hdu=1)			# Read FITS file table and converted to Table
  tsr_not_sorted = Table.read(srcTSList, hdu=1)
# Writing Keywords:
#
# tsr.meta['MY_KEYWD'] = 'my value'
# tsr.meta['COMMENT'] = ['First comment', 'Second comment', 'etc']
# tsr.write('my_table.fits', overwrite=True)
#

# Reading FITS keywords:
obs_id     = tsr.meta['OBS_ID']
instrument = tsr.meta['INSTRUME']
if instrument == 'EPN' or instrument == 'EMOS1' or instrument == 'EMOS2':
    try:
      ra_deg      = tsr.meta['SRC_RA']
      dec_deg     = tsr.meta['SRC_DEC']
      srcnum      = tsr.meta['SRCNUM']
      coord = SkyCoord(ra=ra_deg*u.degree, dec=dec_deg*u.degree)
      coords_hmsdms = coord.to_string('hmsdms')

    except KeyError:
      pass


# Reading Table data and sort by 'BACKV'
tsr.sort('BACKV')				# Sort Table by 'BACKV'

rate  = np.nan_to_num(tsr['RATE']) * np.nan_to_num(tsr['FRACEXP'] )		# rate = RATE * FRACEXP == rate weighted by the ExposureFraction
crate = (np.fmax(rate, rate*0)).cumsum()					# cummulative rate
erate = np.sqrt( ((np.nan_to_num(tsr['ERROR'])**2)*tsr['FRACEXP']).cumsum() )	# cummulative error

crate = np.ma.masked_array(crate, mask=(erate == 0))				# mask with 0 any NaN
erate = np.ma.masked_array(erate, mask=(erate == 0))				# mask with 0 any NaN

s2nc   = np.nan_to_num(crate/erate)						# cummulative S/N
imax   = s2nc.argmax()			# indices of the maximum values
s2nmax = s2nc[imax]			# 
backvmax = tsr['BACKV'][imax]		# Value of 'BACKV' for the maximum S/N

covered_time_fraction = 100.*(imax+1.)/len(s2nc)

s2nlast = s2nc[-1]
diff_max_last = 100*(abs(s2nmax - s2nlast)/s2nmax)

#print('====> ', s2nmax, s2nlast, diff_max_last)

# If diff max-last < 20 % OR max(s2nc) < 10 ==> Include all the the range == no-threshold
if diff_max_last < 20 or s2nmax < 10:
  backvmax = tsr['BACKV'][ isNotNaN(tsr['BACKV']) ][-1]
  #backvmax = tsr['BACKV'][-1]		# Value of 'BACKV' for the last S/N point
  covered_time_fraction = 100.		# Coverage = 100 %


# Remove entries where 'BACKV' = NaN

# Print result to stdout
#
#print ("Max. rate for maximum S/N =", backvmax," ; fraction of time covered =",100.*(imax+1.)/len(s2nc))
print('Bkg rate threshold cut = {0:.4f}; fraction of time covered = {1:.4f} %'.format(backvmax, covered_time_fraction))

# Write Keyword (FLCUTTHR = backvmax) in the original FITS
#
fits.setval(srcTSList, 'FLCUTTHR', value = backvmax, ext=1)


#-------------------------------------------------------------
# Plot
#
if diagPlotTest:

  f = plt.figure(1)	      	# Default figure size = 640x480 = 6.4 x 4.8 inch
  #print(f.get_size_inches())   # Return the figure size in inches
  #f.set_size_inches(8.0,6.0)    # ratio 8/6 = 640/480
  f.set_size_inches(8.0,10.0)
  #print(f.get_size_inches())   # 
  #print(f)

  plt.subplot(3, 1, 1)
  #plt.title('Flaring background filtering\nfor maximum S/N diagnostic plot')
  plt.plot(tsr['BACKV'], s2nc, ms=4, ls='-')
  plt.plot( backvmax, s2nmax, 'rs', fillstyle='none', ms=8)
  plt.axvline( x=backvmax, color='r', ls='--')
  plt.xlabel('{0}'.format('Background rate (BACKV)'))
  plt.ylabel('{0}'.format('Cummulative S/N'))
  plt.annotate('Max. Cumm. S/N: {0:.2f}'.format(s2nmax), xy=(0.4, 0.2), xycoords='axes fraction')
  plt.annotate('Diff. Max. Cumm. S/N vs last value: {0:.2f} %'.format(diff_max_last), xy=(0.4, 0.1), xycoords='axes fraction')
  
  time = tsr_not_sorted['TIME'] - tsr_not_sorted['TIME'][0]	# Ref. Time 0 at TimeSeries start time
  
  plt.subplot(3, 1, 2)
  plt.plot( time ,tsr_not_sorted['RATE'], ms=0.5, ls='-')
  plt.ylabel('{0}'.format('Source rate'))
  
  plt.subplot(3, 1, 3)
  plt.plot( time ,tsr_not_sorted['BACKV'], ms=0.5, ls='-')
  plt.axhline( y=backvmax, color='r', ls='--')
  plt.xlabel('{0}'.format('TIME (sec)'))
  plt.ylabel('{0}'.format('Background rate (BACKV)'))

  f.subplots_adjust(hspace=0.4) # hspace between plots
  
  plt.annotate('Flaring background filtering\nfor maximum S/N diagnostic plot', xy=(0.5, 0.95), xycoords='figure fraction', horizontalalignment='center', size=12)
  plt.annotate('OBS_ID: {0}'.format(obs_id),        xy=(0.13, 0.91), xycoords='figure fraction')
  plt.annotate('EPIC SRCNUM: {0}'.format(srcnum),   xy=(0.50, 0.91), xycoords='figure fraction')
  plt.annotate('RA,Dec: {0}'.format(coords_hmsdms), xy=(0.50, 0.89), xycoords='figure fraction')
  
  #ext='png'
  #figure='figure.' + ext
  plt.savefig( figure, format=ext)

  #plt.show()

#-------------------------------------------------------------