Coulomb's law with two-dimensional sources

Coulomb's law for a two-dimsional source:

$$ {\bf E}({\bf r}) = \frac{1}{4\pi\epsilon_0}\, \iint\frac{{\bf s}}{s^3}\, \sigma({\bf r}^\prime)\, da^\prime. $$

where ${\bf s}$ is Griffiths' script-r vector from source point $dq=\sigma\, da^\prime$ at position ${\bf r}^\prime$ to the observation point at ${\bf r}$, $\sigma$ is the surface charge density of the source object, and $da^\prime$ is the infinitesimal area on the surface.

Two dimensional surfaces can be parametrized in terms of two parameters -- lets call them $p_1$ and $p_2$. These parameters might by $x$ and $y$ for planar area at constant $z$, or $\theta$ and $\phi$ for a spherical surface. (This is just a generalization of the parametrization of a one-dimensional source in terms of the parameter $t$ in the notebook coulomb2.) The eventual integrations will be double integrations over the parameters $p_1$, $p_2$ that are chosen for the surface in question. We will need several functions that will depend on the surface charge:

• $da$ will be of the form $da = \mbox{factor}\times dp_1 dp_2$. We need a function that gives this factor (that will depend on $p_1$ and $p_2$)

• we need a function that will convert the parameters $p_1$ and $p_2$ into cartesian coord-nates

It is convenient to work with dimensionless parameters defined in terms of some characteristic length $R$ and charge $Q$ characterizing the source:

$$ {\bf r}^\ast \equiv \frac{\bf r}{R} \quad\quad {\bf s}^\ast \equiv \frac{\bf s}{R}\quad\quad \sigma^\ast = \frac{\sigma}{Q/R^2}\quad\quad \mbox{and} \quad\quad{\bf E}^\ast \equiv \frac{{\bf E}}{Q/(4\pi \epsilon_0 R^2)}. $$

Coulomb's law for 2-d sources written in terms of these parameters is

$$ {\bf E}^\ast = \iint \frac{{\bf s}^\ast}{{s^\ast}^3}\sigma^\ast\, d{a^\prime}^\ast. $$

This is really three separate double integrals, one for each component. Limits of integration must be determined from analysis of the source.

All calculations in what follow will be in terms of dimensionless variables; asterisks will be implicit. The dimensional electric
field can be recovered by multiplying the numerical results below by the factor $Q/(4\pi\epsilon_0 R^2)$.

This notebook uses the quad method from the integrate submodule of scipy for numerical integration. There are other choices, but this works well (and I like to use scipy tools when they're available).

Marty Ligare, November 2022

import numpy as np
from scipy import integrate
import numdifftools as nd

import matplotlib as mpl  
import matplotlib.pyplot as plt

# Following is an Ipython magic command that puts figures in the  notebook.
# For figures in separate windows, comment out following line and uncomment
# the next line
# Must come before defaults are changed.
%matplotlib notebook
#%matplotlib

# As of Aug. 2017 reverting to 1.x defaults.
# In 2.x text.ustex requires dvipng, texlive-latex-extra, and texlive-fonts-recommended, 
# which don't seem to be universal
# See https://stackoverflow.com/questions/38906356/error-running-matplotlib-in-latex-type1cm?
mpl.style.use('classic')
        
# M.L. modifications of matplotlib defaults using syntax of v.2.0 
# More info at http://matplotlib.org/2.0.0/users/deflt_style_changes.html
# Changes can also be put in matplotlibrc file, or effected using mpl.rcParams[]
plt.rc('figure', figsize = (6, 4.5))            # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True)             # Adjusts supblot parameters for new size

Simple testing: Field of uniformly charged disk (Griffiths Problem 2.25c)

Characterizing the source

Disk with radius $R$ and total charge $Q$ in $x$-$y$ plane centered on the origin. The obvious choice for length and charge scales are $R$ and $Q$. The dimensionless surface charge density is then

\begin{eqnarray*} \sigma^\ast &\equiv& \frac{\sigma}{Q/R^2} \\ &=& \frac{Q}{\pi R^2} \frac{1}{Q/R^2} \\ &=& \frac{1}{\pi} \end{eqnarray*}

In this example, the parameters $p_1$ and $p_2$ could be chosen to be $x$ and $y$, or else $r$ and $\phi$, where $r$ is the radial position in the plane of the disk. In this notebook I'll use polar coordinates:

$$ p_1 \longrightarrow r\quad\quad\mbox{and}\quad\quad p_2\longrightarrow \phi. $$

Note that in the cell Definition of functions describing source the arguments are written as polar variables to make interpretation convenient. The definitions in this cell will depend on the source and its pararmetrization under consideration. The cell below that, Source-independent functions, should not have to be modified from problem to problem, so it is written in terms the generic parameters $p_1$ and $p_2$.

Definition of functions describing source

p1_1 = 0.       # Lower limit for integration over r
p1_2 = 1.       # Upper limit for integration over r
p2_1 = 0.       # Lower limit for integration over phi
p2_2 = 2.*np.pi # Upper limit for integration over r

def sigma(r, phi):
    '''Charge density for disk'''
    return 1/np.pi
    
def da(r, phi):
    '''The factor that multiplies dr*dphi in da = r*dr*dphi'''
    return r

def rp_x(r, phi):
    '''Conversion of (r,phi) to cartesian x-component of r-prime vector'''
    return r*np.cos(phi)

def rp_y(r, phi):
    '''Conversion of (r,phi) to cartesian y-component of r-prime vector'''
    return r*np.sin(phi)

def rp_z(r, phi):
    '''Conversion of (r,phi) to cartesian z-component of r-prime vector'''
    return 0

Source-independent functions

def integrand_x(p1,p2,x,y,z):
    '''x-component of field at (x,y,z) due to infinitesimal source dq at p1, p2
    sx, sy, sz are components of vector separation between observation point and 
    source point'''
    sx = x - rp_x(p1, p2)
    sy = y - rp_y(p1, p2)
    sz = z - rp_z(p1, p2)
    return sx/(sx**2 + sy**2 + sz**2)**(3/2)*sigma(p1,p2)*da(p1,p2)

def integrand_y(p1,p2,x,y,z):
    '''y-component of field at (x,y,z) due to infinitesimal source dq at p1, p2.
    sx, sy, sz are components of vector separation between observation point and 
    source point'''
    sx = x - rp_x(p1, p2)
    sy = y - rp_y(p1, p2)
    sz = z - rp_z(p1, p2)
    return sy/(sx**2 + sy**2 + sz**2)**(3/2)*sigma(p1,p2)*da(p1,p2)

def integrand_z(p1,p2,x,y,z):
    '''z-component of field at (x,y,z) due to infinitesimal source dq at p1, p2.
    sx, sy, sz are components of vector separation between observation point and 
    source point'''
    sx = x - rp_x(p1, p2)
    sy = y - rp_y(p1, p2)
    sz = z - rp_z(p1, p2)
    return sz/(sx**2 + sy**2 + sz**2)**(3/2)*sigma(p1,p2)*da(p1,p2)

def ex(x,y,z):
    '''Perform double integration for x-component'''
    a = integrate.dblquad(integrand_x, p2_1, p2_2, p1_1, p1_2, args=(x,y,z))
    return a[0]

def ey(x,y,z):
    '''Perform double integration for y-component'''
    a = integrate.dblquad(integrand_y, p2_1, p2_2, p1_1, p1_2, args=(x,y,z))
    return a[0]

def ez(x,y,z):
    '''Perform double integration for z-component'''
    a = integrate.dblquad(integrand_z, p2_1, p2_2, p1_1, p1_2, args=(x,y,z))
    return a[0]

Nondimensional field at (0,0,1)

ex(0,0,1), ey(0,0,1), ez(0,0,1)

(-1.634938358441546e-17, -1.3680451329819683e-17, 0.5857864376269051)

$$ {\bf E}^\ast(0,0,1) = 0.5858\, \hat{\bf z} $$

Dimensional field at (0.0,R)

$$ {\bf E} = 0.5858\, \times \frac{𝑄}{4\pi\epsilon_0 R^2}\, \hat{\bf z} $$

Analytical result for points on the $z$-axis:

$$ {\bf E}(0,0,z) =\frac{Q}{2\pi\epsilon_0 R^2} \left(1 - \frac{1}{\sqrt{1+\frac{R^2}{z^2}}}\right)\, \hat{\bf z} $$

giving

\begin{eqnarray*} E_z(0,0,R) &=& 2\left(1 - \frac{1}{\sqrt{2}}\right)\frac{Q}{4\pi\epsilon_0 R^2}\\ &=& 0.5858 \times \frac{𝑄}{4\pi\epsilon_0 R^2} \end{eqnarray*}

We get the same on-axis results as the analytical technique, but it's now trivial to get field at other points.

NEW SOURCE

Hemispherical shell with radius $R$ and charge $Q$ resting on $x$-$y$ plane centered on the origin.

$$ \sigma = \frac{Q}{2\pi R^2} \longrightarrow \sigma^\ast = \frac{1}{2\pi} $$

This surface can be parametrized in terms of the angles $\theta$ and $\phi$

NEW definition of functions describing source

p1_1 = 0.
p1_2 = np.pi/2
p2_1 = 0.
p2_2 = 2.*np.pi

def sigma(theta, phi):
    '''Charge density on hemisphere'''
    return 1/np.pi/2
    
def da(theta, phi):
    '''
    The factor that multiplies dtheta*dphi in da = sin(theta)*dtheta*dphi
    (Remember that r = 1 on surface)
    '''
    return np.sin(theta)

def rp_x(theta, phi):
    '''Conversion of (theta,phi) to cartesian x-component of r-prime vector'''
    return np.sin(theta)*np.cos(phi)

def rp_y(theta, phi):
    '''Conversion of (theta,phi) to cartesian y-component of r-prime vector'''
    return np.sin(theta)*np.sin(phi)

def rp_z(theta, phi):
    '''Conversion of (theta,phi) to cartesian z-component of r-prime vector'''
    return np.cos(theta)

Nondimensional field at (0,0,2)

ex(0,0,2),ey(0,0,2),ez(0,0,2)

(0.0, -4.523305928080087e-18, 0.3618033988749892)

$$ {\bf E}^\ast(0,0,2) = 0.3618\, \hat{\bf z} $$

Is this correct?

I don't have a formula at hand, but I can think of one consistency check. Considering adding a second inverted hemisphere below the $x$-$y$ plane, making a complete hemispherical shell with total charge $2Q$. And note that by symmetry the fields due the bottom and top hemispheres are related:

$$ {\bf E}_{\rm bottom}(0,0,z) = -{\bf E}_{\rm top}(0,0,-z) $$

The total field at $(0,0,2)$ should thus be

$$ {\bf E}_{\rm sphere}(0,0,2) = {\bf E}_{\rm top}(0,0,2) - {\bf E}_{\rm top}(0,0,-2) $$

Nondimensional field of sphere at (0,0,-2)

ez(0,0,2) - ez(0,0,-2)

0.49999999999999967

$$ {\bf E}_{\rm sphere}^\ast(0,0,-2) = \frac{1}{2}\, \hat{\bf z} $$

Dimensional field of sphere at (0,0,-2)

$$ {\bf E}_{\rm sphere}(0,0,-2) = \frac{Q}{4\pi\epsilon_0 R^2} \times \frac{1}{2}\, \hat{\bf z} $$

The field of the spherical shell should be that of a point charge $2Q$ located at the origin, which this is.

Another consistency check: Field anywhere within the spherical shell should be zero.

ez(0,0,0.9) - ez(0,0,-0.9)

1.1102230246251565e-16

Version Information

version_information is from J.R. Johansson (jrjohansson at gmail.com)
See Introduction to scientific computing with Python:
http://nbviewer.jupyter.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-0-Scientific-Computing-with-Python.ipynb
for more information and instructions for package installation.

If version_information has been installed system wide (as it has been on Bucknell linux computers with shared file systems), continue with next cell as written. If not, comment out top line in next cell and uncomment the second line.

%load_ext version_information

#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py

version_information numpy, scipy, numdifftools, matplotlib

Software	Version
Python	3.7.13 64bit [GCC 7.5.0]
IPython	7.31.1
OS	Linux 4.9.0 9 amd64 x86_64 with debian 9.13
numpy	1.21.5
scipy	1.7.3
numdifftools	0.9.39
matplotlib	3.5.2
Mon Nov 14 13:59:10 2022 EST