1.3: Modularizing the Code
1.3.1: Designing a Potential Function
We could continue altering the above code in a straightforward way. For example, we could add more particles by adding variables
x1
,
x2
,
q1
,
q2
, and so forth, and altering our formula for computing
phi
. However, this is not very satisfactory: each time we want to consider a new collection of particle positions or charges, or change the number of particles, we would have to re-write the program's internal "logic"—i.e., the part that computes the potentials. In programming terminology, our program is insufficiently "modular". Ideally, we want to isolate the part of the program that computes the potential from the part that specifies the numerical inputs to the calculation, like the positions and charges.
To modularize the code, let's define a function that computes the potential of an arbitrary set of charged particles, sampled at an arbitrary set of positions. Such a function would need three sets of inputs:
- An array of particle positions \(\vec{x}\equiv [x_{0}, \cdots, x_{N-1}]\). (Don't get confused, by the way: we are using these \(N\) numbers to refer to the positions of \(N\) particles in a 1D space, not the position of a single particle in an \(N\)-dimensional space.)
- An array of particle charges \(\vec{q}\equiv [q_{0}, \cdots, q_{N-1}]\).
- An array of sampling points \(\vec{X}\equiv [X_{0}, \cdots, X_{M-1}]\), which are the points where we want to know \(\phi(X)\).
The number of particles, \(N\), and the number of sampling points, \(M\), should be arbitrary positive integers. Furthermore, \(N\) and \(M\) need not be equal.
The function we intend to write must compute the array
\[\begin{bmatrix}\phi(X_{0})\\ \phi(X_{1}) \\ \vdots \\ \phi(X_{M-1})\end{bmatrix}\]
which contains the value of the total electric potential at each of the sampling points. The total potential can be written as the sum of contributions from all particles. Let us define \(\phi_{j}(x)\) as the potential produced by particle \(j\):
\[\phi_{j(x)}\equiv\frac{q_{j}}{|x-x_{j}|}\]
Then the total potential is
\[\begin{bmatrix}\phi(X_0)\\ \phi(X_1) \\ \vdots \\ \phi(X_{M-1})\end{bmatrix} = \begin{bmatrix}\phi_0(X_0)\\ \phi_0(X_1) \\ \vdots \\ \phi_0(X_{M-1})\end{bmatrix} + \begin{bmatrix}\phi_1(X_0)\\ \phi_1(X_1) \\ \vdots \\ \phi_1(X_{M-1})\end{bmatrix} + \cdots + \begin{bmatrix}\phi_{N-1}(X_0)\\ \phi_{N-1}(X_1) \\ \vdots \\ \phi_{N-1}(X_{M-1})\end{bmatrix}.\]
1.3.2 Writing the Program
Let's code this up. Return to the file
potentials.py
, and
delete the entire contents of the file
. Then replace it with the following:
from scipy import * import matplotlib.pyplot as plt ## Return the potential at measurement points X, due to particles ## at positions xc and charges qc. xc, qc, and X must be 1D arrays, ## with xc and qc of equal length. The return value is an array ## of the same length as X, containing the potentials at each X point. def potential(xc, qc, X): M = len(X) N = len(xc) phi = zeros(M) for j in range(N): phi += qc[j] / abs(X - xc[j]) return phi charges_x = array([0.2, -0.2]) charges_q = array([1.5, -0.1]) xplot = linspace(-3, 3, 500) phi = potential(charges_x, charges_q, xplot) plt.plot(xplot, phi) pmin, pmax = -50, 50 plt.ylim(pmin, pmax) plt.show()
When typing or pasting the above into your file, be sure to preserve the indentation (i.e., the number of spaces at the beginning of each line). Indentation is important in Python; as we'll see, it's used to determine program structure. Now save and run the program again:
-
In the Windows GUI, type
F5
in the editing window showingpotentials.py
. -
On GNU/Linux, type
python -i potentials.py
from the command line. -
Alternatively, from the Python command line, type
import potentials
, which will load and run yourpotentials.py
file.
You should now see a figure like the one on the right, showing the electric potential produced by two particles, one at position \(x_{0}=0.2\) with charge \(q_{0}=1.5\) and the other at position \(x_{1}=-0.2\) with charge \(q_{1}=-0.1\).
There are less than \(20\) lines of actual code in the above program, but they do quite a lot of things. Let's go through them in turn:
Module Imports
The first two lines import the Scipy and Matplotlib modules, for use in our program. We have not yet explained how importing works, so let's do that now.
from scipy import * import matplotlib.pyplot as plt
Each Python module, including Scipy and Matplotlib, defines a variety of functions and variables. If you use multiple modules, you might have a situation where, say, two different modules each define a function with the same name, but doing entirely different things. That would be Very Bad. To help avoid this, Python implements a concept called a namespace . Suppose you import a module (say Scipy) like this:
import scipy
One of the functions defined by Scipy is
linspace
, which we have already seen. This function was defined by the
scipy
module, and lies inside the
scipy
namespace. As a result, when you import the Scipy module using the
import scipy
line, you have to call the
linspace
function like this:
x = scipy.linspace(-3, 3, 500)
The
scipy.
in front says that you're referring to the
linspace
function that was defined in the
scipy
namespace. (Note: the online documentation for
linspace
refers to it as
numpy.linspace
, but the exact same function is also present in the
scipy
namespace. In fact, all
numpy.*
functions are replicated in the
scipy
namespace. So unless stated otherwise, we only have to import
scipy
.)
We will be using a lot of functions that are defined in the
scipy
namespace. Since it would be annoying to have to keep typing
scipy.
all over the place, we opt to use a slightly different import statement:
from scipy import *
This imports all the functions and variables in the
scipy
namespace directly into your program's namespace. Therefore, you can just call
linspace
, without the
scipy.
prefix. Obviously, you don't want to do this for every module you use, otherwise you'll end up with the name-clashing problem we alluded to earlier! The only module we'll use this shortcut with is
scipy
.
Another way to avoid having to type long prefixes is shown by this line:
import matplotlib.pyplot as plt
This imports the
matplotlib.pyplot
module (i.e., the
pyplot
module which is nested inside the
matplotlib
module). That's where
plot
,
show
, and other plotting functions are defined. The
as plt
in the above line says that we will refer to the
matplotlib.pyplot
namespace as the short form
plt
instead. Hence, instead of calling the
plot
function like this:
matplotlib.pyplot.plot(x, y)
we will call it like this:
plt.plot(x, y)
Comments
Let's return to the program we were looking at earlier. The next few lines, beginning with
#
, are "comments". Python ignores the
#
character and everything that follows it, up to the end of the line. Comments are very important, even in simple programs like this.
When you write your own programs, please remember to include comments. You don't need a comment for every line of code—that would be excessive—but at a minimum, each function should have a comment explaining what it does, and what the inputs and return values are.
Function Definition
Now we get to the
function definition
for the function named
potential
, which is the function that computes the potential:
def potential(xc, qc, X): M = len(X) N = len(xc) phi = zeros(M) for j in range(N): phi += qc[j] / abs(X - xc[j]) return phi
The first line, beginning with
def
, is a
function header
. This function header states that the function is named
potential
, and that it has three inputs. In computing terminology, the inputs that a function accepts are called
parameters
. Here, the parameters are named
xc
,
qc
and
X
. As explained by the comments, we intend to use these for the positions of the particles, the charges of the particles, and the positions at which to measure the potential, respectively.
The function definition consists of the function header, together the rest of the indented lines below it. The function definition terminates once we get to a line which is at the same indentation level as the function header. (That terminating line is considered a separate line of code, which is not part of the function definition).
By convention, you should use 4 spaces per indentation level.
The indented lines below the function header are called the function body . This is the code that is run each time the function is called. In this case, the function body consists of six lines of code, which are intended to compute the total electric potential, according to the procedure that we have outlined in the preceding section:
-
The first two lines define two helpful variables,
M
andN
. Their values are set to the lengths of theX
andxc
arrays, respectively. -
The next line calls the
zeros
function. The input tozeros
isM
, the length of theX
array (i.e., our function's third parameter). Therefore,zeros
returns an array, of the same same length asX
, with every element set to 0.0. For now, this represents the electric potential in the absence of any charges. We give this array the namephi
. -
The function then iterates over each of the particles and add up its contribution to the potential, using a construct known as a
for
loop . The codefor j in range(N):
is the loop's "header line", and the next line, indented 4 spaces deeper than the header line, is the "body" of the loop.The header line states that we should run the loop body several times, with the variable
j
set to different values during each run. The values ofj
to loop over are given byrange(N)
. This is a function call to the range function, withN
(the number of electric charges) as the input. Therange(N)
function call returns a sequence specifyingN
successive integers, from 0 toN-1
, inclusive. (Note that the last value in the sequence isN-1
, notN
. Because we start from 0, this means that there is a total ofN
integers in the sequence. Also, callingrange(N)
is the same as callingrange(0,N)
.) -
For each
j
, we computeqc[j] / abs(X - xc[j])
. This is an array whose elements are the values of the electric potential at the set of positionsX
, arising from the individual particle\(j\). In mathematical terms, we are calculating
\[\phi_j(X) \equiv \frac{q_j}{|X - x_j|}\]using the array of positions
X
. We then add this array tophi
. Once this is done for allj
, the arrayphi
will contain the desired total potential,
\[\phi(X) = \sum_{j=0}^{N-1}\phi_j(X).\] -
Finally, we call
return
to specify the function's output, or return value . This is the arrayphi
.
Top-Level Code: Numerical Constants
After the function definition comes the code to use the function:
charges_x = array([0.2, -0.2]) charges_q = array([1.5, -0.1]) xplot = linspace(-3, 3, 500)
Like the import statements at the beginning of the program, these lines of code lie at
top level
, i.e., they are not indented. The function header which defines the
potential
function is also at top level. Running a Python program consists of running its top level code, in sequence.
The above lines define variables to store some numerical constants. In the first two lines,
charges_x
and
charges_q
variables store the numerical values of the positions and charges we are interested in. These are initialized using the
array
function. You may be wondering why the
array
function call has square brackets nested in commas. We'll explain later, in part 2 of the tutorial.
On the third line, the
linspace
function call returns an array, whose contents are initialized to the 500 numbers between -3 and 3 (inclusive).
Next, we call the
potential
function, passing
charges_x
,
charges_q
and
xplot
as the inputs:
phi = potential(charges_x, charges_q, xplot)
These inputs provide the values of the function definition's parameters
xc
,
qc
, and
X
respectively. The return value of the function call is an array containing the total potential, evaluated at each of the positions specified in
xplot
. This return value is stored as the variable named
phi
.
Plotting
Finally, we create the plot:
plt.plot(xplot, phi) pmin, pmax = -50, 50 plt.ylim(pmin, pmax) plt.show()
We have already seen how the
plot
and
show
functions work. Here, prior to calling
plt.show
, we have added two extra lines to make the potential curve is more legible, by adjust the plot's y-axis bounds. The
ylim
function accepts two parameters, the lower and upper bounds of the y-axis. In this case, we set the bounds to -50 and 50 respectively. There is an
xlim
function to do the same for the x-axis.
Notice that in the line
pmin, pmax = -50, 50
, we set two variables (
pmin
and
pmax
) on the same line. This is a little "syntactic sugar" to make the code a little easier to read. It's equivalent to having two separate lines, like this:
pmin = -50 pmax = 50
We'll explain how this construct works in the next part of the tutorial.