#WindowStarPhotometryStats.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 FitSinusoid import *
from histOutline 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')

#Here is the list of the times at which the telescope was moved manually
ReadsMoveRA = [0,        4000,     9844,     13961,    14615,    19436,    22638,    23912,    26380,  \
               27627,    29695,    34480,    37376,    41393,    44659] 
#TimesMoveRA = [03:12:20, 03:30:16, 03:35:32, 03:43:00, 04:08:48, 04:03:12]
FrameNumsRA = [0,        0,        2844,     2961,     3615,     438,      3639,     912,      3380,   \
               627,      2695,     3480,     2375,     2393,     1659]
ReadsCum    = arange(50000)
Moves       = size(ReadsMoveRA)
VMoves      = 12 #Number of valid moves before sun started to rise
m1xs        = zeros(VMoves, dtype=float32)
m2xs        = zeros(VMoves, dtype=float32)
b1xs        = zeros(VMoves, dtype=float32)
b2xs        = zeros(VMoves, dtype=float32)
m1ys        = zeros(VMoves, dtype=float32)
m2ys        = zeros(VMoves, dtype=float32)
b1ys        = zeros(VMoves, dtype=float32)
b2ys        = zeros(VMoves, dtype=float32)
DTK1s_Beg   = []
DTK2s_Beg   = []
DifFitsx_1  = []
DifFitsx_2  = []
DifFitsx_12 = []
DifFitsx_22 = []
DifFitsy_1  = []
DifFitsy_2  = []
 
#Go through and plot
mplot.figure(0)
mplot.clf()
mplot.plot(XCMs_1,'go')
mplot.hold(1)

mplot.figure(1)
mplot.clf()
mplot.plot(XCMs_2,'go')
mplot.hold(1)

mplot.figure(2)
mplot.clf()
mplot.plot(YCMs_1,'go')

mplot.figure(3)
mplot.clf()
mplot.plot(YCMs_2,'go')

for i in arange(VMoves):
  TheseReads   = ReadsCum[ReadsMoveRA[i]:ReadsMoveRA[i+1]]
  #X Component
  TheseXCMs_1  = XCMs_1[ReadsMoveRA[i]:ReadsMoveRA[i+1]]
  TheseXCMs_2  = XCMs_2[ReadsMoveRA[i]:ReadsMoveRA[i+1]]
  m1x, b1x, Error1x, CosE1x = Linefit(TheseReads, TheseXCMs_1, TvFlag=0, RejectCos = 0)
  m2x, b2x, Error2x, CosE2x = Linefit(TheseReads, TheseXCMs_2, TvFlag=0, RejectCos = 0)
  m1xs[i] = m1x ; m2xs[i] = m2x ; b1xs[i] = b1x ; b2xs[i] = b2x
  #Y Component
  TheseYCMs_1  = YCMs_1[ReadsMoveRA[i]:ReadsMoveRA[i+1]]
  TheseYCMs_2  = YCMs_2[ReadsMoveRA[i]:ReadsMoveRA[i+1]]
  m1y, b1y, Error1y, CosE1y = Linefit(TheseReads, TheseYCMs_1, TvFlag=0, RejectCos = 0)
  m2y, b2y, Error2y, CosE2y = Linefit(TheseReads, TheseYCMs_2, TvFlag=0, RejectCos = 0)
  m1ys[i] = m1y ; m2ys[i] = m2y ; b1ys[i] = b1y ; b2ys[i] = b2y

  #Calculate Fits and do Difference
  p = [2.4, 0.002, pi/2, 0.0]
  FirstWinFitx  = m1x*TheseReads+b1x
  DifFitsx_1    = append(DifFitsx_1, TheseXCMs_1-FirstWinFitx)
  FitDiff1      = TheseXCMs_1-FirstWinFitx
  SinFit1       = FitSinusoid(p, TheseReads, FitDiff1)
  Resids_x1     = SinFit1[0]*sin(2*pi*SinFit1[1]*TheseReads-SinFit1[2])-\
                  SinFit1[3]
  DifFitsx_12   = append(DifFitsx_12, FitDiff1-Resids_x1)

  SecWinFitx    = m2x*TheseReads+b2x
  FitDiff2      = TheseXCMs_2-SecWinFitx
  SinFit2       = FitSinusoid(p, TheseReads, FitDiff2)
  DifFitsx_2    = append(DifFitsx_2, TheseXCMs_2-SecWinFitx)
  Resids_x2     = SinFit2[0]*sin(2*pi*SinFit2[1]*TheseReads-SinFit2[2])-\
                  SinFit2[3]
  DifFitsx_22   = append(DifFitsx_22, FitDiff2-Resids_x2)
  FirstWinFity  = m1y*TheseReads+b1y
  DifFitsy_1    = append(DifFitsy_1, TheseYCMs_1-FirstWinFity)
  SecWinFity    = m2y*TheseReads+b2y
  DifFitsy_2    = append(DifFitsy_2, TheseYCMs_2-SecWinFity)

  DTK1s_Beg     = append(DTK1s_Beg, DTKs_1[ReadsMoveRA[i]])
  DTK2s_Beg     = append(DTK2s_Beg, DTKs_2[ReadsMoveRA[i]])

  mplot.figure(0)
  mplot.plot(TheseReads, FirstWinFitx, 'k-')
  mplot.figure(1)
  mplot.plot(TheseReads, SecWinFitx, 'k-')
  mplot.figure(2)
  mplot.plot(TheseReads, FirstWinFity, 'k-')
  mplot.figure(3)
  mplot.plot(TheseReads, SecWinFity, 'k-')

#Histograms
NBins   = 100
BinMax  = 5
BinMin  = -5. 
BinsIn  = BinMin+arange(NBins)/(BinMax-BinMin)  
Bins1x, N1x = histOutline(DifFitsx_12, binsIn = BinsIn)
Bins2x, N2x = histOutline(DifFitsx_22, binsIn = BinsIn)
mplot.figure(10)
mplot.clf()
mplot.plot(Bins1x, N1x, 'k-')
mplot.plot(Bins2x, N2x, 'r-')

BinsIn  = BinMin+arange(NBins)/(BinMax-BinMin)
Bins1x, N1y = histOutline(DifFitsy_1, binsIn = BinsIn)
Bins2x, N2y = histOutline(DifFitsy_2, binsIn = BinsIn) 
mplot.figure(11) 
mplot.clf() 
mplot.plot(Bins1x, N1y, 'k-') 
mplot.plot(Bins2x, N2y, 'r-') 


PlotDifs=1
if PlotDifs == 1:
  mplot.figure(4)
  mplot.plot(DifFitsx_1, 'mo')
  mplot.clf()
  mplot.plot(DifFitsx_12, 'mo')
  mplot.figure(5)

  mplot.plot(DifFitsx_2, 'mo')
  mplot.clf()
  mplot.plot(DifFitsx_22, 'mo')
  mplot.figure(6)
  mplot.plot(DifFitsy_1, 'mo')
  mplot.figure(7)
  mplot.plot(DifFitsy_2, 'mo')

mplot.figure(8)
mplot.plot_date(DTKs_1, YCMs_1, 'bo')
mplot.figure(9)
mplot.plot_date(DTKs_1, XCMs_1, 'bo')


print '--------------------------------------------------------------------'
print 'Drift Errors -------------------------------------------------------'
print 'Times of Exposure : ' 
print DTK1s_Beg
print 'X Slopes for Window 1'
print m1xs*4
print 'Abs Average X Slopes for Window 1'
print mean(abs(m1xs*4))
print 'Min / Max X Slopes for Window 1'
print str(min(abs(m1xs*4))) + ' / ' + str(max(abs(m1xs*4)))
print ''
print 'X Slopes for Window 2'
print m2xs*4
print 'Abs Average X Slopes for Window 1'
print mean(abs(m2xs*4))
print 'Min / Max X Slopes for Window 2'
print str(min(abs(m2xs*4))) + ' / ' + str(max(abs(m2xs*4)))
print 
print ''
print 'Y Slopes for Window 1'
print m1ys*4
print 'Abs Average Y Slopes for Window 1'
print mean(abs(m1ys*4))
print 'Min / Max Y Slopes for Window 1'
print str(min(abs(m1ys*4))) + ' / ' + str(max(abs(m1ys*4)))
print ''
print 'Y Slopes for Window 2'
print m2ys*4
print 'Abs Average Y Slopes for Window 1'
print mean(abs(m2ys*4))
print 'Min / Max Y Slopes for Window 2'
print str(min(abs(m2ys*4))) + ' / ' + str(max(abs(m2ys*4)))