27 Nov 2021 - tsp

Last update 27 Nov 2021

8 mins

This is a simple summary on how one can use `numpy`

and `matplotlib`

to plot simple static electric fields caused by point like charge carriers (not
taking into account for non point like carriers, magnetic fields, currents, time
dependent stuff, etc. - so this will not be a Maxwell equation
solver).

The following python script uses Python 3.6 or higher as well as `matplotlib`

.

For version 3.6 these packages can be simply installed using:

```
pkg install python36
pkg install py36-pip
```

To install `matplotlib`

then in theory it would be sufficient to simply
execute the following command line (don’t execute - read further):

```
python3.6 -m pip install matplotlib
```

Currently on FreeBSD it might be necessary to specify the C and C++ compilers to
be used since `pip`

is a little bit buggy …

```
env CC=clang CXX=clang++ python3.6 -m pip install matplotlib
```

Note that even at time of writing the current and thus recommended Python version would be 3.8

The theory is pretty simple. During classical definition of the electric field we simply assume that the field $\vec{E}$ specifies the force $\vec{F}$ that’s caused on a test charge $q$, the vector $\vec{r}$ just specifying the spatial position.

[ \vec{F}(\vec{r}) = \vec{E}(\vec{r}) * q \\ \to \vec{E}(\vec{r}) = \frac{\vec{F}(\vec{r})}{q} ]Since we know from Coulombs law that the force between two charges $q_1$ and $q_2$ is given by

[ F = \frac{1}{4 \pi \epsilon_0} \frac{q_1 * q_2}{r^2} ]we can simply derive the equation for the field caused by a single point like charge $q$ to be

[ \vec{E}(\vec{r}) = \frac{1}{4 \pi \epsilon_0} \frac{q}{r^2} * \vec{e_r} \\ \vec{E}(\vec{r}) = \frac{1}{4 \pi \epsilon_0} \frac{q}{r^2} * \frac{\vec{r}}{\mid r \mid} \\ \vec{E}(\vec{r}) = \frac{1}{4 \pi \epsilon_0} \frac{q * \vec{r}}{r^3} ]Since we know that the superposition principle applies the fields caused by a number of charge carriers simply sum up. Assuming the charge $q_n$ is located at position $r_{q,n}$ one thus can simple calculate:

[ \vec{E}(\vec{r}) = \frac{1}{4 \pi \epsilon_0} \sum_{n = 0}^{N} \frac{q_n * (\vec{r} - \vec{r_n})}{\mid \vec{r} - \vec{r_n} \mid^3} ]This is then exactly the equation that we’re going to apply for our simple plotting function.

Implementing such a plotting function is really simple. First one has to import
the required packages. We’re going to use `numpy`

for some simple calculations
as well as `matplotlib`

for the actual plotting:

```
import numpy as np
import matplotlib.pyplot as plt
```

To make life more easy we’re going to implement a single function that calculates the field strength caused by a single charge that’s resting at a position $\vec{r_0}$ and is evaluated at the position $(x,y)$:

```
def EFieldSingleCharge(q, r0, x, y):
distance = np.hypot(x - r0[0], y - r0[1])
return q * (x - r0[0]) / (distance**3), q * (y - r0[1]) / (distance**3)
```

Here one can see that we used `np.hypot`

to calculate the distance between
the charge carrier and the position we’re interested at. Then we applied the
formula described above for $\vec{E}(\vec{r}) in both dimensions.

The next thing we do require is a simple array of coordinates for each point of
our grid. First we have to decide what range our coordinates should have and
which subdivision we’re going to use - in the following sample I’m using 128
subdivisions in both axis and then used the `linspace`

function to generate
a simple list in the range from $-2$ to $2$ regularly spaced into $nx$ and $ny$
intervals.

```
nx = 128
ny = 128
x = np.linspace(-2, 2, nx)
y = np.linspace(-2, 2, ny)
```

This yields a simple list like

```
[-2. -1.96850394 -1.93700787 -1.90551181 -1.87401575 -1.84251969
-1.81102362 -1.77952756 -1.7480315 -1.71653543 -1.68503937 -1.65354331
-1.62204724 -1.59055118 -1.55905512 -1.52755906 -1.49606299 -1.46456693
-1.43307087 -1.4015748 -1.37007874 -1.33858268 -1.30708661 -1.27559055
-1.24409449 -1.21259843 -1.18110236 -1.1496063 -1.11811024 -1.08661417
-1.05511811 -1.02362205 -0.99212598 -0.96062992 -0.92913386 -0.8976378
-0.86614173 -0.83464567 -0.80314961 -0.77165354 -0.74015748 -0.70866142
-0.67716535 -0.64566929 -0.61417323 -0.58267717 -0.5511811 -0.51968504
-0.48818898 -0.45669291 -0.42519685 -0.39370079 -0.36220472 -0.33070866
-0.2992126 -0.26771654 -0.23622047 -0.20472441 -0.17322835 -0.14173228
-0.11023622 -0.07874016 -0.04724409 -0.01574803 0.01574803 0.04724409
0.07874016 0.11023622 0.14173228 0.17322835 0.20472441 0.23622047
0.26771654 0.2992126 0.33070866 0.36220472 0.39370079 0.42519685
0.45669291 0.48818898 0.51968504 0.5511811 0.58267717 0.61417323
0.64566929 0.67716535 0.70866142 0.74015748 0.77165354 0.80314961
0.83464567 0.86614173 0.8976378 0.92913386 0.96062992 0.99212598
1.02362205 1.05511811 1.08661417 1.11811024 1.1496063 1.18110236
1.21259843 1.24409449 1.27559055 1.30708661 1.33858268 1.37007874
1.4015748 1.43307087 1.46456693 1.49606299 1.52755906 1.55905512
1.59055118 1.62204724 1.65354331 1.68503937 1.71653543 1.7480315
1.77952756 1.81102362 1.84251969 1.87401575 1.90551181 1.93700787
1.96850394 2. ]
```

This of course would already be sufficient to access the coordinates at each
point if one assumes the mesh is anyways linearly spaces. But since one should always
implement the most general solution - especially when it’s easy - I decided to
use `meshgrid`

to generate two lists of lists that contain the coordinate
components of each point on our grid. Thus the X direction grid contains `ny`

lists with each `nx`

entries that look like the following:

```
[[-2. -1.96850394 -1.93700787 ... 1.93700787 1.96850394
2. ]
[-2. -1.96850394 -1.93700787 ... 1.93700787 1.96850394
2. ]
[-2. -1.96850394 -1.93700787 ... 1.93700787 1.96850394
2. ]
...
[-2. -1.96850394 -1.93700787 ... 1.93700787 1.96850394
2. ]
[-2. -1.96850394 -1.93700787 ... 1.93700787 1.96850394
2. ]
[-2. -1.96850394 -1.93700787 ... 1.93700787 1.96850394
2. ]]
```

This is done by using

```
X, Y = np.meshgrid(x,y)
```

The last missing part are the charges. I’m simply going to create a list of point charges:

```
charges = []
charges.append((1, (-1, 0)))
charges.append((-1, (1, 0)))
charges.append((1, (1, -1)))
charges.append((-1, (-1, -1)))
charges.append((1, (1, 1)))
charges.append((-1, (-1, 1)))
```

Now one can already calculate the fields. We know the fields will have components in $x$ and $y$ direction on all $(nx,ny)$ grid points. Thus one can first initialize zero arrays:

```
Ex = np.zeros((ny, nx))
Ey = np.zeros((ny, nx))
```

The next step is simply calculating the effect of each charge on each grid point and adding up the resulting fields (as implied by the superposition principle)

```
for charge in charges:
ex, ey = EFieldSingleCharge(*charge, x=X, y=Y)
Ex += ex
Ey += ey
```

One can see that I used the `*charge`

argument here - sometimes this
confuses people new to Python. What happens here is that one can recall
the definition of the `EFieldSingleCharge(q, r0, x, y)`

function. As one
can see this requires four parameters - the scalar charge $q$, a vector $r0$
as well as $x$ and $y$ coordinates. The `*`

operator unpacks the `charge`

object that consists of a list of two objects (as defined during the append).
Thus the lines

```
charge = (1, (-1, 0))
EFieldSingleCharge(*charge)
```

would be equal to

```
EFieldSingleCharge(1, (-1,0))
```

since the `*`

operator simply unpacks the list into it’s components.

Now that we’ve calculated the field strengths it’s time to plot the result. So let’s create a figure and create a subplot with 1 column and 1 row:

```
fig = plt.figure()
splot = fig.add_subplot(111)
```

Each point also should have a color proportional to the strength of our field.
The strength of our field is calculated as $\mid \vec{E}(\vec{r}) \mid = \sqrt{E_x^2 + E_y^2}$.
This can be hidden using the `np.hypot`

function again. Then one should also
apply some kind of scaling. For these type of plots it has been nice to scale
using the logarithm:

```
color = np.log(np.hypot(Ex, Ey))
```

This yields a 2 dimensional array containing the logarithm of the magnitude of the field at each grid point.

Now one can perform a simple streamplot:

```
splot.streamplot(x,y,Ex, Ey, color=color, linewidth=0.5, cmap=plt.cm.inferno, density = 2, arrowstyle='->', arrowsize=1)
```

It might also be nice to see the locations of the point like charge carriers so
one might also want to draw circles. This can simply be done using the `Circle`

from `matplotlib.patches`

. So let’s first import this package:

```
from matplotlib.patches import Circle
```

Then iterate over all charge carriers and paint circles at the given positions. One might also add some color to positive and negative carriers:

```
chargeCols = {
True : '#FF0000',
False : '#0000FF'
}
for q, pos in charges:
splot.add_artist(Circle(pos, 0.05, color = chargeCols[q>0]))
```

Then one can simply set axis labels and limits and finally display the plot:

```
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_xlim(-2,2)
ax.set_ylim(-2,2)
plt.show()
```

The code is available in an GitHub GIST

This article is tagged: Tutorial, Physics, Electrodynamics, Programming, Python

Dipl.-Ing. Thomas Spielauer, Wien (webcomplains389t48957@tspi.at)

This webpage is also available via TOR at http://jugujbrirx3irwyx.onion/