You are here: Home Education Computational Plasma Physics GC_TrajMagnDipole

GC_TrajMagnDipole

GC_TrajMagnDipole

GC_TrajMagnDipole.m — Objective-C source code, 1 kB (1232 bytes)

File contents

close all;
clear all;
global B0; global q; global m; global Re; global v_mod; global v_par0;
% 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 % replace m_pr with m_el for electron
q = e
% 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
x0 = 4*Re; 
y0 = 0; 
z0 = 0;

pitch_angle = 30.0 %  initial angle between velocity and mag.field (degrees)
% initial parallel velocity
v_par0 = v_mod*cos(pitch_angle*pi/180);

tfin = 50.0 % in s
time = 0:0.01:tfin;

[x_sol,t] = lsode("Newton_LorenzGC",[x0 ;y0; z0; v_par0],time); % solve GC model *** Octave ***
% [t,x_sol] = ode45(@Newton_LorenzGC,time,[x0 ;y0; z0; v_par0]) % use this for *** Matlab ****
plot3(x_sol(:,1)/Re,x_sol(:,2)/Re,x_sol(:,3)/Re);
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;