用matlab模拟涡旋光经过不同掩膜的衍射图样

程序是根据初始光场在掩膜区域内的傅里叶变化写的，便于修改初始光场和掩膜。程序如下：

close all;clear;clc;

lambda = 5.37e-6;
w0     = 3e-3; 
z      = 100*lambda;                   %传播到z处的距离
k      = 2*pi/lambda;
N      = 500;
x      = linspace(-10*w0,10*w0,N);     %x轴500等分
y      = x;                            %y轴500等分
[X,Y]  = meshgrid(x,y);                %生成直角坐标网格矩阵


ELG = generateLG(lambda, w0, 0, 5, z, X, Y); %初始光场涡旋光


a=100*w0;%长
b=0.2*w0;%宽
t1=zeros(N,N);%透射空间
for i =1:N%在透射空间中画出矩孔，内部透射率为1
    for j =1:N
        if (X(i,j)>=-0.5*a) && (Y(i,j)>=-0.5*b) && (X(i,j)<=0.5*a) && (Y(i,j)<=0.5*b)%判断区域
            t1(i,j)=1;%透射率为1
        end 
    end
end
%d=fftshift(abs(fft2(ELG.*t1))); %快速傅里叶变换，并移动坐标为原点
z1=50*w0;                         %衍射距离
d=dfft(lambda, z1, X, Y, ELG.*t1);%菲涅尔衍射

figure(1)
mesh(X,Y,abs(d));%三维

% figure(2)%将缝和光场放一块对比，方便观察相对大小
% mesh(X,Y,t1)
% hold on
% mesh(X,Y,abs(ELG))

观察高斯光将初始光场改为1即可，其余掩膜修改透过率函数。
下面是菲涅尔衍射子函数：

function u1= dfft(lambda, d, X, Y, u0)
%%
% 两次fft:一次fft和一次ifft，所以叫dfft算法
% 该算法求得的观察平面复振幅对应的空间坐标与初始物光的空间坐标X和Y一致

% input: 波长lambda，衍射距离d，平面网格坐标X和Y，衍射平面(物光)复振幅u0
% output:观察平面复振幅

% 注：会有离散误差，当观察平面的衍射场不符合预期时，请增大采样数目（即加密网格）或改变衍射距离
% 注：约定复振幅的相位为exp(j*k*r)
% 注：输入的平面坐标的原点需在中间（如：[-2 -1 0 1]或者[-2 -1 0 1 2]）

%%
% 波数
k = 2*pi/lambda;

% 采样数目
[Ynum, Xnum] = size(X);

% 空间采样间隔
delx = abs(X(1, 2) - X(1, 1));
dely = abs(Y(2, 1) - Y(1, 1));

% 空间频率域采样间隔
delfx = 1/(Xnum*delx);
delfy = 1/(Ynum*dely);

% 空间频率坐标网格
% fx = (-Xnum/2:Xnum/2-1)*delfx;
% fy = (-Ynum/2:Ynum/2-1)*delfy;
Fx = X/delx * delfx;
Fy = Y/dely * delfy;

% 传递函数
H = exp(1i*k*z*(1-lambda^2/2*(Fx.^2 + Fy.^2)));
u0FFT = fftshift(fft2(u0));
u1 = ifft2(ifftshift(u0FFT.*H));

end