#### Author: kvan637 (page 1 of 9)

Chapter 9 of A guided tour of mathematical methods for the physical sciences introduces the theorem of Stokes. In Chapter 7, you saw a special vector field, for which the divergence and curl are zero. This gave rise to a model for Hurricane winds and the magnetic field around a wire, captured in one of Maxwell’s equations representing the Biot-Savart Law.

### Biot-Savart LawÂ¶

The current in an long wire creates a magnetic field around it. If we align the wire with the z-axis, the field cannot vary with $z$, or $\phi$, so that ${\bf B}(r,\phi,z) = {\bf B}(r)$. The divergence of magnetic fields, and of incompressible flows, is zero, we found in Chapter 7 that

$${\bf B} = \frac{A}{r} \hat{\mathbf \phi}$$

But how do we find the value of $A$? This is where we apply the theorem of Stokes to the Biot-Savart Law:

$$\nabla\times{\bf B} = \mu_0{\bf J}.$$
Integrating this
over a disk $d{\bf S}$ centered on the wire of radius $r$ in the plane perpendicular to the wire, we have:
$$\iint_S \nabla\times{\bf B} \cdot d{\bf S} = \iint_S\mu_0{\bf J}\cdot d{\bf S}.$$

For the left hand side, we apply Stokes Theorem:
$$\iint_S \nabla\times{\bf B} \cdot d{\bf S} = \oint_C {\bf B}\cdot d{\bf r} = 2\pi r{\bf B},$$
because $\bf B$ is constant for each circular path $C$ around the wire.

For the right hand side, we recognize that as long as the radius of the disk $r$ is larger than the thickness of our wire:
$$\iint_S\mu_0{\bf J}\cdot d{\bf S}= \mu_0 I,$$
where $I$ is the current. Putting these two results together, we find that the magnetic field is
$${\bf B}(r) = \frac{\mu_0 I}{2\pi r} \hat{\mathbf \phi}.$$ Breaking the radius $r$ and the unit vector $\hat{\phi}$ down in an x- and y-component (see Problem c of Section 7.2 of the book), we define a python function for the three Cartesian components of the magnetic field:

InÂ [22]:
def B(I,x,y):
mu0 = 1.26 * 10**(-6)
r = np.sqrt((x)**2+(y)**2)
c = mu0*I/(2*np.pi)
Bx = -y*c/r**2
By = x*c/r**2
Bz = z*0
return Bx,By,Bz


Then, we set up a grid of points:

InÂ [62]:
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
%matplotlib notebook

X = np.linspace(-1,1,12)
Y = np.linspace(-1,1,12)
Z = np.linspace(-1,1,12)
x,y,z = np.meshgrid(X,Y,Z)


The magnetic field on the grid points around the wire with current $I$ is:

InÂ [63]:
I = 200000
Bx,By,Bz = B(I,x,y)


And finally, we plot of the magnetic field (our first 3D plot with matplotlib!):

InÂ [64]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.quiver(x,y,z,Bx,By,Bz)
ax.plot([0, 0],[0, 0],[-1,1],linewidth=3,color='r')
plt.xlabel('x')
plt.ylabel('y')
plt.show()


You can now rotate the graph to look at this field from all angles, and/or change the current. From the top, it should look just like the flow in the python notebook for Chapter 7. What happens when we run the current in the other direction? Chapter 9 concludes with applications of Stokes’ Theorem to introduce Lenz’ Law, the mysterious quantum mechanical Aharanov-Bohm effect, and the impact of vortices at the tips of sails and wings. Let us now return to the overview of jupyter notebooks!

Chapter 8 of A guided tour of mathematical methods for the physical sciences introduces the theorem of Gauss. In Chapter 6, we stood at its doorstep already, when we explored the divergence of a vector field ${\bf v}$. All of this was described by the flux, or outward flow, through an elementary volume $dV$, with a surface $d{\bf S}$, defined by its normal vectors, as here:

Section 6.1 showed that the flux — the flow across an infinitesimal surface — is $$d\Phi_v = {\bf v} \cdot d{\bf S}.$$

Section 6.2 says that the same flux is also $$d\Phi_v = \left(\nabla\cdot {\bf v}\right)dV.$$

Chapter 7 equates these two to state the Theorem of Gauss in integral form:
$$\oint_S {\bf v} \cdot d{\bf S} =\int_V \left(\nabla\cdot {\bf v}\right)dV.$$

### Gauss’ Theorem across different dimensionsÂ¶

The image above implies a 3D volume, but the 2D version we already explored in Chapter 6, and returns in the next Chapter about the Theorem of Stokes too. There, dV will be an area, and the surface reduces to a line integral. In 1D, the “dV” is reduced to a line, and the “dS” its two end points. Let’s align the line with the x-axis, and then the 1D Gauss’ Theorem is the familiar definition of integration:
$$v_x(b) – v_x(a) = \int_a^b \frac{\partial v_x}{\partial x} dx.$$

### Gravity in a borehole, or the electric field inside a sphereÂ¶

Let us consider the vector field that represents gravitational acceleration from a mass with ${\bf v}({\bf r}) = {\bf g(\bf{r})} = g(r)\hat{{\bf r}}$. Note that the physics for the electric field from a spherically symmetric charge distribution would lead to the same results!

In Section 8.2, you will have worked your way through the problem of the acceleration outside such a mass:

$${\bf g}(r) = \frac{GM}{r^2}\hat{\bf r}.$$

Actually, Gauss’ Theorem says this result is accurate for any spherically symmetric mass distribution!

Inside a sphere with radius $R$ that has a constant density, we found that

$${\bf g}(r) = \frac{GMr}{R^3}\hat{\bf r}.$$

The combo of these two representations would describe aspects of the gravitational experience down a borehole in the Earth or that of the explorers in the movie The Core:

InÂ [34]:
import numpy as np
import matplotlib.pyplot as plt

# The constants in SI units:
G = 6.67*10**(-11) # The Gravitational constant in m^3 s^(-2) kg^-1
M = 6*10**24 # The Earth's mass in kg
R = 6370*10**3 # The radius the Earth in m

#  A function to compute the magnitude of the (radial) gravitational acceleration, inside or outside the Earth:
def gravity(r):
if r>R:
g = M*G/r**2
else:
g =  M*G*r/R**3
return g

# setting up an array with points from the centre (0) to well beyond the radius of the Earth:
rs= np.arange(0,8*10**6,1000)

# compute the magntiude of the gravitational acceleration for all points in the array rs:
gs = []
for r in rs:
gs.append(gravity(r))

# plotting the results:
plt.plot(gs,rs/1000)
plt.axhline(R/1000,linestyle='dashed',color='r')
plt.xlabel("Gravitational acceleration (m/s$^2$)")
plt.ylabel("Distance to the centre of the Earth (km)")
plt.show()



We see that gravity at the centre of the Earth is zero, as any body there would be pulled in every direction equally, canceling all contributions. As we move to the surface, the gravitational field increases linearly in strength for this mass with constant density, until we reach the Earth’s surface. Then, with no more mass contributions to incorporate, gravitational acceleration decreases as $r^{-2}$.

The observed gravitational accelaration at the Earth’s surface does not vary much from this model: around 9.8 m/s$^2$, or 980 Gals). Slightly more accurate would be an ellipsoidal shape to the Earth, where gravity at the poles is about 9.83 and at the equator 9.78 m/s$^2$. and finally here is a picture from the GRACE mission (figure from wikipedia), showing the variations in $g$ from this best fitting ellipsoid are on the order of 100 mGals:

This means the Earth’s density is pretty close to spherically symmetric. So the g values for r > 6730 km we display here are quite realistic. Remember in this calculation that the values of g inside the Earth assumed that the density is constant. For Earth, density generally increases with depth, according to the PREM model (figure from wikipedia):

We’ll leave extending our calculations here to a more realistic distribution of density (maybe start with a linear increase from 2 to 12 g/cm$^3$) to the reader of this notebook, while we return to the overview of jupyter notebooks!