We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Create classes to represent ionization state distributions
My plan for this PR is to create classes to represent the ionization state distributions of one or more elements.  I am going to add in a bunch of dunder methods like `__getitem__` and maybe `__call__` to help making access to the ionization states more straightfoward and intuitive.  Any suggestions on the naming convention will be helpful so that we can maximize readability. 

Eventually we'll need a way to calculate ionization state distributions assuming collisional ionization equilibrium, but that will be for a different PR.  The purpose of this PR is to set up how to store and access the ionization distributions.  This will be discussed in #352.

This will address some of #352.  It will probably be best to wait until after the `0.1.0` release to merge this, since this PR is only for a partial implementation anyway.



--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `plasmapy/classes/sources/plasma3d.py`
Content:
```
1 """
2 Defines the core Plasma class used by PlasmaPy to represent plasma properties.
3 """
4 
5 import numpy as np
6 import astropy.units as u
7 import itertools
8 
9 from plasmapy.formulary.magnetostatics import MagnetoStatics
10 
11 from astropy.constants import mu0
12 
13 from plasmapy.classes import GenericPlasma
14 
15 __all__ = [
16     "Plasma3D"
17 ]
18 
19 
20 class Plasma3D(GenericPlasma):
21     """
22     Core class for describing and calculating plasma parameters with
23     spatial dimensions.
24 
25     Attributes
26     ----------
27     x : `astropy.units.Quantity`
28         x-coordinates within the plasma domain. Equal to the
29         `domain_x` input parameter.
30     y : `astropy.units.Quantity`
31         y-coordinates within the plasma domain. Equal to the
32         `domain_y` input parameter.
33     z : `astropy.units.Quantity`
34         z-coordinates within the plasma domain. Equal to the
35         `domain_z` input parameter.
36     grid : `astropy.units.Quantity`
37         (3, x, y, z) array containing the values of each coordinate at
38         every point in the domain.
39     domain_shape : tuple
40         Shape of the plasma domain.
41     density : `astropy.units.Quantity`
42         (x, y, z) array of mass density at every point in the domain.
43     momentum : `astropy.units.Quantity`
44         (3, x, y, z) array of the momentum vector at every point in
45         the domain.
46     pressure : `astropy.units.Quantity`
47         (x, y, z) array of pressure at every point in the domain.
48     magnetic_field : `astropy.units.Quantity`
49         (3, x, y, z) array of the magnetic field vector at every point
50         in the domain.
51 
52     Parameters
53     ----------
54     domain_x : `astropy.units.Quantity`
55         1D array of x-coordinates for the plasma domain. Must have
56         units convertable to length.
57     domain_y : `astropy.units.Quantity`
58         1D array of y-coordinates for the plasma domain. Must have
59         units convertable to length.
60     domain_z : `astropy.units.Quantity`
61         1D array of z-coordinates for the plasma domain. Must have
62         units convertable to length.
63 
64     """
65     @u.quantity_input(domain_x=u.m, domain_y=u.m, domain_z=u.m)
66     def __init__(self, domain_x, domain_y, domain_z):
67         # Define domain sizes
68         self.x = domain_x
69         self.y = domain_y
70         self.z = domain_z
71 
72         self.grid = np.array(np.meshgrid(self.x, self.y, self.z,
73                                          indexing='ij'))
74         self.domain_shape = (len(self.x), len(self.y), len(self.z))
75 
76         # Initiate core plasma variables
77         self.density = np.zeros(self.domain_shape) * u.kg / u.m**3
78         self.momentum = np.zeros((3, *self.domain_shape)) * u.kg / (u.m ** 2 * u.s)
79         self.pressure = np.zeros(self.domain_shape) * u.Pa
80         self.magnetic_field = np.zeros((3, *self.domain_shape)) * u.T
81         self.electric_field = np.zeros((3, *self.domain_shape)) * u.V / u.m
82 
83     @property
84     def velocity(self):
85         return self.momentum / self.density
86 
87     @property
88     def magnetic_field_strength(self):
89         B = self.magnetic_field
90         return np.sqrt(np.sum(B * B, axis=0))
91 
92     @property
93     def electric_field_strength(self):
94         E = self.electric_field
95         return np.sqrt(np.sum(E * E, axis=0))
96 
97     @property
98     def alfven_speed(self):
99         B = self.magnetic_field
100         rho = self.density
101         return np.sqrt(np.sum(B * B, axis=0) / (mu0 * rho))
102 
103     @classmethod
104     def is_datasource_for(cls, **kwargs):
105         if len(kwargs) == 3:
106             match = all(f'domain_{direction}' in kwargs.keys() for direction in 'xyz')
107         else:
108             match = False
109         return match
110 
111     def add_magnetostatic(self, *mstats: MagnetoStatics):
112         # for each MagnetoStatic argument
113         for mstat in mstats:
114             # loop over 3D-index (ix,iy,iz)
115             for point_index in itertools.product(*[list(range(n)) for n in self.domain_shape]):
116                 # get coordinate
117                 p = self.grid[(slice(None),)+point_index]  # function as [:, *index]
118                 # calculate magnetic field at this point and add back
119                 self.magnetic_field[(slice(None),)+point_index] += mstat.magnetic_field(p)
120
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/plasmapy/classes/sources/plasma3d.py b/plasmapy/classes/sources/plasma3d.py
--- a/plasmapy/classes/sources/plasma3d.py
+++ b/plasmapy/classes/sources/plasma3d.py
@@ -110,9 +110,10 @@
 
     def add_magnetostatic(self, *mstats: MagnetoStatics):
         # for each MagnetoStatic argument
+        prod = itertools.product(*[list(range(n)) for n in self.domain_shape])
         for mstat in mstats:
             # loop over 3D-index (ix,iy,iz)
-            for point_index in itertools.product(*[list(range(n)) for n in self.domain_shape]):
+            for point_index in prod:
                 # get coordinate
                 p = self.grid[(slice(None),)+point_index]  # function as [:, *index]
                 # calculate magnetic field at this point and add back

{"golden_diff": "diff --git a/plasmapy/classes/sources/plasma3d.py b/plasmapy/classes/sources/plasma3d.py\n--- a/plasmapy/classes/sources/plasma3d.py\n+++ b/plasmapy/classes/sources/plasma3d.py\n@@ -110,9 +110,10 @@\n \n     def add_magnetostatic(self, *mstats: MagnetoStatics):\n         # for each MagnetoStatic argument\n+        prod = itertools.product(*[list(range(n)) for n in self.domain_shape])\n         for mstat in mstats:\n             # loop over 3D-index (ix,iy,iz)\n-            for point_index in itertools.product(*[list(range(n)) for n in self.domain_shape]):\n+            for point_index in prod:\n                 # get coordinate\n                 p = self.grid[(slice(None),)+point_index]  # function as [:, *index]\n                 # calculate magnetic field at this point and add back\n", "issue": "Create classes to represent ionization state distributions\nMy plan for this PR is to create classes to represent the ionization state distributions of one or more elements.  I am going to add in a bunch of dunder methods like `__getitem__` and maybe `__call__` to help making access to the ionization states more straightfoward and intuitive.  Any suggestions on the naming convention will be helpful so that we can maximize readability. \r\n\r\nEventually we'll need a way to calculate ionization state distributions assuming collisional ionization equilibrium, but that will be for a different PR.  The purpose of this PR is to set up how to store and access the ionization distributions.  This will be discussed in #352.\r\n\r\nThis will address some of #352.  It will probably be best to wait until after the `0.1.0` release to merge this, since this PR is only for a partial implementation anyway.\r\n\r\n\n", "before_files": [{"content": "\"\"\"\nDefines the core Plasma class used by PlasmaPy to represent plasma properties.\n\"\"\"\n\nimport numpy as np\nimport astropy.units as u\nimport itertools\n\nfrom plasmapy.formulary.magnetostatics import MagnetoStatics\n\nfrom astropy.constants import mu0\n\nfrom plasmapy.classes import GenericPlasma\n\n__all__ = [\n    \"Plasma3D\"\n]\n\n\nclass Plasma3D(GenericPlasma):\n    \"\"\"\n    Core class for describing and calculating plasma parameters with\n    spatial dimensions.\n\n    Attributes\n    ----------\n    x : `astropy.units.Quantity`\n        x-coordinates within the plasma domain. Equal to the\n        `domain_x` input parameter.\n    y : `astropy.units.Quantity`\n        y-coordinates within the plasma domain. Equal to the\n        `domain_y` input parameter.\n    z : `astropy.units.Quantity`\n        z-coordinates within the plasma domain. Equal to the\n        `domain_z` input parameter.\n    grid : `astropy.units.Quantity`\n        (3, x, y, z) array containing the values of each coordinate at\n        every point in the domain.\n    domain_shape : tuple\n        Shape of the plasma domain.\n    density : `astropy.units.Quantity`\n        (x, y, z) array of mass density at every point in the domain.\n    momentum : `astropy.units.Quantity`\n        (3, x, y, z) array of the momentum vector at every point in\n        the domain.\n    pressure : `astropy.units.Quantity`\n        (x, y, z) array of pressure at every point in the domain.\n    magnetic_field : `astropy.units.Quantity`\n        (3, x, y, z) array of the magnetic field vector at every point\n        in the domain.\n\n    Parameters\n    ----------\n    domain_x : `astropy.units.Quantity`\n        1D array of x-coordinates for the plasma domain. Must have\n        units convertable to length.\n    domain_y : `astropy.units.Quantity`\n        1D array of y-coordinates for the plasma domain. Must have\n        units convertable to length.\n    domain_z : `astropy.units.Quantity`\n        1D array of z-coordinates for the plasma domain. Must have\n        units convertable to length.\n\n    \"\"\"\n    @u.quantity_input(domain_x=u.m, domain_y=u.m, domain_z=u.m)\n    def __init__(self, domain_x, domain_y, domain_z):\n        # Define domain sizes\n        self.x = domain_x\n        self.y = domain_y\n        self.z = domain_z\n\n        self.grid = np.array(np.meshgrid(self.x, self.y, self.z,\n                                         indexing='ij'))\n        self.domain_shape = (len(self.x), len(self.y), len(self.z))\n\n        # Initiate core plasma variables\n        self.density = np.zeros(self.domain_shape) * u.kg / u.m**3\n        self.momentum = np.zeros((3, *self.domain_shape)) * u.kg / (u.m ** 2 * u.s)\n        self.pressure = np.zeros(self.domain_shape) * u.Pa\n        self.magnetic_field = np.zeros((3, *self.domain_shape)) * u.T\n        self.electric_field = np.zeros((3, *self.domain_shape)) * u.V / u.m\n\n    @property\n    def velocity(self):\n        return self.momentum / self.density\n\n    @property\n    def magnetic_field_strength(self):\n        B = self.magnetic_field\n        return np.sqrt(np.sum(B * B, axis=0))\n\n    @property\n    def electric_field_strength(self):\n        E = self.electric_field\n        return np.sqrt(np.sum(E * E, axis=0))\n\n    @property\n    def alfven_speed(self):\n        B = self.magnetic_field\n        rho = self.density\n        return np.sqrt(np.sum(B * B, axis=0) / (mu0 * rho))\n\n    @classmethod\n    def is_datasource_for(cls, **kwargs):\n        if len(kwargs) == 3:\n            match = all(f'domain_{direction}' in kwargs.keys() for direction in 'xyz')\n        else:\n            match = False\n        return match\n\n    def add_magnetostatic(self, *mstats: MagnetoStatics):\n        # for each MagnetoStatic argument\n        for mstat in mstats:\n            # loop over 3D-index (ix,iy,iz)\n            for point_index in itertools.product(*[list(range(n)) for n in self.domain_shape]):\n                # get coordinate\n                p = self.grid[(slice(None),)+point_index]  # function as [:, *index]\n                # calculate magnetic field at this point and add back\n                self.magnetic_field[(slice(None),)+point_index] += mstat.magnetic_field(p)\n", "path": "plasmapy/classes/sources/plasma3d.py"}], "after_files": [{"content": "\"\"\"\nDefines the core Plasma class used by PlasmaPy to represent plasma properties.\n\"\"\"\n\nimport numpy as np\nimport astropy.units as u\nimport itertools\n\nfrom plasmapy.formulary.magnetostatics import MagnetoStatics\n\nfrom astropy.constants import mu0\n\nfrom plasmapy.classes import GenericPlasma\n\n__all__ = [\n    \"Plasma3D\"\n]\n\n\nclass Plasma3D(GenericPlasma):\n    \"\"\"\n    Core class for describing and calculating plasma parameters with\n    spatial dimensions.\n\n    Attributes\n    ----------\n    x : `astropy.units.Quantity`\n        x-coordinates within the plasma domain. Equal to the\n        `domain_x` input parameter.\n    y : `astropy.units.Quantity`\n        y-coordinates within the plasma domain. Equal to the\n        `domain_y` input parameter.\n    z : `astropy.units.Quantity`\n        z-coordinates within the plasma domain. Equal to the\n        `domain_z` input parameter.\n    grid : `astropy.units.Quantity`\n        (3, x, y, z) array containing the values of each coordinate at\n        every point in the domain.\n    domain_shape : tuple\n        Shape of the plasma domain.\n    density : `astropy.units.Quantity`\n        (x, y, z) array of mass density at every point in the domain.\n    momentum : `astropy.units.Quantity`\n        (3, x, y, z) array of the momentum vector at every point in\n        the domain.\n    pressure : `astropy.units.Quantity`\n        (x, y, z) array of pressure at every point in the domain.\n    magnetic_field : `astropy.units.Quantity`\n        (3, x, y, z) array of the magnetic field vector at every point\n        in the domain.\n\n    Parameters\n    ----------\n    domain_x : `astropy.units.Quantity`\n        1D array of x-coordinates for the plasma domain. Must have\n        units convertable to length.\n    domain_y : `astropy.units.Quantity`\n        1D array of y-coordinates for the plasma domain. Must have\n        units convertable to length.\n    domain_z : `astropy.units.Quantity`\n        1D array of z-coordinates for the plasma domain. Must have\n        units convertable to length.\n\n    \"\"\"\n    @u.quantity_input(domain_x=u.m, domain_y=u.m, domain_z=u.m)\n    def __init__(self, domain_x, domain_y, domain_z):\n        # Define domain sizes\n        self.x = domain_x\n        self.y = domain_y\n        self.z = domain_z\n\n        self.grid = np.array(np.meshgrid(self.x, self.y, self.z,\n                                         indexing='ij'))\n        self.domain_shape = (len(self.x), len(self.y), len(self.z))\n\n        # Initiate core plasma variables\n        self.density = np.zeros(self.domain_shape) * u.kg / u.m**3\n        self.momentum = np.zeros((3, *self.domain_shape)) * u.kg / (u.m ** 2 * u.s)\n        self.pressure = np.zeros(self.domain_shape) * u.Pa\n        self.magnetic_field = np.zeros((3, *self.domain_shape)) * u.T\n        self.electric_field = np.zeros((3, *self.domain_shape)) * u.V / u.m\n\n    @property\n    def velocity(self):\n        return self.momentum / self.density\n\n    @property\n    def magnetic_field_strength(self):\n        B = self.magnetic_field\n        return np.sqrt(np.sum(B * B, axis=0))\n\n    @property\n    def electric_field_strength(self):\n        E = self.electric_field\n        return np.sqrt(np.sum(E * E, axis=0))\n\n    @property\n    def alfven_speed(self):\n        B = self.magnetic_field\n        rho = self.density\n        return np.sqrt(np.sum(B * B, axis=0) / (mu0 * rho))\n\n    @classmethod\n    def is_datasource_for(cls, **kwargs):\n        if len(kwargs) == 3:\n            match = all(f'domain_{direction}' in kwargs.keys() for direction in 'xyz')\n        else:\n            match = False\n        return match\n\n    def add_magnetostatic(self, *mstats: MagnetoStatics):\n        # for each MagnetoStatic argument\n        prod = itertools.product(*[list(range(n)) for n in self.domain_shape])\n        for mstat in mstats:\n            # loop over 3D-index (ix,iy,iz)\n            for point_index in prod:\n                # get coordinate\n                p = self.grid[(slice(None),)+point_index]  # function as [:, *index]\n                # calculate magnetic field at this point and add back\n                self.magnetic_field[(slice(None),)+point_index] += mstat.magnetic_field(p)\n", "path": "plasmapy/classes/sources/plasma3d.py"}]}