;function returnmullayphasescreen

;purpose:
;	This function will take in N phase screens along with N wind speeds
;	and N CN2 values, each of which correspond to a thin turbulent layer in
;	the atmosphere.  Assuming the frozen screen hypothesis and ignoring 
;	diffraction effects, it will calculate the resulting screen seen at 
;	the telescope.
;
;
;
;	N    = Number of layers
;	Mmax = Number of Samples from each layer

function return_mullayphasescreen_image_ellip, $
	  Screens, XWinds, YWinds, CN2s, ExpTime,$
	  Nmax, PhaseDim, ApNum, ApDim, DifferenceMeth, $
	  TvFlag

  ;PARAMETERS
  PixSize=0.17/ApNum                    ;Size of a Phase Screen cell in meters
  Diameter=0.17				;Size of the WFS lenselet in meters
  XWindMax=Max(abs(XWinds))	        ;Maximum X Wind Speed in N layers
  YWindMax=Max(abs(YWinds))	 	;Maximum Y Wind Speed in N layers
  WindMax=YWindMax		
  If XWindMax gt YWindMax then WindMax=XWindMax
  Mmax=floor(WindMax*ExpTime/PixSize+1)    ;# times screens will be sampled
  SoarPixSize=ApDim*ApNum+2
  if tvflag eq 3 then SoarPixSize=PhaseDim-40
  SoarLayerS=fltarr(SoarPixSize,SoarPixSize,Nmax)	;Array for N index
  SoarTotalS=fltarr(SoarPixSize,SoarPixSize,Mmax)	;Array for M index
  Center=floor((PhaseDim-SoarPixSize)/2)

  ;ARRAYS FOR THE SLOPES AND SECOND MOMENTS
  XSlopes = Fltarr(ApDim*ApNum, ApDim*ApNum, Mmax)
  YSlopes = Fltarr(ApDim*ApNum, ApDim*ApNum, Mmax)
  IValues = Fltarr(3)					;Ixx, Iyy, Ixy

  ;***********************				******************
  for m=0, Mmax-1 do begin				;Loop over Time

    for n=0, Nmax-1 do begin				;Loop over Layers
      XWindOff=round((Xwinds(N)*ExpTime/PixSize)*(float(m)/Mmax))
      YWindOff=round((Ywinds(N)*ExpTime/PixSize)*(float(m)/Mmax))
      SoarLayers(*,*,n)=CN2s(n)*$
		        Screens(Center-XWindOff:Center-XWindOff+SoarPixSize-1,$
       			        Center-YWindOff:Center-YWindOff+SoarPixSize-1,$
			        n)
      if TvFlag eq 1 then begin
	 set_plot,'x'
         window,xsize=600,ysize=600
	 tvscl,congrid(SoarLayers(*,*,n),300,300),0.5*(n/2),0.5*(n mod 2),/norm
	 stop
      endif  
    endfor

    ;Now we compute all of the Squares of the Angles of Arrival
    If Nmax eq 1 then begin
	SoarTotalS(*,*,m)=SoarLayers(*,*,0)/Total(CN2s)
    endif else begin
        SoarTotalS(*,*,m)=Total(SoarLayers,3)/Total(CN2s)
    endelse

    ;Get the Slopes
    Slopes  = Return_Slopes(SoarTotalS(*,*,m), ApDim, DifferenceMeth, ApNum)
    XSlopes(*,*,m) = Slopes(*,*,0)
    YSlopes(*,*,m) = Slopes(*,*,1)

    if TvFlag eq 3 then begin
      Height=Nmax*100
      Phase3D=fltarr(SoarPixSize,SoarPixSize,Height)
      Phase3D(*,*,0)=SoarLayers(*,*,0)
      Phase3D(*,*,1*float(Height)/4)=SoarLayers(*,*,1)
      Phase3D(*,*,2*float(Height)/4)=SoarLayers(*,*,2)
      Phase3D(*,*,3*float(Height)/4)=SoarLayers(*,*,3)
      scale3,xrange=[0,300],yrange=[0,300],zrange=[0,300]
      shade_volume,Phase3D,5e-8,a,b,/low
      tvscl,polyshade(a,b,/t3d)
    endif
  endfor

  If Mmax eq 1 then begin
    SoarFinalScreen=SoarTotals
  endif else begin
    SoarFinalScreen=avg(SoarTotals,2)
  endelse 
  if TvFlag eq 2 then begin
     set_plot,'x'
     window,xsize=600,ysize=600
     tvscl,congrid(SoarFinalScreen,600,600)
     stop
  endif

  ;Calculate and put in array
  IValues(0) = total(XSlopes^2)
  IValues(1) = total(YSlopes^2)
  IValues(2) = total(2*XSlopes*YSlopes)

  return, IValues

end