请用matlab编写单个阿尔法粒子的卢瑟福散射运动轨迹

在MATLAB中编写单个阿尔法粒子的卢瑟福散射运动轨迹涉及到经典物理学中的卢瑟福散射模型，该模型描述了带正电的阿尔法粒子与原子核相互作用时的散射过程。以下是一个简单的MATLAB脚本示例，用于模拟单个阿尔法粒子在受库仑力作用下的散射轨迹：

% MATLAB代码：单个阿尔法粒子的卢瑟福散射运动轨迹

% 初始化参数
Z = 78; % 原子序数，例如金
alpha_q = 2; % 阿尔法粒子的电荷数，阿尔法粒子有2个质子和2个中子
alpha_m = 4; % 阿尔法粒子的质量，单位为电子质量
alpha_v = 2e7; % 阿尔法粒子的初始速度，单位为m/s
theta = 20 * pi / 180; % 散射角度，转换为弧度
phi = 0; % 散射平面内的方位角，这里假设为0
e = 1.6e-19; % 基本电荷
epsilon_0 = 8.854e-12; % 真空电容率
r0 = 2e-15; % 原子核半径，仅为示例值

% 时间参数
dt = 1e-24; % 时间步长，单位为秒
t_max = 1e-20; % 模拟时间，单位为秒
t = 0:dt:t_max; % 时间数组

% 空间参数
r = zeros(3, length(t)); % 粒子位置数组，初始化为零
v = zeros(3, length(t)); % 粒子速度数组，初始化为零
v(:,1) = [alpha_v*sin(theta)*cos(phi), alpha_v*sin(theta)*sin(phi), alpha_v*cos(theta)]; % 初始速度向量

% 库仑力常数
k = (alpha_q * e^2) / (4 * pi * epsilon_0);

% 运动方程求解
for i = 1:length(t)-1
    r(:,i+1) = r(:,i) + v(:,i) * dt; % 更新位置
    r_norm = norm(r(:,i));
    if r_norm < r0
        error('Alpha particle entered the nucleus, simulation ended.');
    end
    force = k * alpha_q * (r(:,i) / r_norm^3); % 库仑力
    a = force / alpha_m; % 加速度
    v(:,i+1) = v(:,i) + a * dt; % 更新速度
end

% 绘制轨迹
plot3(r(1,:), r(2,:), r(3,:));
xlabel('x (m)');
ylabel('y (m)');
zlabel('z (m)');
title('Alpha Particle Scattering Trajectory');
grid on;

这段代码中，我们首先定义了阿尔法粒子和目标原子核的一些基本参数，包括阿尔法粒子的电荷、质量、初始速度，以及原子核的电荷数（原子序数Z）和质量。我们假设了一个散射角度，并设置了时间步长和总模拟时间。然后，我们初始化了粒子的位置和速度，并通过一个循环来计算每一时间步长下的粒子位置和速度，同时考虑了库仑力作用下的加速度。

请注意，上述代码简化了物理模型，并且假设原子核是静止的。在实际情况中，原子核并非静止，同时还需要考虑相对论效应、量子力学效应等因素，因此实际的模拟会更加复杂。