##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#Simulate_PSF_2_1m_TrackingErrors.py
#
#PURPOSE:
#  To numerically integrate the function believed to represent the 
#  telescope tracking errors.
from scipy.integrate import quad
from ReturnRadialAvg import *
from RadialFits import *
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
Mode    = 0           #Mode to use for radial profile

#GOOD SET FOR E=.18
#m_dec   = 0.00176    #Slope of tracking error in DEC (arcsec/sec)
#m_ra    = 0.00277    #Slope of tracking error in RA  (arcsec/sec)
#eps_see = 6.1         #Seeing disc in pixels          (arcsec)
#A_osc   = 2.2104      #Amplitude of oscillation of RA (arcsec)

#GOOD SET FOR E=0.06
m_dec   = 0.00006     #Slope of tracking error in DEC (arcsec/sec)
m_ra    = 0.00007     #Slope of tracking error in RA  (arcsec/sec)
eps_see = 4.31         #Seeing disc in pixels          (arcsec)
f_osc   = 0.00810     #Frequency of oscillation of RA (1/sec)
A_osc   = 1.3504      #Amplitude of oscillation of RA (arcsec)
t_exp   = 1800.       #Exposure time in seconds       (seconds)
I_o     = 1.          #Maximum intensity of star      (ADU)

PSFFunc = lambda t_exp, x, y, I_o, m_dec, m_ra, f_osc, A_osc, eps_see : \
          I_o*exp(-(x-m_ra*t_exp-A_osc*sin(2*pi*f_osc*t_exp))**2/(2*(eps_see/2.35)**2))\
             *exp(-(y-m_dec*t_exp)**2/(2*(eps_see/2.32)**2))

XSize = 30
XOff  = 15
YSize = 30
YOff  = 15
XVals = frange(XSize)-XOff
YVals = frange(YSize)-YOff

Image = zeros([YSize,XSize], dtype=float)
XCoo  = zeros([YSize,XSize], dtype=int16)
YCoo  = zeros([YSize,XSize], dtype=int16)
for i, Xc in zip(arange(XSize), XVals):
  for j, Yc in zip(arange(YSize), YVals):
    Int = quad(PSFFunc, 0, t_exp, args=(Xc,Yc, I_o, m_dec, m_ra, f_osc, A_osc, eps_see))
    Image[j,i] = Int[0]
    XCoo[j,i]  = i
    YCoo[j,i]  = j

MaxIm         = Image.max()
MaxImLocX     = where(Image == MaxIm)[0] 
MaxImLocY     = where(Image == MaxIm)[1]
BoxSize       = 7
XCStart       = MaxImLocX[0]-BoxSize
XCStop        = MaxImLocX[0]+BoxSize
YCStart       = MaxImLocY[0]-BoxSize
YCStop        = MaxImLocY[0]+BoxSize
#Find the new Centroid
Centroid      = ReturnCentroid(Image, XCStart, XCStop, YCStart, YCStop)  
XStart        = Centroid[0]
XStop         = Centroid[1]
YStart        = Centroid[2]
YStop         = Centroid[3]
IntArr3       = ReturnRadialAvg(Image[YStart:YStop, XStart:XStop], Mode=0)
Radii, IntArr = ReturnRadialAvg(Image[YStart:YStop, XStart:XStop], Mode=1, DoCentroid=0)

#Do some radial fits
Height, MuX, MuY, FWHMX, FWHMY, FWHM, E, EA = \
         moments(IntArr, SubSize=1)
Amp_G,   FWHM_G,   Fit_G = RadialGaussFits(Radii, IntArr, IncFitData=1)
Amp_G_M, FWHM_G_M, Beta_G_M, Fit_G_M = RadialMoffatFits(\
    Radii, IntArr, IncFitData=1)


#Show the image
mplot.figure(0)
mplot.clf()
mplot.matshow(Image, extent=[YVals.min(), YVals.max(), XVals.min(), XVals.max()], \
              fignum=0)
mplot.colorbar()

#Show the radial profile
mplot.figure(1)
mplot.clf()
mplot.plot(Radii, IntArr, 'ro')
mplot.plot(Radii, Fit_G, 'k-') 
print Image.min()
print Image.max()
FWHM_Gauss = 2*sqrt(2*log(2))*FWHM_G
FWHM_Moff  = 2*sqrt(2**(-1/Beta_G_M)-1)*FWHM_G_M
print 'Gauss FWHM  : ' + str(FWHM_Gauss)
print 'Moffat FWHM : ' + str(FWHM_Moff)