MATLAB星空模拟 在二维空间中随机生成n个天体碎片，它们有各自的质量和初始速度。因为万有引力，它们之间相互吸引。如果两个天体距离太近，则会发生碰撞，碰撞的结果有可能合二为一，有可能再分裂为几个碎片。试解释合并和分裂的规则，并模拟这一动态过程。给出MATLAB代码

合并和分裂规则：

合并：当两个天体距离小于它们半径之和时，它们会合并为一个天体。合并后的天体质量等于两个天体质量之和，速度等于两个天体质量加权平均速度。

分裂：当一个天体质量小于一定阈值时，它会分裂成两个或多个碎片。分裂后的碎片质量和速度根据碎片数量等分。

MATLAB代码实现：

% 定义常量
G = 6.67e-11; % 万有引力常量
dt = 0.01; % 时间步长
T = 100; % 模拟时间
n = 50; % 天体碎片数量
m_min = 1e16; % 最小天体质量
m_max = 1e20; % 最大天体质量
v_max = 1e4; % 最大初始速度

% 随机生成天体碎片
m = (m_max - m_min) * rand(n, 1) + m_min; % 质量
r = sqrt(m / pi); % 半径
x = 2 * rand(n, 1) - 1; % 初始位置
y = 2 * rand(n, 1) - 1;
vx = (2 * rand(n, 1) - 1) * v_max; % 初始速度
vy = (2 * rand(n, 1) - 1) * v_max;

% 模拟天体碎片运动
for t = 0:dt:T
    % 计算天体之间的相对距离和引力
    dx = x - x.'; % 行向量减列向量，得到矩阵
    dy = y - y.';
    dr = sqrt(dx.^2 + dy.^2);
    Fx = G * m * m.' .* dx ./ dr.^3; % 引力公式
    Fy = G * m * m.' .* dy ./ dr.^3;
    Fx(isnan(Fx)) = 0; % 处理相同位置天体的引力为NaN的情况
    Fy(isnan(Fy)) = 0;
    Fx = sum(Fx, 2); % 每个天体受到的总引力
    Fy = sum(Fy, 2);
    
    % 计算天体碎片的加速度和速度
    ax = Fx ./ m;
    ay = Fy ./ m;
    vx = vx + ax * dt;
    vy = vy + ay * dt;
    
    % 计算天体碎片的位置
    x = x + vx * dt;
    y = y + vy * dt;
    
    % 处理碰撞
    for i = 1:n
        for j = i+1:n
            if dr(i,j) < r(i) + r(j) % 发生碰撞
                if m(i) < m(j) % 保证m(i) >= m(j)
                    temp = i;
                    i = j;
                    j = temp;
                end
                
                % 合并
                m(i) = m(i) + m(j);
                vx(i) = (m(i) * vx(i) + m(j) * vx(j)) / (m(i) + m(j));
                vy(i) = (m(i) * vy(i) + m(j) * vy(j)) / (m(i) + m(j));
                
                % 分裂
                if m(j) < m_min % 质量过小，分裂
                    n = n + 1;
                    m(n) = m(j) / 2;
                    r(n) = sqrt(m(n) / pi);
                    x(n) = x(j);
                    y(n) = y(j);
                    vx(n) = vx(j) + 2 * v_max * rand(1); % 随机速度
                    vy(n) = vy(j) + 2 * v_max * rand(1);
                    m(j) = m(j) / 2;
                    r(j) = sqrt(m(j) / pi);
                end
            end
        end
    end
    
    % 绘制天体碎片
    clf;
    hold on;
    axis([-1, 1, -1, 1]);
    for i = 1:n
        viscircles([x(i), y(i)], r(i));
    end
    drawnow;
end