帮我写一个模拟不同密度的靶的卢瑟福散射和统计数据并画图的matlab代码

当入射粒子与靶原子发生碰撞时，会发生卢瑟福散射。通过测量散射电子的角度和能量，可以获得有关原子的信息，如原子尺寸和核电荷等。下面是一个简单的模拟代码来模拟不同密度的靶的卢瑟福散射和统计数据并画图的matlab代码。

clear all;
close all;

% Constants
e = 1.602e-19;  % elementary charge
eps0 = 8.854e-12;  % permittivity of free space
me = 9.109e-31;  % electron mass
hbar = 1.055e-34;  % Planck's constant over 2*pi

% Simulation parameters
N = 10000;  % number of particles
theta_min = 0;  % minimum scattering angle
theta_max = pi;  % maximum scattering angle
theta_step = 0.01;  % scattering angle step
E = 10e3 * e;  % kinetic energy of incident particle
Z = 2;  % atomic number of target atom

% Density range
rho_min = 1e-6;  % minimum density (kg/m^3)
rho_max = 10e3;  % maximum density (kg/m^3)
rho_step = 1e-5;  % density step (kg/m^3)

% Calculate differential cross section for each density
dcs = zeros(1, floor((rho_max - rho_min) / rho_step));
for i = 1:length(dcs)
    rho = rho_min + (i-1) * rho_step;  % density (kg/m^3)

    % Calculate the impact parameter for each particle
    b = sqrt((2 * Z * e^2 / (4 * pi * eps0 * me * E))^2 + (hbar^2 / (4 * me * E * sin(theta_min / 2)^2)));
    b_max = sqrt((2 * Z * e^2 / (4 * pi * eps0 * me * E))^2 + (hbar^2 / (4 * me * E * sin(theta_max / 2)^2)));
    b_list = linspace(b, b_max, N);

    % Calculate the scattering angle for each particle
    theta_list = zeros(1, N);
    for j = 1:N
        theta_list(j) = 2 * asin(b * sin(theta_min / 2) / b_list(j));
    end

    % Calculate the differential cross section
    dcs(i) = sum(sin(theta_list).^2) / N / theta_step;
end

% Plot the differential cross section vs density
figure;
plot(linspace(rho_min, rho_max, length(dcs)), dcs);
xlabel('Density (kg/m^3)');
ylabel('Differential Cross Section (m^2/kg)');
title('Differential Cross Section vs Density');

这个代码会计算不同密度的靶的卢瑟福散射的微分截面，并将结果绘制成图表。你可以根据需要修改参数来探索不同的情况。