#PlotVBELynFourierSpectra.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 numpy import *

ThisFile  = 1
if ThisFile == 0:
  FileName    = 'VBELynVBELyn_H2RG_SIPIN_MultipleWindows2000_Reads_Dec18_2007_02_17_43'
  ReadsMoveRA = [0,2000]
elif ThisFile == 1:
  FileName    = 'VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_02_49_25'
  ReadsMoveRA = [0,4000]

#Figure out the relevant numbers for frequency spacing
NumSamp   = ReadsMoveRA[1]-ReadsMoveRA[0]
ReadsCum  = arange(NumSamp)
SampInt   = 0.25                               #Sampling interval is about 4 Hz 
SampTimes = arange(NumSamp)*SampInt            #Times of sampling
FSpac     = 1./(NumSamp*SampInt)               #Frequency Spacing
KVals     = FSpac*(-NumSamp/2+arange(NumSamp)) #KVals for full FFT
KValsReal = FSpac*arange(NumSamp/2+1)          #KVals for real FFT

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 '\
                                      +'FILENAME REGEXP ' + '"' + FileName +'"')

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 '\
                                      +'FILENAME REGEXP ' + '"' + FileName + '"')

#Here is the list of the times at which the telescope was moved manually
CumReadsMoveRA = [0,        2000,     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)
m1s         = zeros(Moves, dtype=float32)
m2s         = zeros(Moves, dtype=float32)
b1s         = zeros(Moves, dtype=float32)
b2s         = zeros(Moves, dtype=float32)
DifFits_1   = []
DifFits_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.figure(3)
mplot.clf()

for i in arange(1):
  TheseReads   = ReadsCum[ReadsMoveRA[i]:ReadsMoveRA[i+1]]
  TheseXCMs_1  = XCMs_1[ReadsMoveRA[i]:ReadsMoveRA[i+1]]
  TheseXCMs_2  = XCMs_2[ReadsMoveRA[i]:ReadsMoveRA[i+1]]
  m1, b1, Error1, CosE1 = Linefit(TheseReads, TheseXCMs_1, TvFlag=0, RejectCos = 0)
  m2, b2, Error2, CosE2 = Linefit(TheseReads, TheseXCMs_2, TvFlag=0, RejectCos = 0)
  m1s[i] = m1 ; m2s[i] = m2 ; b1s[i] = b1 ; b2s[i] = b2

  FirstWinFit = m1*TheseReads+b1
  DifFits_1    = append(DifFits_1, TheseXCMs_1-FirstWinFit)
  SecWinFit   = m2*TheseReads+b2
  DifFits_2    = append(DifFits_2, TheseXCMs_2-SecWinFit)

  mplot.figure(0)
  mplot.plot(TheseReads, FirstWinFit, 'k-')
  mplot.figure(1)
  mplot.plot(TheseReads, SecWinFit, 'k-')

mplot.figure(2)
mplot.plot(DifFits_1, 'mo')
mplot.figure(3)
mplot.plot(DifFits_2, 'mo')

mplot.figure(4)
DifFitsACor1 = mplot.acorr(DifFits_1)
DifFitsACor2 = mplot.acorr(DifFits_2)

mplot.figure(5)
mplot.plot(KValsReal, rfft(DifFits_1).real**2)

mplot.figure(6)
mplot.clf()
SizeDifFitsACor1 = size(DifFitsACor1[1])
DifFitsACor1Real = DifFitsACor1[1][SizeDifFitsACor1/2:]
mplot.plot(KValsReal, rfft(DifFitsACor1Real).real**2, 'ro')

mplot.figure(7)
mplot.plot(YCMs_1, 'y.')