带电粒子在磁镜场中的运动

系统：Windows 10 专业版 1903
matlab版本：Windows版 2018a

使用说明

详细的推导过程见这篇文章，程序需要用到这篇文章中的 magnetic.m 程序，请将它与下面的程序放在一个文件夹下。 orbit.m用来模拟氢粒子在磁镜场中的轨迹，可取变量choice为不同值，模拟不同初始位置和速度下的氢离子轨迹，如果需要画二维图，只需要将程序中的相应地方添加注释和取消注释。

程序

orbit.m

tic;clc;clear;
%choice=1,模拟受到反射的情况，2：逃出磁镜 ，3：飘逸小，
%4：产生飘逸，5：飘逸小。y0为初始位置和速度,3、4、5不适合用动图画
choice=5;
switch choice
    case 1
        y0=[0 0.01 0.8 0 sqrt(2)/2*1e6 sqrt(2)/2*1e6];
    case 2
        y0=[0 0.1 0.75 0 sqrt(2)*1e5 sqrt(2)*1e6]; 
    case 3
        y0=[0 0.00075 -1 0 sqrt(2)/2*1e6 sqrt(2)/2*1e5];
    case 4
        y0=[0 0.01 -1 0 sqrt(2)/2*1e6 sqrt(2)/2*1e5];
    case 5
        y0=[0 0.01 0.6 0 sqrt(2)/2*1e6 sqrt(2)/2*1e5];  
end
if choice<=2
    dt=1e-10;tspan=[0:dt:6e-7];
else
    dt=1e-10;tspan=[0:dt:8e-6];
end
[t,y]=ode45(@odefun,tspan,y0);
%画动图
% for i=1:length(y(:,1))
%     plot3(y(i,1),y(i,2),y(i,3),'.');hold on;
%     pause(0.0000001)
% end
%画三维图
plot3(y(:,1),y(:,2),y(:,3));       % static picture
xlabel('x');ylabel('y');zlabel('z');
%画二维图x-z
%plot(y(:,1),y(:,3));
%xlabel(['\fontname{Arial}x']);
%ylabel(['\fontname{Arial}z']);
% plot(y(:,2),y(:,3));
% xlabel(['\fontname{Arial}y']);
% ylabel(['\fontname{Arial}z']);
set(gca,'FontSize',20);hold on;
%画两个线圈
% h=-1;t=0:0.1:(2*pi);t=[t,0];
% plot3(0+0.2*sin(t),0+0.2*cos(t),h*ones(size(t)),'LineWidth',4)
% h=1;t=0:0.1:(2*pi);t=[t,0];
% plot3(0+0.2*sin(t),0+0.2*cos(t),h*ones(size(t)),'LineWidth',4)
% xlabel('x');ylabel('y');zlabel('z');set(gca,'FontSize',20);
% 写入到文件
time=datestr(now,30);
print(gcf,['orbit_',num2str(choice),'_',time,'.eps'],'-depsc','-r600')
savefig(['orbit_',num2str(choice),'_',time,'.fig'])
print(gcf,['orbit_',num2str(choice),'_',time,'.png'],'-dpng')
disp(['耗时', num2str(toc) 's']);

取choice=1时的结果

取choice=2的结果

取choice=3的结果

取choice=4的结果

取choice=5的结果

源码

欢迎转载，不需注明出处，就说是你写的