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

function- GC trajectory in magnetic dipole

function- GC trajectory in magnetic dipole

Newton_LorenzGC.m — Objective-C source code, 1 kB (1557 bytes)

File contents

function x_res   = Newton_LorenzGC(x_vect,t) % use this definition in *** Octave ***
% function x_res = Newton_LorenzGC(t, x_vect)% use this in *** Matlab ***

global B0; global q; global m; global Re; global v_mod; global v_par0;

x_res = zeros(4,1);

vsq = v_mod^2; % This remains unchanged in present of a static B field


fac1 =  -B0*Re^3/(x_vect(1)^2 + x_vect(2)^2 + x_vect(3)^2)^2.5;
Bx = 3*x_vect(1)*x_vect(3)*fac1;
By = 3*x_vect(2)*x_vect(3)*fac1;
Bz = (2*x_vect(3)^2 -x_vect(1)^2- x_vect(2)^2)*fac1;

B_mod = sqrt(Bx*Bx + By*By + Bz*Bz);

mu = m*(vsq-x_vect(4)^2)/(2*B_mod); % magnetic moment is an adiabatic invariant

d = 0.001*Re;
% calculate gradBx
gradB_x =  (getBmod((x_vect(1)+d),x_vect(2),x_vect(3)) - getBmod((x_vect(1)-d),x_vect(2),x_vect(3)))/(2*d);
gradB_y =  (getBmod(x_vect(1),(x_vect(2)+d),x_vect(3)) - getBmod(x_vect(1),(x_vect(2)-d),x_vect(3)))/(2*d);
gradB_z =  (getBmod(x_vect(1),x_vect(2),(x_vect(3)+d)) - getBmod(x_vect(1),x_vect(2),(x_vect(3)-d)))/(2*d);
% b unit vector
b_unit_x = Bx/B_mod;
b_unit_y = By/B_mod;
b_unit_z = Bz/B_mod;
% b unit vector cross gradB
bxgB_x = b_unit_y*gradB_z - b_unit_z*gradB_y;
bxgB_y = b_unit_z*gradB_x - b_unit_x*gradB_z;
bxgB_z = b_unit_x*gradB_y - b_unit_y*gradB_x;
% b unit vector inner product gradB
dotpr = b_unit_x*gradB_x +  b_unit_y*gradB_y + b_unit_z*gradB_z;

fac = m/(2*q*B_mod^2)*(vsq + x_vect(4)^2);

x_res(1) = fac*bxgB_x + x_vect(4)*b_unit_x;
x_res(2) = fac*bxgB_y + x_vect(4)*b_unit_y;
x_res(3) = fac*bxgB_z + x_vect(4)*b_unit_z;
x_res(4) = -mu/m*dotpr;

endfunction