;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, Screens, XWinds, YWinds, CN2s, ExpTime,$
	  Nmax, PhaseDim, ApNum, ApDim, TvFlag
  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
  ;  SoarPixSize=ApDim*Diameter/PixSize+2*Diameter/PixSize
  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)
  ;***********************				******************
  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
    If Nmax eq 1 then begin
	SoarTotalS(*,*,m)=SoarLayers(*,*,0)/Total(CN2s)
    endif else begin
        SoarTotalS(*,*,m)=Total(SoarLayers,3)/Total(CN2s)
    endelse
    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
  return,SoarFinalScreen

end