#EllipticityVsTrackingError.py
#
#PURPOSE:
#  To plot the expected ellipticity for a variety of tracking error rates
#  based upon the numerical integration performed in Return_PSF_2_1m_TrackingErrors
# 
from Return_PSF_2_1m_TrackingErrors import *
from scipy.integrate import quad
from ReturnRadialAvg import *
from RadialFits import *

#NOTES ON POINTING OF TELESCOPE
#
#HERE ARE THE ZENITH DISTANCES FOR VBELYn
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows1000_Reads_Dec18_2007_02_10_52.fits ,       25.5200
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows100_Reads_Dec18_2007_02_17_07.fits ,       24.5300
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows2000_Reads_Dec18_2007_02_17_43.fits ,       24.4400
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows2000_Reads_Dec18_2007_02_25_58.fits ,       23.1700
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows3000_Reads_Dec18_2007_02_37_04.fits ,       21.5300
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_02_49_25.fits ,       19.8000
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_03_06_25.fits ,       17.6700
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_03_23_08.fits ,       15.9500
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_03_39_29.fits ,       14.7600
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_03_55_50.fits ,       14.1900
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_04_12_11.fits ,       14.3100
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_04_28_32.fits ,       15.1100
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_04_44_53.fits ,       16.4800
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_05_01_14.fits ,       18.3000
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_05_17_34.fits ,       20.4300
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_05_33_55.fits ,       22.8000
#File: VBELynVBELyn_H2RG_SIPIN_MultipleWindows4000_Reads_Dec18_2007_05_50_16.fits ,       25.3300
#
#File: SAO5417_I_NegPerAfterDark_H1RG_SIPIN_2500_Reads_Nov19_2007_23_20_46.fits ,       15.0500
#File: SAO5417_I_NegPerDark_H1RG_SIPIN_2500_Reads_Nov19_2007_22_45_50.fits ,       7.81000
#File: SAO5417_I_NegPerDead_H1RG_SIPIN_2500_Reads_Nov19_2007_22_15_46.fits ,       2.29000
#File: SAO5417_I_NegPerSatRamp_H1RG_SIPIN_500_Reads_Nov19_2007_22_09_16.fits ,       1.90000

execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
Mode    = 0           #Mode to use for radial profile
TvFlag  = 9
BoxSize = 15
XSize   = 51
YSize   = 51
XOff    = 25
YOff    = 25

Guiding = 1
if Guiding == 1:
  #GOOD SET FOR E=.18
  m_dec   =  0.000006   #Slope of tracking error in DEC (arcsec/sec)
  m_ra    =  0.000007   #Slope of tracking error in RA  (arcsec/sec)
  A_osc   =  1.4104     #Amplitude of oscillation of RA (arcsec)
  f_osc   =  0.00810    #Frequency of oscillation of RA (1/sec)
  I_o     =  10.00      #Maximum intensity of star      (ADU)
  
elif Guiding == 0:
  #GOOD SET FOR E=0.06
  m_dec   = -0.0042     #Slope of tracking error in DEC (arcsec/sec)
  m_ra    =  0.0024     #Slope of tracking error in RA  (arcsec/sec)
  f_osc   =  0.00810    #Frequency of oscillation of RA (1/sec)
  A_osc   =  3.1004     #Amplitude of oscillation of RA (arcsec)
  I_o     =  10.00      #Maximum intensity of star      (ADU)


SeeVals      = [4.0, 4.5, 5.0, 5.5, 6.0, 6.5]
m_DECs       = [-0.0042, -0.0010, 0.0005]
m_RAs        = [0.0024,  0.0060,  0.0040]
NumSees      = size(SeeVals)
TotalExpTime = 1800.  #Total exposure time in seconds
NumTimes     = 4.
ExpTimes     = (arange(NumTimes)+1)*TotalExpTime/NumTimes
Heights      = zeros([NumTimes, NumSees], dtype=float32)
MuXs         = zeros([NumTimes, NumSees], dtype=float32)
MuYs         = zeros([NumTimes, NumSees], dtype=float32)
FWHMXs       = zeros([NumTimes, NumSees], dtype=float32)
FWHMYs       = zeros([NumTimes, NumSees], dtype=float32)
FWHMs        = zeros([NumTimes, NumSees], dtype=float32)
Es           = zeros([NumTimes, NumSees], dtype=float32)
EAs          = zeros([NumTimes, NumSees], dtype=float32)
FWHM_Gs      = zeros([NumTimes, NumSees], dtype=float32)
FWHM_G_Ms    = zeros([NumTimes, NumSees], dtype=float32)

for eps_see, s in zip(SeeVals, arange(NumSees)):
  for t_exp, i in zip(ExpTimes, arange(NumTimes)):
    Height, MuX, MuY, FWHMX, FWHMY, FWHM, E, EA, \
    Amp_G, FWHM_G, Fit_G, Amp_G_M, FWHM_G_M, Beta_G_M, Fit_G_M = \
    Return_PSF_2_1m_TrackingErrors(m_dec, m_ra, eps_see, f_osc, A_osc, t_exp, I_o, Mode = 0, \
                               TvFlag=TvFlag, XSize=XSize, YSize=YSize, \
                               XOff=XOff, YOff=YOff, \
                               BoxSize=BoxSize)
    Heights[i,s]      = Height
    MuXs[i,s]         = MuX
    MuYs[i,s]         = MuY
    FWHMXs[i,s]       = FWHMX
    FWHMYs[i,s]       = FWHMY
    FWHMs[i,s]        = FWHM
    Es[i,s]           = E
    EAs[i,s]          = EA
    FWHM_G            = 2*sqrt(2*log(2))*FWHM_G
    FWHM_GStr         = "%.2f"   % FWHM_G
    FWHM_G_M          = 2*sqrt(2**(-1/Beta_G_M)-1)*FWHM_G_M
    FWHM_G_MStr       = "%.2f"   % FWHM_G_M
    print 'FWHM_G_M ' + str(FWHM_G_M)
    print 'FWHM_G_M ' + str(FWHM_G)
    FWHM_Gs[i,s]      = FWHM_G
    FWHM_G_Ms[i,s]    = FWHM_G_M

mplot.figure(0)
mplot.clf()
for i in arange(NumSees):
  mplot.plot(ExpTimes, Es[:,i], 'o', label = r'$\epsilon$ : '+ str(SeeVals[i]))
mplot.legend(numpoints=1)

mplot.figure(1)
mplot.clf()
for i in arange(NumTimes):
  mplot.plot(SeeVals, Es[i,:], 'o', label = r'$t_{exp} = $' + str(ExpTimes[i]))
mplot.legend(numpoints=1)

mplot.figure(2)
mplot.clf()
for i in arange(NumSees):
  mplot.plot(ExpTimes, FWHMXs[:,i], 'o', label = r'$\epsilon$ X: '+ str(SeeVals[i]))
  mplot.plot(ExpTimes, FWHMYs[:,i], 's', label = r'$\epsilon$ Y: '+ str(SeeVals[i]))

mplot.figure(3)
mplot.clf()
for i in arange(NumSees):
  mplot.plot(ExpTimes, FWHM_Gs[:,i], 'o')
  mplot.plot(ExpTimes, FWHM_G_Ms[:,i], 's')