Điện tích chuyển động trong chai từ tính

function Magnetic_bottle
% Author: Tran Hai Cat
% Lecturer in Physics, HCM University of Technology and Education
% - Dai hoc Su pham Ky thuat Tp. Ho Chi Minh
% Created: 2019.08.17
clc
clear variables
close all

global mu0 m q

%% CONST
q_element = 1.6e-19;
mu0 = 4*pi*1e-7;

%% FOR USER:
% Initial velocity:
v0x = 1e6;
v0y = 1*v0x;
v0z = 0;
% Initial position:
x0 = 0;
y0 = 0;
z0 = 0;

%% PARAMETR
m = 9.1095e-31;
q = -q_element;

tmax = 65/v0x;
xmin = -2; xmax = 2;
ymin = -1; ymax = 1;
zmin = -1; zmax = 1;

%% CALCULATION
Nx = 20;
Ny = 10;
Nz = 10;
x_field = linspace(xmin,xmax,Nx);
y_field = linspace(ymin,ymax,Ny);
z_field = linspace(zmin,zmax,Nz);
[X_field,Y_field,Z_field] = meshgrid(x_field,y_field,z_field);
[Bx_matrix,By_matrix,Bz_matrix] = magnetic_field(X_field,Y_field,Z_field);

[T,Y] = ode45(@lorentz,linspace(0,tmax,5000),[x0 y0 z0 v0x v0y v0z]);
x = Y(:,1);
y = Y(:,2);
z = Y(:,3);

%% DISPLAY RESULT
figure('name','CALCULATION PROCESS',...
'color','k','numbertitle','off');
set(gcf,'Units','normalized');
set(gcf,'Position',[0 0.05 1 0.85]);
set(gca,'color','k','xcolor','w','ycolor','w','zcolor','w')
axis equal
axis([zmin zmax xmin xmax ymin ymax]);
rotate3d on

xlabel('x, [m]','fontsize',14);
ylabel('y, [m]','fontsize',14);
zlabel('z, [m]','fontsize',14);
view(30,30);
hold on
streakarrow3d(X_field,Y_field,Z_field,Bx_matrix,By_matrix,Bz_matrix,2,1);
% Draw the ring:
R = 0.5; % Ring's radius
phi = linspace(0,2*pi,180);
x_wire = -1.5*ones(size(phi));
y_wire = R*cos(phi);
z_wire = R*sin(phi);
plot3(x_wire,y_wire,z_wire,'linewidth',3);
x_wire = 1.5*ones(size(phi));
y_wire = R*cos(phi);
z_wire = R*sin(phi);
plot3(x_wire,y_wire,z_wire,'linewidth',3);

plot1 = plot3(x(1),y(1),z(1),'ro','markersize',5,'markerfacecolor','r');
plot2 = animatedline('linewidth',1,'color','g');

for i = 1:length(T)
set(plot1,'xdata',x(i),'ydata',y(i),'zdata',z(i));
addpoints(plot2,x(i),y(i),z(i));
drawnow limitrate
pause(0.01)
end

%% FUNCTIONS
function dy = lorentz(t,y)
%% ODE
global m q
[Bx,By,Bz] = magnetic_field(y(1),y(2),y(3));
dy = zeros(6,1);
dy(1) = y(4);
dy(2) = y(5);
dy(3) = y(6);
dy(4) = q*(y(5)*Bz-y(6)*By)/m;
dy(5) = q*(y(6)*Bx-y(4)*Bz)/m;
dy(6) = q*(y(4)*By-y(5)*Bx)/m;

function [Bx,By,Bz] = magnetic_field(X,Y,Z)
%% Calculating magnetic field
global mu0
Bx = zeros(size(X));
By = zeros(size(Y));
Bz = zeros(size(Z));

I = 1e3; % Ampere
R = 0.5; % Ring's radius
phi = linspace(0,2*pi,180);
x_wire = -1.5*ones(size(phi));
y_wire = R*cos(phi);
z_wire = R*sin(phi);

N = length(x_wire)-1;
kB = mu0/4/pi*I;
for i_wire = 1:N
dx_wire = x_wire(i_wire+1)-x_wire(i_wire);
dy_wire = y_wire(i_wire+1)-y_wire(i_wire);
dz_wire = z_wire(i_wire+1)-z_wire(i_wire);
Delta_x = X-x_wire(i_wire);
Delta_y = Y-y_wire(i_wire);
Delta_z = Z-z_wire(i_wire);
r = sqrt(Delta_x.^2+Delta_y.^2+Delta_z.^2);

% x,y,z components of magnetic vector
Bx = Bx+kB*(dy_wire*Delta_z-dz_wire*Delta_y)./r.^3;
By = By+kB*(dz_wire*Delta_x-dx_wire*Delta_z)./r.^3;
Bz = Bz+kB*(dx_wire*Delta_y-dy_wire*Delta_x)./r.^3;
end

phi = linspace(0,2*pi,180);
x_wire = 1.5*ones(size(phi));
y_wire = R*cos(phi);
z_wire = R*sin(phi);

N = length(x_wire)-1;
kB = mu0/4/pi*I;
for i_wire = 1:N
dx_wire = x_wire(i_wire+1)-x_wire(i_wire);
dy_wire = y_wire(i_wire+1)-y_wire(i_wire);
dz_wire = z_wire(i_wire+1)-z_wire(i_wire);
Delta_x = X-x_wire(i_wire);
Delta_y = Y-y_wire(i_wire);
Delta_z = Z-z_wire(i_wire);
r = sqrt(Delta_x.^2+Delta_y.^2+Delta_z.^2);

% x,y,z components of magnetic vector
Bx = Bx+kB*(dy_wire*Delta_z-dz_wire*Delta_y)./r.^3;
By = By+kB*(dz_wire*Delta_x-dx_wire*Delta_z)./r.^3;
Bz = Bz+kB*(dx_wire*Delta_y-dy_wire*Delta_x)./r.^3;
end

function hh=streakarrow3d(X,Y,Z,U,V,W,np,AC)
%% Bertrand Dano 03-09
% Copyright 1984-2009 The MathWorks, Inc.
DX=abs(X(1,1,1)-X(1,2,1)); DY=abs(Y(1,1,1)-Y(2,1,1)); DZ=abs(Z(1,1,1)-Z(1,1,2));
DD=min([DX DY DZ]);
ks=DD/100; % Size of the "dot" for the tuft graphs
np=np*10;
alpha = 3; % Size of arrow head relative to the length of the vector
beta = .15; % Width of the base of the arrow head relative to the length
XYZ=stream3(X,Y,Z,U,V,W,X,Y,Z);
Vmag=sqrt(U.^2+V.^2+W.^2);
Vmin=min(Vmag(:));Vmax=max(Vmag(:));
Vmag=Vmag(:);
cmap=colormap;
for k=1:length(XYZ)
F=XYZ(k); [L M]=size(F{1});
if L<np
F0{1}=F{1}(1:L,:);
if L==1
F1{1}=F{1}(L,:);
else
F1{1}=F{1}(L-1:L,:);
end

else
F0{1}=F{1}(1:np,:);
F1{1}=F{1}(np-1:np,:);
end
P=F1{1};
if AC==1
vcol=floor((Vmag(k)-Vmin)./(Vmax-Vmin)*64); if vcol==0; vcol=1; end
COL=[cmap(vcol,1) cmap(vcol,2) cmap(vcol,3)];
else
COL='k';
end
hh=streamline(F0);
set(hh,'color',COL,'linewidth',.5);

if L>1
x1=P(1,1); y1=P(1,2); z1=P(1,3);
x2=P(2,1); y2=P(2,2); z2=P(2,3);
u=x2-x1; v=y2-y1; w=z2-z1; uv=sqrt(u.*u + v.*v);

xa1=x2+u-alpha*(u+beta*(v+eps)); xa2=x2+u-alpha*(u-beta*(v+eps)); xa3=x2+u-alpha*u;
ya1=y2+v-alpha*(v-beta.*(u+eps)); ya2=y2+v-alpha*(v+beta.*(u+eps)); ya3=y2+v-alpha*v;
za1=z2+w-alpha*w; za2=z2+w-alpha*(w+beta.*(uv+eps)); za3=z2+w-alpha*(w-beta.*(uv+eps));

plot3([x2 xa1 xa2 x2 xa3 xa3 x2],[y2 ya1 ya2 y2 ya3 ya3 y2],[z2 za1 za1 z2 za2 za3 z2],'color',COL); hold on
end

end
axis image