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

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

```