PlasmaPy Documentation¶
This page was generated by
nbsphinx from
docs/notebooks/diagnostics/proton_radiography_particle_tracing.ipynb.
Interactive online version:
.
Creating Synthetic Proton Radiographs by Particle Tracing¶
Proton radiography is a diagnostic technique often used to interrogate the electric and magnetic fields inside high energy density plasmas. The area of interest is positioned between a bright source of protons and a detector plane. Electric and magnetic fields in the plasma deflect the protons, producing patterns on the detector. Since this represents a nonlinear and lineintegrated measurement of the fields, the interpretation of these “proton radiographs” is complicated.
The SyntheticProtronRadiography class creates a synthetic proton radiographs given a grid of electric and magnetic field (produced either by simulations or analytical models). After the geometry of the problem has been set up, a particle tracing algorithm is run, pushing the protons through the field region. After all of the protons have reached the detector plane, a synthetic proton radiograph is created by making a 2D histogram in that plane.
[1]:
import astropy.constants as const
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from plasmapy.diagnostics import proton_radiography as prad
from plasmapy.plasma.grids import CartesianGrid
To illustrate the use of this package, we’ll first create an example CartesianGrid object and fill it with the analytical electric field produced by a sphere of Gaussian potential.
[2]:
# Create a Cartesian grid
L = 1 * u.mm
grid = CartesianGrid(L, L, num=100)
# Create a spherical potential with a Gaussian radial distribution
radius = np.linalg.norm(grid.grid, axis=3)
arg = (radius / (L / 3)).to(u.dimensionless_unscaled)
potential = 2e5 * np.exp((arg ** 2)) * u.V
# Calculate E from the potential
Ex, Ey, Ez = np.gradient(potential, grid.dax0, grid.dax1, grid.dax2)
Ex = np.where(radius < L / 2, Ex, 0)
Ey = np.where(radius < L / 2, Ey, 0)
Ez = np.where(radius < L / 2, Ez, 0)
# Add those quantities to the grid
grid.add_quantities(E_x=Ex, E_y=Ey, E_z=Ez, phi=potential)
# Plot the Efield
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection="3d")
ax.view_init(30, 30)
# skip some points to make the vector plot intelligable
s = tuple([slice(None, None, 6)] * 3)
ax.quiver(
grid.pts0[s].to(u.mm).value,
grid.pts1[s].to(u.mm).value,
grid.pts2[s].to(u.mm).value,
grid["E_x"][s],
grid["E_y"][s],
grid["E_z"][s],
length=1e6,
)
ax.set_xlabel("X (mm)")
ax.set_ylabel("Y (mm)")
ax.set_zlabel("Z (mm)")
ax.set_xlim(1, 1)
ax.set_ylim(1, 1)
ax.set_zlim(1, 1)
ax.set_title("Gaussian Potential Electric Field")
[2]:
Text(0.5, 0.92, 'Gaussian Potential Electric Field')
Prior to running the particle tracing algorithm, the simulation instance must be instantiated by providing some information about the setup, including the locations of the source and detector relative to the origin of the grid.
The source and detector coordinates are entered as a 3tuple in one of three coordinate systems: Cartesian (\(x\), \(y\), \(z\)), spherical (\(r\), \(\theta\), \(\phi\)) or cylindrical (\(r\), \(\theta\), \(z\)). All values should be astropy.units.Quantity instances with units of either length or angle. The vector from the source to the detector should pass through the origin to maximize the number of particles that pass through the simulated fields.
[3]:
source = (0 * u.mm, 10 * u.mm, 0 * u.mm)
detector = (0 * u.mm, 100 * u.mm, 0 * u.mm)
sim = prad.SyntheticProtonRadiograph(grid, source, detector, verbose=True)
Source: [ 0. 0.01 0. ] m
Detector: [0. 0.1 0. ] m
Magnification: 11.0
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/plasma/grids.py:158: RuntimeWarning: B_x is not specified for the provided grid.This quantity will be assumed to be zero.
RuntimeWarning,
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/plasma/grids.py:158: RuntimeWarning: B_y is not specified for the provided grid.This quantity will be assumed to be zero.
RuntimeWarning,
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/plasma/grids.py:158: RuntimeWarning: B_z is not specified for the provided grid.This quantity will be assumed to be zero.
RuntimeWarning,
Note that, since the example grid did not include a Bfield, the Bfield is assumed to be zero and a warning is printed.
Next, a distribution of nparticles
simulated particles of energy particle_energy
is created using the create_particles() function. Setting the max_theta
parameter eliminates particles with large angles (relative to the sourcedetector axis) which otherwise would likely not hit the detector. Particles with angles less
than \(\theta_{max}\) but greater than \(\theta_{track}\) in the setup figure above will not cross the grid. These particles are retained, but are coasted directly to the detector plane instead of being pushed through the grid.
The particle
keyword sets the type of the particles being traced. The default particle is protons, which is set here explicitly to demonstrate the use of the keyword.
By default, the particle velocities are initialized with random angles (a MonteCarlo approach) with a uniform flux per unit solid angle. However, particles can also be initialized in other ways by setting the distribution
keyword.
[4]:
sim.create_particles(1e5, 3 * u.MeV, max_theta=np.pi / 15 * u.rad, particle="p")
Creating Particles
The simulation is now ready to run. In brief, the steps of the simulation cycle are as follows:
Particles that will never hit the field grid are ignored (until a later step, when they will be automatically advanced to the detector plane).
Particles are advanced to the time when the first particle enters the simulation volume. This is done in one step to save computation time.
While particles are on the grid, the particle pusher advances them each timestep by executing the following steps:
The fields at each particle’s location are interpolated using the interpolators defined in the AbstractGrid subclasses.
The simulation timestep is automatically (and adaptively) calculated based on the proton energy, grid resolution, and field amplitudes. This timestep can be clamped or overridden by setting the
dt
keyword in the run() function.An implementation of the Boris particle push algorithm is used to advance the velocities and positions of the particles in the interpolated fields.
After all of the particles have left the grid, all particles are advanced to the detector plane (again saving time). Particles that are headed away from the detector plane at this point are deleted, as those particles will never be detected.
When the simulation runs, a progress meter will show the number of particles currently on the grid. This bar will start at zero, increase as particles enter the grid, then decrease as they leave it. When almost all particles have left the grid, the simulation ends.
[5]:
sim.run()
Particles on grid: 0% 1.6e+02/5.6e+04 particles
Run completed
Fraction of particles tracked: 55.6%
Fraction of tracked particles that entered the grid: 64.0%
Fraction of tracked particles deflected away from the detector plane: 0.0%
The following plot illustrates that, after the simulation has ended, all particles have been advanced to the detector plane.
[6]:
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection="3d")
ax.view_init(30, 150)
ax.set_xlabel("X (cm)")
ax.set_ylabel("Y (cm)")
ax.set_zlabel("Z (cm)")
# Plot the sourcetodetector axis
ax.quiver(
sim.source[0] * 100,
sim.source[1] * 100,
sim.source[2] * 100,
sim.detector[0] * 100,
sim.detector[1] * 100,
sim.detector[2] * 100,
color="black",
)
# Plot the simulation field grid volume
ax.scatter(0, 0, 0, color="green", marker="s", linewidth=5, label="Simulated Fields")
# Plot the the proton source and detector plane locations
ax.scatter(
sim.source[0] * 100,
sim.source[1] * 100,
sim.source[2] * 100,
color="red",
marker="*",
linewidth=5,
label="Source",
)
ax.scatter(
sim.detector[0] * 100,
sim.detector[1] * 100,
sim.detector[2] * 100,
color="blue",
marker="*",
linewidth=10,
label="Detector",
)
# Plot the final proton positions of some (not all) of the protons
ind = slice(None, None, 200)
ax.scatter(
sim.x[ind, 0] * 100, sim.x[ind, 1] * 100, sim.x[ind, 2] * 100, label="Protons",
)
ax.legend()
[6]:
<matplotlib.legend.Legend at 0x7f024a471750>
A ‘synthetic proton radiograph’ can now be constructed by creating a 2D histogram of proton positions in the image plane. The synthetic radiograph function takes two keywords:
‘size’ gives the locations of the lower left and upper right corners of the detector grid in image plane coordinates.
‘bins’ is the number of histogram bins to be used in the horizontal and vertical directions. Using more bins creates a higher resolution image, but at the cost of more noise.
[7]:
# A function to reduce repetative plotting
def plot_radiograph(hax, vax, intensity):
fig, ax = plt.subplots(figsize=(8, 8))
plot = ax.pcolormesh(
hax.to(u.cm).value,
vax.to(u.cm).value,
intensity.T,
cmap="Blues_r",
shading="auto",
)
cb = fig.colorbar(plot)
cb.ax.set_ylabel("Intensity")
ax.set_aspect("equal")
ax.set_xlabel("X (cm), Image plane")
ax.set_ylabel("Z (cm), Image plane")
ax.set_title("Synthetic Proton Radiograph")
size = np.array([[1, 1], [1, 1]]) * 1.5 * u.cm
bins = [200, 200]
hax, vax, intensity = sim.synthetic_radiograph(size=size, bins=bins)
plot_radiograph(hax, vax, intensity)
As expected, the outwardpointing electric field in the sphere has deflected the protons out of the central region, leaving a dark shadow.
Kugland et al. 2012 and Bott et al. 2017 define the dimensionless “contrast parameter” that separates different regimes of proton radiography:
Where \(l\) is the distance from the source to the grid, \(a\) is the spatial scale of the scattering electromagnetic fields, and \(\alpha\) is the particle deflection angle. The value of \(\mu\) can fall in one of three regimes:
where \(\mu_c \sim 1\) is a characteristic value at which particle paths cross, leading to the formation of bright caustics. Correctly placing a radiograph in the correct regime is necessary to determine which analysis techniques can be applied to it.
The maximum deflection angle can be calculated after the simulation has run by comparing the initial and final velocity vectors of each particle
[8]:
max_deflection = sim.max_deflection
print(f"Maximum deflection α = {np.rad2deg(max_deflection):.2f}")
Maximum deflection α = 2.77 deg
The spatial scale of the field constructed in this example is \(\sim\) 1 mm, and \(l\) is approximately the distance from the source to the grid origin. Therefore, we can calculate the value of \(\mu\)
[9]:
a = 1 * u.mm
l = np.linalg.norm(sim.source * u.m).to(u.mm)
mu = l * max_deflection.value / a
print(f"a = {a}")
print(f"l = {l:.1f}")
print(f"μ = {mu:.2f}")
a = 1.0 mm
l = 10.0 mm
μ = 0.48
which places this example in the nonlinear injective regime.
Options¶
For sake of comparison, here is the result achieved by setting distribution = 'uniform'
in the create_particles() function.
[10]:
sim.create_particles(
1e5, 3 * u.MeV, max_theta=np.pi / 15 * u.rad, distribution="uniform"
)
sim.run()
size = np.array([[1, 1], [1, 1]]) * 1.5 * u.cm
bins = [200, 200]
hax, vax, intensity = sim.synthetic_radiograph(size=size, bins=bins)
plot_radiograph(hax, vax, intensity)
Creating Particles
Particles on grid: 0% 2.2e+02/8.6e+04 particles
Run completed
Fraction of particles tracked: 85.9%
Fraction of tracked particles that entered the grid: 66.0%
Fraction of tracked particles deflected away from the detector plane: 0.0%
This page was generated by
nbsphinx from
docs/notebooks/dispersion/two_fluid_dispersion.ipynb.
Interactive online version:
.
Dispersion: A Full Two Fluid Solution¶
This notebook walks through the functionalty of the two_fluid_dispersion_solution() function. This function computes the wave frequencies for given wavenumbers and plasma parameters based on the analytical solution presented by Bellan 2012 to the Stringer 1963 two fluid dispersion relation. The two fluid dispersion equaiton assumes a uniform magnetic field, a zero D.C. electric field, and lowfrequency waves \(\omega / k c \ll 1\) which equates to
where
\(\omega\) is the wave frequency, \(k\) is the wavenumber, \(v_A\) is the Alfvén velocity, \(c_s\) is the sound speed, \(\omega_{ci}\) is the ion gyrofrequency, and \(\omega_{pe}\) is the electron plasma frequency.
The approach outlined in Section 5 of Bellan 2012 produces exact roots to the above dispersion equation for all three modes (fast, acoustic, and Alfvén) without having to make additional approximations. The following dispersion relation is what the two_fluid_dispersion_solution() function computes.
where \(j = 0\) represents the fast mode, \(j = 1\) represents the Alfvén mode, and \(j = 2\) represents the acoustic mode. Additionally,
Contents:¶
[1]:
%matplotlib inline
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.constants.si import c
from matplotlib import colors
from matplotlib.ticker import MultipleLocator
from mpl_toolkits.axes_grid1 import make_axes_locatable
from plasmapy.dispersion.two_fluid_dispersion import two_fluid_dispersion_solution
from plasmapy.formulary import parameters as pfp
from plasmapy.particles import Particle
plt.rcParams["figure.figsize"] = [10.5, 0.56 * 10.5]
Wave Propagating at 45 Degrees¶
Below we define the required parameters to compute the wave frequencies.
[2]:
# define input parameters
inputs = {
"k": np.linspace(10 ** 7, 10 ** 2, 10000) * u.rad / u.m,
"theta": 45 * u.deg,
"n_i": 5 * u.cm ** 3,
"B": 8.3e9 * u.T,
"T_e": 1.6e6 * u.K,
"T_i": 4.0e5 * u.K,
"ion": Particle("p+"),
}
# a few useful plasma parameters
params = {
"n_e": inputs["n_i"] * abs(inputs["ion"].integer_charge),
"cs": pfp.ion_sound_speed(inputs["T_e"], inputs["T_i"], inputs["ion"],),
"va": pfp.Alfven_speed(inputs["B"], inputs["n_i"], ion=inputs["ion"],),
"wci": pfp.gyrofrequency(inputs["B"], inputs["ion"]),
}
params["lpe"] = pfp.inertial_length(params["n_e"], "e")
params["wpe"] = pfp.plasma_frequency(params["n_e"], "e")
The computed wave frequencies (\(rad/s\)) are returned in a dictionary with keys representing the wave modes and the values being an Astropy Quantity. Since our inputs had a scalar \(\theta\) and a 1D array of \(k\)’s, the computed wave frequencies will be a 1D array of size equal to the size of the \(k\) array.
[3]:
# compute
omegas = two_fluid_dispersion_solution(**inputs)
(list(omegas.keys()), omegas["fast_mode"], omegas["fast_mode"].shape)
[3]:
(['fast_mode', 'alfven_mode', 'acoustic_mode'],
<Quantity [1.63839709e02, 1.80262570e01, 3.44262572e01, ...,
1.52032171e+03, 1.52047365e+03, 1.52062560e+03] rad / s>,
(10000,))
Let’s plot the results of each wave mode.
[4]:
fs = 14 # default font size
figwidth, figheight = plt.rcParams["figure.figsize"]
figheight = 1.6 * figheight
fig = plt.figure(figsize=[figwidth, figheight])
# normalize data
k_prime = inputs["k"] * params["lpe"]
# plot
plt.plot(
k_prime, np.real(omegas["fast_mode"] / params["wpe"]), "r.", ms=1, label="Fast",
)
ax = plt.gca()
ax.plot(
k_prime, np.real(omegas["alfven_mode"] / params["wpe"]), "b.", ms=1, label="Alfvèn",
)
ax.plot(
k_prime,
np.real(omegas["acoustic_mode"] / params["wpe"]),
"g.",
ms=1,
label="Acoustic",
)
# adjust axes
ax.set_xlabel(r"$kc / \omega_{pe}$", fontsize=fs)
ax.set_ylabel(r"$Re(\omega / \omega_{pe})$", fontsize=fs)
ax.set_yscale("log")
ax.set_xscale("log")
ax.set_ylim(1e6, 2e2)
ax.tick_params(
which="both", direction="in", width=1, labelsize=fs, right=True, length=5,
)
# annotate
text = (
f"$v_A/c_s = {params['va'] / params['cs']:.1f} \qquad "
f"c/v_A = 10^{np.log10(c / params['va']):.0f} \qquad "
f"\\theta = {inputs['theta'].value:.0f}"
"^{\\circ}$"
)
ax.text(0.25, 0.95, text, transform=ax.transAxes, fontsize=18)
ax.legend(loc="upper left", markerscale=5, fontsize=fs)
[4]:
<matplotlib.legend.Legend at 0x7f6ad22ca810>
Wave frequencies on the ktheta plane¶
Let us now look at the distribution of \(\omega\) on a \(k\)\(\theta\) plane.
[5]:
# define input parameters
inputs = {
"k": np.linspace(10 ** 7, 10 ** 2, 10000) * u.rad / u.m,
"theta": np.linspace(5, 85, 100) * u.deg,
"n_i": 5 * u.cm ** 3,
"B": 8.3e9 * u.T,
"T_e": 1.6e6 * u.K,
"T_i": 4.0e5 * u.K,
"ion": Particle("p+"),
}
# a few useful plasma parameters
params = {
"n_e": inputs["n_i"] * abs(inputs["ion"].integer_charge),
"cs": pfp.ion_sound_speed(inputs["T_e"], inputs["T_i"], inputs["ion"],),
"va": pfp.Alfven_speed(inputs["B"], inputs["n_i"], ion=inputs["ion"],),
"wci": pfp.gyrofrequency(inputs["B"], inputs["ion"]),
}
params["lpe"] = pfp.inertial_length(params["n_e"], "e")
params["wpe"] = pfp.plasma_frequency(params["n_e"], "e")
Since the \(\theta\) and \(k\) values are now 1D arrays, the returned wave frequencies will be 2D arrays with the first dimension matching the size of \(k\) and the second dimension matching the size of \(\theta\).
[6]:
# compute
omegas = two_fluid_dispersion_solution(**inputs)
(
omegas["fast_mode"].shape,
omegas["fast_mode"].shape[0] == inputs["k"].size,
omegas["fast_mode"].shape[1] == inputs["theta"].size,
)
[6]:
((10000, 100), True, True)
Let’s plot (the fast mode)!
[7]:
fs = 14 # default font size
figwidth, figheight = plt.rcParams["figure.figsize"]
figheight = 1.6 * figheight
fig = plt.figure(figsize=[figwidth, figheight])
# normalize data
k_prime = inputs["k"] * params["lpe"]
zdata = np.transpose(np.real(omegas["fast_mode"].value)) / params["wpe"].value
# plot
im = plt.imshow(
zdata,
aspect="auto",
origin="lower",
extent=[
np.min(k_prime.value),
np.max(k_prime.value),
np.min(inputs["theta"].value),
np.max(inputs["theta"].value),
],
interpolation=None,
cmap=plt.cm.Spectral,
)
ax = plt.gca()
# # adjust axes
ax.set_xscale("linear")
ax.set_xlabel(r"$kc/\omega_{pe}$", fontsize=fs)
ax.set_ylabel(r"$\theta$ [$deg.$]", fontsize=fs)
ax.tick_params(
which="both",
direction="in",
width=2,
labelsize=fs,
right=True,
top=True,
length=10,
)
# Add colorbar
divider = make_axes_locatable(ax)
cax = divider.append_axes("top", size="5%", pad=0.07)
cbar = plt.colorbar(
im, cax=cax, orientation="horizontal", ticks=None, fraction=0.05, pad=0.0,
)
cbar.ax.tick_params(
axis="x",
direction="in",
width=2,
length=10,
top=True,
bottom=False,
labelsize=fs,
pad=0.0,
labeltop=True,
labelbottom=False,
)
cbar.ax.xaxis.set_label_position("top")
cbar.set_label(r"$\omega/\omega_{pe}$", fontsize=fs, labelpad=8)
Reproduce Figure 1 from Bellan 2012¶
Figure 1 of Bellan 2012 chooses parameters such that \(\beta = 0.4\) and \(\Lambda=0.4\). Below we define parameters to approximate Bellan’s assumptions.
[8]:
# define input parameters
inputs = {
"B": 400e4 * u.T,
"ion": Particle("He+"),
"n_i": 6.358e19 * u.m ** 3,
"T_e": 20 * u.eV,
"T_i": 10 * u.eV,
"theta": np.linspace(0, 90) * u.deg,
"k": (2 * np.pi * u.rad) / (0.56547 * u.m),
}
# a few useful plasma parameters
params = {
"n_e": inputs["n_i"] * abs(inputs["ion"].integer_charge),
"cs": pfp.cs_(inputs["T_e"], inputs["T_i"], inputs["ion"]),
"wci": pfp.wc_(inputs["B"], inputs["ion"]),
"va": pfp.va_(inputs["B"], inputs["n_i"], ion=inputs["ion"]),
}
params["beta"] = (params["cs"] / params["va"]).value ** 2
params["wpe"] = pfp.wp_(params["n_e"], "e")
params["Lambda"] = (inputs["k"] * params["va"] / params["wci"]).value ** 2
(params["beta"], params["Lambda"])
[8]:
(0.4000832135717194, 0.4000017351804854)
[9]:
# compute
omegas = two_fluid_dispersion_solution(**inputs)
[10]:
# generate data for plots
plt_vals = {}
for mode, arr in omegas.items():
norm = (np.absolute(arr) / (inputs["k"] * params["va"])).value ** 2
plt_vals[mode] = {
"x": norm * np.sin(inputs["theta"].to(u.rad).value),
"y": norm * np.cos(inputs["theta"].to(u.rad).value),
}
[11]:
fs = 14 # default font size
figwidth, figheight = plt.rcParams["figure.figsize"]
figheight = 1.6 * figheight
fig = plt.figure(figsize=[figwidth, figheight])
# Fast mode
plt.plot(
plt_vals["fast_mode"]["x"], plt_vals["fast_mode"]["y"], linewidth=2, label="Fast",
)
ax = plt.gca()
# adjust axes
ax.set_xlabel(r"$(\omega / k v_A)^2 \, \sin \theta$", fontsize=fs)
ax.set_ylabel(r"$(\omega / k v_A)^2 \, \cos \theta$", fontsize=fs)
ax.set_xlim(0.0, 1.5)
ax.set_ylim(0.0, 2.0)
for spine in ax.spines.values():
spine.set_linewidth(2)
ax.minorticks_on()
ax.tick_params(which="both", labelsize=fs, width=2)
ax.tick_params(which="major", length=10)
ax.tick_params(which="minor", length=5)
ax.xaxis.set_major_locator(MultipleLocator(0.5))
ax.xaxis.set_minor_locator(MultipleLocator(0.1))
ax.yaxis.set_major_locator(MultipleLocator(0.5))
ax.yaxis.set_minor_locator(MultipleLocator(0.1))
# Alfven mode
plt.plot(
plt_vals["alfven_mode"]["x"],
plt_vals["alfven_mode"]["y"],
linewidth=2,
label="Alfv$\`{e}$n",
)
# Acoustic mode
plt.plot(
plt_vals["acoustic_mode"]["x"],
plt_vals["acoustic_mode"]["y"],
linewidth=2,
label="Acoustic",
)
# annotations
plt.legend(fontsize=fs, loc="upper right")
[11]:
<matplotlib.legend.Legend at 0x7f6acfa9cfd0>
This page was generated by
nbsphinx from
docs/notebooks/diagnostics/thomson.ipynb.
Interactive online version:
.
Thomson Scattering: Spectral Density¶
The thomson.spectral_density function calculates the spectral density function S(k,w), which is one of several terms that determine the scattered power spectrum for the Thomson scattering of a probe laser beam by a plasma. In particular, this function calculates \(S(k,w)\) for the case of a plasma consisting of one or more ion species and electron populations under the assumption that all of the ion species and the electron fluid have Maxwellian velocity distribution functions and that the combined plasma is quasineutral. In this regime, the spectral density is given by the equation:
where \(\chi_e\) is the electron component susceptibility of the plasma and \(\epsilon = 1 + \sum_e \chi_e + \sum_i \chi_i\) is the total plasma dielectric function (with \(\chi_i\) being the ion component of the susceptibility), \(Z_i\) is the charge of each ion, \(k\) is the scattering wavenumber, \(\omega\) is the scattering frequency, and the functions \(f_{e0,e}\) and \(f_{i0,i}\) are the Maxwellian velocity distributions for the electrons and ion species respectively.
Thomson scattering can be either noncollective (the scattered spectrum is a linear sum of the light scattered by individual particles) or collective (the scattered spectrum is dominated by scattering off of collective plasma waves). The thomson.spectral_density function can be used in both cases. These regimes are delineated by the dimensionless constant \(\alpha\):
where \(\lambda_{De}\) is the Debye length. \(\alpha > 1\) corresponds to collective scattering, while \(\alpha < 1\) corresponds to noncollective scattering. Depending on which of these regimes applies, fitting the scattered spectrum can provide the electron (and sometimes ion) density and temperature. Doppler shifting of the spectrum can also provide a measurement of the drift velocity of each plasma species.
For a detailed explanation of the underlying physics (and derivations of these expressions), see “Plasma Scattering of Electromagnetic Radiation” by Sheffield et al.
[1]:
%matplotlib inline
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from plasmapy.diagnostics import thomson
Construct parameters that define the Thomson diagnostic setup, the probing beam and scattering collection. These parameters will be used for all examples.
[2]:
# The probe wavelength can in theory be anything, but in practice integer frequency multiples of the Nd:YAG wavelength
# 1064 nm are used (532 corresponds to a frequencydoubled probe beam from such a laser).
probe_wavelength = 532 * u.nm
# Array of wavelengths over which to calcualte the spectral distribution
wavelengths = (
np.arange(probe_wavelength.value  60, probe_wavelength.value + 60, 0.01) * u.nm
)
# The scattering geometry is defined by unit vectors for the orientation of the probe laser beam (probe_n) and
# the path from the scattering volume (where the measurement is made) to the detector (scatter_n).
# These can be setup for any experimental geometry.
probe_vec = np.array([1, 0, 0])
scattering_angle = np.deg2rad(63)
scatter_vec = np.array([np.cos(scattering_angle), np.sin(scattering_angle), 0])
In order to calcluate the scattered spectrum, we must also include some information about the plasma. For this plot we’ll allow the fract
, ion_species
, fluid_vel
, and ion_vel
keywords to keep their default values, describing a singlespecies H+ plasma at rest in the laboratory frame.
[3]:
ne = 2e17 * u.cm ** 3
Te = 12 * u.eV
Ti = 10 * u.eV
alpha, Skw = thomson.spectral_density(
wavelengths,
probe_wavelength,
ne,
Te,
Ti,
probe_vec=probe_vec,
scatter_vec=scatter_vec,
)
fig, ax = plt.subplots()
ax.plot(wavelengths, Skw, lw=2)
ax.set_xlim(probe_wavelength.value  10, probe_wavelength.value + 10)
ax.set_ylim(0, 1e13)
ax.set_xlabel("$\lambda$ (nm)")
ax.set_ylabel("S(k,w)")
ax.set_title("Thomson Scattering Spectral Density")
[3]:
Text(0.5, 1.0, 'Thomson Scattering Spectral Density')
Example Cases in Different Scattering Regimes¶
We will now consider several example cases in different scattering regimes. In order to facilitate this, we’ll set up each example as a dictionary of plasma parameters:
A singlespecies, stationary hydrogen plasma with a density and temperature that results in a scattering spectrum dominated by scattering off of single electrons.
[4]:
non_collective = {
"name": "NonCollective Regime",
"n": 5e15 * u.cm ** 3,
"Te": 40 * u.eV,
"Ti": np.array([10]) * u.eV,
"ion_species": ["H+"],
"electron_vel": np.array([[0, 0, 0]]) * u.km / u.s,
"ion_vel": np.array([[0, 0, 0]]) * u.km / u.s,
}
A singlespecies, stationary hydrogen plasma with a density and temperature that result in weakly collective scattering (scattering paramter \(\alpha\) approaching 1)
[5]:
weakly_collective = {
"name": "Weakly Collective Regime",
"n": 2e17 * u.cm ** 3,
"Te": 20 * u.eV,
"Ti": 10 * u.eV,
"ion_species": ["H+"],
"electron_vel": np.array([[0, 0, 0]]) * u.km / u.s,
"ion_vel": np.array([[0, 0, 0]]) * u.km / u.s,
}
A singlespecies, stationary hydrogen plasma with a density and temperature that result in a spectrum dominated by multiparticle scattering, including scattering off of ions.
[6]:
collective = {
"name": "Collective Regime",
"n": 5e17 * u.cm ** 3,
"Te": 10 * u.eV,
"Ti": 4 * u.eV,
"ion_species": ["H+"],
"electron_vel": np.array([[0, 0, 0]]) * u.km / u.s,
"ion_vel": np.array([[0, 0, 0]]) * u.km / u.s,
}
A case identical to the collective example above, except that now the electron fluid has a substantial drift velocity parallel to the probe laser and the ions have a drift (relative to the electrons) at an angle.
[7]:
drifts = {
"name": "Drift Velocities",
"n": 5e17 * u.cm ** 3,
"Te": 10 * u.eV,
"Ti": 10 * u.eV,
"ion_species": ["H+"],
"electron_vel": np.array([[700, 0, 0]]) * u.km / u.s,
"ion_vel": np.array([[600, 100, 0]]) * u.km / u.s,
}
A case identical to the collective example, except that now the plasma consists 25% He+1 and 75% C+5, and two electron populations exist with different temperatures.
[8]:
two_species = {
"name": "Two Ion and Electron Components",
"n": 5e17 * u.cm ** 3,
"Te": np.array([50, 10]) * u.eV,
"Ti": np.array([10, 50]) * u.eV,
"efract": np.array([0.5, 0.5]),
"ifract": np.array([0.25, 0.75]),
"ion_species": ["He4 1+", "C12 5+"],
"electron_vel": np.array([[0, 0, 0], [0, 0, 0]]) * u.km / u.s,
"ion_vel": np.array([[0, 0, 0], [0, 0, 0]]) * u.km / u.s,
}
[9]:
examples = [non_collective, weakly_collective, collective, drifts, two_species]
For each example, plot the the spectral distribution function over a large range to show the broad electron scattering feature (top row) and a narrow range around the probe wavelength to show the ion scattering feature (bottom row)
[10]:
fig, ax = plt.subplots(ncols=len(examples), nrows=2, figsize=[25, 10])
fig.subplots_adjust(wspace=0.4, hspace=0.4)
lbls = "abcdefg"
for i, x in enumerate(examples):
alpha, Skw = thomson.spectral_density(
wavelengths,
probe_wavelength,
x["n"],
x["Te"],
x["Ti"],
ifract=x.get("ifract"),
efract=x.get("efract"),
ion_species=x["ion_species"],
electron_vel=x["electron_vel"],
probe_vec=probe_vec,
scatter_vec=scatter_vec,
)
ax[0][i].axvline(x=probe_wavelength.value, color="red") # Mark the probe wavelength
ax[0][i].plot(wavelengths, Skw)
ax[0][i].set_xlim(probe_wavelength.value  15, probe_wavelength.value + 15)
ax[0][i].set_ylim(0, 1e13)
ax[0][i].set_xlabel("$\lambda$ (nm)")
ax[0][i].set_title(lbls[i] + ") " + x["name"] + "\n$\\alpha$={:.4f}".format(alpha))
ax[1][i].axvline(x=probe_wavelength.value, color="red") # Mark the probe wavelength
ax[1][i].plot(wavelengths, Skw)
ax[1][i].set_xlim(probe_wavelength.value  1, probe_wavelength.value + 1)
ax[1][i].set_ylim(0, 1.1 * np.max(Skw.value))
ax[1][i].set_xlabel("$\lambda$ (nm)")
Plots of the spectral density function (Skw) which determines the amount of light scattered into different wavelengths.
In the noncollective regime only the electron feature is visible.
In the weakly collective regime (alpha approaches 1) an ion feature starts to appear and the electron feature is distorted
In the collective regime both features split into two peaks, corresponding to scattering off of forward and backwards propagating plasma oscillations.
The introduction of drift velocities introduces several Doppler shifts in the calculations, resulting in a shifted spectrum.
Including multiple ion and electron populations modifies the ion and electron features respectively.
PlasmaPy is an open source communitydeveloped core Python 3.7+ package for plasma physics currently under development.
PlasmaPy’s Vision Statement¶
About PlasmaPy¶
PlasmaPy is a communitydeveloped and communitydriven free and open source Python package that provides common functionality required for plasma physics in a single, reliable codebase.
Motivation¶
In recent years, researchers in many different scientific disciplines have worked together to develop core Python packages such as Astropy, SunPy, and SpacePy. These packages provide core functionality, common frameworks for data visualization and analysis, and educational tools for their respective scientific disciplines. We believe that a similar cooperatively developed package for plasma physics will greatly benefit our field. In this document, we lay out our vision for PlasmaPy: a communitydeveloped and communitydriven open source core Python software package for plasma physics.
There is considerable need in plasma physics for open, general purpose software framework using modern best practices for scientific computing. As most scientific programmers are largely selftaught, software often does not take advantage of these practices and is instead written in a rush to produce results for the next research paper. The resulting code is often difficult to read and maintain, the documentation is usually inadequate, and tests are typically implemented late in the development process if at all. Legacy code is often written in low level languages such as Fortran, which typically makes compiling and installing packages difficult and frustrating, especially if it calls external libraries. It is also unusual to share code, and access to major software is often restricted in some way, resulting in many different programs and packages which do essentially the same thing but with little or no interoperability. These factors lead to research that is difficult to reproduce, and present a significant barrier to entry for new users.
The plasma physics community is slowly moving in the open source direction. Several different types of packages and software have been released under open source licences, including the UCLA PIC codes, PICCANTE, EPOCH, VPIC, PIConGPU, WARP, the FLASH framework, Athena, and PENCIL. These projects are built as individual packages, are written in different programming languages, and often have many dependencies on specific packages. Python packages such as Astropy, SunPy, and SpacePy have had notable success providing open source alternatives to legacy code in related fields. We are grateful to these communities for their hard work, and hope to build upon their accomplishments for the field of plasma physics.
An end user might not always be interested in a complicated powerpack to perform one specific task on supercomputers. She might also be interested in performing some basic plasma physics calculations, running small desktop scale simulations to test preliminary ideas (e.g., 1D MHD/PIC or test particles), or even comparing data from two different sources (simulations vs. spacecraft). Such tasks require a central platform. This is where PlasmaPy comes in.
PlasmaPy Community Code of Conduct¶
Please see the attached Community Code of Conduct.
Organizational Structure¶
The Coordinating Committee (CC) will oversee the PlasmaPy project and code development. The CC will ensure that roles are being filled, facilitate communitywide communication, coordinate and delegate tasks, manage the project repository, oversee the code review process, regulate intercompatibility between different subpackages, seek funding mechanisms, facilitate compromises and cooperation, enforce the code of conduct, and foster a culture of appreciation.
The Community Engagement Committee (CEC) will be responsible for organizing conferences, trainings, and workshops; maintaining and moderating social media groups and accounts; overseeing PlasmaPy’s website; and communicating with the PlasmaPy and plasma physics communities. The CEC will facilitate partnerships with groups such as Software Carpentry.
Each subpackage will have lead and deputy coordinators who will guide and oversee the development of that subpackage.
The Accessibility Coordinator will work to ensure that the PlasmaPy codebase, documentation, and practices are accessible to disabled students and scientists. Additional roles include the Webpage Maintainer, the Release Coordinator, and the Testing Coordinator.
The work undertaken by each of these groups and coordinators should be done openly and transparently, except where confidentiality is needed. We will strive to have multiple subfields from plasma physics in each committee. Major decisions should ideally be made by general consensus among the PlasmaPy community, but when consensus is not possible then the committees may decide via majority vote. Much of this section is following the organizational structure of Astropy.
Development Procedure¶
The initial developers of PlasmaPy will create a flexible development roadmap that outlines and prioritizes subpackages to be developed. The developers will survey existing open source Python packages in plasma physics. Priority will be given to determining how data will be stored and structured. Developers will break up into small groups to work on different subpackages. These small groups will communicate regularly and work towards interoperability and common coding practices.
Because Python is new to many plasma physicists, community engagement is vital. The CEC will arrange occasional informal trainings early in the project that are director towards the initial developers.
New code and edits should be submitted as a pull request to the development branch of the PlasmaPy repository on GitHub. The pull request will undergo a code review by the subpackage maintainers and/or the CC, who will provide suggestions on how the contributor may update the pull request. Subpackage maintainers will generally be responsible for deciding on pull requests with minor changes, while pull requests with major changes should be decided jointly by the subpackage maintainers and the CC. The CC and CEC will develop a friendly guide on how users may contribute new code to PlasmaPy.
New code should conform to the PEP 8 style guide for Python code and the established coding style within PlasmaPy. New code should be submitted with documentation and tests. Documentation should be written primarily in docstrings and follow the numpydoc documentation style guide. Every new module, class and function should have an appropriate docstring. The documentation should describe the interface and the purpose for the method, but generally not the implementation. The code itself should be readable enough to be able to explain how it works. Documentation should be updated when the code is edited. The tests should cover new functionality (especially methods with complex logic), but the tests should also be readable and easy to maintain. Existing tests should be updated when necessary [e.g., during the initial development of a new feature when the application program interface (API) is not yet stable], but with caution since this may imply loss of backwards compatibility.
Members of the PlasmaPy community may submit PlasmaPy Enhancement Proposals (PLEPs) to suggest changes such as major reorganization of a subpackage, creation of a new subpackage, nonbackwards compatible changes to a stable package, or significant changes to policies and procedures related to the organization of this project. The issues list on GitHub will generally be more appropriate for changes that do not require community discussion. The CC shall maintain a GitHub repository of PLEPs. PLEPs will be made openly available for community discussion and transparency for a period of at least four weeks, during which time the proposal may be updated and revised by the proposers. The CC shall approve or decline these proposals after seeking community input. The rationale behind the decision and a summary of the community discussion shall be recorded along with the PLEP.
Programming Guidelines¶
Choice of Languages¶
PlasmaPy shall be written using Python 3. PlasmaPy shall initially guarantee compatibility with Python 3.6 and above. Python 3 is continually growing, so we will proceed on the general principle that future updates to PlasmaPy remain compatible with releases of Python that are up to two years old. Python 2.7 and below will not be supported as these versions will no longer be updated past 2020. The core package will initially be written solely in Python.
Code readability is more important than optimization, except when performance is critical. Code should be optimized only after getting it to work, and primarily for where there is a performance bottleneck. Performancecritical parts of the core package will preferably be written using Numba to achieve compiled speeds while maintaining the significant advantages of using a high level language.
Versioning¶
PlasmaPy will use Semantic Versioning. Releases will be given version numbers of the form MAJOR.MINOR.PATCH, where MAJOR, MINOR, and PATCH are nonnegative integers. Starting with version 1.0, MAJOR will be incremented when backwards incompatible changes are made, MINOR will be incremented when new backwardscompatible functionality is added, and PATCH will be incremented when backwardscompatible bug fixes are made.
Development releases will have MAJOR equal to zero and start at version 0.1. The API should not be considered stable during the development phase. PlasmaPy will release version 1.0 once it has a stable public API that users are depending on for production code.
All releases will be provided with release notes and change log entries, and a table will be provided that describes the stability of the public API for each PlasmaPy subpackage.
Dependencies¶
Dependencies have the advantage of providing capabilities that will enhance PlasmaPy and speed up its development, but the disadvantage that they can make manual installation more difficult and potentially frustrating. Package managers such as Anaconda and Homebrew greatly simplify installation of Python packages, but there will be situations where manual installation is necessary (e.g., on some supercomputers without package managers). The core package should be able to be imported using a minimal number of packages (e.g., NumPy, SciPy, and matplotlib) without getting an import error. Additional packages may be included as dependencies of the core package if there is a strong need for it, and if these packages are easily installed with currently available package managers. Subpackages may use additional dependencies when appropriate.
Affiliated Packages¶
We will follow the practice of Astropy by having a core package and affiliated packages. The core package will contain common tools and base functionality that most plasma physicists will need. The affiliated packages contained in separate repositories will include more specialized functionality that is needed for subfields of plasma physics. This approach will reduce the likelihood of scope creep for the core package while maintaining avenues for broader development.
Units¶
Multiple sets of units are used by plasma physicists. There exist some peculiarities with how units are used within plasma physics, such as how an electron volt is typically used as a measurement of temperature. Code will be most readable and maintainable if written assuming a particular set of units, but there should be enough flexibility for people in different subfields to choose their preferred set of units. As the generally most common accepted international standard, SI base units will be utilized. We will use an existing Python module (e.g., astropy.units or pint) to assign units to variables and allow straightforward conversion between different systems of units.
Installing PlasmaPy¶
Note
If you would like to contribute to PlasmaPy, please refer to the instructions on installing PlasmaPy for development.
Requirements¶
PlasmaPy requires Python version 3.7 or newer. PlasmaPy requires the following packages for installation:
NumPy — 1.18.1 or newer
SciPy — 1.2 or newer
Astropy — 4.0 or newer
pandas — 1.0 or newer
xarray — above 0.14
tqdm — 4.56 or newer
cached_property — 1.5.2 or newer
PlasmaPy also depends on the following packages for optional features:
h5py — 2.8 or newer
lmfit — 1.0.1 or newer
matplotlib — 2.0 or newer
mpmath — 1.0 or newer
Installing PlasmaPy¶
Installation with pip¶
To install the most recent release of PlasmaPy on PyPI with pip into an existing Python environment with both required and optional dependencies, run
pip install plasmapy[all]
PlasmaPy may be installed with the required but not the optional dependencies
with the following command, though this may result in an ImportError
when
using certain specialized functionality.
pip install plasmapy
Note
In some systems, it may be necessary to use pip3
instead of pip
.
Installation with conda¶
We recommend installing PlasmaPy from a Python environment created using Conda. Conda allows us to create and switch between Python environments that are isolated from each other and the system installation (in contrast to this xkcd).
After installing conda, create a PlasmaPy environment with all required and optional dependencies by running:
conda create n env_name python=3.9 plasmapy c condaforge
where env_name
is replaced by the name of the environment.
To activate this environment, run:
conda activate env_name
Building and installing from source code¶
Obtaining source code¶
Stable release¶
The source code for the most recent stable release of PlasmaPy can be downloaded from PyPI or from Zenodo.
Development version on GitHub¶
If you have git installed on your computer, you may clone PlasmaPy’s GitHub repository and access source code from the most recent development version by running:
git clone https://github.com/PlasmaPy/PlasmaPy.git
The repository will be cloned inside a new subdirectory called PlasmaPy
.
If you do not have git installed on your computer, then you may download the most recent source code from PlasmaPy’s GitHub repository by selecting “Clone or Download”, which will give you the option to download a zip file.
Note
Cloning a repository with HTTPS as above is recommended, but you may also clone a repository using SSH as a more secure alternative.
Note
The How to Contribute guide has instructions on how to fork a repository and create branches so that you may make pull requests.
Building and installing¶
In the PlasmaPy
directory, run
pip install e .[all]
where e
makes the installation editable and [all]
will ensure that
all optional dependencies are installed. PlasmaPy could also be installed
by running
python setup.py install
Note, however, that this does not download all the dependencies. Check the
requirements/requirements.txt
file for the current set.
Examples¶
Here we catalog all the example Jupyter notebooks that have been created for
the various functionality contained in plasmapy
.
General¶
This page was generated by
nbsphinx from
docs/notebooks/ExB_drift.ipynb.
Interactive online version:
.
[1]:
%matplotlib inline
import math
import matplotlib.pyplot as plt
import numpy as np
from astropy import units as u
from mpl_toolkits.mplot3d import Axes3D
from plasmapy import formulary, particles
Physics of the ExB drift¶
Consider a single particle of mass \(m\) and charge \(q\) in a constant, uniform magnetic field \(\mathbf{B}=B\ \hat{\mathbf{z}}\). In the absence of external forces, it travels with velocity \(\mathbf{v}\) governed by the equation of motion
which simply equates the net force on the particle to the corresponding Lorentz force. Assuming the particle initially (at time \(t=0\)) has \(\mathbf{v}\) in the \(x,z\) plane (with \(v_y=0\)), solving reveals
while the parallel velocity \(v_z\) is constant. This indicates that the particle gyrates in a circular orbit in the \(x,y\) plane with constant speed \(v_\perp\), angular frequency \(\omega_c = \frac{\lvert q\rvert B}{m}\), and Larmor radius \(r_L=\frac{v_\perp}{\omega_c}\).
As an example, take one proton p+
moving with velocity \(1\ m/s\) in the \(x\)direction at \(t=0\):
[2]:
# Setup proton in uniform B field
B = 5 * u.T
proton = particles.Particle("p+")
omega_c = formulary.parameters.gyrofrequency(B, proton)
v_perp = 1 * u.m / u.s
r_L = formulary.parameters.gyroradius(B, proton, Vperp=v_perp)
We can define a function that evolves the particle’s position according to the relations above describing \(v_x,v_y\), and \(v_z\). The option to add a constant drift velocity \(v_d\) to the solution is included as an argument, though this drift velocity is zero by default:
[3]:
def single_particle_trajectory(v_d=np.array([0, 0, 0])):
# Set time resolution & velocity such that proton goes 1 meter along B per rotation
T = 2 * math.pi / omega_c.value # rotation period
v_parallel = 1 / T * u.m / u.s
dt = T / 1e2 * u.s
# Set initial particle position
x = []
y = []
xt = 0 * u.m
yt = r_L
# Evolve motion
timesteps = np.arange(0, 10 * T, dt.value)
for t in list(timesteps):
v_x = v_perp * math.cos(omega_c.value * t) + v_d[0]
v_y = v_perp * math.sin(omega_c.value * t) + v_d[1]
xt += +v_x * dt
yt += +v_y * dt
x.append(xt.value)
y.append(yt.value)
x = np.array(x)
y = np.array(y)
z = v_parallel.value * timesteps
return x, y, z
Executing with the default argument and plotting the particle trajectory gives the expected helical motion, with a radius equal to the Larmor radius:
[4]:
x, y, z = single_particle_trajectory()
fig = plt.figure()
ax = fig.gca(projection="3d")
ax.plot(x, y, z, label="$\mathbf{F}=0$")
ax.legend()
bound = 3 * r_L.value
ax.set_xlim([bound, bound])
ax.set_ylim([bound, bound])
ax.set_zlim([0, 10])
ax.set_xlabel("x [m]")
ax.set_ylabel("y [m]")
ax.set_zlabel("z [m]")
plt.show()
print(f"r_L = {r_L.value:.2e} m")
print(f"omega_c = {omega_c.value:.2e} rads/s")
r_L = 2.09e09 m
omega_c = 4.79e+08 rads/s
How does this motion change when a constant external force \(\mathbf{F}\) is added? The new equation of motion is
and we can find a solution by considering velocities of the form \(\mathbf{v}=\mathbf{v}_\parallel + \mathbf{v}_L + \mathbf{v}_d\). Here, \(\mathbf{v}_\parallel\) is the velocity parallel to the magnetic field, \(\mathbf{v}_L\) is the Larmor gyration velocity in the absence of \(\mathbf{F}\) found previously, and \(\mathbf{v}_d\) is some constant drift velocity perpendicular to the magnetic field. Then, we find that
and applying the vector triple product yields
In the case where the external force \(\mathbf{F} = q\mathbf{E}\) is due to a constant electric field, this is the constant \(\mathbf{E}\times\mathbf{B}\) drift velocity:
Built in drift functions allow you to account for the new force added to the system in two different ways:
[5]:
E = 0.2 * u.V / u.m # Efield magnitude
ey = np.array([0, 1, 0])
ez = np.array([0, 0, 1])
F = proton.charge * E # force due to Efield
v_d = formulary.drifts.force_drift(F * ey, B * ez, proton.charge)
print("F drift velocity: ", v_d)
v_d = formulary.drifts.ExB_drift(E * ey, B * ez)
print("ExB drift velocity: ", v_d)
F drift velocity: [0.04 0. 0. ] m / s
ExB drift velocity: [0.04 0. 0. ] m / s
The resulting particle trajectory can be compared to the case without drifts by calling our previously defined function with the drift velocity now as an argument. As expected, there is a constant drift in the direction of \(\mathbf{E}\times\mathbf{B}\):
[6]:
x_d, y_d, z_d = single_particle_trajectory(v_d=v_d)
fig = plt.figure()
ax = fig.gca(projection="3d")
ax.plot(x, y, z, label="$\mathbf{F}=0$")
ax.plot(x_d, y_d, z_d, label="$\mathbf{F}=q\mathbf{E}$")
ax.legend()
bound = 3 * r_L.value
ax.set_xlim([bound, bound])
ax.set_ylim([bound, bound])
ax.set_zlim([0, 10])
ax.set_xlabel("x [m]")
ax.set_ylabel("y [m]")
ax.set_zlabel("z [m]")
plt.show()
print(f"r_L = {r_L.value:.2e} m")
print(f"omega_c = {omega_c.value:.2e} rads/s")
r_L = 2.09e09 m
omega_c = 4.79e+08 rads/s
Of course, the implementation in our single_particle_trajectory()
function requires the analytical solution for the velocity \(\mathbf{v}_d\). This solution can be compared with that implemented in the particle stepper notebook. It uses the Boris algorithm to evolve the particle along its trajectory in prescribed \(\mathbf{E}\) and \(\mathbf{B}\) fields, and thus does not require the analytical solution.
This page was generated by
nbsphinx from
docs/notebooks/cold_plasma_tensor_elements.ipynb.
Interactive online version:
.
[1]:
%matplotlib inline
Cold Magnetized Plasma Waves Tensor Elements (S, D, P in Stix’s notation)¶
This example shows how to calculate the values of the cold plasma tensor elements for various electromagnetic wave frequencies.
[2]:
import matplotlib.pyplot as plt
import numpy as np
from astropy import units as u
from plasmapy.formulary import (
cold_plasma_permittivity_LRP,
cold_plasma_permittivity_SDP,
)
Take a look at the docs to cold_plasma_permittivity_SDP()
and cold_plasma_permittivity_LRP()
for more information on these.
Let’s define some parameters, such as the magnetic field magnitude, the plasma species and densities and the frequency band of interest
[3]:
B = 2 * u.T
species = ["e", "D+"]
n = [1e18 * u.m ** 3, 1e18 * u.m ** 3]
f = np.logspace(start=6, stop=11.3, num=3001) # 1 MHz to 200 GHz
omega_RF = f * (2 * np.pi) * (u.rad / u.s)
[4]:
help(cold_plasma_permittivity_SDP)
Help on function cold_plasma_permittivity_SDP in module plasmapy.formulary.dielectric:
cold_plasma_permittivity_SDP(B: Unit("T"), species, n, omega: Unit("rad / s"))
Magnetized cold plasma dielectric permittivity tensor elements.
Elements (S, D, P) are given in the "Stix" frame, i.e. with
:math:`B ∥ \hat{z}`\ .
The :math:`\exp(i ω t)` timeharmonic convention is assumed.
Parameters

B : `~astropy.units.Quantity`
Magnetic field magnitude in units convertible to tesla.
species : `list` of `str`
List of the plasma particle species,
e.g.: ``['e', 'D+']`` or ``['e', 'D+', 'He+']``.
n : `list` of `~astropy.units.Quantity`
`list` of species density in units convertible to per cubic meter
The order of the species densities should follow species.
omega : `~astropy.units.Quantity`
Electromagnetic wave frequency in rad/s.
Returns

sum : `~astropy.units.Quantity`
S ("Sum") dielectric tensor element.
difference : `~astropy.units.Quantity`
D ("Difference") dielectric tensor element.
plasma : `~astropy.units.Quantity`
P ("Plasma") dielectric tensor element.
Notes

The dielectric permittivity tensor is expressed in the Stix frame with
the :math:`\exp(i ω t)` timeharmonic convention as
:math:`ε = ε_0 A`, with :math:`A` being
.. math::
ε = ε_0 \left(\begin{matrix} S & i D & 0 \\
+i D & S & 0 \\
0 & 0 & P \end{matrix}\right)
where:
.. math::
S = 1  \sum_s \frac{ω_{p,s}^2}{ω^2  Ω_{c,s}^2}
D = \sum_s \frac{Ω_{c,s}}{ω}
\frac{ω_{p,s}^2}{ω^2  Ω_{c,s}^2}
P = 1  \sum_s \frac{ω_{p,s}^2}{ω^2}
where :math:`ω_{p,s}` is the plasma frequency and
:math:`Ω_{c,s}` is the signed version of the cyclotron frequency
for the species :math:`s`.
References

 T.H. Stix, Waves in Plasma, 1992.
Examples

>>> from astropy import units as u
>>> from numpy import pi
>>> B = 2*u.T
>>> species = ['e', 'D+']
>>> n = [1e18*u.m**3, 1e18*u.m**3]
>>> omega = 3.7e9*(2*pi)*(u.rad/u.s)
>>> permittivity = S, D, P = cold_plasma_permittivity_SDP(B, species, n, omega)
>>> S
<Quantity 1.02422...>
>>> permittivity.sum # namedtuplestyle access
<Quantity 1.02422...>
>>> D
<Quantity 0.39089...>
>>> P
<Quantity 4.8903...>
[5]:
S, D, P = cold_plasma_permittivity_SDP(B, species, n, omega_RF)
Filter positive and negative values, for display purposes only. Still for display purposes, replace 0 by NaN to NOT plot 0 values
[6]:
S_pos = S * (S > 0)
D_pos = D * (D > 0)
P_pos = P * (P > 0)
S_neg = S * (S < 0)
D_neg = D * (D < 0)
P_neg = P * (P < 0)
S_pos[S_pos == 0] = np.NaN
D_pos[D_pos == 0] = np.NaN
P_pos[P_pos == 0] = np.NaN
S_neg[S_neg == 0] = np.NaN
D_neg[D_neg == 0] = np.NaN
P_neg[P_neg == 0] = np.NaN
[7]:
plt.figure(figsize=(12, 6))
plt.semilogx(f, abs(S_pos), f, abs(D_pos), f, abs(P_pos), lw=2)
plt.semilogx(
f,
abs(S_neg),
"#1f77b4",
f,
abs(D_neg),
"#ff7f0e",
f,
abs(P_neg),
"#2ca02c",
lw=2,
ls="",
)
plt.yscale("log")
plt.grid(True, which="major")
plt.grid(True, which="minor")
plt.ylim(1e4, 1e8)
plt.xlim(1e6, 200e9)
plt.legend(("S > 0", "D > 0", "P > 0", "S < 0", "D < 0", "P < 0"), fontsize=16, ncol=2)
plt.xlabel("RF Frequency [Hz]", size=16)
plt.ylabel("Absolute value", size=16)
plt.tick_params(labelsize=14)
Cold Plasma tensor elements in the rotating basis
[8]:
L, R, P = cold_plasma_permittivity_LRP(B, species, n, omega_RF)
[9]:
L_pos = L * (L > 0)
R_pos = R * (R > 0)
L_neg = L * (L < 0)
R_neg = R * (R < 0)
L_pos[L_pos == 0] = np.NaN
R_pos[R_pos == 0] = np.NaN
L_neg[L_neg == 0] = np.NaN
R_neg[R_neg == 0] = np.NaN
plt.figure(figsize=(12, 6))
plt.semilogx(f, abs(L_pos), f, abs(R_pos), f, abs(P_pos), lw=2)
plt.semilogx(
f,
abs(L_neg),
"#1f77b4",
f,
abs(R_neg),
"#ff7f0e",
f,
abs(P_neg),
"#2ca02c",
lw=2,
ls="",
)
plt.yscale("log")
plt.grid(True, which="major")
plt.grid(True, which="minor")
plt.xlim(1e6, 200e9)
plt.legend(("L > 0", "R > 0", "P > 0", "L < 0", "R < 0", "P < 0"), fontsize=16, ncol=2)
plt.xlabel("RF Frequency [Hz]", size=16)
plt.ylabel("Absolute value", size=16)
plt.tick_params(labelsize=14)
Checks if the values obtained are coherent. They should satisfy S = (R+L)/2 and D = (RL)/2
[10]:
try:
np.testing.assert_allclose(S, (R + L) / 2)
np.testing.assert_allclose(D, (R  L) / 2)
except AssertionError as e:
print(e)
# Checks for R=S+D and L=SD
try:
np.testing.assert_allclose(R, S + D)
np.testing.assert_allclose(L, S  D)
except AssertionError as e:
print(e)
This page was generated by
nbsphinx from
docs/notebooks/distribution.ipynb.
Interactive online version:
.
[1]:
%matplotlib inline
1D Maxwellian distribution function¶
We import the usual modules, and the hero of this notebook, the Maxwellian 1D distribution:
[2]:
import matplotlib.pyplot as plt
import numpy as np
from astropy import units as u
from astropy.constants import k_B, m_e
from plasmapy.formulary import Maxwellian_1D
Take a look at the docs to Maxwellian_1D()
for more information on these.
Given we’ll be plotting, import astropy’s quantity support:
[3]:
from astropy.visualization import quantity_support
quantity_support()
[3]:
<astropy.visualization.units.quantity_support.<locals>.MplQuantityConverter at 0x7f0273636dd0>
As a first example, let’s get the probability density of finding an electron with a speed of 1 m/s if we have a plasma at a temperature of 30 000 K:
[4]:
p_dens = Maxwellian_1D(
v=1 * u.m / u.s, T=30000 * u.K, particle="e", v_drift=0 * u.m / u.s
)
print(p_dens)
5.916328704912824e07 s / m
Note the units! Integrated over speed, this will give us a probability. Let’s test that for a bunch of particles:
[5]:
T = 3e4 * u.K
dv = 10 * u.m / u.s
v = np.arange(5e6, 5e6, 10) * u.m / u.s
Check that the integral over all speeds is 1 (the particle has to be somewhere):
[6]:
for particle in ["p", "e"]:
pdf = Maxwellian_1D(v, T=T, particle=particle)
integral = (pdf).sum() * dv
print(f"Integral value for {particle}: {integral}")
plt.plot(v, pdf, label=particle)
plt.legend()
Integral value for p: 1.0000000000000002
Integral value for e: 0.9999999999998787
[6]:
<matplotlib.legend.Legend at 0x7f0270cd2f10>
The standard deviation of this distribution should give us back the temperature:
[7]:
std = np.sqrt((Maxwellian_1D(v, T=T, particle="e") * v ** 2 * dv).sum())
T_theo = (std ** 2 / k_B * m_e).to(u.K)
print("T from standard deviation:", T_theo)
print("Initial T:", T)
T from standard deviation: 29999.999999792235 K
Initial T: 30000.0 K
This page was generated by
nbsphinx from
docs/notebooks/magnetic_statics.ipynb.
Interactive online version:
.
[1]:
%matplotlib inline
Magnetostatic Fields¶
An example of using PlasmaPy’s Magnetostatic
class in physics
subpackage.
[2]:
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from plasmapy.formulary import magnetostatics
from plasmapy.plasma.sources import Plasma3D
Some common magnetostatic fields can be generated and added to a plasma object. A dipole
[3]:
dipole = magnetostatics.MagneticDipole(
np.array([0, 0, 1]) * u.A * u.m * u.m, np.array([0, 0, 0]) * u.m
)
print(dipole)
MagneticDipole(moment=[0. 0. 1.]A m2, p0=[0. 0. 0.]m)
initialize a a plasma, where the magnetic field will be calculated on
[4]:
plasma = Plasma3D(
domain_x=np.linspace(2, 2, 30) * u.m,
domain_y=np.linspace(0, 0, 1) * u.m,
domain_z=np.linspace(2, 2, 20) * u.m,
)
add the dipole field to it
[5]:
plasma.add_magnetostatic(dipole)
X, Z = plasma.grid[0, :, 0, :], plasma.grid[2, :, 0, :]
U = plasma.magnetic_field[0, :, 0, :].value.T # because grid uses 'ij' indexing
W = plasma.magnetic_field[2, :, 0, :].value.T # because grid uses 'ij' indexing
[6]:
plt.figure()
plt.axis("square")
plt.xlim(2, 2)
plt.ylim(2, 2)
plt.title(
"Dipole field in xz plane, generated by a dipole pointing in the z direction"
)
plt.streamplot(plasma.x.value, plasma.z.value, U, W)
[6]:
<matplotlib.streamplot.StreamplotSet at 0x7fac826a7a10>
A circular currentcarring wire (CircularWire
)
[7]:
cw = magnetostatics.CircularWire(
np.array([0, 0, 1]), np.array([0, 0, 0]) * u.m, 1 * u.m, 1 * u.A
)
print(cw)
CircularWire(normal=[0. 0. 1.], center=[0. 0. 0.]m, radius=1.0m, current=1.0A)
initialize a a plasma, where the magnetic field will be calculated on
[8]:
plasma = Plasma3D(
domain_x=np.linspace(2, 2, 30) * u.m,
domain_y=np.linspace(0, 0, 1) * u.m,
domain_z=np.linspace(2, 2, 20) * u.m,
)
add the circular coil field to it
[9]:
plasma.add_magnetostatic(cw)
X, Z = plasma.grid[0, :, 0, :], plasma.grid[2, :, 0, :]
U = plasma.magnetic_field[0, :, 0, :].value.T # because grid uses 'ij' indexing
W = plasma.magnetic_field[2, :, 0, :].value.T # because grid uses 'ij' indexing
plt.figure()
plt.axis("square")
plt.xlim(2, 2)
plt.ylim(2, 2)
plt.title(
"Circular coil field in xz plane, generated by a circular coil in the xy plane"
)
plt.streamplot(plasma.x.value, plasma.z.value, U, W)
[9]:
<matplotlib.streamplot.StreamplotSet at 0x7fac80d8a950>
a circular wire can be described as parametric equation and converted to GeneralWire
[10]:
gw_cw = cw.to_GeneralWire()
# the calculated magnetic field is close
print(gw_cw.magnetic_field([0, 0, 0])  cw.magnetic_field([0, 0, 0]))
[ 0.00000000e+00 0.00000000e+00 4.13416205e12] T
A infinite straight wire (InfiniteStraightWire
)
[11]:
iw = magnetostatics.InfiniteStraightWire(
np.array([0, 1, 0]), np.array([0, 0, 0]) * u.m, 1 * u.A
)
print(iw)
InfiniteStraightWire(direction=[0. 1. 0.], p0=[0. 0. 0.]m, current=1.0A)
initialize a a plasma, where the magnetic field will be calculated on
[12]:
plasma = Plasma3D(
domain_x=np.linspace(2, 2, 30) * u.m,
domain_y=np.linspace(0, 0, 1) * u.m,
domain_z=np.linspace(2, 2, 20) * u.m,
)
# add the infinite straight wire field to it
plasma.add_magnetostatic(iw)
X, Z = plasma.grid[0, :, 0, :], plasma.grid[2, :, 0, :]
U = plasma.magnetic_field[0, :, 0, :].value.T # because grid uses 'ij' indexing
W = plasma.magnetic_field[2, :, 0, :].value.T # because grid uses 'ij' indexing
plt.figure()
plt.title(
"Dipole field in xz plane, generated by a infinite straight wire "
"pointing in the y direction"
)
plt.axis("square")
plt.xlim(2, 2)
plt.ylim(2, 2)
plt.streamplot(plasma.x.value, plasma.z.value, U, W)
[12]:
<matplotlib.streamplot.StreamplotSet at 0x7fac80a99790>
This page was generated by
nbsphinx from
docs/notebooks/physics.ipynb.
Interactive online version:
.
[1]:
%matplotlib inline
Analysing ITER parameters¶
Let’s try to look at ITER plasma conditions using the physics
subpackage.
[2]:
import matplotlib.pyplot as plt
import numpy as np
from astropy import units as u
from mpl_toolkits.mplot3d import Axes3D
from plasmapy import formulary
The radius of electric field shielding clouds, also known as the Debye_length()
,
would be
[3]:
electron_temperature = 8.8 * u.keV
electron_concentration = 10.1e19 / u.m ** 3
print(formulary.Debye_length(electron_temperature, electron_concentration))
6.939046810942984e05 m
Note that we can also neglect the unit for the concentration, as 1/m^3 is the a standard unit for this kind of Quantity:
[4]:
print(formulary.Debye_length(electron_temperature, 10.1e19))
6.939046810942984e05 m
WARNING: UnitsWarning: The argument 'n_e' to function Debye_length() has no specified units. Assuming units of 1 / m3. To silence this warning, explicitly pass in an astropy Quantity (e.g. 5. * astropy.units.cm) (see http://docs.astropy.org/en/stable/units/) [plasmapy.utils.decorators.validators]
Assuming the magnetic field as 5.3 Teslas (which is the value at the major radius):
[5]:
B = 5.3 * u.T
print(formulary.gyrofrequency(B, particle="e"))
print(formulary.gyroradius(B, T_i=electron_temperature, particle="e"))
932174605709.2465 rad / s
5.968562765183547e05 m
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/utils/decorators/checks.py:1400: RelativityWarning: thermal_speed is yielding a velocity that is 18.559% of the speed of light. Relativistic effects may be important.
RelativityWarning,
The electron inertial_length()
would be
[6]:
print(formulary.inertial_length(electron_concentration, particle="e"))
0.0005287720427518426 m
In these conditions, they should reach thermal velocities of about
[7]:
print(formulary.thermal_speed(T=electron_temperature, particle="e"))
55637426.422858626 m / s
And the Langmuir wave plasma_frequency()
should be on the order of
[8]:
print(formulary.plasma_frequency(electron_concentration, particle="e"))
566959736448.652 rad / s
Let’s try to recreate some plots and get a feel for some of these quantities. We will also compare our data to realworld plasma situations.
[9]:
n_e = np.logspace(4, 30, 100) / u.m ** 3
plt.plot(n_e, formulary.plasma_frequency(n_e, particle="e"))
plt.scatter(
electron_concentration,
formulary.plasma_frequency(electron_concentration, particle="e"),
label="Our Data",
)
# IRT1 Tokamak Data
# http://article.sapub.org/pdf/10.5923.j.jnpp.20110101.03.pdf
n_e = 1.2e19 / u.m ** 3
T_e = 136.8323 * u.eV
B = 0.82 * u.T
plt.scatter(n_e, formulary.plasma_frequency(n_e, particle="e"), label="IRT1 Tokamak")
# Wendelstein 7X Stellerator Data
# https://nucleus.iaea.org/sites/fusionportal/Shared%20Documents/FEC%202016/fec2016preprints/preprint0541.pdf
n_e = 3e19 / u.m ** 3
T_e = 6 * u.keV
B = 3 * u.T
plt.scatter(
n_e, formulary.plasma_frequency(n_e, particle="e"), label="W7X Stellerator"
)
# Solar Corona
# Estimated by Nick Murphy
n_e = 1e15 / u.m ** 3
T_e = 1e6 * u.K
B = 0.005 * u.T
T_e.to(u.eV, equivalencies=u.temperature_energy())
plt.scatter(n_e, formulary.plasma_frequency(n_e, particle="e"), label="Solar Corona")
# Interstellar (warm neutral) Medium
# Estimated by Nick Murphy
n_e = 1e6 / u.m ** 3
T_e = 5e3 * u.K
B = 0.005 * u.T
T_e.to(u.eV, equivalencies=u.temperature_energy())
plt.scatter(
n_e, formulary.plasma_frequency(n_e, particle="e"), label="Interstellar Medium"
)
# Solar Wind at 1 AU
# Estimated by Nick Murphy
n_e = 7e6 / u.m ** 3
T_e = 1e5 * u.K
B = 5e9 * u.T
T_e.to(u.eV, equivalencies=u.temperature_energy())
plt.scatter(
n_e, formulary.plasma_frequency(n_e, particle="e"), label="Solar Wind (1AU)"
)
plt.xlabel("Electron Concentration (m^3)")
plt.ylabel("Langmuir Wave Plasma Frequency (rad/s)")
plt.grid()
plt.xscale("log")
plt.yscale("log")
plt.legend()
plt.title("Logscale plot of plasma frequencies")
plt.show()
Analyses & Diagnostics¶
This page was generated by
nbsphinx from
docs/notebooks/analysis/fit_functions.ipynb.
Interactive online version:
.
Fit Functions¶
Fit functions are a set of callable classes designed to aid in fitting analytical functions to data. A fit function class combines the following functionality:
An analytical function that is callable with given parameters or fitted parameters.
Curve fitting functionality (usually SciPy’s curve_fit() or linregress()), which stores the fit statistics and parameters into the class. This makes the function easily callable with the fitted parameters.
Error propagation calculations.
A root solver that returns either the known analytical solutions or uses SciPy’s fsolve() to calculate the roots.
[1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from pathlib import Path
from plasmapy.analysis import fit_functions as ffuncs
plt.rcParams["figure.figsize"] = [10.5, 0.56 * 10.5]
Contents:¶
Fit function basics¶
There is an ever expanding collection of fit functions, but this notebook will use ExponentialPlusLinear as an example.
A fit function class has no required arguments at time of instantiation.
[2]:
# basic instantiation
explin = ffuncs.ExponentialPlusLinear()
# fit parameters are not set yet
(explin.params, explin.param_errors)
[2]:
(None, None)
Each fit parameter is given a name.
[3]:
explin.param_names
[3]:
('a', 'alpha', 'm', 'b')
These names are used throughout the fit function’s documentation, as well as in its __repr__
, __str__
, and latex_str
methods.
[4]:
(explin, explin.__str__(), explin.latex_str)
[4]:
(f(x) = a exp(alpha x) + m x + b <class 'plasmapy.analysis.fit_functions.ExponentialPlusLinear'>,
'f(x) = a exp(alpha x) + m x + b',
'a \\, \\exp(\\alpha x) + m x + b')
Fitting to data¶
Fit functions provide the curve_fit() method to fit the analytical function to a set of \((x, y)\) data. This is typically done with SciPy’s curve_fit() function, but fitting is done with SciPy’s linregress() for the Linear fit funciton.
Let’s generate some noisy data to fit to…
[5]:
params = (5.0, 0.1, 0.5, 8.0) # (a, alpha, m, b)
xdata = np.linspace(20, 15, num=100)
ydata = explin.func(xdata, *params) + np.random.normal(0.0, 0.6, xdata.size)
plt.plot(xdata, ydata)
plt.xlabel("X", fontsize=14)
plt.ylabel("Y", fontsize=14)
[5]:
Text(0, 0.5, 'Y')
The fit function curve_fit()
shares the same signature as SciPy’s curve_fit(), so any **kwargs
will be passed on. By default, only the \((x, y)\) values are needed.
[6]:
explin.curve_fit(xdata, ydata)
Getting fit results¶
After fitting, the fitted parameters, uncertainties, and coefficient of determination, or \(r^2\), values can be retrieved through their respective properties, params
, parame_errors
, and rsq
.
[7]:
(explin.params, explin.params.a, explin.params.alpha)
[7]:
(FitParamTuple(a=4.7437093495345, alpha=0.10300303124434122, m=0.49262801905579634, b=7.7073863670661975),
4.7437093495345,
0.10300303124434122)
[8]:
(explin.param_errors, explin.param_errors.a, explin.param_errors.alpha)
[8]:
(FitParamTuple(a=0.9294198931258074, alpha=0.008821749639446746, m=0.04296472145442814, b=0.9545154438753005),
0.9294198931258074,
0.008821749639446746)
[9]:
explin.rsq
[9]:
0.951997505392852
Fit function is callable¶
Now that parameters are set, the fit function is callable.
[10]:
explin(0)
[10]:
2.9636770175316975
Associated errors can also be generated.
[11]:
y, y_err = explin(np.linspace(1, 1, num=10), reterr=True)
(y, y_err)
[11]:
(array([2.93534316, 2.94573245, 2.9538276 , 2.9595755 , 2.9629218 ,
2.96381089, 2.96218588, 2.95798856, 2.95115938, 2.94163739]),
array([1.27175981, 1.28415417, 1.2971345 , 1.31071781, 1.32492153,
1.3397636 , 1.35526238, 1.37143675, 1.38830606, 1.40589019]))
Known uncertainties in \(x\) can be specified too.
[12]:
y, y_err = explin(np.linspace(1, 1, num=10), reterr=True, x_err=0.1)
(y, y_err)
[12]:
(array([2.93534316, 2.94573245, 2.9538276 , 2.9595755 , 2.9629218 ,
2.96381089, 2.96218588, 2.95798856, 2.95115938, 2.94163739]),
array([1.27177038, 1.28416092, 1.29713825, 1.31071941, 1.32492188,
1.33976361, 1.35526301, 1.37143899, 1.38831092, 1.40589873]))
Plotting results¶
[13]:
# plot original data
plt.plot(xdata, ydata, marker="o", linestyle=" ", label="Data")
ax = plt.gca()
ax.set_xlabel("X", fontsize=14)
ax.set_ylabel("Y", fontsize=14)
ax.axhline(0.0, color="r", linestyle="")
# plot fitted curve + error
yfit, yfit_err = explin(xdata, reterr=True)
ax.plot(xdata, yfit, color="orange", label="Fit")
ax.fill_between(
xdata,
yfit + yfit_err,
yfit  yfit_err,
color="orange",
alpha=0.12,
zorder=0,
label="Fit Error",
)
# plot annotations
plt.legend(fontsize=14, loc="upper left")
txt = f"$f(x) = {explin.latex_str}$\n" f"$r^2 = {explin.rsq:.3f}$\n"
for name, param, err in zip(explin.param_names, explin.params, explin.param_errors):
txt += f"{name} = {param:.3f} $\\pm$ {err:.3f}\n"
txt_loc = [13.0, ax.get_ylim()[1]]
txt_loc = ax.transAxes.inverted().transform(ax.transData.transform(txt_loc))
txt_loc[0] = 0.02
txt_loc[1] = 0.05
ax.text(
txt_loc[0],
txt_loc[1],
txt,
fontsize="large",
transform=ax.transAxes,
va="top",
linespacing=1.5,
)
[13]:
Text(0.2072727272727273, 0.9500000000000002, '$f(x) = a \\, \\exp(\\alpha x) + m x + b$\n$r^2 = 0.952$\na = 4.744 $\\pm$ 0.929\nalpha = 0.103 $\\pm$ 0.009\nm = 0.493 $\\pm$ 0.043\nb = 7.707 $\\pm$ 0.955\n')
Root solving¶
An exponential plus a linear offset has no analytical solutions for its roots, except for a few specific cases. To get around this, ExponentialPlusLinear().root_solve() uses SciPy’s fsolve() to calculate it’s roots. If a fit function has analytical solutions to its roots (e.g. Linear().root_solve()), then the method is overriden with the known solution.
[14]:
root, err = explin.root_solve(15.0)
(root, err)
[14]:
(13.163847024274588, nan)
Let’s use Linear().root_solve() as an example for a known solution.
[15]:
lin = ffuncs.Linear(params=(1.0, 5.0), param_errors=(0.1, 0.1))
root, err = lin.root_solve()
(root, err)
[15]:
(5.0, 0.5099019513592785)
This page was generated by
nbsphinx from
docs/notebooks/analysis/swept_langmuir/find_floating_potential.ipynb.
Interactive online version:
.
Swept Langmuir Analysis: Floating Potential¶
This notebook covers the use of the **find_floating_potential()** function and how it is used to determine the floating potential from a swept Langmuir trace.
The floating potential, \(V_f\), is defined as the probe bias voltage at which there is no net collected current, \(I=0\). This occurs because the floating potential slows the collected electrons and accelerates the collected ions to a point where the electron and ioncurrents balance each other out.
[1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pprint
from pathlib import Path
from plasmapy.analysis import swept_langmuir as sla
plt.rcParams["figure.figsize"] = [10.5, 0.56 * 10.5]
Contents:¶
How find_floating_potential()
works¶
The passed current array is scanned for points that equal zero and pointpairs that straddle where the current, \(I\), equals zero. This forms an collection of “crossingpoints.”
The crossingpoints are then grouped into “crossingislands” based on the
threshold
keyword.A new island is formed when a successive crossingpoint is more (index) steps away from the previous crossingpoint than defined by
threshold
. For example, ifthreshold=4
then an new island is formed if a crossingpoint candidate is more than 4 steps away from the previous candidate.If multiple crossingislands are identified, then the function will compare the total span of all crossingislands to
min_points
. If the span is greater thanmin_points
, then the function is incapable of identifying \(V_f\) and will returnnumpy.nan
values; otherwise, the span will form one larger crossingisland.
To calculate the floating potential…
If the number of points that make up the crossingisland is less than
min_points
, then each side of the “crossingisland” is equally padded with the nearest neighbor points untilmin_points
is satisfied.If
fit_type="linear"
, then ascipy.stats.linregress
fit is applied to the points that make up the crossingisland.If
fit_type="exponential"
, then ascipy.optimize.curve_fit
fit is applied to the points that make up the crossingisland.
Notes about usage¶
The function provides no signal processing. If needed, the user must smooth, sort, crop, or process the arrays before passing them to the function.
The function requires the voltage array to be monotonically increasing or decreasing.
If the total range spanned by all crossingislands is less than or equal to
min_points
, thenthreshold
is ignored and all crossingislands are grouped into one island.
Knobs to turn¶
fit_type
There are two types of curves that can be fitted to the identified crossing point data:
"linear"
and"exponential"
. The former will fit a line to the data, whereas, the latter will fit an exponential curve with an offset. The default curve is"exponential"
since swept Langmuir data is not typically linear as it passes through \(I=0\).min_points
This variable specifies the minimum number of points that will be used in the curve fitting. As mentioned above, the crossingislands are identified and then padded until
min_points
is satisfied.min_pints = None
(Default) then the larger of 5 andfactor * array_size
is taken, wherefactor = 0.1
for"linear"
and0.2
for"exponential"
.min_points = numpy.inf
then the entire passed array is fitted.min_points >= 1
then this is the minimum number of points used.0 < min_points < 1
then then the minimum number of points is taken asmin_points * array_size
.
threshold
The max allowed index distance between crossingpoints before a new crossingisland is formed.
Calculate the Floating Potential¶
Below we’ll compute the floating potential using the default fitting behavior (fit_type="exponential"
) and a linear fit (fit_type="linear"
).
[2]:
# load data
filename = "Beckers2017_noisy.npy"
filepath = (Path.cwd() / ".." / ".." / "langmuir_samples" / filename).resolve()
voltage, current = np.load(filepath)
# voltage array needs to be monotonically increasing/decreasing
isort = np.argsort(voltage)
voltage = voltage[isort]
current = current[isort]
# get default fit results (exponential fit)
results = sla.find_floating_potential(voltage, current, min_points=0.3)
# get linear fit results
results_lin = sla.find_floating_potential(voltage, current, fit_type="linear")
Interpreting results¶
The find_floating_potential()
function returns a six element named tuple, where…
results[0]
=results.vf
= the determined floating potential (same units as the passvoltage
array)
[3]:
(results[0], results.vf)
[3]:
(5.665729883017972, 5.665729883017972)
results[1]
=results.vf_err
= the associated uncertainty in the \(V_F\) calculation (same units asresults.vf
)
[4]:
(results[1], results.vf_err)
[4]:
(0.45174046877528606, 0.45174046877528606)
results[2]
=results.rsq
= the coefficient of determination (rsquared) value of the fit
[5]:
(results[2], results.rsq)
[5]:
(0.9564766142728285, 0.9564766142728285)
results[3]
=results.func
= the resulting fitted functionresults.func
is a callable representation of the fitted functionI = results.func(V)
.resulst.func
is an instance of a subclass ofAbstractFitFunction
. (FitFuction classes)Since
results.func
is a class instance, there are many other attributes available. For example,results.func.params
is a named tuple of the fitted parametersresults.func.param_errors
is a named tuple of the fitted parameter errorsresults.func.root_solve()
finds the roots of the fitted function. This is how \(V_f\) is calculated.
[6]:
(
results[3],
results.func,
results.func.params,
results.func.params.a,
results.func.param_errors,
results.func.param_errors.a,
results.func(results.vf),
)
[6]:
(f(x) = a exp(alpha x) + b <class 'plasmapy.analysis.fit_functions.ExponentialPlusOffset'>,
f(x) = a exp(alpha x) + b <class 'plasmapy.analysis.fit_functions.ExponentialPlusOffset'>,
FitParamTuple(a=0.01027537751642042, alpha=0.3382781633131013, b=0.001511583538487724),
0.01027537751642042,
FitParamTuple(a=0.00031440591547328146, alpha=0.021631300563905463, b=0.00012999282477305467),
0.00031440591547328146,
4.336808689942018e19)
results[4]
=results.islands
= a list of slice objects representing all the indentified crossingislands
[7]:
(
results[4],
results.islands,
voltage[results.islands[0]],
)
[7]:
([slice(192, 195, None), slice(201, 210, None)],
[slice(192, 195, None), slice(201, 210, None)],
array([7.49600064, 7.22376048, 7.15120308]))
results[5]
=results.indices
= a slice object representing the indices used in the fit
[8]:
(
results[5],
results.indices,
voltage[results.indices],
)
[8]:
(slice(155, 247, None),
slice(155, 247, None),
array([11.65942357, 11.64882654, 11.53579538, 11.38604798,
11.11860447, 11.11179009, 11.09142646, 11.00483855,
10.99442158, 10.57322789, 10.50059601, 10.37949917,
10.315585 , 10.2698363 , 10.14162481, 9.97124346,
9.62883006, 9.5755686 , 9.36686515, 9.28450711,
9.07258268, 9.06560791, 8.91301953, 8.8817602 ,
8.76066758, 8.74794412, 8.74340073, 8.72406653,
8.55441808, 8.34046557, 8.2062713 , 8.13718384,
7.98570514, 7.68529681, 7.65010823, 7.55378266,
7.49960431, 7.49600064, 7.22376048, 7.15120308,
7.02994327, 6.77525529, 6.51507915, 6.39924154,
6.36984977, 6.24817453, 6.20213736, 6.19568807,
6.04787946, 5.76312984, 5.46965878, 5.32623331,
5.119681 , 5.0823194 , 5.0486728 , 4.95381421,
4.88749753, 4.78029289, 4.68512796, 4.51734484,
4.40678539, 4.19565637, 4.17704741, 4.17117347,
4.10283824, 4.00779631, 3.93855879, 3.82827448,
3.78163806, 3.27907775, 3.10329043, 3.09674414,
3.09107256, 2.96068724, 2.80465682, 2.60784067,
2.26018439, 2.17330689, 2.15893471, 1.98233593,
1.94282875, 1.91200309, 1.71890015, 1.60465752,
1.56380995, 1.38993179, 1.32157414, 1.30011675,
1.28811388, 1.25099997, 0.98118767, 0.85042501]))
Plotting results¶
[9]:
figwidth, figheight = plt.rcParams["figure.figsize"]
figheight = 2.0 * figheight
fig, axs = plt.subplots(3, 1, figsize=[figwidth, figheight])
# plot original data
axs[0].set_xlabel("Bias Voltage (V)", fontsize=12)
axs[0].set_ylabel("Current (A)", fontsize=12)
axs[0].plot(voltage, current, zorder=10, label="Sweep Data")
axs[0].axhline(0.0, color="r", linestyle="", label="I = 0")
axs[0].legend(fontsize=12)
# zoom on fit
for ii, label, fit in zip([1, 2], ["Exponential", "Linear"], [results, results_lin]):
# calc island points
isl_pts = np.array([], dtype=np.int64)
for isl in fit.islands:
isl_pts = np.concatenate((isl_pts, np.r_[isl]))
# calc xrange for plot
xlim = [voltage[fit.indices].min(), voltage[fit.indices].max()]
vpad = 0.25 * (xlim[1]  xlim[0])
xlim = [xlim[0]  vpad, xlim[1] + vpad]
# calc data points for fit curve
mask1 = np.where(voltage >= xlim[0], True, False)
mask2 = np.where(voltage <= xlim[1], True, False)
mask = np.logical_and(mask1, mask2)
vfit = np.linspace(xlim[0], xlim[1], 201, endpoint=True)
ifit, ifit_err = fit.func(vfit, reterr=True)
axs[ii].set_xlabel("Bias Voltage (V)", fontsize=12)
axs[ii].set_ylabel("Current (A)", fontsize=12)
axs[ii].set_xlim(xlim)
axs[ii].plot(
voltage[mask], current[mask], marker="o", zorder=10, label="Sweep Data",
)
axs[ii].scatter(
voltage[fit.indices],
current[fit.indices],
linewidth=2,
s=6 ** 2,
facecolors="deepskyblue",
edgecolors="deepskyblue",
zorder=11,
label="Points for Fit",
)
axs[ii].scatter(
voltage[isl_pts],
current[isl_pts],
linewidth=2,
s=8 ** 2,
facecolors="deepskyblue",
edgecolors="black",
zorder=12,
label="Island Points",
)
axs[ii].autoscale(False)
axs[ii].plot(vfit, ifit, color="orange", zorder=13, label=label + " Fit")
axs[ii].fill_between(
vfit,
ifit + ifit_err,
ifit  ifit_err,
color="orange",
alpha=0.12,
zorder=0,
label="Fit Error",
)
axs[ii].axhline(0.0, color="r", linestyle="")
axs[ii].fill_between(
[fit.vf  fit.vf_err, fit.vf + fit.vf_err],
axs[1].get_ylim()[0],
axs[1].get_ylim()[1],
color="grey",
alpha=0.1,
)
axs[ii].axvline(fit.vf, color="grey")
axs[ii].legend(fontsize=12)
# add text
rsq = fit.rsq
txt = f"$V_f = {fit.vf:.2f} \\pm {fit.vf_err:.2f}$ V\n"
txt += f"$r^2 = {rsq:.3f}$"
txt_loc = [fit.vf, axs[ii].get_ylim()[1]]
txt_loc = axs[ii].transData.transform(txt_loc)
txt_loc = axs[ii].transAxes.inverted().transform(txt_loc)
txt_loc[0] = 0.02
txt_loc[1] = 0.26
axs[ii].text(
txt_loc[0],
txt_loc[1],
txt,
fontsize="large",
transform=axs[ii].transAxes,
ha="right",
)
This page was generated by
nbsphinx from
docs/notebooks/diagnostics/langmuir_analysis.ipynb.
Interactive online version:
.
Langmuir probe data analysis¶
Let’s analyze a few Langmuir probe characteristics using the diagnostics.langmuir
subpackage. First we need to import the module and some basics.
[1]:
%matplotlib inline
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import os
from pathlib import Path
from pprint import pprint
from plasmapy.diagnostics.langmuir import Characteristic, swept_probe_analysis
The first characteristic we analyze is a simple singleprobe measurement in
a low (ion) temperature, low density plasma with a cylindrical probe. This
allows us to utilize OML theory implemented in swept_probe_analysis()
.
The data has been preprocessed with some smoothing, which allows us to obtain
a Electron Energy Distribution Function (EEDF) as well.
[2]:
# Load the bias and current values stored in the .p pickle file.
path = (Path.cwd() / ".." / "langmuir_samples" / "Beckers2017.npy").resolve()
bias, current = np.load(path)
# Create the Characteristic object, taking into account the correct units
characteristic = Characteristic(u.Quantity(bias, u.V), u.Quantity(current, u.A))
# Calculate the cylindrical probe surface area
probe_length = 1.145 * u.mm
probe_diameter = 1.57 * u.mm
probe_area = probe_length * np.pi * probe_diameter + np.pi * 0.25 * probe_diameter ** 2
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/diagnostics/langmuir.py:39: FutureWarning: The plasmapy.diagnostics.langmuir module will be deprecated in favor of the plasmapy.analysis.swept_langmuir subpackage and phased out over 2021. The plasmapy.analysis package was released in v0.5.0.
FutureWarning,
Now we can actually perform the analysis. Since the plasma is in Helium an ion mass number of 4 is entered. The results are visualized and the obtained EEDF is also shown.
[3]:
pprint(
swept_probe_analysis(
characteristic, probe_area, "He4+", visualize=True, plot_EEDF=True
)
)
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/diagnostics/langmuir.py:428: MatplotlibDeprecationWarning: The 'nonposy' parameter of __init__() has been renamed 'nonpositive' since Matplotlib 3.3; support for the old name will be dropped two minor releases later.
ax2.set_yscale("log", nonposy="clip")
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/diagnostics/langmuir.py:39: FutureWarning: The plasmapy.diagnostics.langmuir module will be deprecated in favor of the plasmapy.analysis.swept_langmuir subpackage and phased out over 2021. The plasmapy.analysis package was released in v0.5.0.
FutureWarning,
{'I_es': <Quantity 0.01873572 A>,
'I_is': <Quantity 0.00231536 A>,
'T_e': <Quantity 0.4309528 eV>,
'V_F': <Quantity 5.75541887 V>,
'V_P': <Quantity 2.15098227 V>,
'n_e': <Quantity 1.40397606e+17 1 / m3>,
'n_i': <Quantity 9.85345823e+17 1 / m3>,
'n_i_OML': <Quantity 1.69390014e+17 1 / m3>}
The cyan and yellow lines indicate the fitted electron and ion currents, respectively. The green line is the sum of these and agrees nicely with the data. This indicates a successful analysis.
The next sample probe data is provided by David Pace. It is also obtained from a low relatively ion temperature and density plasma, in Argon.
[4]:
# Load the data from a file and create the Characteristic object
path = (Path.cwd() / ".." / "langmuir_samples" / "Pace2015.npy").resolve()
bias, current = np.load(path)
characteristic = Characteristic(u.Quantity(bias, u.V), u.Quantity(current, u.A))
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/diagnostics/langmuir.py:39: FutureWarning: The plasmapy.diagnostics.langmuir module will be deprecated in favor of the plasmapy.analysis.swept_langmuir subpackage and phased out over 2021. The plasmapy.analysis package was released in v0.5.0.
FutureWarning,
Initially the electrons are assumed to be Maxwellian. To check this the fit of the electron growth region will be plotted.
[5]:
swept_probe_analysis(
characteristic,
0.738 * u.cm ** 2,
"Ar40 1+",
bimaxwellian=False,
plot_electron_fit=True,
)
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/diagnostics/langmuir.py:39: FutureWarning: The plasmapy.diagnostics.langmuir module will be deprecated in favor of the plasmapy.analysis.swept_langmuir subpackage and phased out over 2021. The plasmapy.analysis package was released in v0.5.0.
FutureWarning,
[5]:
{'V_P': <Quantity 16.4 V>,
'V_F': <Quantity 35.6 V>,
'I_es': <Quantity 0.00282382 A>,
'I_is': <Quantity 0.000129 A>,
'n_e': <Quantity 2.51287474e+15 1 / m3>,
'n_i': <Quantity 2.06008231e+16 1 / m3>,
'T_e': <Quantity 0.32266985 eV>,
'n_i_OML': <Quantity 3.07522416e+15 1 / m3>}
It can be seen that this plasma is slightly biMaxwellian, as there are two distinct slopes in the exponential section. The analysis is now performed with bimaxwellian set to True, which yields improved results.
[6]:
pprint(
swept_probe_analysis(
characteristic,
0.738 * u.cm ** 2,
"Ar40 1+",
bimaxwellian=True,
visualize=True,
plot_electron_fit=True,
)
)
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/diagnostics/langmuir.py:39: FutureWarning: The plasmapy.diagnostics.langmuir module will be deprecated in favor of the plasmapy.analysis.swept_langmuir subpackage and phased out over 2021. The plasmapy.analysis package was released in v0.5.0.
FutureWarning,
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/diagnostics/langmuir.py:428: MatplotlibDeprecationWarning: The 'nonposy' parameter of __init__() has been renamed 'nonpositive' since Matplotlib 3.3; support for the old name will be dropped two minor releases later.
ax2.set_yscale("log", nonposy="clip")
{'I_es': <Quantity 0.00282382 A>,
'I_is': <Quantity 0.000129 A>,
'T_e': <Quantity [3.146203 , 3.72242199] eV>,
'V_F': <Quantity 35.6 V>,
'V_P': <Quantity 16.4 V>,
'hot_fraction': 0.43911961622308765,
'n_e': <Quantity 7.74211513e+14 1 / m3>,
'n_i': <Quantity 6.34707102e+15 1 / m3>,
'n_i_OML': <Quantity 3.07522416e+15 1 / m3>}
The probe current resolution of the raw data is relatively poor, but the analysis still performs well in the ion current region. The biMaxwellian properties are not significant but do make a difference. Check this analysis without setting bimaxwellian
to True! This is reflected in the results, which indicate that the temperatures of the cold and hot electron population are indeed different, but relatively close.
This Helium plasma is fully biMaxwellian.
[7]:
# Import probe data and calculate probe surface area.
path = (Path.cwd() / ".." / "langmuir_samples" / "Beckers2017b.npy").resolve()
bias, current = np.load(path)
characteristic = Characteristic(u.Quantity(bias, u.V), u.Quantity(current, u.A))
probe_length = 1.145 * u.mm
probe_diameter = 1.57 * u.mm
probe_area = probe_length * np.pi * probe_diameter + np.pi * 0.25 * probe_diameter ** 2
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/diagnostics/langmuir.py:39: FutureWarning: The plasmapy.diagnostics.langmuir module will be deprecated in favor of the plasmapy.analysis.swept_langmuir subpackage and phased out over 2021. The plasmapy.analysis package was released in v0.5.0.
FutureWarning,
plot_electron_fit
is set to True to check the biMaxwellian properties. The fit converges nicely to the two slopes of the electron growth region.
[8]:
pprint(
swept_probe_analysis(
characteristic,
probe_area,
"He4+",
bimaxwellian=True,
plot_electron_fit=True,
visualize=True,
)
)
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/diagnostics/langmuir.py:39: FutureWarning: The plasmapy.diagnostics.langmuir module will be deprecated in favor of the plasmapy.analysis.swept_langmuir subpackage and phased out over 2021. The plasmapy.analysis package was released in v0.5.0.
FutureWarning,
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/diagnostics/langmuir.py:428: MatplotlibDeprecationWarning: The 'nonposy' parameter of __init__() has been renamed 'nonpositive' since Matplotlib 3.3; support for the old name will be dropped two minor releases later.
ax2.set_yscale("log", nonposy="clip")
{'I_es': <Quantity 0.02655063 A>,
'I_is': <Quantity 0.00080287 A>,
'T_e': <Quantity [1.42489166, 7.4511896 ] eV>,
'V_F': <Quantity 21.29773863 V>,
'V_P': <Quantity 2.42370446 V>,
'hot_fraction': 0.2374803685084104,
'n_e': <Quantity 7.72855712e+16 1 / m3>,
'n_i': <Quantity 1.32724647e+17 1 / m3>,
'n_i_OML': <Quantity 6.08613986e+16 1 / m3>}
This page was generated by
nbsphinx from
docs/notebooks/diagnostics/proton_radiography_particle_tracing_custom_source.ipynb.
Interactive online version:
.
Synthetic Radiographs with Custom Source Profiles¶
In real charged particle radiography experiments, the finite size and distribution of the particle source limits the resolution of the radiograph. Some realistic sources produce particles with a nonuniform angular distribution that then superimposes a large scale “source profile” on the radiograph. For these reasons, the SyntheticProtronRadiography particle tracing class allows users to specify their own initial particle positions and velocities. This example will demonstrate how to use this functionality to create a more realistic synthetic radiograph that includes the effects from a nonuniform, finite source profile.
[1]:
import astropy.constants as const
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import warnings
from mpl_toolkits.mplot3d import Axes3D
from plasmapy.diagnostics import proton_radiography as prad
from plasmapy.formulary.mathematics import rot_a_to_b
from plasmapy.particles import Particle
from plasmapy.plasma.grids import CartesianGrid
Contents¶
Creating Particles¶
In this example we will create a source of 1e5 protons with a 5% variance in energy, a nonuniform angular velocity distribution, and a finite size.
[2]:
nparticles = 1e5
particle = Particle("p+")
We will choose a setup in which the sourcedetector axis is parallel to the \(y\)axis.
[3]:
# define location of source and detector plane
source = (0 * u.mm, 10 * u.mm, 0 * u.mm)
detector = (0 * u.mm, 100 * u.mm, 0 * u.mm)
Creating the Initial Particle Velocities¶
We will create the source distribution by utilizing the method of separation of variables,
and separately define the distribution component for each independent variable, \(u(v)\), \(g(\theta)\), and \(h(\phi)\). For geometric convenience, we will generate the velocity vector distribution around the \(z\)axis and then rotate the final velocities to be parallel to the sourcedetector axis (in this case the \(y\)axis).
First we will create the orientation angles polar (\(\theta\)) and azimuthal (\(\phi\)) for each particle. Generating \(\phi\) is simple: we will choose the azimuthal angles to just be uniformly distributed
[4]:
phi = np.random.uniform(high=2 * np.pi, size=int(nparticles))
However, choosing \(\theta\) is more complicated. Since the solid angle \(d\Omega = sin \theta d\theta d\phi\), if we draw a uniform distribution of \(\theta\) we will create a nonuniform distribution of particles in solid angle. This will create a sharp central peak on the detector plane.
[5]:
theta = np.random.uniform(high=np.pi / 2, size=int(nparticles))
fig, ax = plt.subplots(figsize=(6, 6))
theta_per_sa, bins = np.histogram(theta, bins=100, weights=1 / np.sin(theta))
ax.set_xlabel("$\\theta$ (rad)", fontsize=14)
ax.set_ylabel("N/N$_0$ per d$\\Omega$", fontsize=14)
ax.plot(bins[1:], theta_per_sa / np.sum(theta_per_sa))
ax.set_title(f"N$_0$ = {nparticles:.0e}", fontsize=14)
ax.set_yscale("log")
ax.set_xlim(0, np.pi / 2)
ax.set_ylim(None, 1)
[5]:
(0.0005365284988101142, 1)
To create a uniform distribution in solid angle, we need to draw values of \(\theta\) with a probability distribution weighted by \(\sin \theta\). This can be done using the np.random.choice() function, which draws size
elements from a distribution arg
with a probability distribution prob
. Setting the replace
keyword allows the same arguments to be drawn multiple times.
[6]:
arg = np.linspace(0, np.pi / 2, num=int(1e5))
prob = np.sin(arg)
prob *= 1 / np.sum(prob)
theta = np.random.choice(arg, size=int(nparticles), replace=True, p=prob)
fig, ax = plt.subplots(figsize=(6, 6))
theta_per_sa, bins = np.histogram(theta, bins=100, weights=1 / np.sin(theta))
ax.plot(bins[1:], theta_per_sa / np.sum(theta_per_sa))
ax.set_xlabel("$\\theta$ (rad)", fontsize=14)
ax.set_ylabel("N/N$_0$ per d$\\Omega$", fontsize=14)
ax.set_title(f"N$_0$ = {nparticles:.0e}", fontsize=14)
ax.set_yscale("log")
ax.set_xlim(0, np.pi / 2)
ax.set_ylim(None, 0.1)
[6]:
(0.007958407347562273, 0.1)
Now that we have a \(\theta\) distribution that is uniform in solid angle, we can perturb it by adding additional factors to the probability distribution used in np.random.choice(). For this case, let’s create a Gaussian distribution in solid angle.
Since particles moving at large angles will not be seen in the synthetic radiograph, we will set an upper bound \(\theta_{max}\) on the argument here. This is equivalent to setting the max_theta
keyword in create_particles()
[7]:
arg = np.linspace(0, np.pi / 8, num=int(1e5))
prob = np.sin(arg) * np.exp((arg ** 2) / 0.1 ** 2)
prob *= 1 / np.sum(prob)
theta = np.random.choice(arg, size=int(nparticles), replace=True, p=prob)
fig, ax = plt.subplots(figsize=(6, 6))
theta_per_sa, bins = np.histogram(theta, bins=100, weights=1 / np.sin(theta))
ax.plot(bins[1:], theta_per_sa / np.sum(theta_per_sa))
ax.set_title(f"N$_0$ = {nparticles:.0e}", fontsize=14)
ax.set_xlabel("$\\theta$ (rad)", fontsize=14)
ax.set_ylabel("N/N$_0$ per d$\\Omega$", fontsize=14)
ax.set_yscale("log")
ax.set_xlim(0, np.pi / 2)
ax.set_ylim(None, 1)
[7]:
(1.023799843394276e06, 1)
Now that the angular distributions are done, we will determine the energy (speed) for each particle. For this example, we will assume that the particle energy distribution is not a function of angle. We will create a Gaussian distribution of speeds with ~5% variance centered on a particle energy of 15 MeV.
[8]:
v_cent = np.sqrt(2 * 15 * u.MeV / particle.mass).to(u.m / u.s).value
v0 = np.random.normal(loc=v_cent, scale=1e6, size=int(nparticles))
v0 *= u.m / u.s
fig, ax = plt.subplots(figsize=(6, 6))
v_per_bin, bins = np.histogram(v0.si.value, bins=100)
ax.plot(bins[1:], v_per_bin / np.sum(v_per_bin))
ax.set_title(f"N$_0$ = {nparticles:.0e}", fontsize=14)
ax.set_xlabel("v0 (m/s)", fontsize=14)
ax.set_ylabel("N/N$_0$", fontsize=14)
ax.axvline(x=1.05 * v_cent, label="+5%", color="C1")
ax.axvline(x=0.95 * v_cent, label="5%", color="C2")
ax.legend(fontsize=14, loc="upper right")
[8]:
<matplotlib.legend.Legend at 0x7f91a98debd0>
Next, we will construct velocity vectors centered around the zaxis for each particle.
[9]:
vel = np.zeros([int(nparticles), 3]) * u.m / u.s
vel[:, 0] = v0 * np.sin(theta) * np.cos(phi)
vel[:, 1] = v0 * np.sin(theta) * np.sin(phi)
vel[:, 2] = v0 * np.cos(theta)
Finally, we will use the function rot_a_to_b() to create a rotation matrix that will rotate the vel
distribution so the distribution is centered about the \(y\) axis instead of the \(z\) axis.
[10]:
a = np.array([0, 0, 1])
b = np.array([0, 1, 0])
R = rot_a_to_b(a, b)
vel = np.matmul(vel, R)
Since the velocity vector distribution should be symmetric about the \(y\) axis, we can confirm this by checking that the normalized average velocity vector is close to the \(y\) unit vector.
[11]:
avg_v = np.mean(vel, axis=0)
print(avg_v / np.linalg.norm(avg_v))
[2.80752080e04 9.99999960e01 3.81821769e05]
Creating the Initial Particle Positions¶
For this example, we will create an initial position distribution representing a laser spot centered on the source location defined above as source
. The distribution will be cylindrical (oriented along the \(y\)axis) with a uniform distribution in y and a Gaussian distribution in radius (in the xz plane). We therefore need to create distributions in \(y\), \(\theta\), and \(r\), and then transform those into Cartesian positions.
Just as we previously weighted the \(\theta\) distribution with a \(sin \theta\) probability distribution to generate a uniform distribution in solid angle, we need to weight the \(r\) distribution with a \(r\) probability distribution so that the particles are uniformly distributed over the area of the disk.
[12]:
dy = 300 * u.um
y = np.random.uniform(
low=(source[1]  dy).to(u.m).value,
high=(source[1] + dy).to(u.m).value,
size=int(nparticles),
)
arg = np.linspace(1e9, 1e3, num=int(1e5))
prob = arg * np.exp(((arg / 3e4) ** 2))
prob *= 1 / np.sum(prob)
r = np.random.choice(arg, size=int(nparticles), replace=True, p=prob)
theta = np.random.uniform(low=0, high=2 * np.pi, size=int(nparticles))
x = r * np.cos(theta)
z = r * np.sin(theta)
hist, xpos, zpos = np.histogram2d(
x * 1e6, z * 1e6, bins=[100, 100], range=np.array([[5e2, 5e2], [5e2, 5e2]])
)
hist2, xpos2, ypos = np.histogram2d(
x * 1e6,
(y  source[1].to(u.m).value) * 1e6,
bins=[100, 100],
range=np.array([[5e2, 5e2], [5e2, 5e2]]),
)
fig, ax = plt.subplots(ncols=2, figsize=(12, 6))
fig.subplots_adjust(wspace=0.3, right=0.8)
fig.suptitle("Initial Particle Position Distribution", fontsize=14)
vmax = np.max([np.max(hist), np.max(hist2)])
p1 = ax[0].pcolormesh(xpos, zpos, hist.T, vmax=vmax)
ax[0].set_xlabel("x ($\\mu m$)", fontsize=14)
ax[0].set_ylabel("z ($\\mu m$)", fontsize=14)
ax[0].set_aspect("equal")
p2 = ax[1].pcolormesh(xpos2, ypos, hist2.T, vmax=vmax)
ax[1].set_xlabel("x ($\\mu m$)", fontsize=14)
ax[1].set_ylabel("y  $y_0$ ($\\mu m$)", fontsize=14)
ax[1].set_aspect("equal")
cbar_ax = fig.add_axes([0.85, 0.2, 0.03, 0.6])
cbar_ax.set_title("# Particles")
fig.colorbar(p2, cax=cbar_ax)
[12]:
<matplotlib.colorbar.Colorbar at 0x7f91a9733b10>
Finally we will combine these position arrays into an array with units.
[13]:
pos = np.zeros([int(nparticles), 3]) * u.m
pos[:, 0] = x * u.m
pos[:, 1] = y * u.m
pos[:, 2] = z * u.m
Creating a Synthetic Radiograph¶
To create an example synthetic radiograph, we will first create a field grid representing the analytical electric field produced by a sphere of Gaussian potential.
[14]:
# Create a Cartesian grid
L = 1 * u.mm
grid = CartesianGrid(L, L, num=100)
# Create a spherical potential with a Gaussian radial distribution
radius = np.linalg.norm(grid.grid, axis=3)
arg = (radius / (L / 3)).to(u.dimensionless_unscaled)
potential = 6e5 * np.exp((arg ** 2)) * u.V
# Calculate E from the potential
Ex, Ey, Ez = np.gradient(potential, grid.dax0, grid.dax1, grid.dax2)
mask = radius < L / 2
Ex = np.where(mask, Ex, 0)
Ey = np.where(mask, Ey, 0)
Ez = np.where(mask, Ez, 0)
# Add those quantities to the grid
grid.add_quantities(E_x=Ex, E_y=Ey, E_z=Ez, phi=potential)
# Plot the Efield
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection="3d")
ax.view_init(30, 30)
# skip some points to make the vector plot intelligable
s = tuple([slice(None, None, 6)] * 3)
ax.quiver(
grid.pts0[s].to(u.mm).value,
grid.pts1[s].to(u.mm).value,
grid.pts2[s].to(u.mm).value,
grid["E_x"][s],
grid["E_y"][s],
grid["E_z"][s],
length=5e7,
)
ax.set_xlabel("X (mm)", fontsize=14)
ax.set_ylabel("Y (mm)", fontsize=14)
ax.set_zlabel("Z (mm)", fontsize=14)
ax.set_xlim(1, 1)
ax.set_ylim(1, 1)
ax.set_zlim(1, 1)
ax.set_title("Gaussian Potential Electric Field", fontsize=14)
[14]:
Text(0.5, 0.92, 'Gaussian Potential Electric Field')
We will then create the synthetic radiograph object. The warning filter ignores a warning that arises because \(B_x\), \(B_y\), \(B_z\) are not provided in the grid (they will be assumed to be zero).
[15]:
with warnings.catch_warnings():
warnings.simplefilter("ignore")
sim = prad.SyntheticProtonRadiograph(grid, source, detector, verbose=False)
Now, instead of using create_particles() to create the particle distribution, we will use the load_particles() function to use the particles we have created above.
[16]:
sim.load_particles(pos, vel, particle=particle)
Now the particle radiograph simulation can be run as usual.
[17]:
sim.run()
[18]:
size = np.array([[1, 1], [1, 1]]) * 1.5 * u.cm
bins = [200, 200]
hax, vax, intensity = sim.synthetic_radiograph(size=size, bins=bins)
fig, ax = plt.subplots(figsize=(8, 8))
plot = ax.pcolormesh(
hax.to(u.cm).value, vax.to(u.cm).value, intensity.T, cmap="Blues_r", shading="auto",
)
cb = fig.colorbar(plot)
cb.ax.set_ylabel("Intensity", fontsize=14)
ax.set_aspect("equal")
ax.set_xlabel("X (cm), Image plane", fontsize=14)
ax.set_ylabel("Z (cm), Image plane", fontsize=14)
ax.set_title("Synthetic Proton Radiograph", fontsize=14);
[18]:
Text(0.5, 1.0, 'Synthetic Proton Radiograph')
Calling the synthetic_radiograph() function with the ignore_grid
keyword will produce the synthetic radiograph corresponding to the source profile propagated freely through space (i.e. in the absence of any grid fields).
[19]:
hax, vax, intensity = sim.synthetic_radiograph(size=size, bins=bins, ignore_grid=True)
fig, ax = plt.subplots(figsize=(8, 8))
plot = ax.pcolormesh(
hax.to(u.cm).value, vax.to(u.cm).value, intensity.T, cmap="Blues_r", shading="auto",
)
cb = fig.colorbar(plot)
cb.ax.set_ylabel("Intensity", fontsize=14)
ax.set_aspect("equal")
ax.set_xlabel("X (cm), Image plane", fontsize=14)
ax.set_ylabel("Z (cm), Image plane", fontsize=14)
ax.set_title("Source Profile", fontsize=14)
[19]:
Text(0.5, 1.0, 'Source Profile')
[20]:
fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(hax.to(u.cm).value, np.mean(intensity, axis=0))
ax.set_xlabel("X (cm), Image plane", fontsize=14)
ax.set_ylabel("Mean intensity", fontsize=14)
ax.set_title("Mean source profile", fontsize=14)
[20]:
Text(0.5, 1.0, 'Mean source profile')
This page was generated by
nbsphinx from
docs/notebooks/diagnostics/proton_radiography_particle_tracing_wire_mesh.ipynb.
Interactive online version:
.
Synthetic Proton Radiographs with a Wire Mesh¶
In proton radiography experiments a wire mesh grid is often placed between the proton source and the object of interest, leaving a shadow of the grid in the proton fluence. The displacement of these shadow grid lines can then be used to quantitatively extract the lineintegrated force experienced at each grid vertex.
The SyntheticProtronRadiography class includes a method (add_wire_mesh()) that can be used to create synthetic radiographs with a mesh in place. In this example notebook we will illustrate the options available for creating and placing the mesh(s), the demonstrate the use of a mesh grid in a practical example.
[1]:
import astropy.constants as const
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import warnings
from mpl_toolkits.mplot3d import Axes3D
from plasmapy.diagnostics import proton_radiography as prad
from plasmapy.plasma.grids import CartesianGrid
Creating and Placing a Wire Mesh¶
We will begin by creating an empty CartesianGrid object in which the electric and magnetic fields are zero. Particle tracing through this grid will allow us to image just the mesh once we add one in place.
[2]:
empty_grid = CartesianGrid(1 * u.mm, 1 * u.mm, num=50)
SyntheticProtronRadiography will warn us every time we use this grid that the fields are not specified (before assuming that they are zero). The following line will silence this warning.
[3]:
warnings.simplefilter("ignore")
We’ll also define a fixed source and detector that we won’t change for the rest of the example.
[4]:
source = (0 * u.mm, 10 * u.mm, 0 * u.mm)
detector = (0 * u.mm, 200 * u.mm, 0 * u.mm)
Finally, we’ll create an instance of SyntheticProtronRadiography.
[5]:
sim = prad.SyntheticProtonRadiograph(empty_grid, source, detector, verbose=False)
Now it’s time to create the mesh. add_wire_mesh() takes four required parameters:
location
: A vector from the grid origin to the center of the mesh.extent
: The size of the mesh. If two values are given the mesh is assumed to be rectangular (extent is the width, height), but if only one is provided the mesh is assumed to be circular (extent is the diameter).nwires
: The number of wires in each direction. If only one value is given, it’s assumed to be the same for both directions.wire_diameter
: The diameter of each wire.
add_wire_mesh() works by extrapolating the positions of the particles in the mesh plane (based on their initial velocities) and removing those particles that will hit the wires. When add_wire_mesh() is called, the description of the mesh is stored inside the SyntheticProtronRadiography object. Multiple meshes can be added. The particles are then removed when the run() method is called.
[6]:
location = np.array([0, 2, 0]) * u.mm
extent = (1 * u.mm, 1 * u.mm)
nwires = (9, 12)
wire_diameter = 20 * u.um
sim.add_wire_mesh(location, extent, nwires, wire_diameter)
Now that the mesh has been created, we will run the particle tracing simulation and create a synthetic radiograph to visualize the result. We’ll wrap this in a function so we can use it again later.
[7]:
def run_radiograph(sim, vmax=None):
sim.create_particles(1e5, 15 * u.MeV, max_theta=8 * u.deg)
sim.run(field_weighting="nearest neighbor")
h, v, i = sim.synthetic_radiograph(
size=np.array([[1, 1], [1, 1]]) * 2 * u.cm, bins=[200, 200]
)
if vmax is None:
vmax = np.max(i)
fig, ax = plt.subplots(figsize=(8, 8))
ax.pcolormesh(h.to(u.mm).value, v.to(u.mm).value, i.T, cmap="Blues_r", vmax=vmax)
ax.set_aspect("equal")
ax.set_xlabel("X (mm)")
ax.set_ylabel("Y (mm)")
[8]:
run_radiograph(sim)
Notice that the distance from the source to the mesh is \(10  2 = 8\) mm, while the distance from the mesh to the detector is \(200 + 2 = 202\) mm. The magnification is therefore \(M = 1 + 202/8 = 26.25\), so the \(1\) mm wide mesh is \(26.25\) mm wide in the image.
Changing the location
keyword can change both the magnification and shift the mesh center.
[9]:
sim = prad.SyntheticProtonRadiograph(empty_grid, source, detector, verbose=False)
sim.add_wire_mesh(np.array([0.5, 4, 0]) * u.mm, extent, nwires, wire_diameter)
run_radiograph(sim)
Setting the extent
keyword to a single value will create a circular mesh (with a rectangular grid of wires).
[10]:
sim = prad.SyntheticProtonRadiograph(empty_grid, source, detector, verbose=False)
sim.add_wire_mesh(location, (1 * u.mm), nwires, wire_diameter)
run_radiograph(sim)
add_wire_mesh() has two optional keywords that can be used to change the orientation of the mesh. The first, mesh_hdir
is a unit vector that sets the horizontal direction of the mesh plane. This can be used to effectively rotate the mesh. For example the following example will rotate the mesh by \(45^\circ\) (note that
these unit vector inputs are automatically normalized).
[11]:
sim = prad.SyntheticProtonRadiograph(empty_grid, source, detector, verbose=False)
nremoved = sim.add_wire_mesh(
location, extent, nwires, wire_diameter, mesh_hdir=np.array([0.5, 0, 0.5])
)
run_radiograph(sim)
The second keyword argument, mesh_vdir
, overrides the unit vector that defines the vertical direction of the mesh plane. By default this vector is set to be mutually orthogonal to mesh_hdir
and the detector plane normal so that the mesh is parallel to the detector plane. Changing this keyword (alone or in combination with mesh_hdir
) can be used to create a tilted mesh.
[12]:
sim = prad.SyntheticProtonRadiograph(empty_grid, source, detector, verbose=False)
nremoved = sim.add_wire_mesh(
location, extent, nwires, wire_diameter, mesh_vdir=np.array([0, 0.7, 1])
)
run_radiograph(sim)
Using a Wire Mesh in an Example Radiograph¶
To illustrate the use of a mesh in an actual example, we’ll first create an example CartesianGrid object and fill it with the analytical electric field produced by a sphere of Gaussian potential.
[13]:
# Create a Cartesian grid
L = 1 * u.mm
grid = CartesianGrid(L, L, num=150)
# Create a spherical potential with a Gaussian radial distribution
radius = np.linalg.norm(grid.grid, axis=3)
arg = (radius / (L / 3)).to(u.dimensionless_unscaled)
potential = 2e5 * np.exp((arg ** 2)) * u.V
# Calculate E from the potential
Ex, Ey, Ez = np.gradient(potential, grid.dax0, grid.dax1, grid.dax2)
Ex = np.where(radius < L / 2, Ex, 0)
Ey = np.where(radius < L / 2, Ey, 0)
Ez = np.where(radius < L / 2, Ez, 0)
# Add those quantities to the grid
grid.add_quantities(E_x=Ex, E_y=Ey, E_z=Ez, phi=potential)
# Plot the Efield
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection="3d")
ax.view_init(30, 30)
# skip some points to make the vector plot intelligable
s = tuple([slice(None, None, 10)] * 3)
ax.quiver(
grid.pts0[s].to(u.mm).value,
grid.pts1[s].to(u.mm).value,
grid.pts2[s].to(u.mm).value,
grid["E_x"][s],
grid["E_y"][s],
grid["E_z"][s],
length=1e6,
)
ax.set_xlabel("X (mm)")
ax.set_ylabel("Y (mm)")
ax.set_zlabel("Z (mm)")
ax.set_xlim(1, 1)
ax.set_ylim(1, 1)
ax.set_zlim(1, 1)
ax.set_title("Gaussian Potential Electric Field")
[13]:
Text(0.5, 0.92, 'Gaussian Potential Electric Field')
Now we will create a mesh and run the particle tracing simulation
[14]:
sim = prad.SyntheticProtonRadiograph(grid, source, detector, verbose=False)
sim.add_wire_mesh(location, extent, 11, wire_diameter)
sim.create_particles(3e5, 15 * u.MeV, max_theta=8 * u.deg)
run_radiograph(sim, vmax=10)
Notice how the vertices of the grid are displaced by the fields.
Dispersion¶
This page was generated by
nbsphinx from
docs/notebooks/dispersion/dispersion_function.ipynb.
Interactive online version:
.
[1]:
%matplotlib inline
The plasma dispersion function¶
Let’s import some basics (and PlasmaPy
!)
[2]:
import matplotlib.pyplot as plt
import numpy as np
[3]:
import plasmapy.dispersion.dispersionfunction
help(plasmapy.dispersion.dispersionfunction.plasma_dispersion_func)
Help on function plasma_dispersion_func in module plasmapy.dispersion.dispersionfunction:
plasma_dispersion_func(zeta: Union[complex, int, float, numpy.ndarray, astropy.units.quantity.Quantity]) > Union[complex, float, numpy.ndarray, astropy.units.quantity.Quantity]
Calculate the plasma dispersion function.
Parameters

zeta : complex, int, float, ~numpy.ndarray, or ~astropy.units.Quantity
Argument of plasma dispersion function.
Returns

Z : complex, float, or ~numpy.ndarray
Value of plasma dispersion function.
Raises

TypeError
If the argument is of an invalid type.
~astropy.units.UnitsError
If the argument is a `~astropy.units.Quantity` but is not
dimensionless.
ValueError
If the argument is not entirely finite.
See Also

plasma_dispersion_func_deriv
Notes

The plasma dispersion function is defined as:
.. math::
Z(\zeta) = \pi^{0.5} \int_{\infty}^{+\infty}
\frac{e^{x^2}}{x\zeta} dx
where the argument is a complex number [#]_.
In plasma wave theory, the plasma dispersion function appears
frequently when the background medium has a Maxwellian
distribution function. The argument of this function then refers
to the ratio of a wave's phase velocity to a thermal velocity.
References

.. [#] Fried, Burton D. and Samuel D. Conte. 1961.
The Plasma Dispersion Function: The Hilbert Transformation of the
Gaussian. Academic Press (New York and London). ISBN 9781483261737
Examples

>>> plasma_dispersion_func(0)
1.7724538509055159j
>>> plasma_dispersion_func(1j)
0.757872156141312j
>>> plasma_dispersion_func(1.52+0.47j)
(0.6088888957234254+0.33494583882874024j)
Take a look at the docs to plasma_dispersion_func()
for more information on this.
We’ll now make some sample data to visualize the dispersion function:
[4]:
x = np.linspace(1, 1, 1000)
X, Y = np.meshgrid(x, x)
Z = X + 1j * Y
print(Z.shape)
(1000, 1000)
Before we start plotting, let’s make a visualization function first:
[5]:
def plot_complex(X, Y, Z, N=50):
fig, (real_axis, imag_axis) = plt.subplots(1, 2)
real_axis.contourf(X, Y, Z.real, N)
imag_axis.contourf(X, Y, Z.imag, N)
real_axis.set_title("Real values")
imag_axis.set_title("Imaginary values")
for ax in [real_axis, imag_axis]:
ax.set_xlabel("Real values")
ax.set_ylabel("Imaginary values")
fig.tight_layout()
plot_complex(X, Y, Z)
We can now apply our visualization function to our simple dispersion relation
[6]:
# sphinx_gallery_thumbnail_number = 2
F = plasmapy.dispersion.dispersionfunction.plasma_dispersion_func(Z)
plot_complex(X, Y, F)
So this is going to be a hack and I’m not 100% sure the dispersion function is quite what I think it is, but let’s find the area where the dispersion function has a lesser than zero real part because I think it may be important (brb reading Fried and Conte):
[7]:
plot_complex(X, Y, F.real < 0)
We can also visualize the derivative:
[8]:
F = plasmapy.dispersion.dispersionfunction.plasma_dispersion_func_deriv(Z)
plot_complex(X, Y, F)
Plotting the same function on a larger area:
[9]:
x = np.linspace(2, 2, 2000)
X, Y = np.meshgrid(x, x)
Z = X + 1j * Y
print(Z.shape)
(2000, 2000)
[10]:
F = plasmapy.dispersion.dispersionfunction.plasma_dispersion_func(Z)
plot_complex(X, Y, F, 100)
Now we examine the derivative of the dispersion function as a function of the phase velocity of an electromagnetic wave propagating through the plasma. This is recreating figure 5.1 in: J. Sheffield, D. Froula, S. H. Glenzer, and N. C. Luhmann Jr, Plasma scattering of electromagnetic radiation: theory and measurement techniques. Chapter 5 Pg 106 (Academic press, 2010).
[11]:
xs = np.linspace(0, 4, 100)
ws = (1 / 2) * plasmapy.dispersion.dispersionfunction.plasma_dispersion_func_deriv(xs)
wRe = np.real(ws)
wIm = np.imag(ws)
plt.plot(xs, wRe, label="Re")
plt.plot(xs, wIm, label="Im")
plt.axis([0, 4, 0.3, 1])
plt.legend(
loc="upper right", frameon=False, labelspacing=0.001, fontsize=14, borderaxespad=0.1
)
plt.show()
Formulary¶
This page was generated by
nbsphinx from
docs/notebooks/formulary/braginskii.ipynb.
Interactive online version:
.
[1]:
%matplotlib inline
Braginskii coefficients¶
A short example of how to calculate classical transport coefficients from Bragiński’s theory.
[2]:
from astropy import units as u
from plasmapy.formulary import ClassicalTransport
We’ll use some sample ITER data, without much regard for whether the regime is even fit for classical transport theory:
[3]:
thermal_energy_per_electron = 8.8 * u.keV
electron_concentration = 10.1e19 / u.m ** 3
thermal_energy_per_ion = 8.0 * u.keV
ion_concentration = electron_concentration
ion = "D+" # a crude approximation
We now make the default ClassicalTransport
object:
[4]:
braginskii = ClassicalTransport(
thermal_energy_per_electron,
electron_concentration,
thermal_energy_per_ion,
ion_concentration,
ion,
)
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/utils/decorators/checks.py:1400: RelativityWarning: thermal_speed is yielding a velocity that is 18.561% of the speed of light. Relativistic effects may be important.
RelativityWarning,
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/utils/decorators/checks.py:1400: RelativityWarning: V is yielding a velocity that is 18.561% of the speed of light. Relativistic effects may be important.
RelativityWarning,
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/utils/decorators/checks.py:1400: RelativityWarning: thermal_speed is yielding a velocity that is 18.559% of the speed of light. Relativistic effects may be important.
RelativityWarning,
WARNING: UnitsWarning: The argument 'z_mean' to function collision_frequency() has no specified units. Assuming units of . To silence this warning, explicitly pass in an astropy Quantity (e.g. 5. * astropy.units.cm) (see http://docs.astropy.org/en/stable/units/) [plasmapy.utils.decorators.validators]
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/plasmapy0.6.1.dev3+g710abdepy3.7.egg/plasmapy/utils/decorators/checks.py:1400: RelativityWarning: V is yielding a velocity that is 18.559% of the speed of light. Relativistic effects may be important.
RelativityWarning,
These variables are calculated during initialization and can be referred to straight away:
[5]:
print(braginskii.coulomb_log_ei)
print(braginskii.coulomb_log_ii)
print(braginskii.hall_e)
print(braginskii.hall_i)
18.015542112815666
20.41557520752423
0.0
0.0
These quantities are not calculated during initialization and can be referred to via methods. To signify the need to calculate them, we call them via ().
[6]:
print(braginskii.resistivity)
print(braginskii.thermoelectric_conductivity)
print(braginskii.electron_thermal_conductivity)
print(braginskii.ion_thermal_conductivity)
1.1541382845331304e09 m Ohm
0.7110839986207994
1065176076.6192665 W / (K m)
21360464.463766705 W / (K m)
WARNING: UnitsWarning: The argument 'z_mean' to function Coulomb_logarithm() has no specified units. Assuming units of . To silence this warning, explicitly pass in an astropy Quantity (e.g. 5. * astropy.units.cm) (see http://docs.astropy.org/en/stable/units/) [plasmapy.utils.decorators.validators]
They also change with magnetization:
[7]:
mag_braginskii = ClassicalTransport(
thermal_energy_per_electron,
electron_concentration,
thermal_energy_per_ion,
ion_concentration,
ion,
B=0.1 * u.T,
)
print(mag_braginskii.resistivity)
print(mag_braginskii.thermoelectric_conductivity)
print(mag_braginskii.electron_thermal_conductivity)
print(mag_braginskii.ion_thermal_conductivity)
1.1541382845331304e09 m Ohm
0.7110839986207994
1065176076.6192665 W / (K m)
21360464.463766705 W / (K m)
They also change with direction with respect to the magnetic field. Here,
we choose to print out, as arrays, the (parallel, perpendicular,
and cross) directions. Take a look at the docs to
ClassicalTransport
for more information on these.
[8]:
all_direction_braginskii = ClassicalTransport(
thermal_energy_per_electron,
electron_concentration,
thermal_energy_per_ion,
ion_concentration,
ion,
B=0.1 * u.T,
field_orientation="all",
)
print(all_direction_braginskii.resistivity)
print(all_direction_braginskii.thermoelectric_conductivity)
print(all_direction_braginskii.electron_thermal_conductivity)
print(all_direction_braginskii.ion_thermal_conductivity)
[1.15413828e09 2.25078755e09 1.39690568e15] m Ohm
[7.11083999e01 6.76676822e13 5.46328992e07]
[1.06517608e+09 2.08451764e04 3.06777888e+02] W / (K m)
[2.13604645e+07 4.24754851e03 2.69422221e+02] W / (K m)
The viscosities return arrays:
[9]:
print(braginskii.electron_viscosity)
print(mag_braginskii.electron_viscosity)
print(braginskii.ion_viscosity)
print(mag_braginskii.ion_viscosity)
[16.29411376 16.28874805 16.28874805 0. 0. ] Pa s
[1.62941138e+01 2.00480711e25 8.01922844e25 1.47442522e12
2.94885044e12] Pa s
[1271.38945503 1267.52435833 1267.52435833 0. 0. ] Pa s
[1.27138946e+03 5.99222933e17 2.39689173e16 2.57162285e07
5.14324570e07] Pa s
This page was generated by
nbsphinx from
docs/notebooks/formulary/thermal_bremsstrahlung.ipynb.
Interactive online version:
.
Emission of Thermal Bremsstrahlung by a Maxwellian Plasma¶
The radiation.thermal_bremsstrahlung function calculates the bremsstrahlung spectrum emitted by the collision of electrons and ions in a thermal (Maxwellian) plasma. This function calculates this quantity in the RayleighJeans limit where \(\hbar\omega \ll k_B T_e\). In this regime, the power spectrum of the emitted radiation is
where \(w_{pe}\) is the electron plasma frequency and \(E_1\) is the exponential integral
and y is the dimensionless argument
where \(k_{max}\) is a maximum wavenumber arising from binary collisions approximated here as
where \(\lambda_B\) is the electron de Broglie wavelength. In some regimes other values for \(k_{max}\) may be appropriate, so its value may be set using a keyword. Bremsstrahlung emission is greatly reduced below the electron plasma frequency (where the plasma is opaque to EM radiation), so these expressions are only valid in the regime \(w < w_{pe}\).
[1]:
%matplotlib inline
import astropy.constants as const
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from plasmapy.formulary.radiation import thermal_bremsstrahlung
Create an array of frequencies over which to calculate the bremsstrahlung spectrum and convert these frequencies to photon energies for the purpose of plotting the results. Set the plasma density, temperature, and ion species.
[2]:
frequencies = np.arange(15, 16, 0.01)
frequencies = (10 ** frequencies) / u.s
energies = (frequencies * const.h.si).to(u.eV)
ne = 1e22 * u.cm ** 3
Te = 1e2 * u.eV
ion_species = "C12 4+"
Calculate the spectrum, then plot it.
[3]:
spectrum = thermal_bremsstrahlung(frequencies, ne, Te, ion_species=ion_species)
print(spectrum.unit)
lbl = "$T_e$ = {:.1e} eV,\n".format(Te.value) + "$n_e$ = {:.1e} 1/cm^3".format(ne.value)
plt.plot(energies, spectrum, label=lbl)
plt.title(f"Thermal Bremsstrahlung Spectrum")
plt.xlabel("Energy (eV)")
plt.ylabel("Power Spectral Density (W s/m^3)")
plt.legend()
plt.show()
C6 s3 / (F3 J(1/2) kg(3/2) m6)
The power spectrum is the power per angular frequency per volume integrated over \(4\pi\) sr of solid angle, and therefore has units of watts / (rad/s) / m\(^3\) * \(4\pi\) rad = W s/m\(^3\).
[4]:
spectrum = spectrum.to(u.W * u.s / u.m ** 3)
spectrum.unit
[4]:
This means that, for a given volume and time period, the total energy emitted can be determined by integrating the power spectrum
[5]:
t = 5 * u.ns
vol = 0.5 * u.cm ** 3
dw = 2 * np.pi * np.gradient(frequencies) # Frequency step size
total_energy = (np.sum(spectrum * dw) * t * vol).to(u.J)
print("Total Energy: {:.2e} J".format(total_energy.value))
Total Energy: 4.97e+04 J
This page was generated by
nbsphinx from
docs/notebooks/formulary/thermal_speed.ipynb.
Interactive online version:
.
Thermal Speed¶
[1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from astropy import units as u
from plasmapy.formulary import (
Maxwellian_speed_1D,
Maxwellian_speed_2D,
Maxwellian_speed_3D,
)
from plasmapy.formulary.parameters import thermal_speed
The thermal_speed function can be used to calculate the thermal velocity for a Maxwellian velocity distribution. There are three common definitions of the thermal velocity, which can be selected using the “method” keyword, which are defined for a 3D velocity distribution as
‘most_probable’ \(v_{th} = \sqrt{\frac{2 k_B T}{m}}\)
‘rms’ \(v_{th} = \sqrt{\frac{3 k_B T}{m}}\)
‘mean_magnitude’ \(v_{th} = \sqrt{\frac{8 k_B T}{m\pi}}\)
The differences between these velocities can be seen by plotitng them on a 3D Maxwellian speed distribution
[2]:
T = 1e5 * u.K
speeds = np.linspace(0, 8e6, num=600) * u.m / u.s
pdf_3D = Maxwellian_speed_3D(speeds, T=T, particle="e")
fig, ax = plt.subplots(figsize=(4, 3))
v_most_prob = thermal_speed(T=T, particle="e", method="most_probable", ndim=3)
v_rms = thermal_speed(T=T, particle="e", method="rms", ndim=3)
v_mean_magnitude = thermal_speed(T=T, particle="e", method="mean_magnitude", ndim=3)
ax.plot(speeds / v_rms, pdf_3D, color="black", label="Maxwellian")
ax.axvline(x=v_most_prob / v_rms, color="blue", label="Most Probable")
ax.axvline(x=v_rms / v_rms, color="green", label="RMS")
ax.axvline(x=v_mean_magnitude / v_rms, color="red", label="Mean Magnitude")
ax.set_xlim(0.1, 3)
ax.set_ylim(0, None)
ax.set_title("3D")
ax.set_xlabel("v/v$_{rms}$")
ax.set_ylabel("f(v)")
[2]:
Text(0, 0.5, 'f(v)')
Similar speeds are defined for 1D and 2D distributions. The differences between these definitions can be illustrated by plotting them on their respective Maxwellian speed distributions.
[3]:
pdf_1D = Maxwellian_speed_1D(speeds, T=T, particle="e")
pdf_2D = Maxwellian_speed_2D(speeds, T=T, particle="e")
dim = [1, 2, 3]
pdfs = [pdf_1D, pdf_2D, pdf_3D]
plt.tight_layout()
fig, ax = plt.subplots(ncols=3, figsize=(10, 3))
for n, pdf in enumerate(pdfs):
ndim = n + 1
v_most_prob = thermal_speed(T=T, particle="e", method="most_probable", ndim=ndim)
v_rms = thermal_speed(T=T, particle="e", method="rms", ndim=ndim)
v_mean_magnitude = thermal_speed(
T=T, particle="e", method="mean_magnitude", ndim=ndim
)
ax[n].plot(speeds / v_rms, pdf, color="black", label="Maxwellian")
ax[n].axvline(x=v_most_prob / v_rms, color="blue", label="Most Probable")
ax[n].axvline(x=v_rms / v_rms, color="green", label="RMS")
ax[n].axvline(x=v_mean_magnitude / v_rms, color="red", label="Mean Magnitude")
ax[n].set_xlim(0.1, 3)
ax[n].set_ylim(0, None)
ax[n].set_title("{:d}D".format(ndim))
ax[n].set_xlabel("v/v$_{rms}$")
ax[n].set_ylabel("f(v)")
ax[2].legend(bbox_to_anchor=(1.9, 0.8), loc="upper right")
[3]:
<matplotlib.legend.Legend at 0x7f8d5facae90>
<Figure size 432x288 with 0 Axes>
Simulation¶
This page was generated by
nbsphinx from
docs/notebooks/simulation/particle_stepper.ipynb.
Interactive online version:
.
[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]:
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()
We can also take a look at the trajectory in the z direction:
[9]:
trajectory.plot_time_trajectories("z")
or plot the shape of the trajectory in 3D:
[10]:
trajectory.plot_trajectories()
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)
Plasma Objects¶
This page was generated by
nbsphinx from
docs/notebooks/plasma/grids_cartesian.ipynb.
Interactive online version:
.
Grids: UniformlySpaced Cartesian Grids¶
Grids are a datastructure that represent one or more physical quantities that share spatial coordinates. For example, the density or magnetic field in a plasma as specified on a Cartesian grid. In addition to storing data, grids have builtin interpolator functions for estimating the values of quantities in between grid vertices.
Creating a grid¶
[1]:
%matplotlib inline
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from plasmapy.plasma import grids
A grid can be created either by providing three arrays of spatial coordinates for vertices (eg. x,yz positions) or using a np.linspace
like syntax. For example, the two following methods are equivalent:
[2]:
# Method 1
xaxis, yaxis, zaxis = [np.linspace(1 * u.cm, 1 * u.cm, num=20)] * 3
x, y, z = np.meshgrid(xaxis, yaxis, zaxis, indexing="ij")
grid = grids.CartesianGrid(x, y, z)
# Method 2
grid = grids.CartesianGrid(
np.array([1, 1, 1]) * u.cm, np.array([1, 1, 1]) * u.cm, num=(150, 150, 150)
)
The grid object provides access to a number of properties
[3]:
print(f"Is the grid uniformly spaced? {grid.is_uniform}")
print(f"Grid shape: {grid.shape}")
print(f"Grid units: {grid.units}")
print(f"Grid spacing on xaxis: {grid.dax0:.2f}")
Is the grid uniformly spaced? True
Grid shape: (150, 150, 150)
Grid units: [Unit("cm"), Unit("cm"), Unit("cm")]
Grid spacing on xaxis: 0.01 cm
The grid points themselves can be explicitly accessed in one of two forms
[4]:
x, y, z = grid.grids
x.shape
[4]:
(150, 150, 150)
[5]:
xyz = grid.grid
xyz.shape
[5]:
(150, 150, 150, 3)
And the axes can be accessed similarly.
[6]:
xaxis = grid.ax0
xaxis.shape
[6]:
(150,)
Adding Quantities¶
Now that the grid has been initialized, we can add quantities to it that represent physical properties defined on the grid vertices. Each quantity must be a u.Quantity
array of the same shape as the grid.
[7]:
Ex = np.random.rand(*grid.shape) * u.V / u.m
Ey = np.random.rand(*grid.shape) * u.V / u.m
Ez = np.random.rand(*grid.shape) * u.V / u.m
Bz = np.random.rand(*grid.shape) * u.T
Bz.shape
[7]:
(150, 150, 150)
When quantities are added to the grid, they are associated with a key string (just like a dictionary). Any key string can be used, but PlasmaPy functions use a shared set of recognized quantities to automatically interperet quantities. Each entry is stored as a namedtuple with fields (“key”, “description”, “unit”). The full list of recognized quantities can be accessed in the module:
[8]:
for key in grid.recognized_quantities.keys():
rk = grid.recognized_quantities[key]
key, description, unit = rk.key, rk.description, rk.unit
print(f"{key} > {description} ({unit})")
x > x spatial position (m)
y > y spatial position (m)
z > z spatial position (m)
rho > Mass density (kg / m3)
E_x > Electric field (x component) (V / m)
E_y > Electric field (y component) (V / m)
E_z > Electric field (z component) (V / m)
B_x > Magnetic field (x component) (T)
B_y > Magnetic field (y component) (T)
B_z > Magnetic field (z component) (T)
phi > Electric Scalar Potential (V)
Quantities can be added to the grid as keyword arguments. The keyword becomes the key string for the quantity in the dataset.
[9]:
grid.add_quantities(B_z=Bz)
grid.add_quantities(E_x=Ex, E_y=Ey, E_z=Ez)
Adding an unrecognized quantity will lead to a warning, but the quantity will still be added to the grid and can still be accessed by the user later.
[10]:
custom_quantity = np.random.rand(*grid.shape) * u.T * u.mm
grid.add_quantities(int_B=custom_quantity)
/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/v0.6.x/lib/python3.7/sitepackages/ipykernel_launcher.py:2: UserWarning: Warning: int_B is not recognized quantity key
A summary of the grid, including the currentlydefined quantities, can be produced by printing the grid object
[11]:
print(grid)
*** Grid Summary ***
<class 'plasmapy.plasma.grids.CartesianGrid'>
Dimensions: (ax0: 150, ax1: 150, ax2: 150)
Uniformly Spaced: (dax0, dax1, dax2) = (0.013 cm, 0.013 cm, 0.013 cm)

Coordinates:
> ax0 (cm) float64 (150,)
> ax1 (cm) float64 (150,)
> ax2 (cm) float64 (150,)

Recognized Quantities:
> B_z (T) float64 (150, 150, 150)
> E_x (V / m) float64 (150, 150, 150)
> E_y (V / m) float64 (150, 150, 150)
> E_z (V / m) float64 (150, 150, 150)

Unrecognized Quantities:
> int_B (mm T) float64 (150, 150, 150)
A simple list of the defined quantity keys on the grid can also be easily accessed
[12]:
print(grid.quantities)
['B_z', 'E_x', 'E_y', 'E_z', 'int_B']
Methods¶
A number of methods are built into grid objects, and are illustrated here.
The grid.on_grid
method determines which points in an array are within the bounds of the grid. Since our example grid is a cube spanning from 1 to 1 cm on each axis, the first of the following points is on the grid while the second is not.
[13]:
pos = np.array([[0.1, 0.3, 0], [3, 0, 0]]) * u.cm
print(grid.on_grid(pos))
[ True False]
Similarly, the grid.vector_intersects
function determines whether the line between two points passes through the grid.
[14]:
pt0 = np.array([3, 0, 0]) * u.cm
pt1 = np.array([3, 0, 0]) * u.cm
pt2 = np.array([3, 10, 0]) * u.cm
print(f"Line from pt0 to pt1 intersects: {grid.vector_intersects(pt0, pt1)}")
print(f"Line from pt0 to pt2 intersects: {grid.vector_intersects(pt0, pt2)}")
Line from pt0 to pt1 intersects: True
Line from pt0 to pt2 intersects: False
Interpolating Quantities¶
Grid objects contain several interpolator methods to evaluate quantites at positions between the grid vertices. These interpolators use the fact that all of the quantities are defined on the same grid to perform faster interpolations using a nearestneighbor scheme. When an interpolation at a position is requested, the grid indices closest to that position are calculated. Then, the quantity arrays are evaluated at the interpolated indices. Using this method, many quantities can be interpolated at the same positions in almost the same amount of time as is required to interpolate one quantity.
Positions are provided to the interpolator as a u.Quantity
array of shape [N,3] where N is the number of positions and [i,:] represents the x,y,z coordinates of the ith position:
[15]:
pos = np.array([[0.1, 0.3, 0], [0.5, 0.25, 0.8]]) * u.cm
print(f"Pos shape: {pos.shape}")
print(f"Position 1: {pos[0,:]}")
print(f"Position 2: {pos[1,:]}")
Pos shape: (2, 3)
Position 1: [ 0.1 0.3 0. ] cm
Position 2: [0.5 0.25 0.8 ] cm
The simplest interpolator directly returns the nearestneighbor values for each quantity. Positions that are outofbounds return an interpolated value of zero.
[16]:
Ex_vals = grid.nearest_neighbor_interpolator(pos, "E_x")
print(f"Ex at position 1: {Ex_vals[0]:.2f}")
Ex at position 1: 0.03 V / m
Multiple values can be interpolated at the same time by including additional keys as arguments. In this case, the interpolator returns a tuple of arrays as a result.
[17]:
Ex_vals, Ey_vals, Ez_vals = grid.nearest_neighbor_interpolator(pos, "E_x", "E_y", "E_z")
print(f"E at position 1: ({Ex_vals[0]:.2f}, {Ey_vals[0]:.2f}, {Ez_vals[0]:.2f})")
E at position 1: (0.03 V / m, 0.40 V / m, 0.47 V / m)
For a higherorder interpolation, some grids (such as the CartesianGrid subclass) also include a volumeweighted interpolator. This interpolator averages the values on the eight grid vertices surrounding the position (weighted by their distance). The syntax for using this interpolator is the same:
[18]:
Ex_vals, Ey_vals, Ez_vals = grid.volume_averaged_interpolator(pos, "E_x", "E_y", "E_z")
print(f"E at position 1: ({Ex_vals[0]:.2f}, {Ey_vals[0]:.2f}, {Ez_vals[0]:.2f})")
E at position 1: (0.10 V / m, 0.49 V / m, 0.47 V / m)
If repeated identical calls are being made to the same interpolator at different positions (for example, in a simulation loop), setting the persistent
keyword to True will increase performance by not repeatedly reloading the same quantity arrays from the grid. Setting this keyword for a single interpolation will not improve performace, and is not recommended (and is only done here for illustration).
[19]:
Ex_vals, Ey_vals, Ez_vals = grid.volume_averaged_interpolator(
pos, "E_x", "E_y", "E_z", persistent=True
)
print(f"E at position 1: ({Ex_vals[0]:.2f}, {Ey_vals[0]:.2f}, {Ez_vals[0]:.2f})")
E at position 1: (0.10 V / m, 0.49 V / m, 0.47 V / m)
[ ]:
This page was generated by
nbsphinx from
docs/notebooks/plasma/grids_nonuniform.ipynb.
Interactive online version:
.
Grids: NonUniform Grids¶
Some data cannot be easily represented on a grid of uniformly spaced vertices. It is still possible to create a grid object to represent such a dataset.
[1]:
%matplotlib inline
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from plasmapy.plasma import grids
[2]:
grid = grids.NonUniformCartesianGrid(
np.array([1, 1, 1]) * u.cm, np.array([1, 1, 1]) * u.cm, num=(50, 50, 50)
)
Currently, all nonuniform data is stored as an unordered 1D array of points. Therefore, although the dataset created above falls approximately on a Cartesian grid, its treatment is identical to a completely unordered set of points
[3]:
grid.shape
[3]:
(125000,)
Many of the properties defined for uniform grids are inaccessible for nonuniform grids. For example, it is not possible to pull out an axis. However, the following properties still apply
[4]:
print(f"Grid points: {grid.grid.shape}")
print(f"Units: {grid.units}")
Grid points: (125000, 3)
Units: [Unit("cm"), Unit("cm"), Unit("cm")]
Properties can be added in the same way as on uniform grids.
[5]:
Bx = np.random.rand(*grid.shape) * u.T
grid.add_quantities(B_x=Bx)
print(grid)
*** Grid Summary ***
<class 'plasmapy.plasma.grids.NonUniformCartesianGrid'>
Dimensions: (ax: 125000)
NonUniform Spacing

Coordinates:
> ax (cm) object (125000,)

Recognized Quantities:
> B_x (T) float64 (125000,)

Unrecognized Quantities:
None
Methods¶
Many of the methods defined for uniform grids also work for nonuniform grids, however there is usually a substantial performance penalty in the nonuniform case.
For example, grid.on_grid
behaves similarly. In this case, the boundaries of the grid are defined by the furthest point away from the origin in each direction.
[6]:
pos = np.array([[0.1, 0.3, 0], [3, 0, 0]]) * u.cm
print(grid.on_grid(pos))
[ True False]
The same definition is used to define the grid boundaries in grid.vector_intersects
[7]:
pt0 = np.array([3, 0, 0]) * u.cm
pt1 = np.array([3, 0, 0]) * u.cm
pt2 = np.array([3, 10, 0]) * u.cm
print(f"Line from pt0 to pt1 intersects: {grid.vector_intersects(pt0, pt1)}")
print(f"Line from pt0 to pt2 intersects: {grid.vector_intersects(pt0, pt2)}")
Line from pt0 to pt1 intersects: True
Line from pt0 to pt2 intersects: False
Interpolating Quantities¶
Nearestneighbor interpolation also works identically. However, volumeweighted interpolation is not implemented for nonuniform grids.
[8]:
pos = np.array([[0.1, 0.3, 0], [0.5, 0.25, 0.8]]) * u.cm
print(f"Pos shape: {pos.shape}")
print(f"Position 1: {pos[0,:]}")
print(f"Position 2: {pos[1,:]}")
Bx_vals = grid.nearest_neighbor_interpolator(pos, "B_x")
print(f"Bx at position 1: {Bx_vals[0]:.2f}")
Pos shape: (2, 3)
Position 1: [ 0.1 0.3 0. ] cm
Position 2: [0.5 0.25 0.8 ] cm
Bx at position 1: 1.00 T
[ ]:
Feedback and Communication¶
Matrix chat¶
The primary communication channel for PlasmaPy, and the quickest way to ask a question and get a response, is our Matrix room, which is also bridged (mirrored) to Gitter.
GitHub Discussions¶
PlasmaPy’s GitHub Discussions page is a great place to suggest ideas, bring up discussion topics, and ask questions in threaded public discussions.
Mailing list¶
We also have a mailing list that serves as a less volatile discussion forum.
Suggestion box¶
We have a suggestion box if you would like to (optionally anonymously) suggest a feature/topic for consideration. These will be reposted on the mailing list or directly in GitHub issues, as appropriate, for further discussion.
Weekly community meetings¶
We have approximately weekly community meetings in the PlasmaPy room on Jitsi. The schedule of our community meetings is on our calendar, and you may access the minutes and agendas. Any last minute changes will be discussed on Matrix. As of November 2020, our meetings are on Tuesdays at 19:00 UTC.
How to Contribute¶
There are numerous ways to contribute to PlasmaPy, including by providing code and documentation, suggesting and discussing ideas, submitting issues and bug reports, and engaging the broader plasma physics community.
Impostor syndrome disclaimer 1¶
We want your help. No, really.
There may be a little voice inside your head that is telling you that you’re not ready to be an open source contributor; that your skills aren’t nearly good enough to contribute. What could you possibly offer a project like this one?
We assure you  the little voice in your head is wrong. If you can write code at all, you can contribute code to open source. Contributing to open source projects is a fantastic way to advance one’s coding skills. Writing perfect code isn’t the measure of a good developer (that would disqualify all of us!); it’s trying to create something, making mistakes, and learning from those mistakes. That’s how we all improve, and we are happy to help others learn.
Being an open source contributor doesn’t just mean writing code, either. You can help out by writing documentation, tests, or even giving feedback about the project (and yes  that includes giving feedback about the contribution process). Some of these contributions may be the most valuable to the project as a whole, because you’re coming to the project with fresh eyes, so you can see the errors and assumptions that seasoned contributors have glossed over.
Contributing code or documentation to PlasmaPy¶
If you see something you’d like to work on amongst our issues, start hacking away on that! However, please announce your intent first in the relevant issue to make sure there is no work duplication.
Please note that PlasmaPy has a Community Code of Conduct.
Issues marked by the community as help wanted mean just that  either they’re good contributions for outsiders or there’s an issue in the ongoing work that requires a second opinion. Please consider these first!
Work on PlasmaPy is done via GitHub, so you’ll need a (free) account. If you are new to git, helpful resources include documentation on git basics and an interactive git tutorial. You must also install git locally on your computer. We highly recommend getting reasonably familiar with git by going through these tutorials or a Software Carpentry workshop prior to making code contributions. Do note that you can usually find help in the PlasmaPy Matrix chatroom.
For actual guidelines for working on PlasmaPy, please see our Development Guide.
Towncrier changelog entries¶
Every pull request should include a changelog entry. Please see
changelog/README.rst
for instructions.
To summarize, put a file like <PULL REQUEST>.<TYPE>.rst
, where <PULL
REQUEST>
is a pull request number, and <TYPE>
is one of breaking
,
feature
, bugfix
, doc
, removal
, trivial
. If unsure, ask
a maintainer.
Footnotes¶
 1
The imposter syndrome disclaimer was originally written by Adrienne Lowe for a PyCon talk. It was adapted in the README files for MetPy and yt, and was then adapted by PlasmaPy.
Community Code of Conduct¶
The PlasmaPy community strives to follow the best practices in open source software development. New contributors are encouraged to join the team and contribute to the codebase. We anticipate/encourage a global participation from people with diverse backgrounds, skills, interests, and opinions. We believe that such diversity is critical in ensuring a growth of ideas in our community.
Our Pledge¶
In the interest of fostering an open and welcoming environment, we as contributors and maintainers pledge to making participation in our project and our community a harassmentfree experience for everyone, regardless of age, body size, disability, ethnicity, gender identity and expression, level of experience, nationality, personal appearance, race, religion, or sexual identity and orientation.
Our Standards¶
Examples of behavior that contributes to creating a positive environment include:
Using welcoming and inclusive language
Being respectful of differing viewpoints and experiences
Gracefully accepting constructive criticism
Focusing on what is best for the community
Showing empathy towards other community members
We as a community pledge to abide by the following guidelines:
We pledge to treat all people with respect and provide a harassment and bullyingfree environment, regardless of age, sex, sexual orientation and/or gender identity, disability, physical appearance, body size, race, nationality, ethnicity, religion, and level of experience. In particular, sexual language and imagery, sexist, racist, or otherwise exclusionary jokes are not appropriate.
We pledge to respect the work of others by recognizing acknowledgment/citation requests of original authors. As authors, we pledge to be explicit about how we want our own work to be cited or acknowledged.
We pledge to welcome those interested in joining the community, and realize that including people with a variety of opinions and backgrounds will only serve to enrich our community. In particular, discussions relating to pros/cons of various technologies, programming languages, and so on are welcome, but these should be done with respect, taking proactive measure to ensure that all participants are heard and feel confident that they can freely express their opinions.
We pledge to welcome questions and answer them respectfully, paying particular attention to those new to the community. We pledge to provide respectful criticisms and feedback in forums, especially in discussion threads resulting from code contributions.
We pledge to be conscientious of the perceptions of the wider community and to respond to criticism respectfully. We will strive to model behaviors that encourage productive debate and disagreement, both within our community and where we are criticized. We will treat those outside our community with the same respect as people within our community.
We pledge to work from the very beginning of this project to make PlasmaPy accessible to people with disabilities.
We pledge to help the entire community follow these guidelines, and to not remain silent when we see violations of them. We will take action when members of our community violate these guidelines. Members of the PlasmaPy community may contact any member of the Coordinating Committee to report violations. Members of the Coordinating Committee will treat these reports in the strictest confidence. The Coordinating Committee will develop formal procedures for how to handle reported violations.
Our Responsibilities¶
Project maintainers are responsible for clarifying the standards of acceptable behavior and are expected to take appropriate and fair corrective action in response to any instances of unacceptable behavior.
Project maintainers have the right and responsibility to remove, edit, or reject comments, commits, code, wiki edits, issues, and other contributions that are not aligned to this Code of Conduct, or to ban temporarily or permanently any contributor for other behaviors that they deem inappropriate, threatening, offensive, or harmful.
Scope¶
This Code of Conduct applies both within project spaces and in public spaces when an individual is representing the project or its community. Examples of representing a project or community include using an official project email address, posting via an official social media account, or acting as an appointed representative at an online or offline event. Representation of a project may be further defined and clarified by project maintainers.
Enforcement¶
Instances of abusive, harassing, or otherwise unacceptable behavior may be reported by contacting Nick Murphy at namurphy@cfa.harvard.edu or any member of the Coordinating Committee. All complaints will be reviewed and investigated and will result in a response that is deemed necessary and appropriate to the circumstances. The project team is obligated to maintain confidentiality with regard to the reporter of an incident. Project team members should recuse themselves from enforcement of the code of conduct for a given incident if they have an actual or apparent conflict of interest. Further details of specific enforcement policies may be posted separately.
Project maintainers who do not follow or enforce the Code of Conduct in good faith may face temporary or permanent repercussions as determined by other members of the project’s leadership.
Attribution¶
Parts of these guidelines have been adapted from the Contributor Covenant (version 1.4), the Astropy Community Code of Conduct, and the Python Software Foundation code of conduct.
Acknowledging and Citing¶
If you use PlasmaPy for a project that results in a research publication, we ask that you cite the Zenodo record for the specific version of PlasmaPy used in your project. Citing a software package promotes scientific reproducibility, gives credit to its developers, and highlights the importance of software as a vital research product.
Version 0.6.0 of PlasmaPy may be cited with the following reference:
PlasmaPy Community et al. (2021). PlasmaPy, version 0.6.0, Zenodo, http://doi.org/10.5281/zenodo.4602818
This reference may be made, for example, by adding the following line to the methods or acknowledgements section of a paper.
This research made use of PlasmaPy version 0.6.0, a communitydeveloped open source Python package for plasma science (PlasmaPy Community et al. 2021).
A paper written using AASTeX
use the \software{PlasmaPy}
command to cite PlasmaPy.
All public releases of PlasmaPy, as well as many informational resources on PlasmaPy, are openly archived in the PlasmaPy Community on Zenodo.
We encourage authors to acknowledge the packages that PlasmaPy depends on, including but not limited to Astropy, NumPy, and SciPy.
Important
The Analysis and Diagnostic framework is in active development at the moment. For the foreseeable future, the API will be in continuous flux as functionality is added and modified. To follow the package development please visit our GitHub Project ( https://github.com/PlasmaPy/PlasmaPy/projects/19 ) and comment on any of the relevant issues and/or pull requests.
Analysis & Diagnostic Toolkits¶
Analyses and diagnostics go handinhand, but have subtle differences. Thus,
PlasmaPy gives each their own subpackages, plasmapy.analysis
and
plasmapy.diagnostics
respectively.
Think of the plasmapy.analysis
as your toolbox. It has all the tools
(functionality) you need to analyze your data. Functionality is built around
numpy arrays and each
function has a welldefined, focused task. For example, numpy.fft.fft()
does one specific task: compute the onedimensional discrete Fourier Transform.
Similarly, plasmapy.analysis.swept_langmuir.find_floating_potential()
only
finds the floating potential for a given Langmuir trace. It does not have
smoothing. It does not do any filtering. It does not do any signal
conditioning. Its sole task is to find the floating potential of a single
Langmuir trace.
Diagnostics have a much broader scope and leverage the tools
defined in plasmapy.analysis
to give a more integrated user experience when
analyzing data. Diagnostics try to enhance the analysis workflow by focusing
on some of following key areas…
A more humanfriendly way of managing data by building an interface around
xarray
arrays and datasets via custom diagnostic accessors.xarray
provides labelled multidimensional arrays and datasets.Diagnostics selfmanage the computed analysis data within a
xarray
dataset while maintaining the computed data’s relation to the original data.
Quick viewing of analyzed data with default plotting routines.
Fully defining the physical parameters of a diagnostic with purposely designed
Probe
classes that are integrated into the analysis workflow.Adding graphical user interfaces (GUIs) to the analysis workflow via notebook widgets.
Important
The Analysis and Diagnostic framework is in active development at the moment. For the foreseeable future, the API will be in continuous flux as functionality is added and modified. To follow the package development please visit our GitHub Project ( https://github.com/PlasmaPy/PlasmaPy/projects/19 ) and comment on any of the relevant issues and/or pull requests.
plasmapy.analysis
¶
plasmapy.analysis.fit_functions
¶
FitFunction
classes designed to assist in curve fitting of swept Langmuir
traces.
Classes¶

Abstract class for defining fit functions \(f(x)\) and the tools for fitting the function to a set of data. 

A subclass of 

A subclass of 

A subclass of 

A subclass of 
Important
The Analysis and Diagnostic framework is in active development at the moment. For the foreseeable future, the API will be in continuous flux as functionality is added and modified. To follow the package development please visit our GitHub Project ( https://github.com/PlasmaPy/PlasmaPy/projects/19 ) and comment on any of the relevant issues and/or pull requests.
Swept Langmuir Analysis Module¶
Subpackage containing routines for analyzing swept Langmuir probe traces.
Example Notebooks¶
API¶
SubPackages & Modules¶
Functionality for determining the floating potential of a Langmuir sweep. 

Helper functions for analyzing swept Langmuir traces. 
Functions¶

Function for checking that the voltage and current arrays are properly formatted for analysis by 

Determines the floating potential (\(V_f\)) for a given currentvoltage (IV) curve obtained from a swept Langmuir probe. 

Alias to 
Notebooks¶
Important
The Analysis and Diagnostic framework is in active development at the moment. For the foreseeable future, the API will be in continuous flux as functionality is added and modified. To follow the package development please visit our GitHub Project ( https://github.com/PlasmaPy/PlasmaPy/projects/19 ) and comment on any of the relevant issues and/or pull requests.
plasmapy.diagnostics
¶
The diagnostics subpackage contains functionality for defining diagnostic parameters and processing data collected by diagnostics, synthetic or experimental.
Proton radiography analysis¶
Routines for the analysis of proton radiographs. These routines can be broadly classified as either creating synthetic radiographs from prescribed fields or methods of ‘inverting’ experimentally created radiographs to reconstruct the original fields (under some set of assumptions).
Classes¶

Represents a charged particle radiography experiment with simulated or calculated E and B fields given at positions defined by a grid of spatial coordinates. 
Important
The Analysis and Diagnostic framework is in active development at the moment. For the foreseeable future, the API will be in continuous flux as functionality is added and modified. To follow the package development please visit our GitHub Project ( https://github.com/PlasmaPy/PlasmaPy/projects/19 ) and comment on any of the relevant issues and/or pull requests.
Warning
This module will be deprecated in favor of plasmapy.analysis.swept_langmuir
.
Langmuir analysis¶
Defines the Langmuir analysis module as part of the diagnostics package.
Functions¶

Attempt to perform a basic swept probe analysis based on the provided characteristic and probe data. 

Implement the simplest but crudest method for obtaining an estimate of the plasma potential from the probe characteristic. 

Implement the simplest but crudest method for obtaining an estimate of the floating potential from the probe characteristic. 
Obtain an estimate of the electron saturation current corresponding to the obtained plasma potential. 


Implement the simplest but crudest method for obtaining an estimate of the ion saturation current from the probe characteristic. 

Implement the LangmuirMottley (LM) method of obtaining the ion density. 
Implement the LangmuirMottley (LM) method of obtaining the electron density. 


Extract the section of exponential electron current growth from the probe characteristic. 

Extract the section dominated by ion collection from the probe characteristic. 

Obtain the Maxwellian or biMaxwellian electron temperature using the exponential fit method. 

Extrapolate the electron current from the Maxwellian electron temperature obtained in the exponential growth region. 

Reduce a biMaxwellian (dual) temperature to a single mean temperature for a given fraction. 

Implement the Orbital Motion Limit (OML) method of obtaining an estimate of the ion density. 

Extrapolate the ion current from the ion density obtained with the OML method. 

Implement the Druyvesteyn method of obtaining the normalized Electron Energy Distribution Function (EEDF). 
Classes¶

Class representing a single IV probe characteristic for convenient experimental data access and computation. 
Thomson scattering¶
Defines the Thomson scattering analysis module as part of the diagnostics package.
Functions¶

Calculate the spectral density function for Thomson scattering of a probe laser beam by a multispecies Maxwellian plasma. 
Dispersion (plasmapy.dispersion
)¶
The dispersion
subpackage contains functionality associated with
plasma dispersion relations, solvers and analytical solutions.
Example Notebooks¶
SubPackages & Modules¶
Functions¶

Calculate the plasma dispersion function. 
Calculate the derivative of the plasma dispersion function. 


Alias to 

Using the solution provided by Bellan 2012, calculate the analytical solution to the two fluid, lowfrequency (\(\omega/kc \ll 1\)) dispersion relation presented by Stringer 1963. 
Formulary (plasmapy.formulary
)¶
plasmapy.formulary
provides theoretical formulas for calculation of
physical quantities helpful for plasma physics.
Classical transport theory (


Calculate the resistivity. 

Calculate the thermoelectric conductivity. 

Calculate the thermal conductivity for ions. 

Calculate the thermal conductivity for electrons. 

Calculate the ion viscosity. 

Calculate the electron viscosity. 
Classes¶

Classical transport coefficients (e.g. 
Class Inheritance Diagram¶
Collisions (plasmapy.formulary.collisions
)¶
Functions to calculate transport coefficients.
This module includes a number of functions for handling Coulomb collisions spanning weakly coupled (low density) to strongly coupled (high density) regimes.
Coulomb collisions¶
Coulomb collisions are collisions where the interaction force is conveyed via the electric field, instead of any kind of contact force. They usually result in relatively small deflections of particle trajectories. However, given that there are many charged particles in a plasma, one has to take into account the cumulative effects of many such collisions.
Coulomb logarithms¶
Please read the documentation of Coulomb_logarithm
below for an explanation of the
seven PlasmaPysupported methods of computing the Coulomb logarithm.
Collision rates¶
The module gathers a few functions helpful for calculating collision
rates between particles. The most general of these is collision_frequency
,
while if you need average values for a Maxwellian distribution, try
out collision_rate_electron_ion
and collision_rate_ion_ion
. These
use collision_frequency
under the hood.
Macroscopic properties¶
These include:
Functions¶

Compute the Coulomb logarithm. 

Distance of closest approach for a 90° Coulomb collision. 

Impact parameters for classical and quantum Coulomb collision 

Collision frequency of particles in a plasma. 

Cross section for a large angle Coulomb collision. 
Average momentum relaxation rate for a slowly flowing Maxwellian distribution of electrons. 


Average momentum relaxation rate for a slowly flowing Maxwellian distribution of ions. 

Collisional mean free path (m) 

Spitzer resistivity of a plasma 

Return the electrical mobility. 

Knudsen number (dimensionless) 

Ratio of the Coulomb energy to the kinetic (usually thermal) energy. 
Dielectric functions (plasmapy.formulary.dielectric
)¶
Functions to calculate plasma dielectric parameters
Functions¶

Magnetized cold plasma dielectric permittivity tensor elements. 

Magnetized cold plasma dielectric permittivity tensor elements. 

Compute the classical dielectric permittivity for a 1D Maxwellian plasma. 
Dimensionless parameters (plasmapy.formulary.dimensionless
)¶
Module of dimensionless plasma parameters.
These are especially important for determining what regime a plasma is in. (e.g., turbulent, quantum, collisional, etc.).
For example, plasmas at high (much larger than 1) Reynolds numbers are highly turbulent, while turbulence is negligible at low Reynolds numbers.
Functions¶

Compute the ratio of thermal pressure to magnetic pressure. 

Compute the magnetic Reynolds number 

Compare Fermi energy to thermal kinetic energy to check if quantum effects are important. 

Alias to 

Compute the Reynolds number. 

Alias to 
Distribution functions (plasmapy.formulary.distribution
)¶
Common distribution functions for plasmas, such as the Maxwelian or Kappa distributions. Functionality is intended to include generation, fitting and calculation.
Functions¶

Probability distribution function of velocity for a Maxwellian distribution in 1D. 

Probability distribution function of velocity for a Maxwellian distribution in 2D. 

Probability distribution function of velocity for a Maxwellian distribution in 3D. 

Probability distribution function of speed for a Maxwellian distribution in 1D. 

Probability distribution function of speed for a Maxwellian distribution in 2D. 

Probability distribution function of speed for a Maxwellian distribution in 3D. 

Return the probability density at the velocity 

Return the probability density function for finding a particle with velocity components 
Particle drifts (plasmapy.formulary.drifts
)¶
Functions for calculating particle drifts.
Functions¶

Calculate the diamagnetic fluid perpendicular drift. 

Calculate the “electric cross magnetic” particle drift. 

Calculate the general force drift for a particle in a magnetic field. 

Alias to 

Alias to 

Alias to 
Ionization related functionality (plasmapy.formulary.ionization
)¶
Functions related to ionization states and the properties thereof.
Functions¶

Return the average ionization state of ions in a plasma assuming that the numbers of ions in each state are equal. 

Return the ratio of populations of two ionization states. 

Alias for 
Magnetostatics (plasmapy.formulary.magnetostatics
)¶
Define MagneticStatics class to calculate common static magnetic fields as first raised in issue #100.
Classes¶

Circular wire(coil) class. 

Finite length straight wire class. 

General wire class described by its parametric vector equation 

Infinite straight wire class. 

Simple magnetic dipole — two nearby opposite point charges. 
Abstract class for all kinds of magnetic static fields 


Abstract wire class for concrete wires to be inherited from. 
Class Inheritance Diagram¶
Mathematics (plasmapy.formulary.mathematics
)¶
Mathematical formulas relevant to plasma physics.
Functions¶

Calculate the complete FermiDirac integral. 

Calculates the 3D rotation matrix that will rotate vector 
Plasma parameters (plasmapy.formulary.parameters
)¶
Functions to calculate fundamental plasma parameters.
Functions¶

Calculate the Alfvén speed. 

Return the Bohm diffusion coefficient. 

Alias to 

Alias to 

Alias to 

Alias to 

Calculate the characteristic decay length for electric fields, 

Return the number of electrons within a sphere with a radius of the Debye length. 

Calculate the particle gyrofrequency in units of radians per second. 

Return the particle gyroradius. 

Calculate the 

Calculate a charged particle’s inertial length. 

Return the ion sound speed for an electronion plasma. 

Return the most probable speed for a particle within a Kappa distribution. 

Alias to 

Return the lower hybrid frequency. 
Calculate the magnetic energy density. 

Calculate the magnetic pressure. 


Calculate the mass density from a number density. 

Alias to 

Alias to 

Calculate the particle plasma frequency. 

Alias to 

Alias to 

Alias to 

Alias to 

Alias to 

Return the thermal pressure for a Maxwellian distribution. 

Return the most probable speed for a particle within a Maxwellian distribution. 

Alias to 

Return the upper hybrid frequency. 

Alias to 

Alias to 

Alias to 

Alias to 

Alias to 

Alias to 

Alias to 
Quantum physics functions (plasmapy.formulary.quantum
)¶
Functions for quantum parameters, including electron degenerate gases and warm dense matter.
Functions¶

Calculate the ideal chemical potential. 

Return the de Broglie wavelength. 

Alias to 

Calculate the kinetic energy in a degenerate electron gas. 

Alias to 

Alias to 

Calculate the exponential scale length for charge screening for cold and dense plasmas. 
Calculate the thermal de Broglie wavelength for electrons. 

Calculate the WignerSeitz radius, which approximates the interparticle spacing. 
Electromagnetic Radiation Functions (plasmapy.formulary.radiation
)¶
Functions for calculating quantities associated with electromagnetic radiation.
Functions¶

Calculate the bremsstrahlung emission spectrum for a Maxwellian plasma in the RayleighJeans limit \(ℏ ω ≪ k_B T_e\) 
Relativistic functions (plasmapy.formulary.relativity
)¶
Functions for calculating relativistic quantities (\(v \to c\)).
Functions¶
Return the Lorentz factor. 


Calculate the relativistic energy (in joules) of an object of mass 
The subpackage makes heavy use of astropy.units.Quantity
for handling
conversions between different unit systems. This is especially important
for electronvolts, commonly used in plasma physics to denote
temperature, although it is technically a unit of energy.
Most functions expect astropy.units.Quantity
as input, however some
will use the validate_quantities
decorator
to automatically cast arguments to Quantities with appropriate units. If
that happens, you will be notified via an astropy.units.UnitsWarning
.
Please note that well maintained physical constant data with units and
uncertainties can be found in astropy.constants
.
For a general overview of how unitbased input works, take a look at the following example:
Examples¶
Notes for developers¶
Values should be returned as an Astropy Quantity in SI units.
If a quantity has several names, then the function name should be the one that provides the most physical insight into what the quantity represents. For example, ‘gyrofrequency’ indicates gyration, while Larmor frequency indicates that this frequency is somehow related to a human (or perhaps a cat?) named Larmor. Similarly, using omega_ce as a function name for this quantity will make the code less readable to people who are unfamiliar with the notation or use a different symbol.
The docstrings for plasma parameter methods should describe the physics associated with these quantities in ways that are understandable to students who are taking their first course in plasma physics while still being useful to experienced plasma physicists.
Particles (plasmapy.particles
)¶
Introduction¶
The particles
subpackage provides access to information about
atoms, ions, isotopes, and other particles.
Submodules¶
Particle
Class¶
The Particle
class provides an objectoriented
interface to access and represent particle information.
Creating a Particle
object¶
The simplest way to create a Particle
object
is to pass it a str
representing a particle.
>>> from plasmapy.particles import Particle
>>> electron = Particle('e')
The Particle
class accepts a variety of different
str
formats to represent particles. Atomic symbols are casesensitive,
but element names and many aliases are not.
>>> alpha = Particle('alpha')
>>> deuteron = Particle('D+')
>>> triton = Particle('tritium 1+')
>>> iron56 = Particle('Fe56')
>>> helium = Particle('helium')
>>> muon = Particle('mu')
>>> antimuon = Particle('antimuon')
>>> hydride = Particle('H')
An int
may be used as the first positional argument to
Particle
to represent an atomic number. For isotopes
and ions, the mass number may be represented with the mass_numb
keyword and the integer charge may be represented with the Z
keyword.
>>> proton = Particle(1, mass_numb=1, Z=1)
The most frequently used Particle
objects may be
imported directly from the atomic subpackage.
>>> from plasmapy.particles import proton, electron
The Particle
objects that may be imported
directly are: proton
,
electron
, neutron
,
positron
, deuteron
,
triton
, and alpha
.
Accessing particle properties¶
The properties of each particle may be accessed using the attributes of
the corresponding Particle
object.
>>> proton.atomic_number
1
>>> electron.integer_charge
1
>>> triton.mass_number
3
Some of these properties are returned as a Quantity
in
SI units.
>>> alpha.charge
<Quantity 3.20435324e19 C>
>>> deuteron.mass
<Quantity 3.34358372e27 kg>
>>> triton.half_life
<Quantity 3.888e+08 s>
>>> iron56.binding_energy.to('GeV')
<Quantity 0.49225958 GeV>
Strings representing particles may be accessed using the
symbol
,
element
,
isotope
, and
ionic_symbol
attributes.
>>> antimuon.symbol
'mu+'
>>> triton.element
'H'
>>> alpha.isotope
'He4'
>>> deuteron.ionic_symbol
'D 1+'
Categories¶
The categories
attribute returns a set
with the classification categories corresponding to the particle.
>>> sorted(electron.categories)
['charged', 'electron', 'fermion', 'lepton', 'matter', 'stable']
Membership of a particle within a category may be checked using the
is_category
method.
>>> alpha.is_category('lepton')
False
>>> electron.is_category('fermion', 'lepton', 'charged')
True
>>> iron56.is_category(['element', 'isotope'])
True
The particle must be in all of the categories in the require
keyword, at least one of the categories in the any_of
keyword, and
none of the categories in the exclude
in order for it to return
True
.
>>> deuteron.is_category(require={'element', 'isotope', 'ion'})
True
>>> iron56.is_category(any_of=['charged', 'uncharged'])
False
>>> alpha.is_category(exclude='lepton')
True
Calling the is_category
method with no
arguments returns a set containing all of the valid categories for any
particle. Valid categories include: 'actinide'
, 'alkali metal'
,
'alkaline earth metal'
, 'antibaryon'
, 'antilepton'
,
'antimatter'
, 'antineutrino'
, 'baryon'
, 'boson'
,
'charged'
, 'electron'
, 'element'
, 'fermion'
,
'halogen'
, 'ion'
, 'isotope'
, 'lanthanide'
, 'lepton'
,
'matter'
, 'metal'
, 'metalloid'
, 'neutrino'
,
'neutron'
, 'noble gas'
, 'nonmetal'
, 'positron'
,
'posttransition metal'
, 'proton'
, 'stable'
,
'transition metal'
, 'uncharged'
, and 'unstable'
.
Conditionals and equality properties¶
Equality between particles may be tested either between two
Particle
objects, or between a
Particle
object and a str
.
>>> Particle('H1') == Particle('protium 1+')
False
>>> alpha == 'He4 2+'
True
The is_electron
and
is_ion
attributes provide a quick way to
check whether or not a particle is an electron or ion, respectively.
>>> electron.is_electron
True
>>> hydride.is_electron
False
>>> deuteron.is_ion
True
The element
and
isotope
attributes return None
when the
particle does not correspond to an element or isotope. Because
nonempty strings evaluate to True
and None
evaluates to False
when converted to a bool
, these attributes may be used in conditional
statements to test whether or not a particle is in one of these
categories.
particles = [Particle('e'), Particle('Fe56'), Particle('alpha')]
for particle in particles:
if particle.element:
print(f"{particle} corresponds to element {particle.element}")
if particle.isotope:
print(f"{particle} corresponds to isotope {particle.isotope}")
Returning antiparticles¶
The antiparticle of an elementary particle or antiparticle may be found
by either using Python’s unary invert operator (~
) or the
antiparticle
attribute of a
Particle
object.
>>> ~electron
Particle("e+")
>>> antimuon.antiparticle
Particle("mu")
Functions¶
In addition to the Particle
class, the
particles
subpackage has a functional interface.
Symbols and names¶
Several functions in particles
return string representations
of particles, including atomic_symbol
,
isotope_symbol
, ionic_symbol
,
and element_name
.
>>> from plasmapy.particles import *
>>> atomic_symbol('alpha')
'He'
>>> isotope_symbol('alpha')
'He4'
>>> ionic_symbol('alpha')
'He4 2+'
>>> particle_symbol('alpha')
'He4 2+'
>>> element_name('alpha')
'helium'
The full symbol of the particle can be found using
particle_symbol
.
>>> particle_symbol('electron')
'e'
Particle properties¶
The atomic_number
and mass_number
functions are analogous to the corresponding attributes in the
Particle
class.
>>> atomic_number('iron')
26
>>> mass_number('T+')
3
Charge information may be found using integer_charge
and electric_charge
.
>>> integer_charge('H')
1
>>> electric_charge('muon antineutrino')
<Quantity 0. C>
These functions will raise a ChargeError
for
elements and isotopes that lack explicit charge information.
>>> electric_charge('H')
Traceback (most recent call last):
...
plasmapy.particles.exceptions.ChargeError: Charge information is required for electric_charge.
The standard atomic weight for the terrestrial environment may be
accessed using standard_atomic_weight
.
>>> standard_atomic_weight('Pb').to('u')
<Quantity 207.2 u>
The mass of a particle may be accessed through the
particle_mass
function.
>>> particle_mass('deuteron')
<Quantity 3.34358372e27 kg>
Isotopes¶
The relative isotopic abundance of each isotope in the terrestrial
environment may be found using isotopic_abundance
.
>>> isotopic_abundance('H1')
0.999885
>>> isotopic_abundance('D')
0.000115
A list of all discovered isotopes in order of increasing mass number
can be found with known_isotopes
.
>>> known_isotopes('H')
['H1', 'D', 'T', 'H4', 'H5', 'H6', 'H7']
The isotopes of an element with a nonzero isotopic abundance may be
found with common_isotopes
.
>>> common_isotopes('Fe')
['Fe56', 'Fe54', 'Fe57', 'Fe58']
All stable isotopes of an element may be found with
stable_isotopes
.
>>> stable_isotopes('Pb')
['Pb204', 'Pb206', 'Pb207', 'Pb208']
Stability¶
The is_stable
function returns True
for stable
particles and False
for unstable particles.
>>> is_stable('e')
True
>>> is_stable('T')
False
The half_life
function returns the particle’s
halflife as a Quantity
in units of seconds, if known.
>>> half_life('n')
<Quantity 881.5 s>
For stable particles (or particles that have not been discovered to be
unstable), half_life
returns inf
seconds.
>>> half_life('p+')
<Quantity inf s>
If the particle’s halflife is not known to sufficient precision, then
half_life
returns a str
with the estimated value
while issuing a MissingParticleDataWarning
.
Reduced mass¶
The reduced_mass
function is useful in cases of
twobody collisions.
>>> reduced_mass('e', 'p+')
<Quantity 9.10442514e31 kg>
>>> reduced_mass('D+', 'T+')
<Quantity 2.00486597e27 kg>
Nuclear Reactions¶
Binding energy¶
The binding energy of a nuclide may be accessed either as an
attribute of a Particle
object, or by using the
nuclear_binding_energy
function.
>>> from plasmapy.particles import Particle, nuclear_binding_energy
>>> D = Particle('deuterium')
>>> D.binding_energy
<Quantity 3.56414847e13 J>
>>> nuclear_binding_energy('D').to('GeV')
<Quantity 0.00222457 GeV>
Nuclear reaction energy¶
The energy released from a nuclear reaction may be found using the
nuclear_reaction_energy
function. The input may be
a str
representing the reaction.
>>> from plasmapy.particles import nuclear_reaction_energy
>>> nuclear_reaction_energy('Be8 + alpha > carbon12')
<Quantity 1.18025735e12 J>
The reaction may also be inputted using the reactants
and
products
keywords.
>>> nuclear_reaction_energy(reactants=['D', 'T'], products=['alpha', 'n'])
<Quantity 2.81812097e12 J>
Ionization state data structures¶
The ionization state (or charge state) of a plasma refers to the fraction of an element that is at each ionization level. For example, the ionization state of a pure helium plasma could be 5% He⁰⁺, 94% He¹⁺, and 1% He²⁺.
The ionization state of a single element¶
We may use the IonizationState
class
to represent the ionization state of a single element, such as for this
example.
>>> from plasmapy.particles import IonizationState
>>> ionization_state = IonizationState("He", [0.05, 0.94, 0.01])
The ionization state for helium may be accessed using the
ionic_fractions
attribute. These ionic fractions correspond to the
integer_charges
attribute.
>>> ionization_state.ionic_fractions
array([0.05, 0.94, 0.01])
>>> ionization_state.integer_charges
array([0, 1, 2])
The Z_mean
attribute returns the mean integer charge averaged
over all particles in that element.
>>> ionization_state.Z_mean
0.96
The Z_rms
attribute returns the root mean square integer charge.
>>> ionization_state.Z_rms
0.9899...
The Z_most_abundant
attribute returns a list
of the most abundant
ion(s). The list
may contain more than one integer charge in case of
a tie.
>>> ionization_state.Z_most_abundant
[1]
The summarize
method prints out the ionic fraction for the ions with
an abundance of at least 1%.
>>> ionization_state.summarize()
IonizationState instance for He with Z_mean = 0.96

He 0+: 0.050
He 1+: 0.940
He 2+: 0.010

The number density of the element may be specified through the
n_elem
keyword argument.
>>> import astropy.units as u
>>> ionization_state = IonizationState(
... "He", [0.05, 0.94, 0.01], n_elem = 1e19 * u.m ** 3,
... )
The n_e
attribute provides the electron number density as a
Quantity
.
>>> ionization_state.n_e
<Quantity 9.6e+18 1 / m3>
The number_densities
attribute provides the number density of each
ion or neutral.
>>> ionization_state.number_densities
<Quantity [5.0e+17, 9.4e+18, 1.0e+17] 1 / m3>
Ionization states for multiple elements¶
The IonizationStateCollection
class may be used to
represent the ionization state for multiple elements. This can be used,
for example, to describe the various impurities in a fusion plasma or
the charge state distributions of different elements in the solar wind.
>>> from plasmapy.particles import IonizationStateCollection
The minimal input to IonizationStateCollection
is a list
of the elements or isotopes to represent. Integers in the list
will
be treated as atomic numbers.
>>> states = IonizationStateCollection(["H", 2])
To set the ionic fractions for hydrogen, we may do item assignment.
>>> states["H"] = [0.9, 0.1]
We may use indexing to retrieve an IonizationState
instance for an element.
>>> states["H"]
<IonizationState instance for H>
The ionization states for all of the elements may be specified directly as arguments to the class.
>>> states = IonizationStateCollection(
... {"H": [0.01, 0.99], "He": [0.04, 0.95, 0.01]},
... abundances={"H": 1, "He": 0.08},
... n0 = 5e19 * u.m ** 3,
... )
The ionic fractions will be stored as a dict
.
>>> states.ionic_fractions
{'H': array([0.01, 0.99]), 'He': array([0.04, 0.95, 0.01])}
The number density for each element is the product of the number
density scaling factor n0
with that element’s abundance.
The number density for each ion is the product of n0
, the
corresponding element’s abundance, and the ionic fraction.
>>> states.n0
<Quantity 5.e+19 1 / m3>
>>> states.abundances
{'H': 1.0, 'He': 0.08}
>>> states.number_densities["H"]
<Quantity [5.00e+17, 4.95e+19] 1 / m3>
The corresponding element’s abundance, and the ionic fraction.
>>> states.n0
<Quantity 5.e+19 1 / m3>
>>> states.abundances
{'H': 1.0, 'He': 0.08}
>>> states.number_densities["H"]
<Quantity [5.00e+17, 4.95e+19] 1 / m3>
The corresponding element’s abundance, and the ionic fraction.
>>> states.n
<Quantity 5.e+19 1 / m3>
>>> states.abundances
{'H': 1.0, 'He': 0.08}
>>> states.number_densities["H"]
<Quantity [5.00e+17, 4.95e+19] 1 / m3>
The summarize
method may also be
used to get a summary of the ionization states.
>>> states.summarize()

H 1+: 0.990 n_i = 4.95e+19 m**3

He 0+: 0.040 n_i = 1.60e+17 m**3
He 1+: 0.950 n_i = 3.80e+18 m**3

n_e = 5.34e+19 m**3
T_e = 1.30e+04 K

Decorators¶
Passing Particle
objects to functions and methods¶
When calculating plasma parameters, we very frequently need to access
the properties of the particles that make up that plasma. The
particle_input
decorator allows functions
and methods to easily access properties of different particles.
The particle_input
decorator takes valid
representations of particles given in arguments to functions and passes
through the corresponding Particle
object. The
arguments must be annotated with Particle
so that the decorator knows to create the Particle
object. The decorated function can then access particle properties by using
Particle
attributes. This decorator will raise an
InvalidParticleError
if the input does
not correspond to a valid particle.
Here is an example of a decorated function.
from plasmapy.particles import Particle, particle_input
@particle_input
def particle_mass(particle: Particle):
return particle.mass
This function can now accept either Particle
objects
or valid representations of particles.
>>> particle_mass('p+') # string input
<Quantity 1.67262192e27 kg>
>>> proton = Particle("proton")
>>> particle_mass(proton) # Particle object input
<Quantity 1.67262192e27 kg>
If only one positional or keyword argument is annotated with
Particle
, then the keywords mass_numb
and Z
may be used when the decorated function is called.
@particle_input
def integer_charge(particle: Particle, Z: int = None, mass_numb: int = None) > int:
return particle.integer_charge
The above example includes optional type hint annotations for Z
and
mass_numb
and the returned value. The
particle_input
decorator may be used in methods in
classes as well:
class ExampleClass:
@particle_input
def particle_symbol(self, particle: Particle) > str:
return particle.symbol
On occasion it is necessary for a function to accept only certain
categories of particles. The particle_input
decorator enables several ways to allow this.
If an annotated keyword is named element
, isotope
, or ion
;
then particle_input
will raise an
InvalidElementError
,
InvalidIsotopeError
, or
InvalidIonError
if the particle is not
associated with an element, isotope, or ion; respectively.
@particle_input
def capitalized_element_name(element: Particle):
return element.element_name
@particle_input
def number_of_neutrons(isotope: Particle):
return isotope.mass_number  isotope.atomic_number
@particle_input
def number_of_bound_electrons(ion: Particle):
return ion.atomic_number  ion.integer_charge
The keywords require
, any_of
, and exclude
to the
decorator allow further customization of the particle categories
allowed as inputs. These keywords are used as in
is_category
.
@particle_input(require='charged')
def sign_of_charge(charged_particle: Particle):
"""Require a charged particle."""
return '+' if charged_particle.integer_charge > 0 else ''
@particle_input(any_of=['charged', 'uncharged'])
def integer_charge(particle: Particle) > int:
"""Accept only particles with charge information."""
return particle.integer_charge
@particle_input(exclude={'antineutrino', 'neutrino'})
def particle_mass(particle: Particle):
"""
Exclude neutrinos/antineutrinos because these particles have
weakly constrained masses.
"""
return particle.mass
See Also¶
The mendeleev Python package provides access to properties of elements, isotopes, and ions in the periodic table of elements.
Reference/API¶
Functions¶

Return the number of protons in an atom, isotope, or ion. 

Return the atomic symbol. 

Return a list of isotopes of an element with an isotopic abundances greater than zero, or if no input is provided, a list of all such isotopes for every element. 

Return the electric charge (in coulombs) of a particle. 

Return the name of an element. 

Return the halflife in seconds for unstable isotopes and particles, and 

Return the integer charge of a particle. 

Return the ionic symbol of an ion or neutral atom. 

Return 

Return the symbol representing an isotope. 

Return the isotopic abundances if known, and otherwise zero. 

Deserialize a JSON document into the appropriate particle object. 

Deserialize a JSON string into the appropriate particle object. 

Return a list of all known isotopes of an element, or a list of all known isotopes of every element if no input is provided. 

Get the mass number (the number of protons and neutrons) of an isotope. 

Return the nuclear binding energy associated with an isotope. 

Return the released energy from a nuclear reaction. 

Convert arguments to methods and functions to 

Return the mass of a particle. 

Return the symbol of a particle. 

Find the reduced mass between two particles. 

Return a list of all stable isotopes of an element, or if no input is provided, a list of all such isotopes for every element. 

Return the standard (conventional) atomic weight of an element based on the relative abundances of isotopes in terrestrial environments. 
Classes¶
An abstract base class that defines the interface for particles. 

Base class for particles that are defined with physical units. 


A class to represent custom particles. 

A class to represent dimensionless custom particles. 

Representation of the ionic fraction for a single ion. 

Representation of the ionization state distribution of a single element or isotope. 

Describe the ionization state distributions of multiple elements or isotopes. 

A class for an individual particle or antiparticle. 

A custom 

A 
Variables¶
An 

Create an object with taxonomy information for special particles. 

A 

A 

A 

A 

A 

A 

A 
Class Inheritance Diagram¶
Simulation (plasmapy.simulation
)¶
Introduction¶
The simulation
subpackage provides basic, didactic reference
implementations of popular methods of simulating plasmas, and interfaces
to common simulation tools.
Submodules¶
Particle Tracker (plasmapy.simulation.tracker
)¶
Introduction¶
This module contains the ParticleTracker
class, which is a simple particle
stepper implementing the Boris algorithm.
This module is highly unstable and is expected to change a lot in the future.
Reference/API¶
Classes¶

Object representing a species of particles: ions, electrons, or simply a group of particles with a particular initial velocity distribution. 
Class Inheritance Diagram¶
Particle movement integrators¶
Reference/API¶
plasmapy.simulation.particle_integrators Module¶
Particle movement integrators, for particle simulations.
These do not have astropy.units
support, choosing instead to
limit overhead and increase performance.
They act inplace on position and velocity arrays to reduce memory allocation.

The explicit Boris pusher. 
plasmapy.simulation.abstractions
¶
plasmapy.simulation.abstractions Module¶
Abstract classes for numerical simulations.
Classes¶
An abstract base class to represent the normalization constants for systems of equations describing plasmas. 

A prototype abstract interface for numerical simulations. 

A prototype abstract interface for timedependent numerical simulations. 
Class Inheritance Diagram¶
Reference/API¶
Classes¶
A prototype abstract interface for numerical simulations. 

A prototype abstract interface for timedependent numerical simulations. 


Object representing a species of particles: ions, electrons, or simply a group of particles with a particular initial velocity distribution. 
Class Inheritance Diagram¶
PlasmaPy Plasma¶
Overview¶
One of the core classes in PlasmaPy is Plasma
. In order
to make it easy to work with different plasma data in PlasmaPy, the
Plasma
object provides a number of methods for
commonly existing plasmas in nature.
All Plasma objects are created using the Plasma factory
Plasma
.
A number of plasma data structures are supported by subclassing this base object. See Plasma Subclasses to see a list of all of them.
Creating Plasma Objects¶
Plasma objects are constructed using the special factory class
Plasma
:
>>> x = plasmapy.plasma.Plasma(T_e=T_e,
... n_e=n_e,
... Z=Z,
... particle=particle)
The result of a call to Plasma
will be either a
GenericPlasma
object, or a subclass of
GenericPlasma
which deals with a specific type of
data, e.g. PlasmaBlob
or
Plasma3D
(see Plasma Subclasses to see
a list of all of them).

class
plasmapy.plasma.plasma_factory.
PlasmaFactory
(default_widget_type=None, additional_validation_functions=None, registry=None) Plasma factory class. Used to create a variety of Plasma objects. Valid plasma structures are specified by registering them with the factory.
Using Plasma Objects¶
Once a Plasma object has been created using Plasma
it will be a instance or a subclass of the
GenericPlasma
class. The documentation of
GenericPlasma
lists the attributes and methods that
are available on all Plasma objects.
Plasma Classes¶
Defined in plasmapy.plasma.sources
are a set of
GenericPlasma
subclasses which convert the keyword
arguments data to the standard GenericPlasma
interface. These subclasses also provide a method, which describes to
the Plasma
factory
which the data match its plasma data structure.
Classes¶
Registration class for 


A Generic Plasma class. 
Class Inheritance Diagram¶
Plasma Subclasses¶
The Plasma3D
class is a basic structure to contain spatial
information about a plasma. To initialize a Plasma3D system, first
create an instance of the Plasma3D
class and then set the
Plasma3D.density
, Plasma3D.momentum
,
Plasma3D.pressure
and the Plasma3D.magnetic_field
.
This feature is currently under development.
The PlasmaBlob
class is a basic structure to contain just
plasma parameter information about a plasma with no associated
spatial or temporal scales. To initialize a PlasmaBlob system, call
it with arguments: electron temperature PlasmaBlob.T_e
,
and electron density PlasmaBlob.n_e
. You may also optionally
define the ionization, PlasmaBlob.Z
, and relevant plasma
particle, PlasmaBlob.particle
.
This feature is currently under development.
Classes¶



Core class for describing and calculating plasma parameters with spatial dimensions. 

Class for describing and calculating plasma parameters without spatial/temporal description. 
Class Inheritance Diagram¶
Writing a new Plasma subclass¶
Any subclass of GenericPlasma
which defines a method
named is_datasource_for
will
automatically be registered with the
Plasma
factory. The
is_datasource_for
method describes the form of the data for which
the GenericPlasma
subclass is valid. For example,
it might check the number and types of keyword arguments. This makes it
straightforward to define your own GenericPlasma
subclass for a new data structure or a custom data source like simulated
data. These classes only have to be imported for this to work, as
demonstrated by the following example.
import plasmapy.plasma
import astropy.units as u
class FuturePlasma(plasmapy.plasma.GenericPlasma):
def __init__(self, **kwargs):
super(FuturePlasma, self).__init__(**kwargs)
# Specify a classmethod that determines if the input data matches
# this new subclass
@classmethod
def is_datasource_for(cls, **kwargs):
"""
Determines if any of keyword arguments have a dimensionless value.
"""
for _, value in kwargs.items():
try:
if value.unit == u.dimensionless_unscaled:
return True
except AttributeError:
pass
return False
This class will now be available through the
Plasma
factory as long
as this class has been defined, i.e. imported into the current session.
If you do not want to create a method named is_datasource_for
you
can manually register your class and matching method using the following
method.
import plasmapy.plasma
plasmapy.plasma.Plasma.register(FuturePlasma, FuturePlasma.some_matching_method)
API¶
SubPackages & Modules¶
Exceptions and warnings for functionality defined in 

Defines the AbstractGrid class and child classes 

Module for defining the base framework of the plasma classes. 

Module for defining the framework around the plasma factory. 

Classes¶
Registration class for 


A Generic Plasma class. 
Class Inheritance Diagram¶
Core package utilities (plasmapy.utils
)¶
Introduction¶
The utils
subpackage contains functionality that is needed
across multiple subpackages or does not fit nicely in any other subpackage.
Functionality contained in utils
includes:
Warnings and exceptions used in PlasmaPy, such as
RelativityWarning
orPhysicsError
.Decorators we use for reusable physical
Quantity
computation and checking, such asvalidate_quantities
andcheck_relativistic
.Helper utilities for importing and testing packages.
Reference/API¶
plasmapy.utils.decorators Package¶
A module to contain various decorators used to build readable and useful code.
Functions¶
A decorator that adds to a function the ability to convert the function’s return from angular frequency (rad/s) to frequency (Hz). 


Warns or raises an exception when the output of the decorated function is greater than 

A decorator to ‘check’ – limit/control – the values of input and return arguments to a function or method. 

A decorator to ‘check’ – limit/control – the units of input and return arguments to a function or method. 

A decorator which programmatically prepends and/or appends the docstring of the decorated method/function. 
A decorator for decorators, which preserves the signature of the function being wrapped. 


A decorator to ‘validate’ – control and convert – the units and values of input and return arguments to a function or method. 
Classes¶

Base class for ‘Check’ decorator classes. 

A decorator class to ‘check’ – limit/control – the units of input and return arguments to a function or method. 

A decorator class to ‘check’ – limit/control – the values of input and return arguments to a function or method. 

A decorator class to ‘validate’ – control and convert – the units and values of input and return arguments to a function or method. 
Class Inheritance Diagram¶
plasmapy.utils.exceptions Module¶
Exceptions and warnings specific to PlasmaPy.
Classes¶
Base class of PlasmaPy custom errors. 

The base exception for physicsrelated errors. 

An exception to be raised when the input is not a valid Roman numeral. 

An exception to be raised for integers that outside of the range that can be converted to Roman numerals. 

An exception for speeds greater than the speed of light. 

A base exception for errors from 

Base class of PlasmaPy custom warnings. 

A warning for functions that rely on a particular coupling regime to be valid. 

The base warning for warnings related to nonphysical situations. 

A warning for when relativistic velocities are being used in or are returned by nonrelativistic functionality. 
Class Inheritance Diagram¶
plasmapy.utils.code_repr Module¶
Tools for formatting strings, including for error messages.
Functions¶

Approximate a call of a function or class with positional and keyword arguments. 

Approximate the command to instantiate a class, and access an attribute of the resulting class instance. 

Approximate the command to instantiate a class, and then call a method in the resulting class instance. 
plasmapy.utils.pytest_helpers Package¶
Functions¶

Test for ability to handle numpy array quantities. 

Test that a function or class returns the expected result, raises the expected exception, or issues an expected warning for the supplied positional and keyword arguments. 

Test that different functions/inputs return equivalent results. 
Development Guide¶
Installing PlasmaPy for Development¶
Obtaining PlasmaPy source code¶
After creating your GitHub account, go to the PlasmaPy repository on GitHub and fork a copy of PlasmaPy to your account.
To access Git commands on Windows, try Git Bash.
Next you must clone your fork to your computer. Go to the directory that will host your PlasmaPy directory, and run one of the following commands (after changing yourusername to your username). If you would like to use HTTPS (which is the default and easier to set up), then run:
git clone https://github.com/yourusername/PlasmaPy.git
SSH is a more secure option, but requires you to set up an SSH key beforehand. The equivalent SSH command is:
git clone git@github.com:yourusername/PlasmaPy.git
After cloning, we must tell git where the development version of PlasmaPy is by running:
cd PlasmaPy
git remote add upstream git://github.com/PlasmaPy/PlasmaPy.git
To check on which remotes exist, run git remote v
. You should get
something like this:
origin git@github.com:namurphy/PlasmaPy.git (fetch)
origin git@github.com:namurphy/PlasmaPy.git (push)
upstream git@github.com:PlasmaPy/PlasmaPy.git (fetch)
upstream git@github.com:PlasmaPy/PlasmaPy.git (push)
Setting up an environment for development¶
Setup procedures for the two most popular virtual environments, conda and virtualenv, are listed below.
Conda¶
To set up a development environment for PlasmaPy, we strongly recommend the Anaconda distribution.
Activate Anaconda¶
After installing Anaconda, launch any conda environment. By default,
conda installs a root
environment, which you should be able to
activate via
source /home/user/anaconda3/bin/activate root
where /home/user/anaconda3/
can be swapped to wherever your anaconda
installation resides.
On newer versions of Anaconda the recommended activation process has changed to:
. /home/user/anaconda3/etc/profile.d/conda.sh
conda activate
Note
On Windows, the way to do this is via running Anaconda Prompt
from the Start Menu. Git Bash
may also work if you have added
Anaconda to PATH
.
Create your environment¶
Having activated Anaconda, enter PlasmaPy’s repository root directory and create an environment with our suggested packages by executing the following:
conda env create f requirements/environment.yml
You may now enter the environment via
source activate plasmapy
Note
On Windows, skip the source
part of the previous command.
In newer Conda versions, the command to run is
conda activate plasmapy
Virtualenv¶
Create a directory for holding the PlasmaPy repository, move into it and create the virtual environment
virtualenv p python3 .
You may need to make sure that this directory’s path doesn’t contain any spaces, otherwise virtualenv may throw an error.
Your virtual environment should now be created. If you run ls
you
will notice that virtualenv has created a number of subdirectories:
bin/
, lib/
, and include/
. This is why we’re not creating the
virtualenv within the repository itself  so as to not pollute it. To
activate the virtualenv you will run:
source ./bin/activate
You should now see that your shell session is prepended with (plasmapy), like so:
(plasmapy) user@name:~/programming/plasmapy$
This indicates that the virtualenv is running. Congratulations! When your’re done working on PlasmaPy, you can deactivate the virtualenv by running
source deactivate
Now that you have plasmapy on your local computer and you have a virtual environment, you will want to “install” this development version of PlasmaPy along with its dependencies. Start by activating your virtual environment. Next you want install the PlasmaPy dependencies. One way to do this is to do
(plasmapy) user@name:~/programming/plasmapy$ pip install r requirements/environment.txt
Next, setup the development version of PlasmaPy which you just cloned by moving into the root directory of the cloned repo and running the setup.py script there:
(plasmapy) user@name:~/programming/plasmapy/PlasmaPy$ pip install e .
You should now be all set to run development versions of PlasmaPy
modules via import PlasmaPy
in your test scripts!
Running anaconda with virtualenv¶
If you are running the Anaconda suite and want to use virtualenv to
setup your virtual environment, you will have to let the system know
where the Python interpreter can be found. On Linux this is done with
(for example, assuming having installed Anaconda into ~/anaconda3
):
export LD_LIBRARY_PATH="$HOME/anaconda3/lib/"
Exporting the library path to the dynamic linker will only last for the duration of the current shell session.
You will have to add the python library directory to LD_LIBRARY_PATH, as described in a previous step, prior to activating the virtualenv for every new shell session.
Installing your own dev version¶
To be able to import PlasmaPy from your source version, enter the repository root and use one of
python setup.py develop
pip install e .
Note
If you are not working within a virtual environment, this may end in
a permission error  this can be avoided via also adding the
user
flag. But seriously, use a virtual environment and spare
yourself the trouble.
Either one of these commands will create a soft link to your cloned
repository. Any changes in Python code you make there will be there
when you import plasmapy
from an interactive session.
Code Development Guidelines¶
This document describes the coding requirements and guidelines to be followed during the development of PlasmaPy and affiliated packages.
Code written for PlasmaPy must be compatible with Python 3.7 and later.
Coding Style¶
TL;DR: use precommit¶
PlasmaPy has a configuration for the precommit framework that takes care of style mostly automatically.
Install it with pip install precommit
, then use precommit install
within
the repository.
This will cause precommit to download the right versions of linters we use,
then run an automated style checking suite on every commit. Do note that this
works better with a git add
, then git commit
workflow than a git commit
a
workflow — that way, you can check via git diff
what the automated
changes actually did.
Note that the “Style linters / precommit (pull_request)” part of our
Continuous Integration system can and will (metaphorically) shout at you if it
finds you didn’t apply the linters. Also note that the linters’ output may vary
with version, so, rather than apply black
and isort
manually, let
precommit do the version management for you instead!
Our precommit suite can be found in .precommitconfig.yaml. It includes
PlasmaPy Code Style Guide, codified¶
PlasmaPy follows the PEP8 Style Guide for Python Code. This style choice helps ensure that the code will be consistent and readable.
Line lengths should be chosen to maximize the readability and elegance of the code. The maximum line length for Python code in PlasmaPy is 88 characters.
Docstrings and comments should generally be limited to about 72 characters.
During code development, use black to automatically format code and ensure a consistent code style throughout the package and isort to automatically sort imports.
Follow the existing coding style within a subpackage. This includes, for example, variable naming conventions.
Use standard abbreviations for imported packages when possible, such as
import numpy as np
,import matplotlib as mpl
,import matplotlib.pyplot as plt
, and