#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')
PlotDir = '/nfs/slac/g/ki/ki04/lances/LSST/KPNO/Latex/Thesis/Photometry/Figures/'
import WindowStarPhotometry_Extract_Database
reload(WindowStarPhotometry_Extract_Database)
from WindowStarPhotometry_Extract_Database import *
from sg_filter import *
from MedFilter import *
from scipy.signal import *
import datetime

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')

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

###############################################################RAW DATA PLOT 1
##FIRST PLOT THE RAW MAGNITUDES:
##Correct the Magnitudes
RefAdj       = 3.8
Mags_1_Adj   = Mags_1-RefAdj
Mags_2_Adj   = Mags_2-RefAdj 
mplot.figure(1)
mplot.clf()

##Deal with subplot 1
mplot.subplot(2,1,1)
mplot.plot_date(DTKs_1, Mags_1_Adj)
ax1=gca()
Xlabs=ax1.xaxis.get_majorticklabels()
Xlines = ax1.xaxis.get_majorticklines()
StartInd = 5000
xlim(DTKs_1[StartInd], DTKs_1[size(DTKs_1)-1])
SizeLabs  = size(Xlabs)
SizeLines = size(Xlines)
for i, lab, in zip(arange(SizeLabs), ax1.xaxis.get_majorticklabels()):
    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)
ylim(8.6,9.1)

mplot.ylabel(r'\textbf{Magnitude}')
##Deal with subplot 2
mplot.subplot(2,1,2)
mplot.plot_date(DTKs_1, Mags_2_Adj)
ax1=gca()
Xlabs=ax1.xaxis.get_majorticklabels()
Xlines = ax1.xaxis.get_majorticklines()
StartInd = 5000
xlim(DTKs_1[StartInd], DTKs_1[size(DTKs_1)-1])
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)
ylim(10,12)
mplot.xlabel(r'\textbf{Time}')
mplot.ylabel(r'\textbf{Magnitude}')

mplot.savefig(PlotDir+'BELyn_Raw_Photometry.eps')
mplot.savefig(PlotDir+'BELyn_Raw_Photometry.png')

########################################################CORRECTED DATA PLOT 2
########################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
FindFluxes_1c[(where(isfinite(FluxRatios) == 0))[0]]=0
FindFluxes_2c[(where(isfinite(FluxRatios) == 0))[0]]=0
FluxRatios[(where(isfinite(FluxRatios) == 0))[0]]=0

##Do the adjustment after taking the ratio
FluxRatios_Adj = FluxRatios*mean(FindFluxes_2c) #Scale back to normal

##Filter the Adjusted Fluxes
FluxRatios_Med, FluxRatios_Mean  = \
  MedFilter(FluxRatios_Adj, FilterSize=2*FilterSizeRef, RejVal=0, UseDJS=UseDJS)

##Eliminate the zeros
GoodInds_Med   = (where(FluxRatios_Mean != 0))
FluxRatios_Med = FluxRatios_Med[GoodInds_Med]
DTKs_1_Med     = DTKs_1[GoodInds_Med]
GoodInds_Adj   = (where(FluxRatios_Adj != 0))
FluxRatios_Adj = FluxRatios_Adj[GoodInds_Adj]
DTKs_1_Adj     = DTKs_1[GoodInds_Adj]

##Convert to Magnitudes
##ADJUST TO STANDARD MAGS 
##BE LYN                   V = 8.8 
##REF STAR GSC 03425-00544 V = 11.1
RefAdj   = 21.2
Mags_Med = -2.5*log10(FluxRatios_Med)+RefAdj
Mags_Adj = -2.5*log10(FluxRatios_Adj)+RefAdj

##PERIOD INFO
##P=0.0958697 days measured by OJA 1985
##
##According to formula in Liu 1994
##P=
##In [14]:run ../functions/jd2gd 2449018.2684
##
## Task to convert a list of julian dates to gregorian dates.
## Written by Enno Middelberg 2002
##
##2449018.2684 = January 30, 1993,  18:26:29.75999504327774 UTC
##
##Observations made on December 17, 2007:
##To get the amount of elapsed periods (E)
##IDL> print, (15*365.25+180)/0.0958696
##      59025.5
##And now to get the change in period Delta P
##IDL> print, 59025.5^2*8.8e-12*.5+0.09586963
##     0.111199
##So Delta P 

##2.3008728 hours
##138.0523 minutes
##8283.14  seconds
## 02:18:03 

##
##Plot : Period changes based upon length of time between files
Period   = 24980
StartInd = 5000
MarkX1   = 10000
MarkX2   = 35226
mplot.figure(0)
mplot.clf()
mplot.plot(DTKs_1_Adj[StartInd:], Mags_Adj[StartInd:], 'co')
mplot.plot(DTKs_1_Med[StartInd:], Mags_Med[StartInd:], 'r-', linewidth=4)
mplot.plot_date([DTKs_1_Adj[MarkX1],DTKs_1_Adj[MarkX1]], [8.5, 9.3], 'k--', linewidth=2)
mplot.plot_date([DTKs_1_Adj[MarkX2],DTKs_1_Adj[MarkX2]], [8.5, 9.3], 'k--', linewidth=2)
xlim(DTKs_1[StartInd], DTKs_1[size(DTKs_1)-1])
#xlim(DTKs_1_Adj[StartInd], DTKs_1_Adj[size(DTKs_1_Adj)-1])
ylim(8.5,9.3)
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}')
mplot.ylabel(r'\textbf{Magnitude}')
draw()

mplot.savefig(PlotDir+'BELyn_ReferenceCorrected_Photometry.eps')
mplot.savefig(PlotDir+'BELyn_ReferenceCorrected_Photometry.png')