;phase_screen_power_spectrum
;
;PURPOSE:
;       This script will take in phase screens generated by turb2d.c (written
;       by Garret Jernigan) and compute their power spectrum
;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

pro phase_screen_power_spectrum,$
                   xWindSpeed=xWindSpeed,yWindSpeed=yWindSpeed,$
                   PhaseDim=PhaseDim,$
                   low=low,lcut=lcut,$
                   MultiWindSpeeds=MultiWindSpeeds, AvgSub=AvgSub,$
                   UseFid=UseFid,Covariance=Covariance, TvFlag=TvFlag

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(low)) then low=1000.0
If (not keyword_set(lcut)) then lcut=100000.0
If (not keyword_set(MultiWindSpeeds)) then MultiWindSpeeds=2
If (not keyword_set(x_off)) then x_off=0
If (not keyword_set(y_off)) then y_off=0
If (not keyword_set(ExpTime)) then ExpTime=.03
If (not keyword_set(TvFlag)) then TvFlag = 0
;***************************************************************************
;Calculate The Inner and Outer scales of the turbulence
lowstr=strtrim(string(low,format='(%"%10.1f")'),2)
lcutstr=strtrim(string(lcut,format='(%"%10.1f")'),2)
LowInMeters=PhaseDim*.17/low
HighInMeters=PhaseDim*.17/lcut

;FIDUCIAL VALUES
nonfiducial=return_fiducial()
complement,findgen(626),nonfiducial,fiducial
;*************************************************************************
;FILES
AtmosphereDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/'
PostscriptDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/phase_screen_figs/'
ScreenDirectory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/'

PhaseFiles=file_search(ScreenDirectory+'phase_'+strtrim(lowstr,2)+'_*.out',$
                       count=NumPhaseFiles)
oldway=0
If oldway eq 1 then begin
low=1.7
lcut=100000.0
lowstr=strtrim(string(low,format='(%"%10.1f")'),2)
lcutstr=strtrim(string(lcut,format='(%"%10.1f")'),2)
PhaseFiles=file_search(AtmosphereDir+'phase_'+strtrim(lowstr,2)+'_'+$
           strtrim(lcutstr,2)+'*.out',$
             count=NumPhaseFiles)
endif
KMax=PhaseDim/2
Num=PhaseDim
xwinds=[0.]
ywinds=[0.]
CN2s=[1]
NumLayers=N_elements(CN2s)
TotalPhaseScreens=NumPhaseFiles/NumLayers
PhiPowerArr=Fltarr(PhaseDim,PhaseDim,TotalPhaseScreens)

;***************************************************************************
for GlobalScreenNum=0, TotalPhaseScreens-1 do begin    ;Final Phase at pupil
  Screens=fltarr(PhaseDim,PhaseDim,NumLayers)	       ;Array to Hold Screens
  for LayerScreenNum=0, NumLayers -1 do begin          ;Ind. Phase Layers
    ScreenNum=GlobalScreenNum*NumLayers+LayerScreenNum
    ;Open the phase screen binary file and read floating point (data_type=4)
    openr,Phase,PhaseFiles(ScreenNum),/get_lun
    PhaseScreenSoar=read_binary(phase,$
                           data_dims=[PhaseDim,PhaseDim],data_type=4)
    close,phase
    free_lun,phase    
  endfor
  PhiPower=return_entirepowerspectrum(PhaseScreenSoar,FittedOnly,Num)
  PhiPowerArr(*,*,GlobalScreenNum)=PhiPower
  print,"Screen Number", GlobalScreenNum

endfor

TotalPhasePower=avg(PhiPowerArr,2)
stop

end