入射振幅1的初始相位为0的平面光，以标准的双曲相位为调制，并没有聚焦的光强分布，算法优化超透镜时每次迭代用lumericalFDTD软件计算，如何解决？（相关搜索：matlab计算|公式计算）

Question

我在做超透镜仿真，算法优化超透镜时每次迭代用lumericalFDTD软件计算，特别耗时。想转用matlab计算试试看，应该会更省时间，但是我在写代码时遇到了困难。入射振幅1的初始相位为0的平面光，以标准的双曲相位为调制，采用菲涅尔衍射公式计算，并没有聚焦的光强分布，结果如下

这是我的代码

%% --------------------------------------------------------------------------%%
%% 清除空间
clc;
clear ;
%%
%%以双曲相位为例的相位屏
f=100e-6;
lambda=1.064*1e-6; 
M=50;
resolution=M;
p=0.5e-6;  
r=(p*1e6*M/2)*1e-6;
x0=linspace(-r+p/2,r-p/2,M);
y0=linspace(-r+p/2,r-p/2,M);
[X0,Y0]=meshgrid(x0, y0);
k0=2*pi*(1/lambda);
nt=1;
phase_contour = -k0*nt*(sqrt(X0.^2+Y0.^2+f^2)-f);
phase_contour=wrapTo2Pi(phase_contour);
%细化相位屏
new_size = 100;
[X1, Y1] = meshgrid(linspace(1, 50, new_size), linspace(1, 50, new_size));
smoothed_phase_contour = interp2(phase_contour, X1, Y1, 'spline');
figure;
subplot(1, 2, 1);
imagesc(phase_contour);
title('原始相位轮廓');
subplot(1, 2, 2);
imagesc(smoothed_phase_contour);
title('插值后的相位轮廓');
%%
%%计算传输
phase = smoothed_phase_contour;%
E1=exp(1i* phase );%振幅为1初始相位为0的平面光穿过相位屏幕后的复振幅
z_f =f;%要观察的远场位置
E2=zeros(new_size,new_size);%远场光强
e2=zeros(new_size,new_size);%计算的一部分
x1=linspace(-r+p/2,r-p/2,new_size);
y1=linspace(-r+p/2,r-p/2,new_size);
[X1, Y1] = meshgrid(x1, y1);

e1=exp(1i*k0*z_f)/(1i*lambda*z_f);
for a=1:new_size
    for b=1:new_size
        e2=1i*k0*(X1(a)-X1).^2/2/z_f+1i*k0*(Y1(b)-Y1).^2/2/z_f;
        e3=sum(sum(E1.*exp(e2)));
        E2(a,b)=e1*e3;
    end
end
intensity=E2.*conj(E2);
figure;
imagesc(x0,y0, intensity);
title('远场光强分布');

这是菲涅尔衍射公式

Answer 1

引用文心一言及其思考回答：
从您的代码中可以看出，您正在尝试使用Matlab模拟并计算光的衍射过程，然后计算并显示光强的分布。

首先，我注意到您的代码中可能存在一些问题，我会逐一解释并提供可能的解决方案。

1. 在计算相位轮廓的部分，您使用了wrapTo2Pi函数，这个函数会把相位映射到[0, 2π)的区间内。然而在您的代码中，相位被映射到了一个更小的区间，这可能会导致错误。如果需要将相位映射到[0, 2π)，您应该使用mod函数而不是wrapTo2Pi。
2. 插值后相位轮廓图像显示部分可能会出现错误，因为插值后的相位轮廓是一个二维数组，不能直接用imagesc来显示。您可以使用imshow函数来显示插值后的相位轮廓。
3. 在计算传输的部分，我注意到您的代码可能有一些混淆的地方。在计算e2时，您使用了X1(a)-X1和Y1(b)-Y1，但这些变量应该使用小写的x1和y1代替，因为在后面的循环中您使用了大写的X1和Y1来定义网格。
4. 在计算e3时，您使用了sum(sum(E1.*exp(e2)))，这里应该使用小写的e2代替大写的E2，因为您在前面的循环中已经定义了小写的e2。
5. 在计算远场光强分布的部分，我注意到您的代码可能存在一些错误。首先，您在计算远场光强分布时使用了变量x0和y0，但这些变量并没有在前面的代码中定义。其次，您应该使用小写的x1和y1来代替大写的X1和Y1，因为在前面的循环中您已经定义了小写的x1和y1。

根据以上建议，我为您的代码进行了修改并重新排列了部分代码：

```matlab
% 清除空间
clc;
clear ;

% 以双曲相位为例的相位屏
f = 100e-6;
lambda = 1.0641e-6;
M = 50;
resolution = M;
p = 0.5e-6;
r = (p1e6M/2)1e-6;
x0 = linspace(-r+p/2, r-p/2, M);
y0 = linspace(-r+p/2, r-p/2, M);
[X0, Y0] = meshgrid(x0, y0);
k0 = 2pi(1/lambda);
nt = 1;
phase_contour = -k0nt(sqrt(X0.^2+Y0.^2+f^2)-f);
phase_contour = mod(phase_contour, 2*pi); % 映射到[0, 2π)区间
% 细化相位屏
new_size = 100;
[X1, Y1] = meshgrid(linspace(1, 50, new_size), linspace(1, 50, new_size));
smoothed_phase_contour = interp2(phase_contour, X1, Y1, 'spline');
figure;
subplot(1, 2, 1);
imagesc(phase_contour);
title('原始相位轮廓');
subplot(1, 2, 2);
imshow(smoothed_phase_contour); % 使用imshow来显示插值后的相位轮廓
title('插值后的相位轮廓');

% 计算传输
phase = smoothed_phase_contour; %
E1 = exp(1i*phase); % 振幅为1初始相位为0的平面光穿过相位屏幕后的复振幅
z_f = f; % 要观察的远场位置
E2 = zeros(new_size, new_size); % 远场光强
e2 = zeros(new_size, new_size); % 计算的一部分
x1 = linspace(-r+p/2, r-p/2, new_size);
y1 = linspace(-r+p/2, r-p/2, new_size);
[X1, Y1] = meshgrid(x1, y1); % 使用小写的x1和y1代替大写的X1和Y1

e1 = exp(1ik0