用matlab仿真 光线在不同折射率的介质中传播时光场内任一点处的点扩散函数图像分布

在MATLAB中，可以使用ray tracing方法模拟光线在不同折射率的介质中传播，并计算光场内任一点处的点扩散函数图像分布。

以下是一个简单的光线追踪模拟程序：

% 定义折射率
n1 = 1; % 空气
n2 = 1.5; % 玻璃

% 定义光线和起点
theta = linspace(0, 2*pi, 20);
x0 = 0;
y0 = 0;

% 计算光线路径
x = zeros(length(theta), 2);
y = zeros(length(theta), 2);
for i = 1:length(theta)
    [x(i,:), y(i,:)] = raytrace(n1, n2, x0, y0, theta(i));
end

% 绘制光线路径
figure; hold on;
xlim([-10 10]); ylim([-10 10]);
plot(x', y', 'b');
scatter(x0, y0, 'r', 'filled');

% 计算光场内任一点处的点扩散函数
[X, Y] = meshgrid(-10:0.1:10);
Z = zeros(size(X));
for i = 1:length(theta)
    [Ztmp, ~] = pointspread(n1, n2, x0, y0, theta(i), X, Y);
    Z = Z + Ztmp;
end

% 绘制点扩散函数图像
figure; imagesc(Z); axis image; colormap hot;

其中，raytrace函数计算光线路径，pointspread函数计算点扩散函数。这两个函数的实现可以参考以下代码：

function [x, y] = raytrace(n1, n2, x0, y0, theta)
% 计算光线路径

% 计算入射角
theta1 = atan2(y0, x0);
phi1 = pi/2 - theta1;
phi2 = asin(n1/n2*sin(phi1));
theta2 = pi/2 - phi2;

% 计算折射后的光线方向
theta3 = theta1 - theta2 + theta;
x3 = cos(theta3);
y3 = sin(theta3);

% 计算折射后的光线起点
x = [x0 x0 + 10*x3];
y = [y0 y0 + 10*y3];

end

function [Z, R] = pointspread(n1, n2, x0, y0, theta, X, Y)
% 计算点扩散函数

% 计算入射角
theta1 = atan2(y0, x0);
phi1 = pi/2 - theta1;
phi2 = asin(n1/n2*sin(phi1));
theta2 = pi/2 - phi2;

% 计算折射后的光线方向
theta3 = theta1 - theta2 + theta;
x3 = cos(theta3);
y3 = sin(theta3);

% 计算折射后的光线起点
x = X - x0;
y = Y - y0;
r = sqrt(x.^2 + y.^2);
x1 = x./r;
y1 = y./r;
x2 = n1/n2*x1;
y2 = n1/n2*y1;
x4 = cos(theta3);
y4 = sin(theta3);

% 计算点扩散函数
Z = zeros(size(X));
for i = 1:size(X, 1)
    for j = 1:size(X, 2)
        if r(i,j) > 0
            [xtmp, ytmp] = lineintersect([x1(i,j) y1(i,j)], [x2(i,j) y2(i,j)], [x3 y3], [x4 y4]);
            if ~isempty(xtmp)
                d = sqrt((xtmp-x0)^2 + (ytmp-y0)^2);
                Z(i,j) = Z(i,j) + sinc(d*pi/(n2*500e-9))^2;
            end
        end
    end
end
R = [];
end

function [x, y] = lineintersect(a, b, c, d)
% 计算两条直线的交点

s1 = a(2) - b(2);
t1 = b(1) - a(1);
s2 = c(2) - d(2);
t2 = d(1) - c(1);
denom = s1*t2 - s2*t1;
if denom == 0
    x = [];
    y = [];
else
    num1 = a(1)*b(2) - a(2)*b(1);
    num2 = c(1)*d(2) - c(2)*d(1);
    x = (num1*t2 - num2*t1) / denom;
    y = (num2*s1 - num1*s2) / denom;
end

end

注意，这段代码仅为示例，没有考虑多次反射和折射等情况，仅仅是一个简单的光路追踪程序。如果需要更加精确的模拟结果，需要考虑更多的物理因素。