;return_screenellipticity
;
;PURPOSE:
;	This script takes 3 25x25 arrays of floating point numbers that 
;	represent the tilts of a wavefront across the ShackHartmann wavefront 
;	sensor and the phase across
;	the SOAR telescope and returns the ellipticity of a star that would 
;	be formed by the wavefront in the case of geometric optics.
;
;       It returns the Ellipticity (filtered and unfiltered) scalar number and 
;       direction, and power spectra along the kx=0 and ky=0 directions,
;
;Return:  EllipVals(0,*)=High Passed Ellipticities
;	  EllipVals(1,*)=High Passed Pointing Angles
;	  EllipVals(2,*)=Low Passed Ellipticities
;	  EllipVals(3,*)-Low Passed Pointing Angles
;         EllipVals(4,*)=Power Spectra of XSlopes along ky=0
;         EllipVals(5,*)=Power Spectra of YSlopes along kx=0
;         EllipVals(4,*)=Power Spectra of Phase along ky=0
;         EllipVals(5,*)=Power Spectra of Phase along kx=0

function return_screenellipticity,XSlopes,YSlopes,PhaseScreenSoar,FittedOnly,$
	                        Num,Fit
KMax=Num/2
EllipVals=fltarr(8,KMax)     ;Array to be returned

XSlopesLow=fltarr(num,num)
XSlopesHigh=fltarr(num,num)
XSlopesFft=complexarr(num,num)
XSlopesFftRe=complexarr(num,num)
XSlopesFftHigh=complexarr(num,num)
XSlopesFftLow=complexarr(num,num)
YSlopesFftHighRe=complexarr(num,num)
XSlopesFftLowRe=complexarr(num,num)
XSlopesFft(*,*)=fft(XSlopes(*,*))

YSlopesLow=fltarr(num,num) 
YSlopesHigh=fltarr(num,num)
YSlopesFft=complexarr(num,num)
YSlopesFftRe=complexarr(num,num)
YSlopesFftHigh=complexarr(num,num)
YSlopesFftLow=complexarr(num,num)
YSlopesFftHighRe=complexarr(num,num)
YSlopesFftLowRe=complexarr(num,num)
YSlopesFft(*,*)=fft(YSlopes(*,*))

PhaseScreenFft=complexarr(num,num)
PhaseScreenFftRe=complexarr(num,num)
PhaseScreenFft(*,*)=fft(PhaseScreenSoar)

XSlopesFftRe=Shift(XSlopesFft,KMax,KMax)
YSlopesFftRe=Shift(YSlopesFft,KMax,KMax)
PhaseScreenFftRe=Shift(PhaseScreenFft,KMax,KMax)
 
EllipVals(4,*)=(abs(XSlopesFftRe(KMax+1:2*Kmax,KMax)))^2
EllipVals(5,*)=(abs(YSlopesFftRe(KMax,Kmax+1:2*KMax)))^2
EllipVals(6,*)=(abs(PhaseScreenFftRe(KMax,KMax+1:2*KMax)))^2
EllipVals(7,*)=(abs(PhaseScreenFftRe(KMax+1:2*KMax,KMax)))^2

;BEGIN FILTERING
for kval=0,kmax-1 do begin
  XSlopesFftHighRe=XSlopesFftRe
  YSlopesFftHighRe=YSlopesFftRe
  XSlopesFftHighRe((KMax-kval):(KMax+kval),(KMax-kval):(KMax+kval))=0
  YSlopesFftHighRe((KMax-kval):(KMax+kval),(KMax-kval):(KMax+kval))=0
  YSlopesFftLowRe=complexarr(num,num)
  XSlopesFftLowRe=complexarr(num,num)
  XSlopesFftLowRe((KMax-kval):(KMax+kval),(KMax-kval):(KMax+kval))=$
         XSlopesFftRe((KMax-kval):(KMax+kval),(KMax-kval):(KMax+kval))
  YSlopesFftLowRe((KMax-kval):(KMax+kval),(KMax-kval):(KMax+kval))=$
         YSlopesFftRe((KMax-kval):(KMax+kval),(KMax-kval):(KMax+kval))

  XSlopesFftLow=Shift(XSlopesFftLowRe,-KMax,-KMax)
  YSlopesFftLow=Shift(YSlopesFftLowRe,-KMax,-KMax)
  XSlopesFftHigh=Shift(XSlopesFftHighRe,-KMax,-KMax)
  YSlopesFftHigh=Shift(YSlopesFftHighRe,-KMax,-KMax)

  XSlopesLow(*,*)=real_part(fft(XSlopesFftLow,/inverse))
  YSlopesLow(*,*)=real_part(fft(YSlopesFftLow,/inverse))
  XSlopesHigh(*,*)=real_part(fft(XSlopesFftHigh,/inverse))
  YSlopesHigh(*,*)=real_part(fft(YSlopesFftHigh,/inverse))
  
  ;Calculate High Pass Values   
  TimeE1NumHigh=XSlopesHigh*XSlopesHigh-YSlopesHigh*YSlopesHigh
  TimeE2NumHigh=2*XSlopesHigh*YSlopesHigh
  TimeEDenHigh=XSlopesHigh^2+YSlopesHigh^2

  if FittedOnly eq 1 then begin
     FrameEDenHigh=mean(TimeEDenHigh)
     FrameE1High=mean(TimeE1NumHigh)/mean(TimeEDenHigh)
     FrameE2High=mean(TimeE2NumHigh)/mean(TimeEDenHigh)
  endif else if FittedOnly eq 2 then begin
     FrameEDenHigh=mean(TimeEDenHigh(Fit))
     FrameE1High=mean(TimeE1NumHigh(Fit))/mean(TimeEDenHigh(Fit))
     FrameE2High=mean(TimeE2NumHigh(Fit))/mean(TimeEDenHigh(Fit))
  endif
  FrameEHigh=sqrt(FrameE1High^2+FrameE2High^2)
  FrameThetaHigh=.5*atan(FrameE2High,FrameE1High)

  ;Calculate Low Pass Values
  TimeE1NumLow=XSlopesLow*XSlopesLow-YSlopesLow*YSlopesLow
  TimeE2NumLow=2*XSlopesLow*YSlopesLow
  TimeEDenLow=XSlopesLow^2+YSlopesLow^2

  if FittedOnly eq 1 then begin
     FrameEDenLow=mean(TimeEDenLow)
     FrameE1Low=mean(TimeE1NumLow)/mean(TimeEDenLow)
     FrameE2Low=mean(TimeE2NumLow)/mean(TimeEDenLow)
  endif else if FittedOnly eq 2 then begin
     FrameEDenLow=mean(TimeEDenLow(Fit))
     FrameE1Low=mean(TimeE1NumLow(Fit))/mean(TimeEDenLow(Fit))
     FrameE2Low=mean(TimeE2NumLow(Fit))/mean(TimeEDenLow(Fit))
  endif
  FrameELow=sqrt(FrameE1Low^2+FrameE2Low^2)
  FrameThetaLow=.5*atan(FrameE2Low,FrameE1Low)

  EllipVals(0,kval)=FrameEHigh
  EllipVals(1,kval)=FrameThetaHigh
  EllipVals(2,kval)=FrameELow
  EllipVals(3,kval)=FrameThetaLow 
endfor

return,EllipVals

end