Refactoring of ParticleTracker, Plasma class

docs/conf.py

@@ -44,7 +44,7 @@
     'python': ('https://docs.python.org/3', None),
     'numpy': ('https://docs.scipy.org/doc/numpy', None),
     'scipy': ('https://docs.scipy.org/doc/scipy/reference/', None),
-    'pandas': ('http://pandas.pydata.org/pandas-docs/stable/', None),
+    'pandas': ('https://pandas.pydata.org/pandas-docs/stable/', None),
     'astropy': ('http://docs.astropy.org/en/stable/', None)}
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']

docs/index.rst

@@ -8,13 +8,13 @@ PlasmaPy Documentation
 
 .. raw:: html
 
-    <div style="width: 30%;" class="sphx-glr-thumbcontainer" tooltip="Let&#x27;s try to look at ITER plasma conditions using the physics subpackage. ">
+    <div style="width: 30%;" class="sphx-glr-thumbcontainer">
 
 .. only:: html
 
-    .. figure:: /auto_examples/images/thumb/sphx_glr_plot_physics_thumb.png
+    .. figure:: /auto_examples/images/thumb/sphx_glr_plot_particle_stepper_thumb.png
 
-        :ref:`sphx_glr_auto_examples_plot_physics.py`
+        :ref:`sphx_glr_auto_examples_plot_particle_stepper.py`
 
 .. raw:: html
 

plasmapy/classes/sources/__init__.py

@@ -1,3 +1,4 @@
 from .plasma3d import Plasma3D
 from .plasmablob import PlasmaBlob
+from .analyticalplasma import AnalyticalPlasma
 from .openpmd_hdf5 import HDF5Reader

plasmapy/classes/sources/analyticalplasma.py

@@ -0,0 +1,43 @@
+"""
+Defines the core Plasma class used by PlasmaPy to represent plasma properties.
+"""
+import warnings
+
+import astropy.units as u
+
+from plasmapy.physics.parameters import _grab_charge
+from plasmapy.physics.dimensionless import (quantum_theta,
+                                            )
+from plasmapy.transport import (coupling_parameter,
+                                )
+from plasmapy.atomic import particle_mass
+
+from plasmapy.utils import CouplingWarning
+from plasmapy.utils.pytest_helpers import call_string
+
+from plasmapy.classes import GenericPlasma
+
+__all__ = [
+    "PlasmaBlob"
+]
+
+
+class AnalyticalPlasma(GenericPlasma):
+    """
+    Class for describing and calculating plasma parameters without
+    spatial/temporal description.
+    """
+
+    def __init__(self, magnetic_field, electric_field):
+        """
+        Initialize plasma paramters.
+        The most basic description is composition (ion), temperature,
+        density, and ionization.
+        """
+        self.interpolate_B = magnetic_field
+        self.interpolate_E = electric_field
+
+    # @classmethod
+    # def is_datasource_for(cls, **kwargs):
+    #     match = 'E' in kwargs.keys() and 'n_e' in kwargs.keys()
+    #     return match

plasmapy/classes/sources/plasma3d.py

@@ -7,6 +7,7 @@
 import itertools
 
 from plasmapy.physics.magnetostatics import MagnetoStatics
+import scipy.interpolate as interp
 
 from astropy.constants import (m_p,
                                m_e,
@@ -87,6 +88,30 @@ def __init__(self, domain_x, domain_y, domain_z):
         self.magnetic_field = np.zeros((3, *self.domain_shape)) * u.T
         self.electric_field = np.zeros((3, *self.domain_shape)) * u.V / u.m
 
+        # create intermediate array of dimension (nx,ny,nz,3) in order to allow
+        # interpolation on non-equal spatial domain dimensions
+        self._B_interpolator = interp.RegularGridInterpolator(
+            (self.x.si.value,
+             self.y.si.value,
+             self.z.si.value),
+            self.magnetic_field.si.value.reshape((*self.domain_shape, 3)),
+            method="linear",
+            bounds_error=True)
+
+        self._E_interpolator = interp.RegularGridInterpolator(
+            (self.x.si.value,
+             self.y.si.value,
+             self.z.si.value),
+            self.electric_field.si.value.reshape((*self.domain_shape, 3)),
+            method="linear",
+            bounds_error=True)
+
+    def interpolate_E(self, x):
+        return self._E_interpolator(x.si.value) * u.V / u.m
+
+    def interpolate_B(self, x):
+        return self._B_interpolator(x.si.value) * u.T
+
     @property
     def velocity(self):
         return self.momentum / self.density
@@ -115,6 +140,7 @@ def is_datasource_for(cls, **kwargs):
             match = False
         return match
 
+    # TODO this needs more numba and less for loop
     def add_magnetostatic(self, *mstats: MagnetoStatics):
         # for each MagnetoStatic argument
         for mstat in mstats:

plasmapy/conftest.py

@@ -4,4 +4,4 @@
 except ImportError:
     pass
 else:
-    matplotlib.use('Agg')
+    matplotlib.use("Agg")

plasmapy/examples/plot_particle_stepper.py

@@ -2,49 +2,50 @@
 Particle stepper
 ================
 
-An example of PlasmaPy's particle stepper class, currently in need of a rewrite
-for speed.
+An example of using PlasmaPy's particle stepper class.
 """
 
 
 import numpy as np
 from astropy import units as u
-from plasmapy.classes import Plasma
+from plasmapy.classes.sources.analyticalplasma import AnalyticalPlasma
 from plasmapy.simulation import ParticleTracker
 from plasmapy.physics.parameters import gyrofrequency
 
 ############################################################
 # Initialize a plasma. This will be a source of electric and magnetic
 # fields for our particles to move in.
-
-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
+# We'll take B in the x direction
 # and E in the y direction, which gets us an E cross B drift
 # in the z direction.
 
-B0 = 4 * u.T
-plasma.magnetic_field[0, :, :, :] = np.ones((10, 10, 10)) * B0
 
-E0 = 2 * u.V / u.m
-plasma.electric_field[1, :, :, :] = np.ones((10, 10, 10)) * E0
+def magnetic_field(r):
+    return u.Quantity([[4, 0, 0]]*len(r), u.T)
+
+
+E_unit = u.V / u.m
+
+
+def electric_field(r):
+    return u.Quantity([[0, 2, 0]]*len(r), E_unit)
+
+
+plasma = AnalyticalPlasma(magnetic_field, electric_field)
 
 ############################################################
 # 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.
 
-freq = gyrofrequency(B0, 'p').to(u.Hz, equivalencies=u.dimensionless_angles())
+freq = gyrofrequency(4 * u.T, 'p').to(u.Hz, equivalencies=u.dimensionless_angles())
 gyroperiod = (1/freq).to(u.s)
 steps_to_gyroperiod = 10
 timestep = gyroperiod / steps_to_gyroperiod
 
 ############################################################
 # Initialize the trajectory calculation.
 
-number_steps = steps_to_gyroperiod * int(2 * np.pi)
+number_steps = 5 * steps_to_gyroperiod * int(2 * np.pi)
 trajectory = ParticleTracker(plasma, 'p', 1, 1, timestep, number_steps)
 
 ############################################################
@@ -58,7 +59,7 @@
 # Run the pusher and plot the trajectory versus time.
 
 trajectory.run()
-trajectory.plot_time_trajectories()
+trajectory.plot_time_trajectories('yz')
 
 ############################################################
 # Plot the shape of the trajectory in 3D.
@@ -71,9 +72,24 @@
 
 vmean = trajectory.velocity_history[:, :, 2].mean()
 print(f"The calculated drift velocity is {vmean:.4f} to compare with the"
-      f"theoretical E0/B0 = {E0/B0:.4f}")
+      f"theoretical E0/B0 = {0.5 * u.m / u.s}")
 
 ############################################################
 # and from position:
 Vdrift = trajectory.position_history[-1, 0, 2] / (trajectory.NT * trajectory.dt)
 print(f"The calculated drift velocity from position is {Vdrift:.4f}")
+
+############################################################
+# Supposing we wanted to examine the effect of the initial velocity in the x-y plane on the trajectory:
+N = 20
+np.random.seed(0)
+trajectory = ParticleTracker(plasma, 'p', N, 1, timestep/10, number_steps*10)
+trajectory._v[:, :2] = np.random.normal(size=(N, 2))
+# we choose this as our example's thumbnail:
+# sphinx_gallery_thumbnail_number = 3
+trajectory.run()
+trajectory.plot_trajectories(alpha=0.8)
+
+############################################################
+# Note how while each trajectory fans out in a different way,
+# each one traverses the z direction in about the same time!