You are here: Home Education Computational Plasma Physics main - trajectory in magnetic dipole

main - trajectory in magnetic dipole

main - trajectory in magnetic dipole

TrajMagnDipole.m — Objective-C source code, 1 kB (1477 bytes)

File contents

close all;
clear all;

global B0; global q; global m; global Re; 

% parameters
e = 1.602176565e-19; % Elementary charge (Coulomb)
m_pr = 1.672621777e-27; % Proton mass (kg)
m_el = 9.10938291e-31; % Electron mass (kg)
c = 299792458; % speed of light (m/s)
% Earth parameters
B0 = 3.07e-5 % Tesla
Re = 6378137 % meter (Earth radius)

% use electron/proton
m = m_pr % we choose proton, change to m_el when simulating electron
q = e % positive charge, change to -e when simulating electron

% Trajectory of a proton with 10MeV kinetic energy in dipole field
K = 1e7    % kinetic energy in eV
K = K*e;   % convert to Joule
% Find corresponding speed:
v_mod = c/sqrt(1+(m*c^2)/K) 
% initial position: equatorial plane 4Re from Earth
x0 = 4*Re; 
y0 = 0; 
z0 = 0;

pitch_angle = 30.0 %  initial angle between velocity and mag.field (degrees)

% initial velocity
u0 = 0.0;
v0 = v_mod*sin(pitch_angle*pi/180);
w0 = v_mod*cos(pitch_angle*pi/180);

tfin = 50 % in s
time = 0:0.01:tfin;
[x_sol,t] = lsode("Newton_Lorenz",[x0 ;y0; z0; u0; v0; w0],time); % solve equation of motion *** Octave ***
% [t,x_sol] = ode23(@Newton_Lorenz,[0 tfin],[x0 ;y0; z0; u0; v0; w0]) % use this for *** Matlab ****

% plots
plot3(x_sol(:,1)/Re,x_sol(:,2)/Re,x_sol(:,3)/Re,'r'); % plot trajectory in 3D
xlabel('x[Re]')
ylabel('y[Re]')
zlabel('z[Re]')
title('Proton 10 MeV starting from x = 4Re, y = 0, z = 0')
hold on;
sphere(40);
hold off;

%plot(time,0.5*m*(x_sol(:,4).^2 + x_sol(:,5).^2 + x_sol(:,6).^2))