#PlotVBELynPhotometry.py
#
#PURPOSE:
#  To extract the quantities from the WindowStarPhotometryStars database and analyze them 
#  quantitatively
#
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
execfile('/afs/slac/u/ki/lances/python/ThesisPlotSettings.py')

import WindowStarPhotometry_Extract_Database
reload(WindowStarPhotometry_Extract_Database)
from WindowStarPhotometry_Extract_Database import *
from sg_filter import *
from MedFilter import *
from scipy.signal import *

FrameNums_1, WindowNums_1, XCs_1, YCs_1, XCMs_1, YCMs_1, \
   Sharps_1, Rounds_1, Residuals_1, Mags_1, FindFluxes_1, AperFluxes_1, \
   Peaks_1, G1s_1, G2s_1, SigXs_1, SigYs_1, FWHMXs_1, FWHMYs_1, \
   DTKs_1, FracSecs_1=\
WindowStarPhotometry_Extract_Database('H2RG-32-147','ASIC','VBELynVBELyn', QueryString=' WHERE WINDOWNUM = 0')

FrameNums_2, WindowNums_2, XCs_2, YCs_2, XCMs_2, YCMs_2, \
   Sharps_2, Rounds_2, Residuals_2, Mags_2, FindFluxes_2, AperFluxes_2, \
   Peaks_2, G1s_2, G2s_2, SigXs_2, SigYs_2, FWHMXs_2, FWHMYs_2, \
   DTKs_2, FracSecs_2=\
WindowStarPhotometry_Extract_Database('H2RG-32-147','ASIC','VBELynVBELyn', QueryString=' WHERE WINDOWNUM = 1')

########################FLUXES 1###################################
##DO SOME FILTERING
FindFluxes_1c = FindFluxes_1.copy()
FindFluxes_1c[(where(FindFluxes_1c < 0))[0]]=0
FindFluxes_2c = FindFluxes_2.copy()
FindFluxes_2c[(where(FindFluxes_2c < 0))[0]]=0
FluxRatios = FindFluxes_1c/FindFluxes_2c
FluxRatios[(where(isfinite(FluxRatios) == 0))[0]]=0


UseDJS = 1
'''FilterSize     = 401
FindFluxes_1Med, FindFluxes_1Mean  = MedFilter(FindFluxes_1c, FilterSize=2*FilterSize, RejVal=0, UseDJS=UseDJS)
##FindFluxes_1mf = medfilt(FindFluxes_1c, FilterSize)

FilterSizeRef  = 401
FindFluxes_2Med, FindFluxes_2Mean  = MedFilter(FindFluxes_2c, FilterSize=2*FilterSizeRef, RejVal=0, UseDJS=UseDJS)
##FindFluxes_2mf = medfilt(FindFluxes_2c, FilterSizeRef)

mplot.figure(0)
mplot.clf()
mplot.subplot(2,1,1)
mplot.plot(FindFluxes_1Mean, 'go')
xlim(0,50000)
ylim(75000,105000)
mplot.subplot(2,1,2)
mplot.plot(FindFluxes_2Mean, 'bo')
xlim(0,50000)
ylim(11000,12500)
'''

##Plot the ratios filtered
mplot.figure(1)
mplot.clf()
FilterSizeRef  = 3301
FluxRatios = FindFluxes_1c/FindFluxes_2c
FluxRatios[(where(isfinite(FluxRatios) == 0))[0]]=0
FluxRatios_Med, FluxRatios_Mean  = \
  MedFilter(FluxRatios, FilterSize=2*FilterSizeRef, RejVal=0, UseDJS=UseDJS)
mplot.plot_date(DTKs_1, FluxRatios_Med, 'co') 

mplot.figure(2)
mplot.clf()
mplot.plot_date(DTKs_1, FluxRatios_Mean, 'bo') 

c=date2num(DTKs_1)-date2num(DTKs_1[0])
mplot.figure(3)
mplot.clf()
b=mplot.acorr(FluxRatios_Mean-FluxRatios_Mean.mean())

'''
mplot.plot(FluxRatios, 'go')
FluxRatios_mf = FindFluxes_1mf/FindFluxes_2mf
FluxRatios_mf[(where(isfinite(FluxRatios_mf) == 0))[0]]=0
mplot.plot(FluxRatios_mf, 'bo')
ylim(6,10)

########################FLUXES 2###################################
##DO SOME FILTERING
Mags_1c = Mags_1.copy()
Mags_1c[(where(Mags_1c < 0))[0]]=0
Mags_2c = Mags_2.copy()
Mags_2c[(where(Mags_2c < 0))[0]]=0

#Form the Coefficients for the Savgol filter
Mags_1m  = MedFilter(Mags_1c, FilterSize=2*FilterSize, RejVal=0, UseDJS=1)
Mags_1mf = medfilt(Mags_1c, FilterSize)
Mags_2m  = MedFilter(Mags_2c, FilterSize=2*FilterSize, RejVal=0, UseDJS=1)
Mags_2mf = medfilt(Mags_2c, FilterSize)

mplot.figure(2) 
mplot.clf() 
mplot.subplot(2,1,1) 
mplot.plot(Mags_1m, 'go')
mplot.plot(Mags_1mf, 'bo')
mplot.subplot(2,1,2) 
mplot.plot(Mags_2m, 'go')
mplot.plot(Mags_2mf, 'bo')
 
mplot.figure(3) 
mplot.clf() 
FluxRatios_m = Mags_1m/Mags_2m
FluxRatios_m[(where(isfinite(FluxRatios_m) == 0))[0]]=0
mplot.subplot(2,1,1)
mplot.plot(FluxRatios_m, 'go')
mplot.subplot(2,1,2)
mplot.plot(FluxRatios, 'bo')

FluxRatios_mf = Mags_1mf/Mags_2mf
FluxRatios_mf[(where(isfinite(FluxRatios_mf) == 0))[0]]=0
mplot.plot(FluxRatios_mf, 'bo')

mplot.figure(5)
mplot.subplot(2,1,1)
mplot.plot(FindFluxes_1, 'go')
mplot.subplot(2,1,2)
mplot.plot(FindFluxes_2, 'bo')

mplot.figure(6)
ScaleAdj = mean(FindFluxes_2)
FindFluxes_1_wRef = ScaleAdj*FluxRatios
mplot.plot(FindFluxes_1_wRef, 'go')

mplot.figure(7)
mplot.clf()
FindFluxes_1_wRefSmo = FindFluxes_1/FindFluxes_2m
FindFluxes_1_wRefSmo[(where(isfinite(FindFluxes_1_wRefSmo) == 0))[0]]=0
mplot.plot(FindFluxes_1_wRefSmo,'go')

##OTHER PLOTS OF RAW DATA
##FIRST MATCH THE TIMES SINCE THE LAST READ OF THE SECOND WINDOW GOT CUT OFF

mplot.figure(4)
mplot.clf()
mplot.plot(Mags_1,'ro')
mplot.plot(Mags_2,'bo')

mplot.figure(5)
mplot.clf()
mplot.plot(FindFluxes_1,'ro')
mplot.plot(FindFluxes_2,'bo')

mplot.figure(2)
mplot.clf()
mplot.plot_date(DTKs_1, FindFluxes_1,'ro')
mplot.plot_date(DTKs_2, FindFluxes_2,'bo')

mplot.figure(3)
mplot.clf()
mplot.plot(FindFluxes_1[0:50000]/FindFluxes_2[0:50000],'co')
ylim([0,12])

mplot.figure(7)
mplot.clf()
mplot.plot(FWHMXs_1, 'go')
mplot.plot(FWHMYs_1, 'ro')

mplot.figure(8)
mplot.clf()
mplot.plot(FWHMYs_1/FWHMXs_1, 'co')
mplot.plot(FWHMYs_2/FWHMXs_2, 'mo')
'''