

Code chương trình Matlab
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 | function charged_particle_UNIFORM_mag_field % 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.16 clc clear variables close all global m q Bx By Bz %% CONST q_element = 1.6e-19; %% PARAMETR B0 = 840e-6; % [Tesla] Bx = B0; By = 0; Bz = 0; m = 9.1095e-31; q = -q_element; v0 = 1.11e7; % [m/s] % Initial position: x0 = 0; y0 = 0; z0 = 0; v0x = v0/10; v0y = 0; v0z = v0; dt = 7e-11; % [second] tmax = 5e3*dt; xmin = 0; xmax = 50e-2; ymin = -16e-2; ymax = 1e-2; zmin = -10e-2; zmax = 10e-2; %% 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 = B0*ones(size(X_field)); By_matrix = zeros(size(Y_field)); Bz_matrix = zeros(size(Z_field)); [T,Y] = ode45(@lorentz,linspace(0,tmax,500),[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([xmin xmax ymin ymax zmin zmax]); rotate3d on xlabel('x, [m]','fontsize',14); ylabel('y, [m]','fontsize',14); zlabel('z, [m]','fontsize',14); view(30,30); hold on quiver3(X_field,Y_field,Z_field,Bx_matrix,By_matrix,Bz_matrix,... 'AutoScaleFactor',1,'color','g'); plot1 = plot3(x(1),y(1),z(1),'ro','markersize',5,'markerfacecolor','r'); plot2 = animatedline('linewidth',2,'color','b'); 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 function dy = lorentz(t,y) %% ODE global m q Bx By Bz 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; |