{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Coulomb's law with two-dimensional sources\n", "\n", "\n", "### Coulomb's law for a two-dimsional source:\n", "\n", "$$\n", "{\\bf E}({\\bf r}) = \\frac{1}{4\\pi\\epsilon_0}\\, \\iint\\frac{{\\bf s}}{s^3}\\, \\sigma({\\bf r}^\\prime)\\, da^\\prime.\n", "$$\n", "\n", "where ${\\bf s}$ is Griffiths' script-r vector from source point $dq=\\sigma\\, da^\\prime$ at position \n", "${\\bf r}^\\prime$ to the observation point at ${\\bf r}$, \n", "$\\sigma$ is the surface charge density of the source object, and $da^\\prime$ is the \n", "infinitesimal area on the surface. \n", "\n", "Two dimensional surfaces can be parametrized in terms of two parameters --\n", "lets call them $p_1$ and $p_2$. These parameters might by $x$ and $y$ for planar area at constant $z$, or \n", "$\\theta$ and $\\phi$ for a spherical surface. (This is just a generalization of the parametrization \n", "of a one-dimensional source in terms of the parameter $t$ in the notebook `coulomb2`.) The eventual \n", "integrations will be double integrations over the parameters $p_1$, $p_2$ that are chosen for the \n", "surface in question. We will need several functions that will depend on the surface charge:\n", "\n", "+ $da$ will be of the form $da = \\mbox{factor}\\times dp_1 dp_2$. We need a function that gives\n", "this factor (that will depend on $p_1$ and $p_2$)\n", "\n", "+ we need a function that will convert the parameters $p_1$ and $p_2$ into cartesian coord-nates \n", "\n", "\n", "\n", "It is convenient to work with dimensionless parameters defined in terms of some characteristic length\n", "$R$ and charge $Q$ characterizing the source:\n", "\n", "$$\n", "{\\bf r}^\\ast \\equiv \\frac{\\bf r}{R} \\quad\\quad \n", "{\\bf s}^\\ast \\equiv \\frac{\\bf s}{R}\\quad\\quad \\sigma^\\ast = \\frac{\\sigma}{Q/R^2}\\quad\\quad \n", "\\mbox{and} \\quad\\quad{\\bf E}^\\ast \\equiv \\frac{{\\bf E}}{Q/(4\\pi \\epsilon_0 R^2)}.\n", "$$\n", "\n", "Coulomb's law for 2-d sources written in terms of these parameters is \n", "\n", "$$\n", "{\\bf E}^\\ast = \\iint \\frac{{\\bf s}^\\ast}{{s^\\ast}^3}\\sigma^\\ast\\, d{a^\\prime}^\\ast. \n", "$$\n", "\n", "This is really three separate double integrals, one for each component. Limits of integration must \n", "be determined from analysis of the source.\n", "\n", "All calculations in what follow will be in terms of dimensionless variables; asterisks will \n", "be implicit. The dimensional electric \n", "field can be recovered by multiplying the numerical results below by the factor $Q/(4\\pi\\epsilon_0 R^2)$.\n", "\n", "This notebook uses the quad method from the integrate submodule of scipy\n", "for numerical integration. There are other choices, but this works well (and I like to use scipy\n", "tools when they're available).\n", "\n", "\n", "Marty Ligare, November 2022" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import integrate\n", "import numdifftools as nd\n", "\n", "import matplotlib as mpl \n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Following is an Ipython magic command that puts figures in the notebook.\n", "# For figures in separate windows, comment out following line and uncomment\n", "# the next line\n", "# Must come before defaults are changed.\n", "%matplotlib notebook\n", "#%matplotlib\n", "\n", "# As of Aug. 2017 reverting to 1.x defaults.\n", "# In 2.x text.ustex requires dvipng, texlive-latex-extra, and texlive-fonts-recommended, \n", "# which don't seem to be universal\n", "# See https://stackoverflow.com/questions/38906356/error-running-matplotlib-in-latex-type1cm?\n", "mpl.style.use('classic')\n", " \n", "# M.L. modifications of matplotlib defaults using syntax of v.2.0 \n", "# More info at http://matplotlib.org/2.0.0/users/deflt_style_changes.html\n", "# Changes can also be put in matplotlibrc file, or effected using mpl.rcParams[]\n", "plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures\n", "plt.rc('axes', labelsize=16, titlesize=14)\n", "plt.rc('figure', autolayout = True) # Adjusts supblot parameters for new size" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simple testing: Field of uniformly charged disk (Griffiths Problem 2.25c)\n", "\n", "#### Characterizing the source\n", "Disk with radius $R$ and total charge $Q$ in $x$-$y$ plane centered on the origin. \n", "The obvious choice for length and charge scales are $R$ and $Q$. \n", "The dimensionless surface charge density is then\n", "\n", "\\begin{eqnarray*}\n", "\\sigma^\\ast &\\equiv& \\frac{\\sigma}{Q/R^2} \\\\\n", " &=& \\frac{Q}{\\pi R^2} \\frac{1}{Q/R^2} \\\\\n", " &=& \\frac{1}{\\pi}\n", "\\end{eqnarray*}\n", "\n", "In this example, the parameters $p_1$ and $p_2$ could be chosen to be $x$ and $y$, or else $r$ and $\\phi$,\n", "where $r$ is the radial position in the plane of the disk. In this notebook I'll use polar coordinates:\n", "\n", "$$\n", "p_1 \\longrightarrow r\\quad\\quad\\mbox{and}\\quad\\quad p_2\\longrightarrow \\phi.\n", "$$\n", "\n", "Note that in the cell Definition of functions describing source the arguments are written as\n", "polar variables to make interpretation convenient. The definitions in this cell will depend on the source and \n", "its pararmetrization under consideration. The cell below that, Source-independent functions, should not have to \n", "be modified from problem to problem, so it is written in terms the generic parameters $p_1$ and $p_2$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Definition of functions describing source" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "p1_1 = 0. # Lower limit for integration over r\n", "p1_2 = 1. # Upper limit for integration over r\n", "p2_1 = 0. # Lower limit for integration over phi\n", "p2_2 = 2.*np.pi # Upper limit for integration over r\n", "\n", "def sigma(r, phi):\n", " '''Charge density for disk'''\n", " return 1/np.pi\n", " \n", "def da(r, phi):\n", " '''The factor that multiplies dr*dphi in da = r*dr*dphi'''\n", " return r\n", "\n", "def rp_x(r, phi):\n", " '''Conversion of (r,phi) to cartesian x-component of r-prime vector'''\n", " return r*np.cos(phi)\n", "\n", "def rp_y(r, phi):\n", " '''Conversion of (r,phi) to cartesian y-component of r-prime vector'''\n", " return r*np.sin(phi)\n", "\n", "def rp_z(r, phi):\n", " '''Conversion of (r,phi) to cartesian z-component of r-prime vector'''\n", " return 0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Source-independent functions" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def integrand_x(p1,p2,x,y,z):\n", " '''x-component of field at (x,y,z) due to infinitesimal source dq at p1, p2\n", " sx, sy, sz are components of vector separation between observation point and \n", " source point'''\n", " sx = x - rp_x(p1, p2)\n", " sy = y - rp_y(p1, p2)\n", " sz = z - rp_z(p1, p2)\n", " return sx/(sx**2 + sy**2 + sz**2)**(3/2)*sigma(p1,p2)*da(p1,p2)\n", "\n", "def integrand_y(p1,p2,x,y,z):\n", " '''y-component of field at (x,y,z) due to infinitesimal source dq at p1, p2.\n", " sx, sy, sz are components of vector separation between observation point and \n", " source point'''\n", " sx = x - rp_x(p1, p2)\n", " sy = y - rp_y(p1, p2)\n", " sz = z - rp_z(p1, p2)\n", " return sy/(sx**2 + sy**2 + sz**2)**(3/2)*sigma(p1,p2)*da(p1,p2)\n", "\n", "def integrand_z(p1,p2,x,y,z):\n", " '''z-component of field at (x,y,z) due to infinitesimal source dq at p1, p2.\n", " sx, sy, sz are components of vector separation between observation point and \n", " source point'''\n", " sx = x - rp_x(p1, p2)\n", " sy = y - rp_y(p1, p2)\n", " sz = z - rp_z(p1, p2)\n", " return sz/(sx**2 + sy**2 + sz**2)**(3/2)*sigma(p1,p2)*da(p1,p2)\n", "\n", "def ex(x,y,z):\n", " '''Perform double integration for x-component'''\n", " a = integrate.dblquad(integrand_x, p2_1, p2_2, p1_1, p1_2, args=(x,y,z))\n", " return a[0]\n", "\n", "def ey(x,y,z):\n", " '''Perform double integration for y-component'''\n", " a = integrate.dblquad(integrand_y, p2_1, p2_2, p1_1, p1_2, args=(x,y,z))\n", " return a[0]\n", "\n", "def ez(x,y,z):\n", " '''Perform double integration for z-component'''\n", " a = integrate.dblquad(integrand_z, p2_1, p2_2, p1_1, p1_2, args=(x,y,z))\n", " return a[0]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Nondimensional field at (0,0,1)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-1.634938358441546e-17, -1.3680451329819683e-17, 0.5857864376269051)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ex(0,0,1), ey(0,0,1), ez(0,0,1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "{\\bf E}^\\ast(0,0,1) = 0.5858\\, \\hat{\\bf z}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Dimensional field at (0.0,R)\n", "\n", "$$ \n", "{\\bf E} = 0.5858\\, \\times \\frac{𝑄}{4\\pi\\epsilon_0 R^2}\\, \\hat{\\bf z}\n", "$$\n", "\n", "Analytical result for points on the $z$-axis:\n", "\n", "$$\n", "{\\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}\n", "$$\n", "\n", "giving\n", "\n", "\\begin{eqnarray*}\n", "E_z(0,0,R) &=& 2\\left(1 - \\frac{1}{\\sqrt{2}}\\right)\\frac{Q}{4\\pi\\epsilon_0 R^2}\\\\\n", " &=& 0.5858 \\times \\frac{𝑄}{4\\pi\\epsilon_0 R^2}\n", "\\end{eqnarray*}\n", "\n", "We get the same on-axis results as the analytical technique, but it's now trivial to get field at other points." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### NEW SOURCE\n", "\n", "#### Hemispherical shell with radius $R$ and charge $Q$ resting on $x$-$y$ plane centered on the origin.\n", "\n", "$$\n", "\\sigma = \\frac{Q}{2\\pi R^2} \\longrightarrow \\sigma^\\ast = \\frac{1}{2\\pi}\n", "$$\n", "\n", "This surface can be parametrized in terms of the angles $\\theta$ and $\\phi$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### NEW definition of functions describing source" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "p1_1 = 0.\n", "p1_2 = np.pi/2\n", "p2_1 = 0.\n", "p2_2 = 2.*np.pi\n", "\n", "def sigma(theta, phi):\n", " '''Charge density on hemisphere'''\n", " return 1/np.pi/2\n", " \n", "def da(theta, phi):\n", " '''\n", " The factor that multiplies dtheta*dphi in da = sin(theta)*dtheta*dphi\n", " (Remember that r = 1 on surface)\n", " '''\n", " return np.sin(theta)\n", "\n", "def rp_x(theta, phi):\n", " '''Conversion of (theta,phi) to cartesian x-component of r-prime vector'''\n", " return np.sin(theta)*np.cos(phi)\n", "\n", "def rp_y(theta, phi):\n", " '''Conversion of (theta,phi) to cartesian y-component of r-prime vector'''\n", " return np.sin(theta)*np.sin(phi)\n", "\n", "def rp_z(theta, phi):\n", " '''Conversion of (theta,phi) to cartesian z-component of r-prime vector'''\n", " return np.cos(theta)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Nondimensional field at (0,0,2)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.0, -4.523305928080087e-18, 0.3618033988749892)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ex(0,0,2),ey(0,0,2),ez(0,0,2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "{\\bf E}^\\ast(0,0,2) = 0.3618\\, \\hat{\\bf z}\n", "$$\n", "\n", "#### Is this correct? \n", "I don't have a formula at hand, but I can think of one consistency check. Considering adding a second \n", "inverted hemisphere below the $x$-$y$ plane, making a complete hemispherical shell with total charge $2Q$.\n", "And note that by symmetry the fields due the bottom and top hemispheres are related:\n", "\n", "$$ \n", "{\\bf E}_{\\rm bottom}(0,0,z) = -{\\bf E}_{\\rm top}(0,0,-z)\n", "$$\n", "\n", "The total field at $(0,0,2)$ should thus be\n", "\n", "$$\n", "{\\bf E}_{\\rm sphere}(0,0,2) = {\\bf E}_{\\rm top}(0,0,2) - {\\bf E}_{\\rm top}(0,0,-2)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Nondimensional field of sphere at (0,0,-2)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.49999999999999967" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ez(0,0,2) - ez(0,0,-2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "{\\bf E}_{\\rm sphere}^\\ast(0,0,-2) = \\frac{1}{2}\\, \\hat{\\bf z}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Dimensional field of sphere at (0,0,-2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "{\\bf E}_{\\rm sphere}(0,0,-2) = \\frac{Q}{4\\pi\\epsilon_0 R^2} \\times \\frac{1}{2}\\, \\hat{\\bf z}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The field of the spherical shell should be that of a point charge $2Q$ located at the origin, which this is." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Another consistency check: Field anywhere within the spherical shell should be zero. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.1102230246251565e-16" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ez(0,0,0.9) - ez(0,0,-0.9)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Version Information \n", "`version_information` is from J.R. Johansson (jrjohansson at gmail.com)
\n", "See Introduction to scientific computing with Python:
\n", "http://nbviewer.jupyter.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-0-Scientific-Computing-with-Python.ipynb
\n", "for more information and instructions for package installation.
\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "%load_ext version_information\n", "\n", "#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "application/json": { "Software versions": [ { "module": "Python", "version": "3.7.13 64bit [GCC 7.5.0]" }, { "module": "IPython", "version": "7.31.1" }, { "module": "OS", "version": "Linux 4.9.0 9 amd64 x86_64 with debian 9.13" }, { "module": "numpy", "version": "1.21.5" }, { "module": "scipy", "version": "1.7.3" }, { "module": "numdifftools", "version": "0.9.39" }, { "module": "matplotlib", "version": "3.5.2" } ] }, "text/html": [ "
SoftwareVersion
Python3.7.13 64bit [GCC 7.5.0]
IPython7.31.1
OSLinux 4.9.0 9 amd64 x86_64 with debian 9.13
numpy1.21.5
scipy1.7.3
numdifftools0.9.39
matplotlib3.5.2
Mon Nov 14 13:59:10 2022 EST
" ], "text/latex": [ "\\begin{tabular}{|l|l|}\\hline\n", "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n", "Python & 3.7.13 64bit [GCC 7.5.0] \\\\ \\hline\n", "IPython & 7.31.1 \\\\ \\hline\n", "OS & Linux 4.9.0 9 amd64 x86\\_64 with debian 9.13 \\\\ \\hline\n", "numpy & 1.21.5 \\\\ \\hline\n", "scipy & 1.7.3 \\\\ \\hline\n", "numdifftools & 0.9.39 \\\\ \\hline\n", "matplotlib & 3.5.2 \\\\ \\hline\n", "\\hline \\multicolumn{2}{|l|}{Mon Nov 14 13:59:10 2022 EST} \\\\ \\hline\n", "\\end{tabular}\n" ], "text/plain": [ "Software versions\n", "Python 3.7.13 64bit [GCC 7.5.0]\n", "IPython 7.31.1\n", "OS Linux 4.9.0 9 amd64 x86_64 with debian 9.13\n", "numpy 1.21.5\n", "scipy 1.7.3\n", "numdifftools 0.9.39\n", "matplotlib 3.5.2\n", "Mon Nov 14 13:59:10 2022 EST" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "version_information numpy, scipy, numdifftools, matplotlib" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.13" } }, "nbformat": 4, "nbformat_minor": 1 }