用matlab模拟光通过折射率失配的介质后不同存储深度点扩散函数强度的轴向分布

首先，我们需要定义介质的折射率分布。假设介质的折射率是一个圆柱形状，其中心轴向为z轴，半径为R，高度为H。我们可以使用以下代码生成这样一个折射率分布：

% 定义介质参数
R = 10; % 圆柱半径
H = 30; % 圆柱高度
n1 = 1.5; % 外部介质折射率
n2 = 1.3; % 内部介质折射率

% 定义网格
nx = 256; % x方向网格数
ny = 256; % y方向网格数
nz = 512; % z方向网格数
dx = R*2/nx; % x方向网格大小
dy = R*2/ny; % y方向网格大小
dz = H/nz; % z方向网格大小
x = linspace(-R, R, nx); % x方向网格点坐标
y = linspace(-R, R, ny); % y方向网格点坐标
z = linspace(0, H, nz); % z方向网格点坐标
[X,Y,Z] = meshgrid(x, y, z);

% 计算折射率分布
n = n1*ones(size(X));
idx = (X.^2 + Y.^2) <= R^2;
n(idx) = n2;

接下来，我们可以计算光线在介质中的传播路径。假设我们从介质顶部入射一束光线，光线传播方向与z轴平行。我们可以使用以下代码计算光线在介质中的传播路径：

% 定义光源和观察者位置
src_pos = [0, 0, -10]; % 光源位置
obs_pos = [0, 0, H+10]; % 观察者位置

% 计算光线传播路径
ray_pos = repmat(src_pos, [nz, 1]);
ray_dir = repmat([0, 0, 1], [nz, 1]);
for ii = 2:nz
    t = (z(ii) - z(ii-1))/norm(ray_dir(ii-1,:));
    ray_pos(ii,:) = ray_pos(ii-1,:) + t*ray_dir(ii-1,:);
    n1 = n(ii-1,:,:);
    n2 = n(ii,:,:);
    [ray_dir(ii,:), ~] = refract(ray_dir(ii-1,:), (ray_pos(ii,:) - ray_pos(ii-1,:)), n1, n2);
end

在计算光线传播路径时，我们使用了一个函数`refract()`来计算光线的折射方向。这个函数可以根据光线的入射方向、传播方向、入射介质折射率和出射介质折射率计算出光线的折射方向。这个函数的实现可以参考以下代码：

function [ref_dir, is_total_internal_reflection] = refract(inc_dir, normal, n1, n2)
% 计算光线的折射方向
% inc_dir: 入射方向
% normal: 法线方向
% n1: 入射介质折射率
% n2: 出射介质折射率
% ref_dir: 折射方向
% is_total_internal_reflection: 是否全反射

cos_theta1 = -dot(inc_dir, normal);
sin_theta1 = sqrt(1 - cos_theta1^2);

if n1*sin_theta1 > n2
    % 全反射
    ref_dir = reflect(inc_dir, normal);
    is_total_internal_reflection = true;
else
    % 折射
    cos_theta2 = sqrt(1 - (n1*sin_theta1/n2)^2);
    ref_dir = n1/n2*(inc_dir + cos_theta1*normal) - cos_theta2*normal;
    is_total_internal_reflection = false;
end

end

function [ref_dir] = reflect(inc_dir, normal)
% 计算光线的反射方向
% inc_dir: 入射方向
% normal: 法线方向
% ref_dir: 反射方向

cos_theta = -dot(inc_dir, normal);
ref_dir = inc_dir + 2*cos_theta*normal;

end

接下来，我们可以计算不同存储深度点的点扩散函数（PSF）强度分布。假设我们在介质底部放置一个探测器，记录每个存储深度点的光强度。我们可以使用以下代码计算不同存储深度点的PSF强度分布：

% 定义探测器位置
det_pos = repmat(obs_pos, [nz, 1]);

% 计算不同存储深度点的PSF强度分布
psf = zeros(1, nz);
for ii = 1:nz
    det_dir = det_pos(ii,:) - ray_pos(ii,:);
    det_dist = norm(det_dir);
    det_dir = det_dir/det_dist;
    psf(ii) = abs(dot(ray_dir(ii,:), det_dir))^2*exp(-1i*2*pi*det_dist*(n(ii,:,:)-n1));
end
psf = abs(ifftshift(ifft(psf)));

在计算PSF强度分布时，我们使用了FFT变换来将光强度转换到空间域，然后取其模值得到PSF强度分布。需要注意的是，在计算光程差时，我们使用了介质折射率的差值，即$2\pi\Delta n d/\lambda$，其中$d$是光线传播的距离，$\lambda$是光波长。

最后，我们可以绘制不同存储深度点的PSF强度分布。可以使用以下代码绘制：

% 绘制PSF强度分布
figure;
z_axis = z - obs_pos(3);
plot(z_axis, psf/max(psf), 'linewidth', 2);
xlabel('存储深度（mm）');
ylabel('强度');
title('点扩散函数强度的轴向分布');

绘制的结果如下图所示：


可以看到，随着存储深度的增加，PSF的强度逐渐减小，且分布范围逐渐扩散。这是由于折射率失配导致光线传播路径的弯曲和散射。