;phase_screen_multi_layer
;
;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

pro phase_screen_multi_layer,$
		   xWindSpeed=xWindSpeed,yWindSpeed=yWindSpeed,$
	           PhaseDim=PhaseDim,$
		   low=low,lcut=lcut,$
		   MultiWindSpeeds=MultiWindSpeeds, AvgSub=AvgSub,$
	           UseFid=UseFid,Covariance=Covariance, TvFlag=TvFlag,$
		   OldWay=OldWay
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
If (not keyword_set(OldWay)) then OldWay = 0
If (not keyword_set(Spacing)) then Spacing = 1
;***************************************************************************
;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)
If oldway eq 1 then begin
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

EvalsFileName=ScreenDirectory+'ellipticity_AvgSub_low.txt'
get_lun, Evals
openw,Evals,EvalsFileName

;*************************************************************************
if MultiWindSpeeds eq 1 then begin
   xwinds=[-18.3 , -17.9 , -17.6, -17.3, -17.0, $
	   -17.7,  -17.3 , -16.9, -16.6, -16.3, $
	   -17.0,  -17.3 , -17.6, -17.9, -18.1]

   ywinds=[ 6.6 , 2.9, -0.8, -4.5, -8.2, $
	    1.7 ,-0.2, -2.0, -3.8, -5.6, $
	    5.6, 3.3,  1.0,  -1.3, -3.5]
endif else if MultiWindSpeeds eq 2 then begin
   xwinds=[0]
   ywinds=[0]
endif

;*******************************************************VALUES FOR LOOP
KMax=12
Num=25
xwinds=[0.,0,0]
ywinds=[0.,0,0]
CN2s=[1,2,1]
xwinds=[0.]
ywinds=[0.]
CN2s=[1]
NumLayers=N_elements(CN2s)
TotalPhaseScreens=NumPhaseFiles/NumLayers
EllipticityArr=fltarr(8,KMax,TotalPhaseScreens)
PhaseScreens=dblarr(num,num,TotalPhaseScreens)
CovarianceArr=fltarr(4,25,TotalPhaseScreens)
PhiPowerArr=Fltarr(Num,Num,TotalPhaseScreens)

;***************************************************************************
for GlobalScreenNum=0, TotalPhaseScreens-1 do begin      ;Final Phase at pupil
  Screens=fltarr(PhaseDim,PhaseDim,NumLayers)
  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
    Screens(*,*,LayerScreenNum)=read_binary(phase,$
			   data_dims=[PhaseDim,PhaseDim],data_type=4)
    close,phase 
    free_lun,phase

  endfor

  PhaseScreen=return_mullayphasescreen(Screens, XWinds, YWinds, CN2s, ExpTime,$
          NumLayers, PhaseDim, TvFlag) 
  PhaseScreenSoar=PhaseScreen(1:25,1:25)
  XSlopes=(PhaseScreen(2:26,1:25)-PhaseScreen(0:24,1:25))/2
  YSlopes=(PhaseScreen(1:25,2:26)-PhaseScreen(1:25,0:24))/2
  RMS=sqrt(mean(XSlopes^2+YSlopes^2))
  XSlopes=XSlopes/RMS
  YSlopes=YSlopes/RMS
  EllipVals=return_SoarEllipticity(XSlopes,YSlopes,PhaseScreenSoar,FittedOnly,$
		                   Num)
  EllipticityArr(*,*,GlobalScreenNum)=EllipVals
  PhiPower=return_entirepowerspectrum(PhaseScreenSoar,FittedOnly,Num)
  PhiPowerArr(*,*,GlobalScreenNum)=PhiPower
  print,EllipVals(0,0)

  CovarianceArr(*,*,GlobalScreenNum)=return_covariance(XSlopes,YSlopes)
endfor

PlotStyle=0
TotalPhasePower=avg(PhiPowerArr,2)
PlotE=plot_ellipticities(EllipticityArr,PlotStyle)
PlotC=plot_covariance(CovarianceArr,PlotStyle)
PlotF=plot_powerspectrum(EllipticityArr,PlotStyle,Spacing)

stop
 
end