clc
clear
close all
%% 物理与模拟参数（已优化为星系参数）
G = 1.0; % 引力常数（调大让引力更明显）
rho = 0.001; % 密度
Np = 150; % 星体数量
dt = 0.1; % 时间步长（更小更稳定）
totalSteps = 50000; % 总步数
drawStep =200; % 绘图间隔
%% 边界（固定为可视范围）
lim = 200;
liml = -lim;
limu = lim;
%% 1. 初始化：星系核心+旋转速度（关键！星系必须有切向旋转）
Xa = zeros(Np, 2); % 位置
U = zeros(Np, 2); % 速度
R = rand(Np, 1)*1.5 + 2; % 星体半径
M = rho * pi * R.^2; % 质量
% 生成圆盘状初始分布（星系是圆盘，不是随机方块）
for i = 1:Np
r = (rand()+0.5) * 100; % 径向距离
theta = rand() * 2*pi; % 角度
Xa(i,1) = r * cos(theta);
Xa(i,2) = r * sin(theta);
% 关键：给初始旋转速度（开普勒旋转，星系形成的核心）
v = sqrt(G * sum(M) / (r + 10)); % 轨道速度
U(i,1) = -v * sin(theta); % 切向速度
U(i,2) = v * cos(theta);
end
%% 2. 模拟主循环
for tStep = 1:totalSteps
Ac = zeros(Np, 2); % 加速度归零
% N体引力计算（优化了距离、软ening、力计算）
for i = 1:Np
for j = i+1:Np
% 位置差
dx = Xa(i,1) - Xa(j,1);
dy = Xa(i,2) - Xa(j,2);
distSq = dx^2 + dy^2;
dist = sqrt(distSq);
% 软ening长度：防止距离为0导致爆炸（关键修复）
softening = 2;
dist = max(dist, softening);
disR=R(i)+R(j)+0.5;
% 万有引力公式
if dist>disR
F = G * M(i) * M(j) / dist^2;
else
F=-0.0001/abs((dist-R(i)-R(j)-1));
end
% 单位向量
fx = dx / dist;
fy = dy / dist;
% 加速度更新
Ac(i,1) = Ac(i,1) - F/M(i) * fx;
Ac(i,2) = Ac(i,2) - F/M(i) * fy;
Ac(j,1) = Ac(j,1) + F/M(j) * fx;
Ac(j,2) = Ac(j,2) + F/M(j) * fy;
end
end
% 3. 速度 + 位置更新（欧拉积分，加极弱阻尼保持稳定）
damping = 0.9999;
U = U * damping + Ac * dt;
Xa = Xa + U * dt;
% 4. 边界：不循环！直接固定可视范围（星系不会乱跑）
Xa(Xa > limu) = limu;
Xa(Xa < liml) = liml;
% 5. 实时绘图
if mod(tStep, drawStep) == 0
clf; hold on; axis equal; axis off;
set(gcf,'Color',[0.1 0.1 0.2]); % 黑色背景更像宇宙
set(gca, ...
'XLim', [liml,limu], ... % X 轴显示范围
'YLim', [liml,limu]); % Y 轴显示范围
% 绘制所有星体
for i = 1:Np
ang = 0:0.1:2*pi;
xCircle = Xa(i,1) + R(i)*cos(ang);
yCircle = Xa(i,2) + R(i)*sin(ang);
fill(xCircle, yCircle, [1,1,0.8], 'EdgeColor','none');
end
hold off;
export_fig(gcf,'-png', '-r300', '-opengl', ['PNG/',num2str(100000+tStep),'.png'],'-nocrop');
% drawnow;
end
end