;Generate_PhaseScreen
;
;PURPOSE:
;	This code will generate an NDimxNDim two-dimensional phase screen with 
;	Kolmogorov properties.
;INDimPUTS:
;  NDim:
;	The number of pixels along a dimension.  Total screen is NDimxNDim
;  DeltaX:
;	The number of meters per pixel
;  Lo:
;	The length of the outer scale in meters
;  ro:
;	The Fried Parameter.  Determines the strength of turbulence, i.e. 
;	the rms of the phase about zero mean.
;	
function return_kolmogorovscreen, NDim, DeltaX, Lo, ro, Seed

Low=Double(Lo)	;For some reason, IDL treats function arguments as global
NDim=long(NDim)
TotalLength=NDim*DeltaX	;Total Length of Screen in meters
KMax=floor(float(NDim)/2+.5)
FreqSpacing=1/((NDim)*DeltaX)

;Generate KValues for the phase screen
Ks=FreqSpacing*(-NDim/2+findgen(NDim))
KInd=dindgen(NDim*NDim)
KArr=dblarr(NDim,NDim)
KArr(Kind mod NDim,Kind/NDim)=Ks(Kind mod NDim)^2+Ks(Kind/NDim)^2

;Phase Screen is Scaled by Ro Fried Paramater
ScreenStr=0.1517*(TotalLength/ro)^(5./6)*.5e-6/(2*!pi)
PhaseFFT=Complexarr(NDim,NDim)
PhaseFFT(*,*)=randomn(Seed,NDim,NDim)+complex(0,1)*randomn(Seed,NDim,NDim)
PhaseFFT(*,*)=ScreenStr*PhaseFFT(*,*)

;Eliminate the Alternating Sum if N is Even by setting these to zero
If NDim mod 2 eq 0 then begin
 PhaseFFT(0,*)=0
 PhaseFFT(*,0)=0
endif

;Scale Random Numbers by k^(-11./3)
Spectrum=PhaseFFT/(KArr+(TotalLength/Low)^2)^(11./12)

;INDDEX REARRANDGINDG
KRe=shift(Karr,KMax,KMax)
SpectrumRe=shift(Spectrum,KMax,KMax)
SpectrumRe(0,0)=0	;Elminate average offset
Phase=fft(SpectrumRe,/inverse)
return,real_part(shift(Phase,KMax,KMax))
end