power_sprecum (2)_oam_OAM谱_OAM相位_OAM相位_OAMpowersprecum_

%此程序是任意一个高斯光束在通过大气湍流之前的螺旋普 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);