%此程序是任意一个高斯光束在通过大气湍流之前的螺旋普
clc
clear all
N=300;
w0=0.03;
s=3;
z=1000;
lamda=1.550e-6;
k=2*pi/lamda;
zr = pi*w0.^2/lamda;
wz = w0*sqrt(1+(z/zr).^2);
b=2;%图像坐标范围
dx=b/N;
Cn2=1e-14;
dz =100;
Numz =z/dz;
m2=-N/2:N/2-1;
pp =linspace(0,b,N);
qq= linspace(0,2*pi,N);
[r,theta]=meshgrid(pp,qq);
L=0.7; %窗口宽度
df=1/L;
[fx,fy]=meshgrid(m2*df);%空间频率
fr=sqrt(fx.^2+fy.^2);
kx=2*pi*fx;
ky=2*pi*fy;
kr=2*pi*fr;%空间波数
u0 =sqrt(2./(pi*factorial(abs(s))))*(1/wz).*(r.*sqrt(2)/wz).^abs(s).*exp(-r.^2/(wz.^2)).*exp((-1i*k*r.^2*z)/(2*(z.^2+zr.^2))).*exp(1i*(abs(s)+1).*atan(z/zr)).*exp(1i*s*theta);
u=u0;
% phi=phase_screen(lamda,dz,Cn2);
phi=importdata('cn=2.5e-14.mat');
for j=1:Numz
u1=fft2(exp(1i*phi).*u); %1.先对u进行傅立叶变换，转换到空域上
u1=fftshift(u1);
u=ifft2(ifftshift(exp(1i*k.*dz.*sqrt(1-lamda.^2.*(fx.^2+fy.^2))).*u1));%2.并在空域上乘上变换算子,3.再逆变换回时域
end
m =-10:10;
for j = 1:length(m)
for i = 1:length(r)
am(i,j) = trapz(qq,u(:,i).*exp((-1j)*m(j).*qq'))/sqrt(2*pi);
end
end
for i = 1:length(m)
cm(i) = trapz(pp,(abs(am(:,i)).^2).*pp');
end
p = abs(cm)/sum(abs(cm));
bar(m, p);
评论4