;analyse_kphasescreen
;
;PURPOSE:
;       This script will take in phase screens generated by turb2d.c (written
;       by Garret Jernigan) and simulate N layers with 
;       structure constants CN2(n) drifting across a telescope aperture.
;       
;       Assuming geometrical optics, the ellipticity of a star that would be
;       formed by the wavefront is calculated.  High and low pass filters are
;       applied and the ellipticity of these filtered wavefronts is calculated.
;       In addition, if the Covariance keyword is set, the covariance of each
;        wavefront will be calculated.
;
;CALLING SEQUENCE:
;       phase_screen_ellipticity_filter,[date,[FrameDrift,[depth]]]...
;
;KEYWORDS:
;       Windspeed: 
;           The speed at which the frozen screen will be dragged across 
;                  the screen.  The units 
;       PhasDim: 
;           The dimension of the phase screen (default is 1024)
;       AvgSub: 
;            0 = don't subtract any offset
;            1 = subtract the average offset of a given frame
;            2 = set the kx=0,ky=0 componenet of FFT to zero
;            3 = subtract a scaled offset from the distorted grid
;       FittedOnly:
;            1 = use only the fitted points on each frame
;            2 = use only the non-fitted points one each frame
;            3 = use only points that have non-zero value
;       UseFid:
;            1 = apply the fiducial cut
;            0 = don't apply the fiducial cut
;       Covariance:
;            0 = don't do covariance analysis
;            1 = do covariance analaysis
;       WindSpeed:
;            1 = Use multiple wind speeds
;            2 = Use only one wind speed
;	DifferenceMeth
;	     1 = Finite Centered Difference
;	     2 = Finite Forward Difference
;	     3 = Finite Backward Difference

pro analyse_kphasescreens,$
                   xWindSpeed=xWindSpeed,yWindSpeed=yWindSpeed,$
                   PhaseDim=PhaseDim,$
                   Lo=Lo,lcut=lcut,$
                   MultiWindSpeeds=MultiWindSpeeds, AvgSub=AvgSub,$
                   UseFid=UseFid,Covariance=Covariance, TvFlag=TvFlag,$
                   OldWay=OldWay,Seed=Seed,$
		   TotalPhaseScreens=TotalPhaseScreens,$
		   ApDim=ApDim,DifferenceMeth=DifferenceMeth,$
		   ApNum=ApNum
If (not keyword_set(AvgSub)) then AvgSub =1
If (not keyword_set(FittedOnly)) then FittedOnly = 1
If (not keyword_set(UseFid)) then UseFid =0
If (not keyword_set(Covariance)) then Covariance = 1
If (not keyword_set(PhaseDim)) then PhaseDim=long(1024)
If (not keyword_set(ApDim)) then ApDim=PhaseDim-2
If (not keyword_set(Lo)) then Lo=1000000000.0
If (not keyword_set(DeltaX)) then DeltaX = 1.
If (not keyword_set(Ro)) then Ro=0.1
If (not keyword_set(MultiWindSpeeds)) then MultiWindSpeeds=2
If (not keyword_set(ExpTime)) then ExpTime=.03
If (not keyword_set(TvFlag)) then TvFlag = 0
If (not keyword_set(Spacing)) then Spacing = 1
If (not keyword_set(Seed)) then Seed= 20
If (not keyword_set(TotalPhaseScreens)) then TotalPhaseScreens=200
If (not keyword_set(DifferenceMeth)) then DifferenceMeth=1
If (not keyword_set(ApNum)) then ApNum = 1
;***************************************************************************
Kmax=PhaseDim/2
Border=(PhaseDim-ApDim)/2
SubSec=30
PhasePowerArr=fltarr(PhaseDim,PhaseDim,TotalPhaseScreens)
PhaseApPowerArr=Fltarr(ApDim,ApDim,TotalPhaseScreens)
XSlopesPowerArr=Fltarr(ApDim,ApDim,TotalPhaseScreens)
YSlopesPowerArr=Fltarr(ApDim,ApDim,TotalPhaseScreens)

for GlobalScreenNum=0, TotalPhaseScreens-1 do begin      ;Final Phase at pupil
  Seed=Seed+1
  PhaseScreen=return_kolmogorovscreen(PhaseDim, DeltaX, Lo, Ro, Seed)
  PhaseFFT=complexarr(PhaseDim,PhaseDim)
  PhaseFFT(*,*)=FFT(PhaseScreen(*,*))
  PhasePowerArr(*,*,GlobalScreenNum)=abs(shift(PhaseFFT,KMax,KMax))^2
  XSlopes=fltarr(ApDim,ApDim)
  YSlopes=fltarr(ApDim,ApDim)
  PhaseScreenAp=fltarr(ApDim,ApDim)
  Slopes=return_slopes(PhaseScreen,ApDim,DifferenceMeth, ApNum)
  XSlopes(*,*)=Slopes(*,*,0)
  YSlopes(*,*)=Slopes(*,*,1)
  PhaseScreenAp(*,*)=PhaseScreen(Border:ApDim+Border-1,Border:ApDim+Border-1)
  PowerSpectra=return_entirepowerspectrum(PhaseScreenAp,XSlopes,YSlopes,$
                                          UseFid,ApDIm)
  if TvFlag eq 1 then begin
    PhaseScreenTv=indgen(PhaseDim,PhaseDim,3)
    PhaseScreenTvSub=indgen(SubSec,SubSec,3)
    loadct,0
    tvlct,R,G,B,/get
    MaxScreen=max(PhaseScreen)
    MinScreen=min(PhaseScreen)
    ScreenVal=floor(255*((PhaseScreen-MinScreen)/(MaxScreen-MinScreen)))
    PhaseScreenTV(*,*,0)=R(ScreenVal)
    PhaseScreenTV(*,*,1)=G(ScreenVal)
    PhaseScreenTV(*,*,2)=B(ScreenVal)
    for Segment=0,200 do begin
      PhaseScreenTvSub[*,*,0]=PhaseScreenTV[Segment:Segment+SubSec-1,$
         Segment:Segment+SubSec-1,0]
      PhaseScreenTvSub[*,*,1]=PhaseScreenTV[Segment:Segment+SubSec-1,$
         Segment:Segment+SubSec-1,1]
      PhaseScreenTvSub[*,*,2]=PhaseScreenTV[Segment:Segment+SubSec-1,$
         Segment:Segment+SubSec-1,2]
    wait,.001
    tv,congrid(PhaseScreenTvSub,768,768,3),true=3
   endfor
  endif
   
  PhaseApPowerArr(*,*,GlobalScreenNum)=PowerSpectra(*,*,0)
  XSlopesPowerArr(*,*,GlobalScreenNum)=PowerSpectra(*,*,1)
  YSlopesPowerArr(*,*,GlobalScreenNum)=PowerSpectra(*,*,2)
endfor

PlotStyle=0
TotalPhasePower=avg(PhasePowerArr,2)
TotalApPhasePower=avg(PhaseApPowerArr,2)
TotalXPower=avg(XSlopesPowerArr,2)
TotalYPower=avg(YSlopesPowerArr,2)

PlotF=plot_powerspectrum(TotalPhasePower,TotalApPhasePower,$
			 TotalXPower,TotalYPower,$
                         PlotStyle,DeltaX,PhaseDim,ApDim,Ro)

stop
end