This page was generated by nbsphinx from docs/notebooks/simulation/particle_tracker.ipynb.
Interactive online version: .

%matplotlib inline

Particle Tracker

An example of PlasmaPy’s particle tracker class.

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

from plasmapy.formulary import ExB_drift, gyrofrequency
from plasmapy.particles import Particle
from plasmapy.plasma.grids import CartesianGrid
from plasmapy.simulation.particle_tracker import (
    IntervalSaveRoutine,
    ParticleTracker,
    TimeElapsedTerminationCondition,
)

Take a look at the docs to gyrofrequency() and ParticleTracker for more information

Initialize a :class:~plasmapy.plasma.grids.CartesianGrid object. This will be the source of electric and magnetic fields for our particles to move in.

grid_length = 10
grid = CartesianGrid(-1 * u.m, 1 * u.m, num=grid_length)

Initialize the fields. We’ll take B in the x direction and E in the y direction, which gets us an E cross B drift in the negative z direction.

Bx_fill = 4 * u.T
Bx = np.full(grid.shape, Bx_fill.value) * u.T

Ey_fill = 2 * u.V / u.m
Ey = np.full(grid.shape, Ey_fill.value) * u.V / u.m

grid.add_quantities(B_x=Bx, E_y=Ey)
ExB_drift(np.asarray([0, Ey_fill.value, 0]) * u.V / u.m, np.asarray([Bx_fill.value, 0, 0]) * u.T)

$[0,~0,~-0.5] \; \mathrm{\frac{m}{s}}$

|ParticleTracker| takes arrays of particle positions and velocities of the shape [nparticles, 3], so these arrays represent one particle starting at the origin.

x0 = [[0, 0, 0]] * u.m
v0 = [[1, 0, 0]] * u.m / u.s
particle = Particle("p+")

Initialize our stop condition and save routine. We can determine a relevant duration for the experiment by calculating the gyroperiod for the particle.

particle_gyroperiod = 1 / gyrofrequency(Bx_fill, particle).to(u.Hz, equivalencies=u.dimensionless_angles())

simulation_duration = 100 * particle_gyroperiod
save_interval = particle_gyroperiod / 10

termination_condition = TimeElapsedTerminationCondition(simulation_duration)
save_routine = IntervalSaveRoutine(save_interval)

Initialize the trajectory calculation.

simulation = ParticleTracker(grid, save_routine=save_routine, termination_condition=termination_condition, verbose=False)

/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/latest/lib/python3.12/site-packages/plasmapy/simulation/particle_tracker.py:633: RuntimeWarning: Quantities should go to zero at edges of grid to avoid non-physical effects, but a value of 2.00E+00 V / m was found on the edge of the E_y array. Consider applying a envelope function to force the quantities at the edge to go to zero.
  warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/latest/lib/python3.12/site-packages/plasmapy/simulation/particle_tracker.py:633: RuntimeWarning: Quantities should go to zero at edges of grid to avoid non-physical effects, but a value of 4.00E+00 T was found on the edge of the B_x array. Consider applying a envelope function to force the quantities at the edge to go to zero.
  warnings.warn(

We still have to initialize the particle’s velocity. We’ll limit ourselves to one in the x direction, parallel to the magnetic field B - that way, it won’t turn in the z direction.

simulation.load_particles(x0, v0, particle)

Run the simulation.

simulation.run()

We can take a look at the trajectory in the z direction:

t, x, v = save_routine.results()
particle_trajectory = x[:, 0]
particle_position_z = particle_trajectory[:, 2]

plt.scatter(t, particle_position_z)

<matplotlib.collections.PathCollection at 0x7f92a7079d00>

or plot the shape of the trajectory in 3D:

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

ax.plot(*particle_trajectory.T)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")

Text(0.5, 0, 'Z')

As a test, we calculate the mean velocity in the z direction from the velocity:

v_mean = v[:, :, 2].mean()
print(
    f"The calculated drift velocity is {v_mean:.4f} to compare with the "
    f"expected E0/B0 = {-(Ey_fill/Bx_fill).value:.4f}"
)

The calculated drift velocity is -0.5076 m / s to compare with the expected E0/B0 = -0.5000