;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_1d_kolmogorov_region, PhaseDim=PhaseDim, DeltaX=DeltaX, Lo=Lo, $
				     ro=ro, Seed=Seed,ApDim=ApDim

If (not keyword_set(PhaseDim)) then PhaseDim=long(1024)
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=1000
BorderSize=(PhaseDim-ApDim)/2
!p.multi=[0,2,3]
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))
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)
 Spectrum=PhaseFFT/(KSqrd+(TotalLength/Lo)^2)^(11./12)
 SpectrumRe=shift(Spectrum,KMax)
 SpectrumRe(KMax+1:PhaseDim-1)=conj(reverse(SpectrumRe(1:KMax-1)))
 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(BorderSize:PhaseDim-BorderSize-1)
 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(BorderSize:PhaseDim-BorderSize-1)=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

endfor

AvgPowerSpectrumFull=avg(PowerSpectraFullArr,0)
AvgPowerSpectrumSamp=avg(PowerSpectraSampArr,0)
AvgPowerSpectrumSampZeros=avg(PowerSpectraSampZerosArr,0)
AvgWindowedSpectrum=avg(WindowedSpectraArr,0)
AvgSlowSpectrum=avg(SlowSpectraArr,0)

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)
        
stop

;Here is the Signal
x=DeltaX*findgen(SignalN)
y1=Strength1*sin((2*!pi/Period1)*x)+Strength2*cos((2*!pi/Period2)*x)
plot,x,y1

FreqSpacing1=1/(SignalLength)
Kvals1=FreqSpacing1*(-KMax1+findgen(SignalN))
FftY1=fft(y1)
plot,KVals1,abs(shift(ffty1,KMax1))^2


;Now Sample the function and see what we get Sample=y2
SampleSignalN=N_elements(y1)/2
SampleSignalLength=SampleSignalN*DeltaX
FreqSpacing2=1/(SampleSignalLength)
KMax2=SampleSignalN/2
Kvals2=FreqSpacing2*(-KMax2+findgen(SampleSignalN))

;Here is the sampled signal
x=DeltaX*findgen(SampleSignalN)
y2=y1(0:SampleSignalN)
plot,x,y2

FftY2=fft(Y2)
plot,KVals2,abs(shift(ffty2,KMax2))^2

end