Arbitrary screen diffraction

clear（用实时脚本打开）
%%%%%%%衍射模板%%%%%%%
N=500;                                %采样点 500*500像素
m=N;
n=N;
R=zeros(m,n);                         %预设矩阵存放图像数据
A = imread('picture.bmp');            %衍射屏图片 黑色为不透光 白色为透光（用photoshop制作，文件注意为bmp格式，放在同一文件夹即可）
A=A(:,:,1);
R=im2double(A);                       %屏函数
%%%%%%%光波数据%%%%%%%
lamda=632.2e-9;                       %波长（nm）
U=ones(N);                            %正入射平行光
dfx=1./m;
dfy=1./n;
xd=N;                                 %x方向采样点
yd=N;                                 %y方向采样点
F1=R.*U;                              %经过衍射屏得光波
F2=fft2(F1);
%%%%%%%衍射图案%%%%%%%
Ui=fftshift(F2);                      %计算光波通过衍射屏后在频率域的复振幅分布
d=1.5;                                %光波传播距离
T=exp(j*d*2*pi/lamda)/(j*lamda*d);    %传递函数
Ud=Ui.*T;                             %计算光波传播距离d后在频率域的复振幅分布
U0=ifft2(Ud);                         %计算光波在观察屏的复振幅分布
U01=(abs(Ud)).^2;                     %复振幅转化为光强
U01=U01.*10^-16;                      %放大，使得衍射可视化
%%%%%%%画图%%%%%%%
subplot(2,2,1);
imshow(R);
subplot(2,2,2);
imshow(U01);
subplot(2,2,3);
plot(U01(251,:));
subplot(2,2,4);
meshz(U01);
（引流一下---B站up主：学物理的兔八哥）