function f = earthMagneticForce(q,r,v)

% Lorentz force due to Earth's magnetic field
% Input:
%   q      : particle charge (C)
%   r(3,1) : particle position vector (m)
%   v(3,1) : particle velocity (m/s)
% Output:
%   f(3,1) : particle magnetic force (N)
% J.M.Soler, Feb.2017

  % Set parameters
  c = 8.06e15;     % Earth's magnetic dipole moment, times mu0/(4*pi) (T*m^3)

  % Find magnetic field, assuming dipole model
  mu = [0;0;1];    % unit vector along Earth's axis
  ru = r/norm(r);  % unit vector along r
  B = c*(mu-(mu'*ru)*ru)/norm(r)^3;
  
  % Find Lorentz force
  f = q*cross(v,B);
  
end

