[物理应用] 基于 Matlab 图形用户界面的三体运动仿真 [含 Matlab 源代码 第 871 期]。

% %% function varargout = Threebody(varargin) % THREEBODY MATLAB code for Threebody.fig % THREEBODY, by itself, creates a new THREEBODY or raises the existing % singleton*. % % H = THREEBODY returns the handle to a new THREEBODY or the handle to % the existing singleton*. % % THREEBODY('CALLBACK',hObject,eventData,handles,...) calls the local % function named CALLBACK in THREEBODY.M with the given input arguments. % % THREEBODY('Property','Value',...) creates a new THREEBODY or raises the % existing singleton*. Starting from the left, property value pairs are % applied to the GUI before Threebody_OpeningFcn gets called. An % unrecognized property name or invalid value makes property application % pause. All inputs are passed to Threebody_OpeningFcn via varargin. % % *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one % instance to run (singleton)". % % See also: GUIDE, GUIDATA, GUIHANDLES % Edit the above text to modify the response to menu_help Threebody % Last Modified by GUIDE v2.5 28-Oct-2017 15:50:30 % Begin initialization code - DO NOT EDIT gui_Singleton = 1; gui_State = struct('gui_Name', mfilename, ... 'gui_Singleton', gui_Singleton, ... 'gui_OpeningFcn', @Threebody_OpeningFcn, ... 'gui_OutputFcn', @Threebody_OutputFcn, ... 'gui_LayoutFcn', [] , ... 'gui_Callback', []); if nargin && ischar(varargin{1}) gui_State.gui_Callback = str2func(varargin{1}); end if nargout [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); else gui_mainfcn(gui_State, varargin{:}); end % End initialization code - DO NOT EDIT %% % --- Executes just before Threebody is made visible. function Threebody_OpeningFcn(hObject, eventdata, handles, varargin) %% % This function has no output args, see OutputFcn. % hObject handle to figure % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % varargin command line arguments to Threebody (see VARARGIN) %初始化天体的质量坐标和速度 handles.Star=Model_Init_1(); handles.type=1; handles.pretype=handles.type; %运行标志位，判断是否在正在运行 global Is_Running; %初始化标志位，判断是否在已初始化 global Is_Init; %默认状态为已初始化且正在运行 Is_Running=true; Is_Init=true; handles.output = hObject; % Update handles structure guidata(hObject, handles); % UIWAIT makes Threebody wait for user response (see UIRESUME) % uiwait(handles.figure1); %% % --- Outputs from this function are returned to the command line. function varargout = Threebody_OutputFcn(hObject, eventdata, handles) %% % varargout cell array for returning output args (see VARARGOUT); % hObject handle to figure % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Get default command line output from handles structure varargout{1} = handles.output; %持续运行run()函数，等待中断 run(); %% %星体运动函数 function run() %% global Is_Running; global Is_Init; %万有引力常数 G=6.67*10^(-11); %以0.00001的时间步长迭代 T=0.00001; %注意死循环 while (1) handles=guidata(gcf); %退出 if(isempty(handles)) close all; return; end %已初始化意为一个阶段内的运算结束，将结果显示 if(Is_Init==true) %获取参数 star=handles.Star; %--------------------一号星体-------------------- M1=star(1).M; X1=star(1).X; Y1=star(1).Y; Z1=star(1).Z; U1=star(1).U; V1=star(1).V; W1=star(1).W; %--------------------二号星体-------------------- M2=star(2).M; X2=star(2).X; Y2=star(2).Y; Z2=star(2).Z; U2=star(2).U; V2=star(2).V; W2=star(2).W; %--------------------三号星体-------------------- M3=star(3).M; X3=star(3).X; Y3=star(3).Y; Z3=star(3).Z; U3=star(3).U; V3=star(3).V; W3=star(3).W; %清空 cla; axis([-1000 1010 -1000 1000 -1000 1000]); %不采用默认坐标样式，后面自己画线 axis off; %默认视角下的右手系 set(gca,'XDir','reverse'); set(gca,'YDir','reverse'); %画网格 x=-1000:125:1000; y=-1000:125:1000; for i=1:length(x) line([x(i) x(i)],[1000 -1000],[0 0],'Color',[0.1 0.1 0.1],'LineWidth',0.5); line([1000 -1000],[y(i) y(i)],[0 0],'Color',[0.1 0.1 0.1],'LineWidth',0.5); end %画箭头 line([1000 0],[0 0],[0 0],'Color',[1 0 1],'LineWidth',2); line([0 0],[1000 0],[0 0],'Color',[1 0 1],'LineWidth',2); line([0 0],[0 0],[1000 0],'Color',[1 0 1],'LineWidth',2); %标记坐标轴 text(1050,0,0,'X','Color',[1 0 1],'FontSize',12); text(0,1050,0,'Y','Color',[1 0 1],'FontSize',12); text(0,0,1050,'Z','Color',[1 0 1],'FontSize',12); %标记单位长度 text(125,0,0,'125','Color',[1 0 1],'FontSize',12); text(0,125,0,'125','Color',[1 0 1],'FontSize',12); %初始设置天体颜色、点型、大小等参数，每次循环都重新设置整个画面 if(M1~=0) h=line('Color',[1 0 0],'Marker','.','MarkerSize',25,'erasemode','normal'); line(X1,Y1,Z1,'Color',[1 0 0],'Marker','.','MarkerSize',10); end if(M2~=0) i=line('Color',[0 1 0],'Marker','.','MarkerSize',25,'erasemode','normal'); line(X2,Y2,Z2,'Color',[0 1 0],'Marker','.','MarkerSize',10); end if(M3~=0) l=line('Color',[0 0 1],'Marker','.','MarkerSize',25,'erasemode','normal'); line(X3,Y3,Z3,'Color',[0 0 1],'Marker','.','MarkerSize',10); end %刷新画面 drawnow; Is_Init=false; else [y,fs] = audioread('Windows Background.wav'); %每计算200次重绘一次 for k=1:200; %计算距离 R12=sqrt((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2); R13=sqrt((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2); R23=sqrt((X2-X3)^2+(Y2-Y3)^2+(Z2-Z3)^2); %先判断星球会不会碰撞爆炸，相对距离以8位限 if ((R12<=8)&&(M1~=0)&&(M2~=0)) M1=0; M2=0; delete(h); delete(i); %碰撞音效 sound(y,fs); end if ((R23<=8)&&(M3~=0)&&(M3~=0)) M2=0; M3=0; delete(i); delete(l); sound(y,fs); end if ((R13<=8)&&(M1~=0)&&(M3~=0)) M1=0; M3=0; delete(h); delete(l); sound(y,fs); end %万有引力定律 f12=G*M1*M2/R12^3*[X1-X2,Y1-Y2,Z1-Z2]; f23=G*M2*M3/R23^3*[X2-X3,Y2-Y3,Z2-Z3]; f31=G*M3*M1/R13^3*[X3-X1,Y3-Y1,Z3-Z1]; if(Is_Running==true) t=T; %停止运动 else t=0; end %星体1 if (M1~=0) Ax1=(f31(1)-f12(1))/M1; Ay1=(f31(2)-f12(2))/M1; Az1=(f31(3)-f12(3))/M1; %注意，需先计算坐标，后计算速度 X1=X1+U1*t+1/2*Ax1*t^2; Y1=Y1+V1*t+1/2*Ay1*t^2; Z1=Z1+W1*t+1/2*Az1*t^2; U1=U1+Ax1*t; V1=V1+Ay1*t; W1=W1+Az1*t; end %星体2 if (M2~=0) Ax2=(f12(1)-f23(1))/M2; Ay2=(f12(2)-f23(2))/M2; Az2=(f12(3)-f23(3))/M2; X2=X2+U2*t+1/2*Ax2*t^2; Y2=Y2+V2*t+1/2*Ay2*t^2; Z2=Z2+W2*t+1/2*Az2*t^2; U2=U2+Ax2*t; V2=V2+Ay2*t; W2=W2+Az2*t; end %星体3 if (M3~=0) Ax3=(f23(1)-f31(1))/M3; Ay3=(f23(2)-f31(2))/M3; Az3=(f23(3)-f31(3))/M3; X3=X3+U3*t+1/2*Ax3*t^2; Y3=Y3+V3*t+1/2*Ay3*t^2; Z3=Z3+W3*t+1/2*Az3*t^2; U3=U3+Ax3*t; V3=V3+Ay3*t; W3=W3+Az3*t; end end %重置星体的位置 if(M1~=0) set(h,'XData',X1,'YData',Y1,'ZData',Z1); line(X1,Y1,Z1,'Color',[1 0 0],'Marker','.','MarkerSize',10); end if(M2~=0) set(i,'XData',X2,'YData',Y2,'ZData',Z2); line(X2,Y2,Z2,'Color',[0 1 0],'Marker','.','MarkerSize',10); end if(M3~=0) set(l,'XData',X3,'YData',Y3,'ZData',Z3); line(X3,Y3,Z3,'Color',[0 0 1],'Marker','.','MarkerSize',10); end %显示数据 %显示第一颗星体的数据 set(handles.Text_m1,'String',num2str(M1*10^(-19))); set(handles.Text_x1,'String',num2str(X1)); set(handles.Text_y1,'String',num2str(Y1)); set(handles.Text_z1,'String',num2str(Z1)); set(handles.Text_v1,'String',num2str(V1*10^(-3))); set(handles.Text_u1,'String',num2str(U1*10^(-3))); set(handles.Text_w1,'String',num2str(W1*10^(-3))); %显示第二颗星体的数据 set(handles.Text_m2,'String',num2str(M2*10^(-19))); set(handles.Text_x2,'String',num2str(X2)); set(handles.Text_y2,'String',num2str(Y2)); set(handles.Text_z2,'String',num2str(Z2)); set(handles.Text_u2,'String',num2str(U2*10^(-3))); set(handles.Text_v2,'String',num2str(V2*10^(-3))); set(handles.Text_w2,'String',num2str(W2*10^(-3))); %显示第三颗星体的数据 set(handles.Text_m3,'String',num2str(M3*10^(-19))); set(handles.Text_x3,'String',num2str(X3)); set(handles.Text_y3,'String',num2str(Y3)); set(handles.Text_z3,'String',num2str(Z3)); set(handles.Text_u3,'String',num2str(U3*10^(-3))); set(handles.Text_v3,'String',num2str(V3*10^(-3))); set(handles.Text_w3,'String',num2str(W3*10^(-3))); %显示星体之间的距离 set(handles.value_r12,'String',num2str(R12)); set(handles.value_r23,'String',num2str(R23)); set(handles.value_r13,'String',num2str(R13)); %设置字体颜色 %第一颗星体字体对应为红色 set(handles.Text_m1,'ForegroundColor',[1 0 0]); set(handles.Text_x1,'ForegroundColor',[1 0 0]); set(handles.Text_y1,'ForegroundColor',[1 0 0]); set(handles.Text_z1,'ForegroundColor',[1 0 0]); set(handles.Text_u1,'ForegroundColor',[1 0 0]); set(handles.Text_v1,'ForegroundColor',[1 0 0]); set(handles.Text_w1,'ForegroundColor',[1 0 0]); %第二颗星体字体对应为绿色 set(handles.Text_m2,'ForegroundColor',[0 1 0]); set(handles.Text_x2,'ForegroundColor',[0 1 0]); set(handles.Text_y2,'ForegroundColor',[0 1 0]); set(handles.Text_z2,'ForegroundColor',[0 1 0]); set(handles.Text_u2,'ForegroundColor',[0 1 0]); set(handles.Text_v2,'ForegroundColor',[0 1 0]); set(handles.Text_w2,'ForegroundColor',[0 1 0]);