#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/ThesisLongPlotSettings.py')
PlotDir1 = '/nfs/slac/g/ki/ki03/lances/KPNO/Latex/Thesis/GuideMode/Figures/'
PlotDir2 = '/nfs/slac/g/ki/ki03/lances/KPNO/Latex/H2RG_GuideMode/Figures/'
SaveImg  = 1

import datetime
import WindowStarPhotometry_Extract_Database
reload(WindowStarPhotometry_Extract_Database)
from WindowStarPhotometry_Extract_Database import *

ThisFile  = 3
UseSinFit = 1
A         = 2.2        #Amplitude of oscillation for sinusoidal fit
f_osc     = 0.00810/4 #Frequency of oscillation in units of reads (4 reads ~ 1 second)

if ThisFile == 0:
  FileNameQS    = 'FILENAME REGEXP "VBELynVBELyn_H2RG_SIPIN_MultipleWindows2000_Reads_Dec18_2007_02_17_43"'
  ReadsMoveRA = [0,2000]
  SinOff = 0
  XCMLims_Y   = [18,32]
  YCMLims_Y   = [28,34]
elif ThisFile == 1:
  FileNameQS    = 'FILENAME REGEXP "VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_02_49_25"'
  ReadsMoveRA = [0,4000]
  SinOff = 0
  XCMLims_Y   = [18,32]
  YCMLims_Y   = [28,34]
elif ThisFile == 2:
  FileNameQS    = '(FILENAME REGEXP "VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_04_12_11" or FILENAME REGEXP "VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_04_28_32")'
  ReadsMoveRA = [3000,7000]
  SinOff      = 0
  XCMLims_Y   = [18,32]
  YCMLims_Y   = [28,34] 
elif ThisFile == 3:
  FileNameQS = 'FILENAME REGEXP "VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_03_23_08"' 
  ReadsMoveRA = [0,4000]
  SinOff = 80
  XCMLims_Y   = [18,32]
  YCMLims_Y   = [26,34]
elif ThisFile == 4:
  FileNameQS = 'FILENAME REGEXP "VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_03_06_25"'
  ReadsMoveRA = [0,4000]
  SinOff = 0
  XCMLims_Y   = [18,32]
  YCMLims_Y   = [28,34]
elif ThisFile == 5:
  FileNameQS = ' FILENAME REGEXP "VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_03_39_29"'
  ReadsMoveRA = [0,4000]
  SinOff = 0
  XCMLims_Y   = [18,32]
  YCMLims_Y   = [28,34]

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 AND '\
                                      + FileNameQS)

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 AND '\
                                      + FileNameQS) 

for i in arange(size(DTKs_1)):
  DTKs_1[i]=DTKs_1[i]+datetime.timedelta(hours=7)  
  #Convert the times from Mountain to UTC

#Here is the list of the times at which the telescope was moved manually
CumReadsMoveRA = [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)
DifFitsx_1  = []
DifFitsx_2  = []
DifFitsy_1  = []
DifFitsy_2  = []
 
for i in arange(1):
  TheseReads   = ReadsCum[ReadsMoveRA[i]:ReadsMoveRA[i+1]]
  TheseDTKs_1  = DTKs_1[ReadsMoveRA[i]:ReadsMoveRA[i+1]]
  TheseDTKs_2  = DTKs_2[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
  FirstWinFitx  = m1x*TheseReads+b1x
  if UseSinFit  == 1:
    FirstWinFitx = FirstWinFitx+A*sin(2*pi*f_osc*(-SinOff+TheseReads))
  DifFitsx_1     = append(DifFitsx_1, TheseXCMs_1-FirstWinFitx)
  SecWinFitx     = m2x*TheseReads+b2x
  if UseSinFit  == 1:
    SecWinFitx   = SecWinFitx+A*sin(2*pi*f_osc*(-SinOff+TheseReads))
  DifFitsx_2     = append(DifFitsx_2, TheseXCMs_2-SecWinFitx)
  FirstWinFity   = m1y*TheseReads+b1y
  DifFitsy_1     = append(DifFitsy_1, TheseYCMs_1-FirstWinFity)
  SecWinFity     = m2y*TheseReads+b2y
  DifFitsy_2     = append(DifFitsy_2, TheseYCMs_2-SecWinFity)

Fig=mplot.figure(0)
mplot.clf()
#mplot.subplot(2,1,1)
Fig.add_axes([0.075,0.56,0.9,0.41])
mplot.plot_date(TheseDTKs_1, TheseXCMs_1, 'r+')
mplot.plot_date(TheseDTKs_1, FirstWinFitx, 'k-',linewidth=2)
mplot.ylabel(r'\textbf{X Centroids (pixels)}')
mplot.ylim(XCMLims_Y)
mplot.xlim(DTKs_1.min(),DTKs_1.max())
ax0=gca()
Xlabs=ax0.get_xmajorticklabels()
Xlines = ax0.xaxis.get_majorticklines()
SizeLabs  = size(Xlabs)
SizeLines = size(Xlines)

for i in (arange(SizeLabs)):
  Xlabs[i].set_visible(False)
for i, line in zip(arange(SizeLines),ax0.xaxis.get_majorticklines()):
  if mod(i-2,8) == 0 or mod(i-3,8) == 0:
    line.set_markeredgewidth(3)
    line.set_markersize(8)
  else:
    line.set_markeredgewidth(1)
    line.set_markersize(4)
draw()

#mplot.subplot(2,1,2)
Fig.add_axes([0.075,0.09,0.9,0.41])
mplot.plot_date(TheseDTKs_1, TheseYCMs_1, 'r+', tz='UTC')
mplot.plot_date(TheseDTKs_1, FirstWinFity, 'k-', tz='UTC',linewidth=2)
mplot.ylabel(r'\textbf{Y Centroids (pixels)}')
mplot.xlabel(r'\textbf{Time}')
mplot.ylim(YCMLims_Y)
mplot.xlim(DTKs_1.min(),DTKs_1.max())
ax1=gca()
Xlabs=ax1.xaxis.get_majorticklabels()
Xlines = ax1.xaxis.get_majorticklines()
SizeLabs  = size(Xlabs)
SizeLines = size(Xlines)
for i, lab, in zip(arange(SizeLabs), ax1.xaxis.get_majorticklabels()):
  if mod(i+3,4) == 0:
    lab.set_visible(True)
  else: 
    lab.set_visible(False)

for i, line in zip(arange(SizeLines),ax1.xaxis.get_majorticklines()):
  if mod(i-2,8) == 0 or mod(i-3,8) == 0:
    line.set_markeredgewidth(3)
    line.set_markersize(8)
  else:
    line.set_markeredgewidth(1)
    line.set_markersize(4)
mplot.xlabel(r'\textbf{Time}')
draw()

if SaveImg == 1:
  mplot.savefig(PlotDir1+'TelescopeDrift_2_1m_FileNum'+str(ThisFile)+'.png')
  mplot.savefig(PlotDir1+'TelescopeDrift_2_1m_FileNum'+str(ThisFile)+'.eps')
  mplot.savefig(PlotDir2+'TelescopeDrift_2_1m_FileNum'+str(ThisFile)+'.png')
  mplot.savefig(PlotDir2+'TelescopeDrift_2_1m_FileNum'+str(ThisFile)+'.eps')