This page was generated by
nbsphinx from
docs/notebooks/simulation/particle_tracker.ipynb.

Interactive online version:
.

```
[1]:
```

```
%matplotlib inline
```

# Particle Tracker

An example of PlasmaPy’s particle tracker class.

```
[2]:
```

```
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.

```
[3]:
```

```
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.

```
[4]:
```

```
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)
```

```
[4]:
```

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

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

```
[6]:
```

```
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.

```
[7]:
```

```
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 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(
/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(
```

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.

```
[8]:
```

```
simulation.load_particles(x0, v0, particle)
```

Run the simulation.

```
[9]:
```

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

```
[10]:
```

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

```
[10]:
```

```
<matplotlib.collections.PathCollection at 0x7f29d43c7590>
```

or plot the shape of the trajectory in 3D:

```
[11]:
```

```
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")
```

```
[11]:
```

```
Text(0.5, 0, 'Z')
```

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

```
[12]:
```

```
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
```