##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#Return_PSF_2_1m_TrackingErrors.py
#
#PURPOSE:
#  To numerically integrate the function believed to represent the 
#  telescope tracking errors and produce a PSF for the system.
#
#INPUTS:
#  m_dec: float    
#    Slope of tracking error in DEC (arcsec/sec)
#  m_ra: float
#    Slope of tracking error in RA  (arcsec/sec)
#  eps_see: float 
#    Seeing disc in pixels          (arcsec)
#  f_osc: float
#    Frequency of oscillation of RA (1/sec)
#  A_osc: float
#    Amplitude of oscillation of RA (arcsec)
#  t_exp: float 
#    Exposure time in seconds       (seconds)
#  I_o: float
#    Maximum intensity of star      (ADU)
#KEYWORDS:
#  XSize, YSize: ints
#    The x and y size of the box used to contain the stars
#  XOff, YOff, ints
#    The offset from 0 for each box
#  Mode: int
#    The mode to use for the radial average
#
#
from scipy.integrate import quad
from ReturnRadialAvg import *
from RadialFits import *
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')

def Return_PSF_2_1m_TrackingErrors(m_dec, m_ra, eps_see, A_osc, \
    f_osc, A_osc, t_exp, I_o, XSize=30, YSize=30, XOff=15, YOff=15, Mode = 0)
 
  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))
  
  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)
  
  if TvFlag == 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)
  
  return 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