clc;clear;close;
figure('name','三体运动');
axis equal
grid on
hold on
colordef black
view(3)
xlabel('x')
ylabel('y')
zlabel('z')
hold on
%初始化参数设置(质量、初始位置、初始速度、万有引力常数、时间步长)
m1=rand()*10^23;
m2=rand()*10^23;
m3=rand()*10^23;
pos1=[rand()*10^7,rand()*10^7,rand()*10^7];
pos2=[rand()*10^7,rand()*10^7,rand()*10^7];
pos3=[rand()*10^7,rand()*10^7,rand()*10^7];
v1=[rand()*1000,rand()*1000,rand()*1000];
v2=[rand()*1000,rand()*1000,rand()*1000];
v3=[rand()*1000,rand()*1000,rand()*1000];
G=6.67*10^(-11);
dt=0.005;
%行星和轨迹（定义行星位置和轨迹变量）
planet1=plot3(pos1(1),pos1(2),pos1(3),'b:.','markersize',20);
planet2=plot3(pos2(1),pos2(2),pos2(3),'r:.','markersize',20);
planet3=plot3(pos3(1),pos3(2),pos3(3),'y:.','markersize',20);
h1=animatedline('color','b');
h2=animatedline('color','r');
h3=animatedline('color','y');
frame=0;
while true
%不断循环计算距离、相互作用力、速度、加速度、位置
r12=normest(pos2-pos1);
r23=normest(pos3-pos2);
r13=normest(pos3-pos1);
F12_val=G*m1*m2/(r12^2);
F23_val=G*m2*m3/(r23^2);
F13_val=G*m1*m3/(r13^2);
F12_dir=(pos2-pos1)/normest(pos2-pos1);
F23_dir=(pos3-pos2)/normest(pos3-pos2);
F13_dir=(pos3-pos1)/normest(pos3-pos1);
a1=((F12_val)*(F12_dir)+(F13_val)*(F13_dir))/m1;
a2=((-F12_val)*(F12_dir)+(-F23_val)*(F23_dir))/m2;
a3=((-F13_val)*(F13_dir)+(-F23_val)*(F23_dir))/m3;
v1=v1+a1*dt;
v2=v2+a2*dt;
v3=v3+a3*dt;
pos1=pos1+v1*dt+1/2*a1*dt^2;
pos2=pos2+v2*dt+1/2*a2*dt^2;
pos3=pos3+v3*dt+1/2*a3*dt^2;
frame=frame+1;
if frame==10000
frame=0;
% set(planet1,'position',[pos1(1),pos1(2),pos1(3)]);
%set(planet3,'position',pos3(1),'YData',pos3(2),'ZData',pos3(3));
addpoints(h1,pos1(1),pos1(2),pos1(3));
addpoints(h2,pos2(1),pos2(2),pos2(3));
addpoints(h3,pos3(1),pos3(2),pos3(3));
drawnow
end
if r12==0||r23==0||r13==0 %判断是否相撞
break
end
end
