Matlab模拟三体系统

clc
clear
clf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% G=6.754*10^-11;
% M1=1.989*10^30;
% M2=1.989*10^30;
% M3=1.989*10^30;
% s1=[01*10^9.46*10^12];
% s2=[1*10^9.46*10^120];
% s3=[-1*10^9.46*10^120];
% v1=[03000];
% v2=[00];
% v3=[00];

%等边
G=1;
M1=25000;
M2=25000;
M3=25000;
T=3;
r=（（（T/2/pi）^2*10）^（1/3））*14.0005;
v=2*pi*r/T;
s1=[-r/2*（3^（1/2））0];
s2=[r/2*（3^（1/2））0];
s3=[0r];
v1=[1/2*v-3^0.5/2*v];
v2=[1/2*v3^0.5/2*v];
v3=[-v0];

range=200;%监视范围
dt=0.01;%次数间隔dt秒，dt时间内认为物理量均为线性变化
time=20000;%最大次数
pause_time=0.01;%控制运行速度
show=time;%150;%尾迹长度

% G=3;
% M1=2500;
% M2=2500;
% M3=2500;
% T=3;
% r=（（（T/2/pi）^2*10）^（1/3））;
% v=2*pi*r/T;
% s1=[-30*r0];
% s2=[30*r0];
% s3=[030*r];
% v1=[0v];
% v2=[0-2*v];
% v3=[v0];   
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%函数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
point_number=time;

co_x1=[];
co_y1=[];
co_x2=[];
co_y2=[];
co_x3=[];
co_y3=[];

co_x1=[co_x1 s1（1）];
co_y1=[co_y1 s1（2）];
co_x2=[co_x2 s2（1）];
co_y2=[co_y2 s2（2）];
co_x3=[co_x3 s3（1）];
co_y3=[co_y3 s3（2）];

for calucate_count = 1:point_number;
    
    a1=G*（M2*（s2-s1）/（norm（s2-s1））^3+M3*（s3-s1）/（norm（s3-s1））^3）;
    a2=G*（M1*（s1-s2）/（norm（s1-s2））^3+M3*（s3-s2）/（norm（s3-s2））^3）;
    a3=G*（M1*（s1-s3）/（norm（s1-s3））^3+M2*（s2-s3）/（norm（s2-s3））^3）;
    
    s1=s1+v1*dt;
    s2=s2+v2*dt;
    s3=s3+v3*dt;
   
    v1=v1+a1*dt;
    v2=v2+a2*dt;
    v3=v3+a3*dt;
    
co_x1=[co_x1 s1（1）];
co_y1=[co_y1 s1（2）];
co_x2=[co_x2 s2（1）];
co_y2=[co_y2 s2（2）];
co_x3=[co_x3 s3（1）];
co_y3=[co_y3 s3（2）];
    
end
figure（2）
axis（[-rangerange-range range]）%确定坐标轴范围
hold on
plot（co_x1co_y1）
plot（co_x2co_y2）
plot（co_x3co_y3）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%环境及端点%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure（‘Renderer‘‘opengl‘）
hold on
fill_x=[-1010]*0.05;
fill_y=[0-101]*0.05;
pen_1 = plot（00‘XDataSource‘‘drew_x1‘‘