matlab模拟带电粒子在恒定电磁场中运动 (Boris算法)

激情三角

编辑于 2023年12月01日 10:21

收录于文集

共7篇

matlab 代码:

 代码块
Python
自动换行
复制代码
%% parameters 
q=1; m=1; g=0;%比荷 重力加速度
E=[	1 1 0 ]; %Ex; Ey; Ez;
B=[2 0 0 ]; %Bx; By; Bz; 
p0=[0 0 0 ] ; %x0=0; y0=0; z0=0; 粒子位置
v0=[0 1 0 ]; %vx0=0; vy0=1; vz0=0;
t=10*pi;dt=0.02;
colordef black;
figure(&amp;#39;name&amp;#39;,&amp;#39;move&amp;#39;); 
colordef black;
grid on;
plane=scatter3(p0(1),p0(2),p0(3),&amp;#39;c&amp;#39;,&amp;#39;o&amp;#39;);
planev=text(p0(1)+0.4,p0(2)+0.4,p0(3)+0.31,[&amp;#39;v=&amp;#39;,&amp;#39;(&amp;#39;,num2str(single(v0)),&amp;#39;)&amp;#39;]);
planep=text(p0(1)+0.4,p0(2)+0.4,p0(3)+0.205,[&amp;#39;p=&amp;#39;,&amp;#39;(&amp;#39;,num2str(single(p0)),&amp;#39;)&amp;#39;]);
planet=text(p0(1)+0.4,p0(2)+0.4,p0(3)+0.4,[&amp;#39;t=&amp;#39;,&amp;#39;0&amp;#39;]);
h=animatedline(&amp;#39;color&amp;#39;,&amp;#39;b&amp;#39;);
grid on
hold on 
view(3) 
xlabel(&amp;#39;x&amp;#39;) 
ylabel(&amp;#39;y&amp;#39;) 
zlabel(&amp;#39;z&amp;#39;) 
hold on
tic
%% motion equation 
for t=dt:t/dt+1
[v]=loca(p0,v0,E,B,g,dt,q/m);
p=p0+v*dt;
addpoints(h,p(1),p(2),p(3));
set(plane,&amp;#39;Xdata&amp;#39;,p(1),&amp;#39;Ydata&amp;#39;,p(2),&amp;#39;Zdata&amp;#39;,p(3));
planep.String=[&amp;#39;p=&amp;#39;,&amp;#39;(&amp;#39;,num2str(single(p)),&amp;#39;)&amp;#39;];
planev.String=[&amp;#39;v=&amp;#39;,&amp;#39;(&amp;#39;,num2str(single(v)),&amp;#39;)&amp;#39;];
planet.String=[&amp;#39;t=&amp;#39;,num2str(t*dt)];
drawnow
%pause(0.005);
v0=v;
p0=p;
end
% drawnow
toc
%%% 
function [v]=loca(p,v,E,B,g,dt,qm)
%此处加速度于位置p无关
A=[0,B(3),-B(2);-B(3),0,B(1);B(2),-B(1),0];
C=(qm*E*dt/2)&amp;#39;+[0,0,-g*dt/2]&amp;#39;;
P=(eye(3)+qm*dt*A/2)/(eye(3)-qm*dt*A/2);
v=P*(v&amp;#39;+C)+C;%速度迭代
v=v&amp;#39;;%转置
end复制成功