数值模拟结果不对，求修改

Question

matlab程序。用拉盖尔高斯光束和平面波干涉的干涉条纹做光栅，发生菲涅尔衍射，得到拉盖尔高斯光束。现程序无报错，但结果和推导不符。

```bash
clc
close all

N = 1024*3;            % 取样点数
lambda = 632.8e-9;   % 波长632.8nm
k = 2*pi/lambda;     % 波数
w0 = 5e-2;
x = linspace(-5e-6, 5e-6, N); y = x;
[X, Y] = meshgrid(x, y);
[theta, r] = cart2pol(X, Y);
p = 0;                 % 径向指数
Z_R = pi*w0^2/lambda;  % 瑞利长度
z = 0;
w_z = w0*sqrt(1 + (z/Z_R)^2);  % 光束在z位置的半径
L=1;

% 计算光束的复振幅
E1 = sqrt(2*factorial(p)/pi/(p+factorial(abs(L))))*(1/w_z)*(sqrt(2)*r/w_z).^abs(L)...
        .*exp(-r.^2/w_z^2).*laguerre(p, abs(L), 2*r.^2/w_z^2).*exp(-1i*L*theta).*exp(-1i*k*z)...
        .*exp(-1i*k*r.^2*z/2/(z^2+Z_R^2))*exp(-1i*(2*p+abs(L)+1)*atan(z/Z_R));
E2 = exp(-1i*k*X);       % 平面波
C1 = E1 + E2;           % 干涉场的复振幅

% 计算干涉图样的强度分布
I = C1.*conj(C1);
% 菲涅尔衍射
z = 4e-3 ;  % 菲涅尔衍射距离
T = fftshift(fft2(I));% 光栅的频域表示
H = exp(1i*k*z)*exp(-1i*pi*lambda*z*(X.^2 + Y.^2));  % 菲涅尔传播函数
U = H.*T;  % 菲涅尔衍射
u = ifftshift(ifft2(U));
uu= u.*conj(u);
% 显示衍射图像
figure ;
imagesc(x, y, uu);
axis square;
colormap gray;
title('Fresnel Diffraction');
xlabel('x');
ylabel('y');
%拉盖尔多项式
function result = laguerre(p,L,x)
result = 0;
if p == 0    
    result = 1;
elseif p == 1    
result = 1+abs(L)-x;
else   
  result = (1/p)*((2*p+L-1-x).*laguerre(p-1,abs(L),x)-(p+L-1)*laguerre(p-2,abs(L),x));
end
end

```

Answer 1

不知道你这个问题是否已经解决, 如果还没有解决的话:

• 这个问题的回答你可以参考下: https://ask.csdn.net/questions/7509232
• 这篇博客你也可以参考下：Matlab：涡旋光束与平面光波、球面光波干涉
• 除此之外, 这篇博客: 利用MATLAB绘制箱线图—箱线图在高光谱图像处理中的应用中的 一、箱线图在高光谱图像处理中的应用 部分也许能够解决你的问题, 你可以仔细阅读以下内容或者直接跳转源博客中阅读:
  

  箱线图也称箱须图、箱形图、盒图，用于反映一组或多组连续型定量数据分布的中心位置和散布范围。箱形图包含数学统计量，不仅能够分析不同类别数据各层次水平差异，还能揭示数据间离散程度、异常值、分布差异等等。

  目前已经在异常检测，目标探测，红外小目标检测、分类等多个方向上有广泛的应用，具体箱线图的含义，读者可以百度，或者查看第三部分的相关参考文献，此处不再赘述。

如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^