MATLAB Answers

3D Plot of particle trajectory in E and B field

66 views (last 30 days)
I have the following code that solves Lorents equation using Boris method time integrator and it calculates the gyro radius and frequency asw well. The code works for me now and it provides a 2D plot, however I am trying to plot a 3D plot of the particle trajectory. So an example would be like this:
Do you have sugesstions to do this in MATLAB? Thanks
Code:
%HW1: Problem 2: Boris Solver for 1) E = 0 and B =/ 0 2) E =/0 ande B =/0
%Start setting up constants
dt = 1e-2; % time step
mass = 1.0; % mass of particle
charge = 1.0; % charge of particle
n = 500; % number of time steps
%Initial parameters:
v = [0, 1, 0]; % initial velocity
x = [0, 0, 0]; % initial position
B = [0, 0, 10]; % initial mag. field along z directions
E = [0, 0, 0]; % initial E field, for case 1) E = 0 and B =/ 0
% Compute the cyclotron frequency and radius (Larmor radius)
omega = (charge .*B(3)) ./mass;
rL = v(2) ./omega; % try v(2) as well since initial v is the perp velocity to B
fprintf('Cyclotron frequency is %4.2f \n',omega); %Test
fprintf('Larmor radius is %4.2f \n',rL);
X = zeros(n,3); % initialize an array of zeros with size nx3 for positions
V = zeros(n,3);
for time = 1:1:n
t = (charge./mass) .*(B .*0.5 .*dt); %compute the t vector
s = (2 .*t) ./(1+ t.*t); %compute the s vector
% start the Boris algorithm:
v_minus = v + (charge./mass) .*(E .*0.5 .*dt); % First half E acceleration: generally for E =/0 ande B =/0
v_prime = v_minus + cross(v_minus,t); % First half B rotation
v_plus = v_minus + cross(v_prime,s); % Second hald B rotation
v = v_plus + (charge./mass) .*(E .*0.5 .*dt); % Second half E acceleration
% Finally update position
x = x + v .*dt;
X(time,:) = x;
V(time,:) = v;
end
% plotting
figure;
plot(X(:,1),X(:,2),'k','Linewidth',2); hold on;
set(gca,'TickLabelInterpreter','Latex','Fontsize',14)
ylabel('$ y / d_{\rm p} $','Interpreter','Latex','Fontsize',14);
xlabel('$ x / d_{\rm p} $','Interpreter','Latex','Fontsize',14);

Accepted Answer

David Goodmanson
David Goodmanson on 3 Feb 2021
Edited: David Goodmanson on 3 Feb 2021
Hi Lujain,
After setting v = [0 1 1] to provide some initial z velocity,
figure(2)
plot3(X(:,1),X(:,2),X(:,3))
grid on
and note that if you put the mouse cursor inside the window and click on the symbol with the circle enclosing the cube, you can rotate the plot with a click-hold mouse.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!