#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

##PERIOD OF BE LYN
##138.052368 minutes
##2.3008728  hours
##
##This amounts to about 34300 indices
##
UseDJS = 1
FilterSizeRef  = 3301
FluxRatios = FindFluxes_1c/FindFluxes_2c
FluxRatios[(where(isfinite(FluxRatios) == 0))[0]]=0

##Just do the bloody plot.  The zeros are too much to handle.
##Plot the ratios filtered
mplot.figure(0)
mplot.clf()
FluxRatiosAdj = FluxRatios*mean(FindFluxes_2c) #Scale back to normal
mplot.plot_date(DTKs_1, FluxRatiosAdj, 'co') 
FluxRatios_Med, FluxRatios_Mean  = \
  MedFilter(FluxRatios, FilterSize=2*FilterSizeRef, RejVal=0, UseDJS=UseDJS)
mplot.plot_date(DTKs_1, FluxRatios_Med*mean(FindFluxes_2c), 'co') 

##Now convert to magnitudes
Magnitude = 

mplot.figure(1)
mplot.clf()
mplot.plot_date(DTKs_1, FluxRatios_Mean, 'bo') 
NumTimes             = 53000
Times                = date2num(DTKs_1)-date2num(DTKs_1[0])
Times                = Times+FracSecs_1/(24*60.*60.)
TimesSec             = Times*24.*3600
FullTimeArr          = frange(NumTimes)
FullFluxRatioArr     = zeros(NumTimes, dtype=double)
FullFluxRatioMedArr  = zeros(NumTimes, dtype=double)
FullFluxRatioMeanArr = zeros(NumTimes, dtype=double)
TimeInt              = 0.24174
for i, TimeEl in zip(arange(NumTimes), FullTimeArr):
  ThisTime     = TimeEl*TimeInt
  MatchTimeInd = (where(abs(ThisTime - TimesSec) < 0.2))[0]
  if size(MatchTimeInd) == 2:
    FullFluxRatioArr[i]     = FluxRatios[MatchTimeInd[0]]
    FullFluxRatioMedArr[i]  = FluxRatios_Med[MatchTimeInd[0]]
    FullFluxRatioMeanArr[i] = FluxRatios_Mean[MatchTimeInd[0]]
  if size(MatchTimeInd) == 1:
    FullFluxRatioArr[i]     = FluxRatios[MatchTimeInd]
    FullFluxRatioMedArr[i]  = FluxRatios_Med[MatchTimeInd]
    FullFluxRatioMeanArr[i] = FluxRatios_Mean[MatchTimeInd]
  elif size(MatchTimeInd) == 0:
    FullFluxRatioArr[i]     = 0
    FullFluxRatioMedArr[i]  = 0
    FullFluxRatioMeanArr[i] = 0

##Plot the autocorrelation
mplot.figure(2)
b=mplot.acorr(FullFluxRatioArr-FullFluxRatioArr.mean())
c=mplot.acorr(FullFluxRatioMeanArr-FullFluxRatioMeanArr.mean())

##MatchTime
FullFluxRatioMeanArrSub = FullFluxRatioMeanArr[5000:52299]-mean(FullFluxRatioMeanArr[5000:52299])
FullFluxRatioArrSub     = FullFluxRatioArr[5000:52299]-mean(FullFluxRatioArr[5000:52299])

#The period should be about 34400 indices long, so search around this
NumPeriods = 5000
MatchInt   = 3000
BeginInt   = 0000
Periods    = arange(NumPeriods)+31000
PeriodCo1   = zeros(NumPeriods, dtype=float32)
PeriodCo2   = zeros(NumPeriods, dtype=float32)

for i, Period in zip(arange(NumPeriods), Periods):
  Power = FullFluxRatioMeanArrSub[BeginInt:(BeginInt+MatchInt)]*\
          FullFluxRatioMeanArrSub[(BeginInt+Period):(BeginInt+Period+MatchInt)]
  PeriodCo1[i] = Power.sum()
  Power = FullFluxRatioArrSub[BeginInt:(BeginInt+MatchInt)]*\
          FullFluxRatioArrSub[(BeginInt+Period):(BeginInt+Period+MatchInt)]
  PeriodCo2[i] = Power.sum()

##This gives us 34350

##SHOULD ACTUALLY BE 34264.