Atmospheric Turbulence Screen
Abu Bucker Siddik
New Mexico State University
Presentation Objectives �
Random Phase Screen with Linear Filtering
From Numerical Construction of Random Phase Screen, D. Voelz
Random Phase Screen with Linear Filtering
From Numerical Construction of Random Phase Screen, D. Voelz
Random Phase Screen with Linear Filtering
From Numerical Construction of Random Phase Screen, D. Voelz
Random Phase Screen with Linear Filtering
From Numerical Construction of Random Phase Screen, D. Voelz
Atmospheric Turbulence Screen
Structure Function
Propagation of Plane Wave through Atmospheric Turbulence Screen
Matlab Code for Atmospheric Turbulence Screen�
N=500; %# samples 1D
L=1; %spatial grid side length
dx=L/N; %sample interval
x=linspace(-L/2, L/2-dx, N); %linear coordinates
[X, Y]= meshgrid(x,x); %2D cordinates
B=1/dx; %spatial freq full bandwidth (2-sided)
dv=B/N; %spatial freq interval
vx=linspace(-B/2, B/2-dv, N);%spatial freq coordinates
[Vx, Vy]= meshgrid(vx, vx); %2D freq coordinates
V=sqrt(Vx.^2 + Vy.^2); %radial symmetric freq coordinates
Matlab Code for Atmospheric Turbulence Screen�
var_n =1; std_n=sqrt(var_n); %noise variance & std
N02= var_n/B^2; %white noise PSD constant
Cn2 = 1e-15; %measure of strength of atmospheric turbulence
z=100; %thickness of atmospheric turbulence screen
lambda= 0.25e-6; %wavelength
k= (2*pi)/lambda; %wavenumber
r0= (0.423*(k^2)*Cn2 * z)^(-3/5);%fried diameter
% Filter TF for given Kolmogorov power spectrum (shifted)
H=fftshift(sqrt(0.023 * (r0^(-5/3)) * (V.^(-11/3)) *(1/N02))); H(1,1)=0;
Matlab Code for Atmospheric Turbulence Screen�
nr=std_n*randn(N); %white noise array
Nr=fft2(fftshift(nr))*dx^2;%white noise transform (shifted)
phi0= ifft2(H.* Nr)* N^2 * dv^2; %phase screen (shifted)
phi= fftshift(phi0);
figure
imagesc(x, x, phi);
axis square; axis xy; colorbar; xlabel('x'); ylabel('y');
Matlab Code for Atmospheric Turbulence Screen�
%loop to generate average phase screen PSD and autocorrelation
Nrelz= 100;
G=zeros(N);
for i=1:Nrelz
nr=std_n*randn(N);
Nr=fft2(fftshift(nr))*dx^2;
phi0= ifft2(H.* Nr)* N^2 * dv^2;
G0=(abs(fft2(phi0)*dx^2).^2)/L^2; %power spectrum
G= G+G0; %accumulate PSD
end
% power spectrum
G = G/Nrelz; %computational average PSD estimate (shifted)
Matlab Code for Atmospheric Turbulence Screen�
A=ifft2(G)*N^2*dv^2; %computational autocorrelation estimate (shifted)
D=2*A(1,1)-2*A(1,:); %structure function
D_A=6.88*(abs(x)./r0).^(5/3); %theoretical structure function
z1=0;
A=ones(N);
U=A*exp(-1i*k*z1);%plane wave field
I=abs(U).^2;%intensity of plane wave before atmospheric turbulence screen
Matlab Code for Atmospheric Turbulence Screen�
figure
imagesc(x,x,I);
axis xy; axis square; colorbar; %display intensity
xlabel('x (m)'); ylabel('y (m)');
title('Intensity of plane wave before atmospheric turbulence screen');
figure
plot(x,I(N/2+1,:));
xlabel('x (m)'); ylabel('mag U’);
title('slice of Intensity of plane wave before atmospheric turbulence screen’);
u2=U.*exp(1i*phi);%plane wave field expression after the screen
z2=[0 5 10 20 50 100]; %propagation distance for TF propagator
Matlab Code for Atmospheric Turbulence Screen�
for p=1:6
u3=propTF(u2,L,lambda,z2(p));
x2=x; %obs coords
y2=x; %obs coords
I2=abs(u3.^2); %obs intensity
figure %display obs intensity
imagesc(x2,y2,I2);
axis square; axis xy;
label('x (m)'); ylabel('y (m)');colorbar;
title(['z= ',num2str(z2(p)),' m']);
figure %intensity profile
plot(x2,I2(N/2+1,:));
xlabel('x (m)'); ylabel('Irradiance’);
title(['z= ',num2str(z2(p)),' m']);
end
Atmospheric Turbulence Screen�
Figure: Atmospheric turbulence screen
Structure Function�
Figure: Structure Function for atmospheric turbulence screen
Propagation of Plane Wave through Atmospheric Turbulence Screen
�Intensity of Plane Wave before Atmospheric Turbulence screen�
Figure: Intensity of plane wave before atmospheric turbulence screen (a) pattern (b) profile
(a)
(b)
Intensity at Observation Plane
Figure: Observation plane intensity at z= 0 m (a) pattern (b) profile
(a)
(b)
Intensity at Observation Plane
Figure: Observation plane intensity at z= 20 m (a) pattern (b) profile
(a)
(b)
Intensity at Observation Plane
Figure: Observation plane intensity at z= 50 m (a) pattern (b) profile
(a)
(b)
Intensity at Observation Plane
Figure: Observation plane intensity at z= 100 m (a) pattern (b) profile
(a)
(b)
Conclusions �