36 views (last 30 days)

Show older comments

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);

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.

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

Start Hunting!