Rendering static two dimensional electric fields caused by point charges using Python

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).

Software requirements

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

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.

Since we know from Coulombs law that the force between two charges $q_1$ and $q_2$ is given by

we can simply derive the equation for the field caused by a single point like charge $q$ to be

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:

This is then exactly the equation that we’re going to apply for our simple plotting function.

The 2D implementation

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

The code is available in an GitHub GIST

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