;PURPOSE:
;	The purpose of this script is to show the effect of sampling a function
;	on its fft
;
;	3 Power Spectra Are calculated
;	------------------------------
;	1-Original -11/3 power spectrum generated with random numbers scaled
;	  by power law with length=PhaseDim
;	2-Sampled region of inverted power spectrum with length=ApDim
;	3-Sampled region of inverted power spectrum with length=PhaseDim
;	  and ApDim nonzero elements.
pro fft_sample_power_spectrum, PhaseDim=PhaseDim, DeltaX=DeltaX, Lo=Lo, $
				     ro=ro, Seed=Seed,ApDim=ApDim

If (not keyword_set(PhaseDim)) then PhaseDim=long(64)
If (not keyword_set(ApDim)) then ApDim=24
If (not keyword_set(Lo)) then Lo=100000000.
If (not keyword_set(DeltaX)) then DeltaX = 0.17
If (not keyword_set(Ro)) then Ro=0.1
If (not keyword_set(Seed)) then Seed= 10
If (not keyword_set(TotalPhaseScreens)) then TotalPhaseScreens=1
BorderSize=(PhaseDim-ApDim)/2
if (ApDim+PhaseDim) mod 2 eq 1 then begin
  BorderLow=BorderSize+1
  BorderHigh=PhaseDim-BorderSize-1
endif else begin
  BorderLow=BorderSize
  BorderHigh=PhaseDim-BorderSize-1
endelse

!p.multi=[0,2,4]
window,xsize=600,ysize=600
!p.charsize=2


Low=Double(Lo)  ;For some reason, IDL treats function arguments as global
PhaseDim=long(PhaseDim)
TotalLength=PhaseDim*DeltaX ;Total Length of Screen in meters
;KMax=floor(float(PhaseDim)/2+.5)
KMax=PhaseDim/2
FreqSpacing=1/((PhaseDim)*DeltaX)
FreqSpacingSamp=1/((ApDim*DeltaX))
;Generate KValues for the phase screen
Ks=FreqSpacing*(-PhaseDim/2+findgen(PhaseDim))
KsSamp=FreqSpacingSamp*(-ApDim/2+findgen(ApDim))
KsSlope=FreqSpacingSamp*(-(ApDim-2)/2+findgen(ApDim-2))
Ksqrd=Ks^2
ScreenStr=0.1517*(TotalLength/ro)^(5./6)*.5e-6/(2*!pi)
;And Frequencies
NumFrequencies=PhaseDim/2
SamplingFreq=1./DeltaX
NyquistFreq=SamplingFreq/2
Frequencies=(Findgen(NumFrequencies)/NumFrequencies)*NyquistFreq

PowerSpectraFullArr=fltarr(TotalPhaseScreens,PhaseDim)
PowerSpectraSampZerosArr=fltarr(TotalPhaseScreens,PhaseDim)
PowerSpectraSampArr=fltarr(TotalPhaseScreens,ApDim)
WindowedSpectraArr=fltarr(TotalPhaseScreens,ApDim)
SlowSpectraArr=fltarr(TotalPhaseScreens,NumFrequencies)
;*********************************************************************
for ScreenNum=0, TotalPhaseScreens-1 do begin
 seed=seed+2
 PhaseFFT=Complexarr(PhaseDim)
; PhaseFFT(*)=randomn(Seed,PhaseDim)+complex(0,1)*randomn(Seed,PhaseDim)
 PhaseFFT(*)=3
 Spectrum=PhaseFFT/(KSqrd+(TotalLength/Lo)^2)^(11./12)
 SpectrumRe=shift(Spectrum,KMax)
 if PhaseDim mod 2 eq 0 then begin
   SpectrumRe(KMax+1:PhaseDim-1)=conj(reverse(SpectrumRe(1:KMax-1)))
 endif else begin
   SpectrumRe(KMax+1:PhaseDim-1)=conj(reverse(SpectrumRe(1:KMax)))
 endelse
 SpectrumRe(0)=0       ;Elminate average offset
 SpectrumRe(KMax)=0
 Spectrum(KMax)=0
 PowerSpectraFullArr(ScreenNum,*)=abs(Spectrum)^2
 PhaseScreen=real_part(fft(SpectrumRe,/inverse))
 PhaseScreen=PhaseScreen-Mean(PhaseScreen)

 ;Sampled Phase Screen
 PhaseScreenSamp=PhaseScreen(BorderLow:BorderHigh)
 PhaseScreenSamp=PhaseScreenSamp-Mean(PhaseScreenSamp)
 PowerSpectrumSamp=abs(fft(PhaseScreenSamp))^2
 PowerSpectrumSamp(ApDim/2)=0
 PowerSpectraSampArr(ScreenNum,*)=abs(shift(PowerSpectrumSamp,ApDim/2))
 ;Sampled Phase Screen, Zero on Borders
 PhaseScreenSampZeros=fltarr(PhaseDim)
 PhaseScreenSampZeros(BorderLow:BorderHigh)=PhaseScreenSamp
 PowerSpectrumSampZeros=abs(fft(PhaseScreenSampZeros))^2
 PowerSpectrumSampZeros(PhaseDim/2)=0
 PowerSpectraSampZerosArr(ScreenNum,*)=$
			  abs(shift(PowerSpectrumSampZeros,PhaseDim/2))
 ;Sampled Phase Screen, Manning Window Function Applied
 WindowFunction=.5*(1-cos(2*!pi*findgen(ApDim)/(ApDim-1)))
 WindowedPhaseScreen=PhaseScreenSamp*WindowFunction
 WindowedPowerSpectrum=abs(fft(WindowedPhaseScreen))^2
 WindowedPowerSpectrum(ApDim/2)=0
 WindowedSpectraArr(ScreenNum,*)=abs(shift(WindowedPowerSpectrum,ApDim/2))

 ;Slow Fourier Transform
 SlowFFT=Return_slow_ft(PhaseScreenSamp,PhaseDim,DeltaX)
 SlowPowerSpectrum=(abs(SlowFFT))^2
 SlowPowerSpectrum(0)=0
 SlowSpectraArr(ScreenNum,*)=SlowPowerSpectrum

 ;Get the Slope of the phase
 Slopes=PhaseScreenSamp(2:ApDim-1)-PhaseScreenSamp(0:ApDim-3)
 SlopesFFT=shift(fft(Slopes),ApDim/2-1)
 SlopesPower=abs(SlopesFFT)^2
endfor

AvgPowerSpectrumFull=PowerSpectraFullArr
AvgPowerSpectrumSamp=PowerSpectraSampArr
AvgPowerSpectrumSampZeros=PowerSpectraSampZerosArr
AvgWindowedSpectrum=WindowedSpectraArr
AvgSlowSpectrum=SlowSpectraArr

plot,PhaseScreen
plot,PhaseScreenSamp

plot,ks(PhaseDim/2+1:PhaseDim-1),$	
	AvgPowerSpectrumFull(PhaseDim/2+1:PhaseDim-1),/ylog,/xlog
oplot,ks(PhaseDim/2+1:PhaseDim-1),ks(PhaseDim/2+1:PhaseDim-1)^(-11./3)

plot,ksSamp(ApDim/2+1:ApDim-1),$      
       AvgPowerSpectrumSamp(ApDim/2+1:ApDim-1),/ylog,/xlog
oplot,ksSamp(ApDim/2+1:ApDim-1),ksSamp(ApDim/2+1:ApDim-1)^(-11./3)

plot,KsSamp(ApDim/2+1:ApDim-1),WindowedPowerSpectrum(ApDim/2+1:ApDim-1),/ylog

plot,ks(PhaseDim/2+1:PhaseDim-1),$      
        AvgPowerSpectrumSampZeros(PhaseDim/2+1:PhaseDim-1),/ylog
oplot,ks(PhaseDim/2+1:PhaseDim-1),ks(PhaseDim/2+1:PhaseDim-1)^(-11./3)

plot,Frequencies,AvgSlowSpectrum,/ylog
oplot,Frequencies,Frequencies^(-11./3)

plot,ksSlope(ApDim/2:ApDim-3),SlopesPower(ApDim/2:ApDim-3),/xlog,/ylog
oplot,KsSlope(ApDim/2:ApDim-3),10*KsSlope(ApDim/2:ApDim-3)^(-5./3)
stop

end