module icon

Moving Charge/Cyclotron

This simulation shows how a charged particle moves in a constant magnetic field. The equation describing the motion of such a particle is the Lorentz equation. For a constant magnetic field and zero electric field, the Lorentz equation can be written \[m\frac{\text{d}}{\text{d}t} \mathbf{v}(t) = q\mathbf{v}(t)\mathbf{B}\] where \(\mathbf{v}\) is the velocity of the particle, \(q\) is the charge of particle, \(m\) is the mass of the particle and \(\mathbf{B}\) is the magnetic field. The calculation is non-relativistic.

When the particle is moving perpendicular to the magnetic field the solution to the equation of motion is a circular orbit of the charged particle. The frequency of this orbit is given by \[f = \frac{qB}{2\pi m}.\] This relationship is utilised within mass spectrometers to identify particles. The simulation below illustrates this motion.

This principle is also the basis for motion in a Cyclotron which is a device used to accelerate charges to high speeds. An oscillating electric field between two "Dees" is used to accelerate the charges while the magnetic field controls the direction of motion. Activating the maximum electric potential listed in the simulation allows modelling of this effect. The electric field oscillates sinusoidally at the orbit frequency to ensure that the particle's energy increases each time it crosses the gap between the dees. While in the gap, the equation of motion becomes \[m\frac{\text{d}}{\text{d}t} \mathbf{v}(t) = q\mathbf{v}(t)\mathbf{B} + q\mathbf{E}(t)\] where \(\mathbf{E}\) is the electric field at the location of the charge when between the dees.

Adjust the initial speed, charge, mass, magnetic field strength, and electric potential using the sliders or test boxes below the simulation. Keep in mind that each variable has a maximum and minimum value. Also shown is the time elapsed, the position of the last click, the cyclotron frequency, the particle speed (as a percentage of the speed of light) and the number of rotations. Click in a few places on the grid, the position appears in the `x' and `y' boxes. After the simulation has started, the settings are fixed. The simulation must be reset to change the settings. The reset button makes the simulation return to \(t=0\), but does not reset the variable settings.

Code samples of circular motion in a magnetic field are also provided in Matlab, Python and C++/Gnuplot. These solve the ODE for motion in a constant magnetic field using a first order method. These code samples can be modified to include more complex field configurations. Note that the first order method used to solve the ODE in the code samples is not numerically stable (this will lead to non-circular orbits for long times and other non-physical results).

Charge (proton):
Mass (proton):
Init. Speed (% of c):
Magnetic Field (T):
D-Potential (kV):
Time (μs):
x (m):
y (m):
Freq. (MHz):
Speed (% of c):
Rotations:

Open this simulation in a separate window.

Tim McIntyre, Elise Kenny, Isaac Lenton
Updated 2023

Code Samples

Matlab
xx = [0.0, 0.0];                    % Initial position (m)
vv = [5.0, 0.0]*1.0e7;              % Initial velocity (m/s)
charge = 1.60217662e-19;            % Charge of electron (C)
field = 5.0e-3;                     % Field strength (T)
mass = 9.10938356e-31;              % Mass of electron (kg)
dt = 1.0e-12;                       % Time step size (s)
outputRate = 200;                   % Number of evolutions between outputs

% Arrays to track the particle position
x = xx(1);
y = xx(2);

figure;
h = plot(x, y);
d = 0.03;   % circle diameter
pxc = @(xx) xx(1) - d/2;
pyc = @(xx) xx(2) - d/2;
ch = rectangle('Position',[pxc(xx) pyc(xx) d d],'Curvature',[1,1]);
grid on;
axis square;
axis([-0.3 0.3 -0.3 0.3]);

% Evolve the system and show the animation
while ishandle(h)
    % Evolve position (1st order, inacurate for long times)
    for n = 1:outputRate
        vv = vv + charge/mass*field*[vv(2), -vv(1)]*dt;
        xx = xx + vv*dt;
    end
    
    % Update the graph
    x = [x xx(1)]; set(h, 'xdata', x);
    y = [y xx(2)]; set(h, 'ydata', y);
    set(ch, 'Position', [pxc(xx) pyc(xx) d d]);
    
    pause(0.05);
    drawnow;
end
Python
from numpy import array as npa
import matplotlib.pyplot as plt
import matplotlib.animation as animation

xx = npa([0.0, 0.0]);               # Initial position (m)
vv = npa([5.0, 0.0])*1.0e7;         # Initial velocity (m/s)
charge = 1.60217662e-19;            # Charge of electron (C)
field = 5.0e-3;                     # Field strength (T)
mass = 9.10938356e-31;              # Mass of electron (kg)
dt = 1.0e-12;                       # Time step size (s)
outputRate = 200;                   # Number of evolutions between outputs
outputLength = 400;                 # Number of output frames

# Arrays to track the particle position
x = [xx[0]];
y = [xx[1]];

# Setup plot
fig, ax = plt.subplots()
line, = ax.plot(x, y, lw=1)
circle = plt.Circle((xx[0], xx[1]), 0.015)
ax.add_artist(circle)
ax.grid('on')
ax.axis('equal')
ax.axis([-0.3, 0.3, -0.3, 0.3])

def animate(step_number):

    global xx
    global vv

    # Evolve position (1st order, inacurate for long times)
    for i in xrange(0, outputRate):
        vv = vv + charge/mass*field*npa([vv[1], -vv[0]])*dt;
        xx = xx + vv*dt;

    # Update line/circle
    x.append(xx[0])
    y.append(xx[1])
    line.set_data(x, y)
    circle.center = xx[0], xx[1]

    return line,circle

ani = animation.FuncAnimation(fig, animate, xrange(1, outputLength),
    interval=25, blit=True, repeat=False)
#ani.save('output.mp4')
plt.show()
Cpp Gnuplot
#include <fstream>

const double x0 = 0.0, y0 = 0.0;        // Initial position (m)
const double vx0 = 5.0e7, vy0 = 0.0;    // Initial velocity (m/s)
const double charge = 1.60217662e-19;   // Charge of electron (C)
const double field = 5.0e-3;            // Field strength (T)
const double mass = 9.10938356e-31;     // Mass of electron (kg)
const double dt = 1.0e-12;              // Time step size (s)
const unsigned outputRate = 200;        // Number of evolutions between outputs
const unsigned outputLength = 400;      // Number of output frames

int main(int argc, char** argv) {

  double x = x0, y = y0, vx = vx0, vy = vy0;

  std::ofstream fp("output.dat");
  fp << "# t\tx\ty\n";
  fp << 0.0 << '\t' << x << '\t' << y << '\n';

  for (unsigned i = 0; i < outputLength; ++i) {

    // Evolve position (1st order, inaccurate for long times)
    for (unsigned j = 0; j < outputRate; ++j) {
      double ovx = vx, ovy = vy;
      vx += charge/mass*field*ovy*dt;
      vy -= charge/mass*field*ovx*dt;
      x += vx*dt; y += vy*dt;
    }

    // Write output
    double t = dt*outputRate*(i+1);
    fp << t << '\t' << x << '\t' << y << '\n';
  }

  return 0;
}
set terminal gif animate 10
set output 'output.gif'
set xrange [-0.3:0.3]
set yrange [-0.3:0.3]
set size ratio -1
set grid

radius = 0.015

stats "output.dat" nooutput
do for [i=0:int(STATS_records)-1] {
  plot "output.dat" u 2:3 every ::::i w lines notitle, \
    "output.dat" u 2:3:(radius) every ::i::i with circles notitle
}