;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_image,$
		   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']
SOARPowerFileNames=FinalStatsDir+'PowerEllip'+Date+'.txt'
SOAREHistFileNames=FinalStatsDir+'EHistogram'+Date+'.txt'

;==========================================================================
;*************************************************************************
;***********************************FILEEXISTS=0  CALCULATE VALUES
If ComputeAnyway eq 1 then FileExists=0
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)
 PhaseScreens=dblarr(ApDim,ApDim,TotalPhaseScreens)
 IValArr = dblarr(3, TotalPhaseScreens)
 ;=================================================================================
 ;*********************************************************************************
 ;Loop For Each Segment of the Night	NumFile(FNo) Screens -- NumFile(FNo)*NumLayers
 Seed=MySeed
 GlobalScreenNum=0
 for FNo = 0, 4 do begin; 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, 2 do begin;NumFiles(FNo)-1 do begin       ;Final Phase at pupil
   Screens=fltarr(PhaseDim,PhaseDim,NumLayers)

   ;Generate NumLayers Phase Screens
   for LayerScreenNum=0, NumLayers-1 do begin	       ;Ind. Phase Layers
    Screens(*,*,LayerScreenNum)=return_kolmogorovscreen(PhaseDim,Deltax,$
			          Lo,Ro,Seed)
   endfor
  
   ;Loop through the time increment and layers and get one value for Ixx, Iyy, Ixy
   IValArr(*,GlobalScreenNum) = return_mullayphasescreen_image_ellip(Screens, $
		XWinds, YWinds, CN2s, $
 		ExpTime, NumLayers, PhaseDim, ApNum, ApDim, DifferenceMeth, TvFlag) 
   
   GlobalScreenNum=GlobalScreenNum+1
   print, GlobalScreenNum
  endfor
  print,FNo
 endfor

 ;Form Ellipticity Array to Pass on to the Hist_Ellipticities function
 EllipticityArr = Fltarr(4, TotalPhaseScreens)					   
 EllipticityArr(0,*)  = (IValArr(0,*) - IValArr(1,*))/(IValArr(0,*) + IValArr(1,*))
 EllipticityArr(1,*)  = (IValArr(2,*))/(IValArr(0,*) + IValArr(1,*))
 EllipticityArr(2,*)  = sqrt(EllipticityArr(0,*)^2+EllipticityArr(1,*)^2)
 EllipticityArr(3,*)  = .5*atan(E2, E1) 
 stop

 ;=================================================================================
 ;********************WRITE VALUES TO FILE
 
 HistEllip = Hist_Ellipticities(EllipticityArr)
 EHistFileName=AtmosphereDir+ParamStr+'EHistorgramImage.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'
 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)
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