|
% 参数设置
lambda = 840e-9; % 光波长 (单位：米)
k = 2 * pi / lambda; % 波数
n = 1.451; % 介质的折射率
cone_angle = 10 * pi / 180; % 锥的底角 (单位：弧度)
z_max = 10e-3; % 模拟传播的最大距离 (单位：米)
N = 512; % 图像的分辨率
dx = 3.45e-6; % 每个像素的尺寸 (单位：米)
dy = dx; % 假设为正方形像素 % 空间坐标系
x = (-N/2:N/2 - 1) * dx;
y = (-N/2:N/2 - 1) * dy;
= meshgrid(x, y);
Rxy = sqrt(X.^2 + Y.^2); % 径向坐标 % 锥透镜的相位调制
phi_cone = k * Rxy * cone_angle * (n-1); % 锥透镜的相位因子 % 初始平面波场
U0 = exp(1i * zeros(size(X))); % 平面波，初始相位为0 % 通过锥透镜调制光波
U_coned = U0 .* exp(-1i * phi_cone); % 平面波经过锥透镜调制 % 计算频率空间 (角谱传播)
fx = (-N/2:N/2-1) / (N * dx); % 频率坐标
fy = fx;
= meshgrid(fx, fy);
F = sqrt(Fx.^2 + Fy.^2); % 频率空间的径向坐标 % 传播传递函数，考虑传播距离 z_max
H = @(z) exp(1i * k * z * sqrt(1 - (lambda * F).^2)); % 角谱传播传递函数，动态传递 % 设置传播距离的采样点（例如，模拟从 0 到 z_max 的多个传播距离）
z_samples = linspace(0, z_max, 50); % 50 个传播距离 % 对每个传播距离进行计算并可视化
intensity_xz = zeros(N, length(z_samples)); % 预分配空间
for i = 1:length(z_samples)
z = z_samples(i);
% 计算传播到 z 处的场
U_freq = fftshift(fft2(U_coned)); % 2D 傅里叶变换
U_freq = U_freq .* H(z); % 应用传播传递函数
U_propagated = ifft2(ifftshift(U_freq)); % 逆傅里叶变换回到空间域 % 提取 y = 0 处的 xz 平面数据
y_idx = ceil(N / 2); % 选择 y = 0 处，即 Ny 中心
U_xz = U_propagated(:, y_idx); % 提取 xz 平面的数据 % 计算强度分布 |U(x, 0, z)|^2并保存
intensity_xz(:, i) = abs(U_xz).^2; % 直接存入预分配的矩阵
end % 绘制强度分布（xz平面的切片）
figure(1); % 5 行 10 列的子图 % 图像左右翻转使用 flip
intensity_xz_flipped = flip(intensity_xz, 2); % 对列进行翻转，相当于左右翻转
imagesc(z_samples, x, intensity_xz_flipped); colormap(‘parula’);
axis on; % 打开坐标轴
xlabel(‘z (m)’); % 添加 x 轴标签
ylabel(‘x (m)’); % 添加 y 轴标签
colorbar;
% 获取当前图像的坐标轴范围
xlim_val = xlim; % x 轴范围
ylim_val = ylim; % y 轴范围
% 在图像上添加文字注释
text(mean(xlim_val), mean(ylim_val)*0.2, ‘10°’, …
‘HorizontalAlignment’, ‘center’, ‘FontSize’, 10, ‘Color’, ‘yellow’); % 在图像上添加文字
% 调整坐标轴位置，去除不必要的空白
% set(gca, ‘Position’,
); % 设置坐标轴位置为图形窗口的整个区域 % 手动设置图形窗口为正方形
set(gcf, ‘Position’,
); 其中，z_max为锥透镜最大无衍射距离，现在代码这个部分不太正确，请结合锥透镜无衍射距离公式帮我修改上面代码
| 
