;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)
;	ApNum:
;	    The number of pixels to average over in order to get the wavefront
;	    slopes
;	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
;	PlotStyle:
;	     1 = Use an X Window
;	     2 = Plot each set of data separately
;	     3 = Plot the figures for the paper

pro phase_screen_multi_layer,$
		   xWindSpeed=xWindSpeed,yWindSpeed=yWindSpeed,$
	           PhaseDim=PhaseDim,$
		   Lo=Lo,lcut=lcut,$
		   MultiWindSpeeds=MultiWindSpeeds, AvgSub=AvgSub,$
	           UseFid=UseFid,Covariance=Covariance, TvFlag=TvFlag,$
		   OldWay=OldWay,MySeed=MySeed,ApDim=ApDim,$
		   PlotStyle=PlotStyle,ComputeAnyway=ComputeAnyway,$
		   DifferenceMeth=DifferenceMeth,$
		   ApNum=ApNum
If (not keyword_set(AvgSub)) then AvgSub =1 
If (not keyword_set(FittedOnly)) then FittedOnly = 2
If (not keyword_set(UseFid)) then UseFid =2
If (not keyword_set(ApNum)) then ApNum = 6
If (not keyword_set(Covariance)) then Covariance = 1
If (not keyword_set(PhaseDim)) then PhaseDim=long(250)
If (not keyword_set(ApDim)) then ApDim=25
If (not keyword_set(Lo)) then Lo=100.
If (not keyword_set(DeltaX)) then DeltaX = 0.17/ApNum
If (not keyword_set(Ro)) then Ro=0.1
If (not keyword_set(MultiWindSpeeds)) then MultiWindSpeeds=1
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(MySeed)) then MySeed= 20
If (not keyword_set(PlotStyle)) then PlotStyle=0
If (not keyword_set(ComputeAnyway)) then ComputeAnyway=0
If (not keyword_set(DifferenceMeth)) then DifferenceMeth=1

;***************************************************************************
;FIDUCIAL VALUES
NonFiducial=return_fiducial()
complement,findgen(626),NonFiducial,Fiducial
Border=1					;Only Need 1 unit for Fi Diff
;*************************************************************************
;FILES
;********************INPUT
AtmosphereDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/'
PostscriptDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/final_paper_figs/'
ScreenDirectory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/'
FinalStatsDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/FinalNumbers/'
FileExists=0
;********************OUTPUT

If DifferenceMeth eq 1 then begin
   DifferenceMethSt="FCD"
endif else if DifferenceMeth eq 2 then begin
   DifferenceMethSt="FFD"
endif else if DifferenceMeth eq 3 then begin
   DifferenceMethSt="FBD"
endif else if DifferenceMeth eq 4 then begin
   DifferenceMethSt="ComPade"
endif
If UseFid eq 1 then begin
   FidStr="_Fid_"
endif else if UseFid eq 0 then begin
   FidStr="_NoFid_"
endif else begin
   FidStr='_'
endelse

ParamStr='MultiLayerPhase'+FidStr+DifferenceMethSt+"_"+'Avg'+strtrim(ApNum,2)+$
	'Aps_'+strtrim(MySeed,2)+'Seed_'+strtrim(Lo,2)+'Low_'+$
	strtrim(PhaseDim,2)
FileName=AtmosphereDir+ParamStr+'DimCovar.txt'
If File_Test(FileName) then begin 
   FileExists=1
   print,"File: "+FileName+"Already exists...."
   print,"Taking Data from File"
endif
If PlotStyle eq 3 then begin
  PSFileName=PostscriptDir+ParamStr+'.ps'
endif else if PlotStyle eq 2 then begin
  PSFileName=PostscriptDir+ParamStr+'IndPlots.ps'
endif

Date=['10','11','12']
SOARCovarFileNames=FinalStatsDir+'Covariances'+Date+'.txt'
SOARPowerFileNames=FinalStatsDir+'PowerEllip'+Date+'.txt'
SOAREHistFileNames=FinalStatsDir+'EHistogram'+Date+'.txt'
;==========================================================================
;*************************************************************************
;***********************************FILEEXISTS=0  CALCULATE VALUES
If ComputeAnyway eq 1 then FileExists=0
FileExists = 1
If FileExists eq 0 then begin
 if MultiWindSpeeds eq 1 then begin
   readcol,'/nfs/slac/g/ki/ki08/lsst/CPanalysis/winds/PhaseScreenWinds.txt',$
	   NumFiles,Wx1,Wy1,Wgt1,Wx2,Wy2,Wgt2,Wx3,Wy3,Wgt3,Wx4,Wy4,Wgt4,$
	   format='d,f,f,f,f,f,f,f,f,f,f,f,f'
 endif
 ;*******************************************************VALUES FOR LOOP
 KMax=12
 ApDim=25
 NumLayers=4
 TotalPhaseScreens=Total(NumFiles)
 EllipticityArr=fltarr(8,KMax,TotalPhaseScreens)
 PhaseScreens=dblarr(ApDim,ApDim,TotalPhaseScreens)
 CovarianceArr=fltarr(4,25,TotalPhaseScreens)
 PhasePowerArr=Fltarr(ApDim,ApDim,TotalPhaseScreens)
 XSlopesPowerArr=Fltarr(ApDim,ApDim,TotalPhaseScreens)
 YSlopesPowerArr=Fltarr(ApDim,ApDim,TotalPhaseScreens)
 ;***************************************************************************
 ;Loop For Each Segment of the Night	NumFile(FNo) Screens -- NumFile(FNo)*NumLayers
 Seed=MySeed
 GlobalScreenNum=0
 for FNo = 0, N_elements(NumFiles)-1 do begin
  if NumLayers eq 4 then begin
    XWinds=[Wx1(FNo),Wx2(FNo),Wx3(FNo),Wx4(FNo)]
    YWinds=[Wy1(FNo),Wy2(FNo),Wy3(FNo),Wy4(FNo)]
    CN2s=[Wgt1(FNo),Wgt2(FNo),Wgt3(FNo),Wgt4(FNo)]
  endif else if NumLayers eq 1 then begin
    XWinds=[Wx4(FNo)]
    YWinds=[Wy4(FNo)]
    CN2s=[Wgt1(FNo)]
  endif
  for LocalScreenNum=0, NumFiles(FNo)-1 do begin       ;Final Phase at pupil
   Screens=fltarr(PhaseDim,PhaseDim,NumLayers)
   for LayerScreenNum=0, NumLayers-1 do begin	       ;Ind. Phase Layers
    Screens(*,*,LayerScreenNum)=return_kolmogorovscreen(PhaseDim,Deltax,$
			          Lo,Ro,Seed)
   endfor
   PhaseScreen=return_mullayphasescreen(Screens, XWinds, YWinds, CN2s, $
 		ExpTime, NumLayers, PhaseDim, ApNum, ApDim, TvFlag) 
   PhaseScreenAp=fltarr(ApDim,ApDim)
   XSlopes=fltarr(ApDim,ApDim)
   YSlopes=fltarr(ApDim,ApDim)
   Slopes=return_slopes(PhaseScreen, ApDim, DifferenceMeth, ApNum, /ShackHartman)
   XSlopes(*,*)=Slopes(*,*,0)
   YSlopes(*,*)=Slopes(*,*,1)

   if UseFid eq 1 then begin
    XSlopes(NonFiducial)=0
    YSlopes(NonFiducial)=0
   endif
   if AvgSub eq 1 then begin
    Xfit=where(XSlopes ne 0)
    Yfit=where(YSlopes ne 0)
    MeanXSlopes=mean(XSlopes(Xfit))
    MeanYSlopes=mean(YSlopes(Yfit))
    XSlopes(Xfit)=XSlopes(Xfit)-MeanXSlopes
    YSlopes(Yfit)=YSlopes(Yfit)-MeanYSlopes
   endif
   RMS=sqrt(mean(XSlopes^2+YSlopes^2))
   XSlopes=XSlopes/RMS
   YSlopes=YSlopes/RMS
   PhaseScreenAp(*,*)=PhaseScreen(Border:ApDim+Border-1,Border:ApDim+Border-1)
   PowerSpectra=return_entirepowerspectrum(PhaseScreenAp,XSlopes,YSlopes,$
                                          FittedOnly,ApDIm)
   PhasePowerArr(*,*,GlobalScreenNum)=PowerSpectra(*,*,0)
   XSlopesPowerArr(*,*,GlobalScreenNum)=PowerSpectra(*,*,1)
   YSlopesPowerArr(*,*,GlobalScreenNum)=PowerSpectra(*,*,2)
   EllipVals=return_ScreenEllipticity(XSlopes,YSlopes,PhaseScreenAp,$
				      FittedOnly,ApDim,Fiducial)
   EllipticityArr(*,*,GlobalScreenNum)=EllipVals
   CovarianceArr(*,*,GlobalScreenNum)=return_covariance(XSlopes,YSlopes)
   GlobalScreenNum=GlobalScreenNum+1
  endfor
  print,FNo
 endfor
 ;OBTAIN AVERAGE VALUES 
 TotalPhasePower=avg(PhasePowerArr,2)
 TotalXPower=avg(XSlopesPowerArr,2)
 TotalYPower=avg(YSlopesPowerArr,2)
 AvgCovarianceArr=avg(CovarianceArr,2)
 Cxl=AvgCovarianceArr(0,*)	;X Longitudinal
 Cxt=AvgCovarianceArr(1,*)	;X Transverse
 Cyl=AvgCovarianceArr(2,*)	;Y Longitudinal
 Cyt=AvgCovarianceArr(3,*)	;Y Transverse
 AvgEllipticityArr=avg(EllipticityArr,2)
 FrameEHigh=AvgEllipticityArr(0,*)
 FrameThetaHigh=AvgEllipticityArr(1,*)
 FrameELow=AvgEllipticityArr(2,*)
 FrameThetaLow=AvgEllipticityArr(3,*)
 XPower=AvgEllipticityArr(4,*)
 YPower=AvgEllipticityArr(5,*)
 XPhasePower=AvgEllipticityArr(6,*)
 YPhasePower=AvgEllipticityArr(7,*)
 ;******************************************FILE FOR COVARIANCES
 ScreenCovarFileName=AtmosphereDir+ParamStr+'DimCovar.txt'   
 get_lun,Quantities
 openw,Quantities,ScreenCovarFileName
 printf,Quantities,format='(%"%s\t%s\t%s\t%s")',$
         'XLong','XTrans','YLong','YTrans'
 for i=0, N_elements(Cxl)-1 do begin
   printf,Quantities,$
    format='(%"%f\t%f\t%f\t%f")',$
    Cxl(i),Cxt(i),Cyl(i),Cyt(i)
 endfor
 close,Quantities
 free_lun,Quantities
 ;**********************************************FILE FOR ELLIPTICITIES &POWER
 PowerFileName=AtmosphereDir+ParamStr+'DimPower.txt'
 get_lun,Quantities
 openw,Quantities,PowerFileName
 printf,Quantities,format='(%"%s\t%s\t%s\t%s\t%s\t%s")',$
         'XPower','YPower','XPhaseP',$
         'YPhaseP','ELowPass','EHighPass'
 for i=0, N_elements(XPower)-1 do begin
   printf,Quantities,$
    format='(%"%11.4e\t%11.4e\t%11.4e\t%11.4e\t%f\t%f")',$
    XPower(i),YPower(i),XPhasePower(i),$
    YPhasePower(i),FrameELow(i),FrameEHigh(i)
 endfor
 close,Quantities
 free_lun,Quantities
 ;**********************************************FILE FOR ELLIPTICITY HISTOGRAM
 HistEllip=Hist_Ellipticities(EllipticityArr)
 EHistFileName=AtmosphereDir+ParamStr+'DimEHistorgram.txt'
 get_lun,E
 openw,E,EHistFileName
 printf,E,format='(%"%s\t%s\t%s\t%s")',$
     'E Locs','E Hist','Theta Locs', 'ThetaHist'
 for i=0, N_elements(HistEllip(0,*))-1 do begin
    printf,E,$
     format='(%"%f\t%f\t%f\t%f")',$
     HistEllip(0,i),HistEllip(1,i),HistEllip(2,i),HistEllip(3,i)
 endfor
 close,E
 free_lun,E
 SOARDataFileName='/nfs/slac/g/ki/ki08/lsst/CPanalysis/Ellip_Wind_Pointing.txt'

endif else begin

;===========================================================================
;***************************************************************************
;**********************************************READ FROM FILE FILEEXISTS=1
 PowerFileName=AtmosphereDir+ParamStr+'DimPower.txt'
 ScreenCovarFileName=AtmosphereDir+ParamStr+'DimCovar.txt'
 EHistFileName=AtmosphereDir+ParamStr+'DimEHistorgram.txt'
 SOARDataFileName='/nfs/slac/g/ki/ki08/lsst/CPanalysis/Ellip_Wind_Pointing.txt'
endelse

;==========================================================================
;**********************************************PLOTTING
If PlotStyle eq 1 then begin
   Set_plot,'z'
   device,set_resolution=[1400,1200]
   !p.charsize=3
   !p.charthick=5
   !x.thick=4
   !y.thick=4
   !p.thick=3.5
   !p.noerase=0

   loadct, 39
   white='FFFFFF'x
   black='000000'x
   red='FF0000'x
   !P.multi=[0,1,1]

endif else if PlotStyle eq 2 or PlotStyle eq 3 then begin
   Set_plot,'ps'
   device,filename=PSFileName
   loadct,39                       ;load color table 39
   device,/color                   ;allow color on the postscript
   device,ysize=8.5,/inches        ;Height of plot in y
   device,xsize=12.0,/inches        ;Width of plot in x
   device,yoffset=0.0,/inches      ;Y position of lower left corner
   device,xoffset=0.0,/inches

   !p.charsize=1
   !p.thick=4.5
   !p.charthick=4
   !x.thick=4
   !y.thick=4
   white='FFFFFF'x
   black='000000'x
   red='FF0000'x
   green='00FF00'x
   blue='0000FF'x

endif else begin
   set_plot,'X'
endelse

PlotE=plot_elliphistogram_ff(EHistFileName,SOARDataFileName,PlotStyle,ParamStr,$
			PostScriptDir, /IncludeSoarData)
PlotC=plot_covariance_ff(ScreenCovarFileName,SOARCovarFileNames,PlotStyle,$
			ParamStr, PostScriptDir)
PlotF=plot_powerspectrum_ff(PowerFileName,SOARPowerFileNames,PlotStyle,$
		            DeltaX,ApDim,ApDim,ApNum,Ro, ParamStr, PostScriptDir)

if PlotStyle eq 2 or PlotStyle eq 3 then begin
   device,/close
endif
set_plot,'x'

stop
 
end