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

[1]:
%matplotlib inline

Particle stepper

An example of PlasmaPy’s particle stepper class, currently in need of a rewrite for speed.

[2]:
import numpy as np

from astropy import units as u

from plasmapy.formulary import ExB_drift, gyrofrequency
from plasmapy.plasma import Plasma
from plasmapy.simulation import ParticleTracker

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

Initialize a plasma. This will be a source of electric and magnetic fields for our particles to move in.

[3]:
plasma = Plasma(
    domain_x=np.linspace(-1, 1, 10) * u.m,
    domain_y=np.linspace(-1, 1, 10) * u.m,
    domain_z=np.linspace(-1, 1, 10) * u.m,
)

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]:
B0 = 4 * u.T
plasma.magnetic_field[0, ...] = B0

E0 = 2 * u.V / u.m
plasma.electric_field[1, ...] = E0

ExB_drift(plasma.electric_field[:, 0, 0, 0], plasma.magnetic_field[:, 0, 0, 0])
[4]:
$[0,~0,~-0.5] \; \mathrm{\frac{m}{s}}$

Calculate the timestep. We’ll take one proton p, take its gyrofrequency, invert that to get to the gyroperiod, and resolve that into 10 steps for higher accuracy.

[5]:
freq = gyrofrequency(B0, "p").to(u.Hz, equivalencies=u.dimensionless_angles())
gyroperiod = (1 / freq).to(u.s)
steps_to_gyroperiod = 50
timestep = 20 * gyroperiod / steps_to_gyroperiod

Initialize the trajectory calculation.

[6]:
number_steps = steps_to_gyroperiod * int(2 * np.pi)
trajectory = ParticleTracker(plasma, "p", 1, 1, timestep, number_steps)

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.

[7]:
trajectory.v[0][0] = 1 * (u.m / u.s)

Run the pusher and plot the trajectory versus time.

[8]:
trajectory.run()
trajectory.plot_time_trajectories()
../../_images/notebooks_simulation_particle_stepper_15_0.png

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

[9]:
trajectory.plot_time_trajectories("z")
../../_images/notebooks_simulation_particle_stepper_17_0.png

or plot the shape of the trajectory in 3D:

[10]:
trajectory.plot_trajectories()
../../_images/notebooks_simulation_particle_stepper_19_0.png

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

[11]:
vmean = trajectory.velocity_history[:, :, 2].mean()
print(
    f"The calculated drift velocity is {vmean:.4f} to compare with the "
    f"expected E0/B0 = {-E0/B0:.4f}"
)
The calculated drift velocity is -0.5025 m / s to compare with the expected E0/B0 = -0.5000 V / (m T)