三体系统在北太天元上的模拟

（文章灵感来自卢朓老师的B站视频）最近电视剧《三体》的大热，引起了大家对三体系统的注意力，今天就让我们在北太天元上面模拟一下三体系统的运动轨迹首先，什么是三体系统呢？

三体（three-body problem）
天体力学中的基本力学模型。研究三个可视为质点的天体在相互之间万有引力作用下的运动规律问题。
这三个天体的质量、初始位置和初始速度都是任意的。
----------------------------摘自百度百科

简单的说，就是我们需要模拟三个恒星组成的系统的三体运动。详细代码如下，代码摘自互联网

%模拟三个恒星组成的系统的三体运动
clear
load_plugin("time"); 
%为了使用北太天元软件的pause插件函数
close all
% 三个恒星的质量都是1
ms = 1 ;
mt = 1 ;
mj = 1 ;
% 无量纲后万有引力常数设置为1
G = 1 ;
%初始条件 [xs,ys,xt,yt,xj,yj,vxs,vys,vxt,vyt,vjx,vjt]
CI = [0 -0.1 2 2 5 0 0 0 0 0 0 0];
%初始时刻to = 0;
%计算终止时刻t
f = 120; 
%由位置的导数速度，速度的导数是加速，牛顿第二定律
% 以及万有引力定律得到常微分方程组
fxy = @(ps, pt, pj,ms,mt,mj) ...G*( mt.*(pt-ps)./norm(pt-ps).^3 ...+ mj.*(pj-ps)./norm(pj-ps).^3 );
F = @(t,Y) [Y(7);Y(8);Y(9);Y(10);Y(11);Y(12);...
      fxy(Y([1,2]),Y([3,4]),Y([5,6]),ms,mt,mj); ...
      fxy(Y([3,4]),Y([1,2]),Y([5,6]),mt,ms,mj); ...
      fxy(Y([5,6]),Y([3,4]),Y([1,2]),mj,mt,ms); ...];
%使用ode45求解常微分方程组的初值问题
[t,Y]=ode45(F,[to,tf],CI); 
%plot(Y(:,1),Y(:,2),'r',Y(:,3),Y(:,4),'g',Y(:,5),Y(:,6),'b')
yo = Y(1) ;
dto = 0.3 ;
plotmax = 100 ;
T=to ; 
xmin = min(min(Y(:,[1,3,5]))); 
%三个质点的x坐标（在所有时刻)的最小值
xmax = max(max(Y(:,[1,3,5])));
ymin = min(min(Y(:,[2,4,6]))); 
%三个质点的y坐标（在所有时刻)的最小值
ymax = max(max(Y(:,[2,4,6]))); 
clf
close all
figure('Position',[0 0 1550 800])
hold off
told = 0;
for i = 1:length(Y(:,1))
    dt = abs(Y(i,1)-yo)/abs(Y(i,7));
    if dt >= dto
        if i>plotmax
            shift = plotmax;
        else
            shift = i-1;
        end
        plot(...
            [xmin,xmax],[ymin,ymax], 'w', ... %画一个白色的斜线代替axis([xmin,xmax,ymin,ymax])设置画图范围
            Y(i-shift:i,1),Y(i-shift:i,2),'r','LineWidth',2, ... %画第一个恒星在i-shift个时刻和第i个时刻件的轨迹
            Y(i,1),Y(i,2),'-or','LineWidth',4, ... %画第一个恒星在第i个时刻所在的位置
            Y(i-shift:i,3),Y(i-shift:i,4),'g','LineWidth',2, ...
            Y(i,3),Y(i,4),'-og','LineWidth',4, ...
            Y(i-shift:i,5),Y(i-shift:i,6),'b','LineWidth',2, ...
            Y(i,5),Y(i,6),'-ob','LineWidth',4)
        title(sprintf('时间=%f',t(i)))
        T=[T;t(i)];
        yo = Y(i,1) ;
        vo = Y(i,7) ;
        end
 pause(0.01)
 end
 X=[0:1:length(T)-1];
 figure(2)
 plot(X,T)
 plot(Y(:,1),Y(:,2),'r', 'LineWidth',2, ...
 Y(:,3),Y(:,4),'g','LineWidth',2, ...
 Y(:,5),Y(:,6),'b', 'LineWidth',2)
 unload_plugin("time")

【春节篇】使用北太天元画灯笼

佳节临近，今天给大家分享一下如何使用北太天元绘制一个灯笼灯笼中有很多大面积的颜色区域，是需要使用北太天元的 fill 函数来完成，fill 函数的帮助如下：

>> help fill
填充的二维多边形
语法:    fill(X,Y)
示例:    
    % 函数创建红色多边形。    
    % X 是顶点的x坐标，Y是顶点的y坐标    
    % 例如 一个三角形的三个顶点的坐标是p1(0,0), p2(1,0), p3(0.5,0.5)    
    % 画出红色的三角形    
    x = [0, 1, 0.5];    
    y = [0, 0, 0.5];    
    fill(x,y,'r')

为了调用 fill 函数，我们需要先创建椭圆数据点生成函数，用来生成椭圆和圆形的边界的 X 和 Y 坐标

% 椭圆数据点生成函数
function [X,Y]=getEllipse(Mu,XR,YR,theta,pntNum)
  % Mu      |  中心点
  % XR,YR   |  旋转前X,Y半轴长度
  % theta   |  旋转角度
  % pntNum |  生成数据点个数
  tList = linspace(0,2*pi,pntNum);
  X = cos(tList).*XR;
  Y = sin(tList).*YR;
  rotateMat = [cos(theta),-sin(theta);sin(theta),cos(theta)];
  XY = rotateMat*[X;Y]+Mu(:);
  X = XY(1,:);
  Y = XY(2,:);
end

下面我们就可以开始灯笼的绘制了

clf
hold on
n = 0.4;
x = 1;
a = 0.1;
w = 2.7;h = 2.5; %灯笼主体椭圆的宽和高

% 绘制灯笼主体
[X0,Y0]=getEllipse([0,0],2.9,h,0,200);
fill(X0,Y0,[184,20,25]./255,'EdgeColor',[153,12,40]./255,'LineWidth',1.2)

% 绘制辐线
[X,Y]=getEllipse([0,0],w,h,0,200);plot(X,Y,'Color',[236,136,74]./255,'LineWidth',1.2)
[X,Y]=getEllipse([0,0],w-n,h,0,200);plot(X,Y,'Color',[236,136,74]./255,'LineWidth',1.2)
[X,Y]=getEllipse([0,0],w-2*n,h,0,200);plot(X,Y,'Color',[236,136,74]./255,'LineWidth',1.2)
[X,Y]=getEllipse([0,0],w-3*n,h,0,200);plot(X,Y,'Color',[236,136,74]./255,'LineWidth',1.2)
[X,Y]=getEllipse([0,0],w-4*n+0.05,h,0,200);plot(X,Y,'Color',[236,136,74]./255,'LineWidth',1.2)
[X,Y]=getEllipse([0,0],w-5*n+0.1,h,0,200);plot(X,Y,'Color',[236,136,74]./255,'LineWidth',1.2)
[X,Y]=getEllipse([0,0],w-6*n+0.1,h,0,200);plot(X,Y,'Color',[236,136,74]./255,'LineWidth',1.2)
plot([0,0],[-h,h],'Color',[236,136,74]./255,'LineWidth',1.2)

% 计算其它部件需要的边界坐标
X1 = X0(abs(X0)<=x);
X2 = X1(1:size(X1,2)/2);
X3 = X1(size(X1,2)/2+1:end);
Y1 = Y0(abs(X0)<=x);
Y2 = Y1(1:size(Y1,2)/2);
Y3 = Y1(size(Y1,2)/2+1:end);
XX1 = X0(abs(X0)<=a);
XX2 = XX1(1:size(XX1,2)/2);
XX3 = XX1(size(XX1,2)/2+1:end);	
YY1 = ones(1,size(XX2,2)).*(-w-0.2);
YY2 = ones(1,size(XX2,2)).*(-h-2.5);
YY3 = ones(1,size(XX2,2)).*(-h-2);


% 绘制其它部分
Y4 = ones(1,size(Y2,2)).*(w+0.2);
fill(X1,[Y2,Y4],[76,24,38]./255,'EdgeColor',[236,136,74]./255,'LineWidth',1.2)

Y5 = ones(1,size(Y2,2)).*-(w+0.2);
fill(X1,[Y3,Y5],[76,24,38]./255,'EdgeColor',[236,136,74]./255,'LineWidth',1.2)

fill(XX1,[YY1,YY3],[76,24,38]./255,'EdgeColor',[236,136,74]./255,'LineWidth',1.2)
fill(XX1.+0.7,[YY1,YY3],[76,24,38]./255,'EdgeColor',[236,136,74]./255,'LineWidth',1.2)
fill(XX1.-0.7,[YY1,YY3],[76,24,38]./255,'EdgeColor',[236,136,74]./255,'LineWidth',1.2)
fill(XX1.+0.35,[YY1,YY3],[76,24,38]./255,'EdgeColor',[236,136,74]./255,'LineWidth',1.2)
fill(XX1.-0.35,[YY1,YY3],[76,24,38]./255,'EdgeColor',[236,136,74]./255,'LineWidth',1.2)
fill(XX1,[-YY1,-YY2],[76,24,38]./255,'EdgeColor',[236,136,74]./255,'LineWidth',1.2)

[X6,Y6]=getEllipse([0,h+1],0.2,0.2,0,200);
fill(X6,Y6,[184,20,25]./255,'EdgeColor')

% 调整坐标轴
axis([-5 5 -5 5])
hold off
title("新春快乐！")