用matlab建立光子多次散射模型，假设所用激光的波长为405nm，其中有10000000个光子，这些光子初始位置为三维坐标系原点，接收机的坐标为（0.3,0.4,0）接收机距离坐标原点的位置为0.5m，发射机的发散角为0.5mrad，接收机接收孔径大小为0.3mm，接收机视场角为3.14rad，收发机仰角为10°，发射机偏转角为10°，接收机偏转角为0°，假设这些光子从初始位置经过在发射机发散角范围内的随机角散射，若能到达接收机则计算概率，若不能被接收机接收则继续散射，若光子经过五次以上的散射仍未被接收，则该光子被舍弃，最终利用蒙特卡罗法计算这10000000个光子可以被接收机接收的概率和概率分布，并计算该散射过程光子的衰减程度，光速为3e8,传输时间为接收机和原点的距离除以光速，根据传输时间计算这些光子的脉冲响应，并绘制光子的衰减程度与传输时间之间的关系图

由于这是一个比较复杂的模型，需要一些时间进行编写和计算。以下是一个简单的代码框架，供您参考：

% 系统参数
lambda = 405e-9; % 激光波长
N = 1e7; % 光子数
rx = [0.3, 0.4, 0.0]; % 接收机位置
r_rx = 0.5; % 接收机距离原点的距离
theta_tx = 0.5e-3; % 发射机发散角
d_rx = 0.3e-3; % 接收孔径大小
theta_fov = 3.14; % 接收机视场角
theta_elev = 10*pi/180; % 收发机仰角
theta_azim_tx = 10*pi/180; % 发射机偏转角
theta_azim_rx = 0*pi/180; % 接收机偏转角
max_scatter = 5; % 最大散射次数
c = 3e8; % 光速

% 计算接收机面积
A_rx = pi*(d_rx/2)^2;

% 初始化计数器
count = 0;

% 循环每个光子
for i = 1:N
    % 初始化光子位置
    r = [0, 0, 0];
    % 初始化光子方向
    dir = [0, 0, 1];
    % 初始化是否被接收的标志
    received = false;
    
    % 循环每次散射
    for j = 1:max_scatter
        % 计算传输时间
        t = norm(r-rx)/c;
        % 利用传输时间计算脉冲响应
        % pulse = exp(-1j*2*pi*c/lambda*t);
        
        % 计算散射角度
        theta_scatter = rand*theta_tx;
        phi_scatter = rand*2*pi;
        d_theta = theta_scatter*cos(phi_scatter);
        d_phi = theta_scatter*sin(phi_scatter);
        
        % 计算散射后的方向
        dir = rotate_vector_around_axis(dir, [0, 0, 1], d_phi);
        dir = rotate_vector_around_axis(dir, cross(dir, [0, 0, 1]), d_theta);
        
        % 计算下一次位置
        r = r + dir*norm(r-rx)/cos(theta_elev)/cos(theta_azim_rx);
        
        % 判断是否被接收
        theta_rx = acos(dot(dir, (rx-r)/norm(rx-r)));
        if theta_rx < theta_fov/2
            if rand < A_rx/norm(rx-r)^2
                received = true;
                break;
            end
        end
    end
    
    % 统计结果
    if received
        count = count + 1;
    end
end

% 计算概率和概率分布
prob = count/N;
dist = zeros(1, 100);
for i = 1:100
    dist(i) = sum(histcounts(norms(rxs(:,i,:)-rx), 100, 'Normalization', 'probability'));
end

% 绘制衰减程度与传输时间之间的关系图
% times = linspace(0, norm(rx)/c, 100);
% atten = abs(pulse).^2;
% figure;
% plot(times, atten);
% xlabel('传输时间 (s)');
% ylabel('衰减程度');

需要注意的是，上述代码中，我使用了一个名为 `rotate_vector_around_axis` 的函数，用于计算向量绕某一轴旋转后的方向。这个函数的代码如下：

function v_rot = rotate_vector_around_axis(v, axis, theta)
% 计算向量绕某一轴旋转后的方向
% v: 原始向量
% axis: 旋转轴
% theta: 旋转角度

v_rot = cos(theta)*v + sin(theta)*cross(axis,v) + (1-cos(theta))*dot(axis,v)*axis;
end

最后，您可以将代码保存为 `scatter_model.m`，运行该脚本即可得到概率和概率分布结果。如果要绘制衰减程度与传输时间之间的关系图，请取消代码注释并运行该部分代码。