---
jupyter:
  description: 'The Helmholtz coil with planex2d and NGSolve: two coaxial windings in the plane through their axis, the
    center field against (4/5)^(3/2) μ0 n I / R and against the exact field of the windings on the axis, and the
    field along the axis for three spacings, flattest to the fourth order when the spacing is the radius.'
---

# The Helmholtz coil

Why are the two coils of a Helmholtz pair exactly one radius apart? Two equal coils on one axis, carrying the same
current the same way, add their fields. In the middle between them, the odd derivatives of the field along the axis
vanish by symmetry, and the second derivative vanishes when the spacing equals the radius. The field is then uniform
to the fourth order,

$$\frac{B(z)}{B_0} \approx 1 - \frac{144}{125} \left(\frac{z}{R}\right)^4, \qquad B_0 = \left(\frac{4}{5}\right)^{3/2} \frac{\mu_0 n I}{R},$$

for coils of radius $R$ with $n$ turns carrying the current $I$. Such pairs cancel the Earth's field around an
experiment or set up a known field to calibrate magnetic sensors.

The example models two windings of 100 turns each, with a square section of 1 cm on a radius of 15 cm, in the plane
through their axis. NGSolve solves for their field, which is checked against the exact field of the windings on the
axis, and the field along the axis for three spacings shows why the radius is the right one.

This page is part of a collection of [field and simulation examples](fields-and-simulation.md) using planex2d.

![The field of a Helmholtz pair in the plane through its axis, uniform within a percent in the cream region](../_static/gallery/helmholtz-coil.png){.result}

```python
import math

import matplotlib.pyplot as plt
import ngsolve as ngs
import numpy as np
import planex as px
from scipy.integrate import dblquad
```

```python nbsphinx="hidden"
from docs_figures import MAP, Chart, Row
```

## The coil pair

A coil's field is the same in every plane through its axis, so the half plane of the radius $r \ge 0$ and the height
$z$ is enough. The windings are two squares of side 1 cm at $r = R$, the spacing $s$ apart. The air around them reaches
80 radii out and is meshed finely within two radii of the middle:

```python
mu0 = 4e-7 * math.pi
R, c, turns, I = 0.15, 0.01, 100, 1.0
J = turns * I / c**2


def pair(s):
    """the half plane r >= 0 around two windings at z = ±s/2, with the air near them and far away"""
    half = px.QuickSketch((0, -80 * R)).Rectangle(80 * R, 160 * R).Face()
    outer, inner = (px.Circle((0, 0), size).Face() * half for size in (80 * R, 2 * R))
    windings = [px.QuickSketch().RectangleC((R, z), c, c).Face() for z in (-s / 2, s / 2)]
    for winding in windings:
        winding.name, winding.maxh = "winding", c / 4
    near = inner - windings[0] - windings[1]
    near.name, near.maxh = "air", R / 40
    rest = outer - inner
    rest.name = "air"
    shape = px.Glue([rest, near] + windings)
    for edge in shape.edges:
        (r0, z0), (r1, z1) = edge.bounding_box
        if r1 > 1e-9 and max(r1, -z0, z1) > 40 * R:
            edge.name = "far"
    return shape


print("faces", [face.name for face in pair(R).faces])
# px.Plot(pair(R))
```

```python nbsphinx="hide-input"
Row((pair(R), "The half plane, 80 radii out"),
    (pair(R) * px.QuickSketch((0, -0.32)).Rectangle(0.33, 0.64).Face(), "The windings and the air near them"),
    caption="The coil pair in the half plane through its axis")
```

## The field

The unknown is $u = A_\varphi / r$, the vector potential around the axis divided by the distance from it. Unlike
$A_\varphi$ itself, it is smooth on the axis, and the flux density follows from it as

$$B_r = -r \frac{\partial u}{\partial z}, \qquad B_z = 2u + r \frac{\partial u}{\partial r}.$$

In weak form, $\int r^3 \, \nabla u \cdot \nabla v \, dr \, dz = \int \mu_0 J r^2 v \, dr \, dz$ for every test
function $v$, with the current density $J$ in the windings. The axis needs no boundary condition, and $u$ = 0 far away:

```python
def field(s, order=5):
    """the mesh, the flux density (B_r, B_z) and u of the pair with the spacing s"""
    mesh = ngs.Mesh(px.OCCGeometry(pair(s)).GenerateMesh(maxh=80 * R / 6))
    mesh.Curve(order)
    V = ngs.H1(mesh, order=order, dirichlet="far")
    trial, test = V.TnT()
    r = ngs.x
    a = ngs.BilinearForm(r**3 * ngs.grad(trial) * ngs.grad(test) * ngs.dx).Assemble()
    f = ngs.LinearForm(mu0 * J * r**2 * test * ngs.dx("winding")).Assemble()
    u = ngs.GridFunction(V)
    u.vec.data = a.mat.Inverse(V.FreeDofs(), inverse="sparsecholesky") * f.vec
    return mesh, ngs.CF((-r * ngs.grad(u)[1], 2 * u + r * ngs.grad(u)[0])), u


mesh, B, u = field(R)
print(mesh.ne, "triangles | field at the center", B(mesh(0, 0))[1], "T")
```

The formula for $B_0$ holds for thin loops. On the axis, one loop of radius $a$ at the height $h$ gives the field
$\mu_0 I a^2 / 2 \left(a^2 + (z - h)^2\right)^{3/2}$, and the field of a winding with a section is this loop field
summed over the section with SciPy's `dblquad`. The computed field at the center lies within a few millionths of it, and
the windings' section moves it from the thin-loop formula by less than a ten-thousandth:

```python
def on_axis(z, s):
    """the field on the axis at the height z, from the loop formula summed over the sections of both windings"""
    def loop(a, h):
        return mu0 * J * a**2 / (2 * (a**2 + (z - h) ** 2) ** 1.5)
    return sum(dblquad(loop, h - c / 2, h + c / 2, R - c / 2, R + c / 2, epsabs=1e-16, epsrel=1e-13)[0]
               for h in (-s / 2, s / 2))


computed, windings, thin = B(mesh(0, 0))[1], on_axis(0, R), (4 / 5) ** 1.5 * mu0 * turns * I / R
print("center: computed", computed, "| windings", windings, "| thin loops", thin)
print("computed / windings - 1", computed / windings - 1, "| windings / thin loops - 1", windings / thin - 1)
```

## The spacing

Along the axis, a spacing of 0.8 radii gives a field that is strongest in the middle, 1.2 radii one that is weakest
there, and one radius a field that stays flat the longest. Its drop at a tenth of the radius from the middle is close
to the fourth-order term of thin loops, $144/125 \times 10^{-4}$ (Figure 2):

```python
heights = np.linspace(0, 0.5 * R, 11)
profiles, solved = {}, {}
for spacing in (0.8, 1.0, 1.2):
    mesh, B, u = field(spacing * R)
    computed = np.array([B(mesh(0, z))[1] for z in heights])
    exact = np.array([on_axis(z, spacing * R) for z in heights])
    profiles[spacing] = (computed / computed[0] - 1, exact / exact[0] - 1)
    solved[spacing] = (mesh, B, computed[0])
    print("spacing", spacing, "R | at 0.1 R:", profiles[spacing][0][2], "| exact", profiles[spacing][1][2],
          "| at 0.3 R:", profiles[spacing][0][6], "| exact", profiles[spacing][1][6])
print("fourth-order term at 0.1 R", -144 / 125 * 0.1**4)
# plt.plot(heights / R, 100 * profiles[1.0][0]); plt.show()
# px.PlotField(ngs.Norm(B), mesh, pair(1.2 * R))
```

```python nbsphinx="hide-input"
fine = np.linspace(0, 0.5 * R, 101)
Chart("Height above the middle, in radii", "Field against the middle in %",
      lines={f"Spacing {s} R": (fine / R, [100 * (on_axis(z, s * R) / on_axis(0, s * R) - 1) for z in fine])
             for s in (0.8, 1.0, 1.2)},
      dots={"NGSolve": (np.tile(heights / R, 3), 100 * np.concatenate([p[0] for p in profiles.values()]))},
      title="The field along the axis for three spacings, the lines from the loop formula")
```

Away from the axis, the field is uniform in a region around the middle that is largest at the Helmholtz spacing.
Mirrored about the axis, the plane through the coils shows where the field lies within 2 % of its value in the middle
(Figure 3).

```python nbsphinx="hide-input"
def section(s):
    """the plane through the axis around the middle, mirrored, without the windings"""
    plane = px.QuickSketch().RectangleC((0, 0), 2.6 * R, 1.8 * R).Face()
    for r in (-R, R):
        for z in (-s / 2, s / 2):
            plane = plane - px.QuickSketch().RectangleC((r, z), c, c).Face()
    return plane


def sampled(mesh, B, b0):
    """|B| against the field in the middle, at points of the mirrored section"""
    def strength(face, plot=True, **_):
        (r0, z0), (r1, z1) = face.bounding_box
        gr, gz = np.meshgrid(np.linspace(r0, r1, 521), np.linspace(z0, z1, 361))
        rz = np.c_[gr.ravel(), gz.ravel()]
        rz = rz[face.Distance(rz, signed=True) < -1e-4]
        values = np.asarray(ngs.Norm(B)(mesh(np.abs(rz[:, 0]), rz[:, 1]))).ravel() / b0
        return {"points": [tuple(p) for p in rz], "values": values, "label": "|B| / B0"}
    return strength


Row(*[(section(s * R), f"Spacing {s} R", {"field": sampled(*solved[s]), "clim": (0.98, 1.02)}) for s in solved],
    caption="The field strength against the middle's, between 0.98 and 1.02")
```

See [Identifications](../integrations/netgen/identifications.md) for rotational symmetry with periodic edges, and
[The Halbach array](halbach-array.md) for a field made uniform on one side by the magnets' layout.

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
from docs_figures import INK
from gallery import save_figure

mesh, B, u = field(R)
fig, ax = plt.subplots(figsize=(6, 4))
px.PlotField(sampled(mesh, B, B(mesh(0, 0))[1]), section(R), ax=ax, colormap=MAP, min=0.98, max=1.02)
# field lines: the level lines of the flux r^2 u, at the levels that cross the middle plane at 0.1, 0.25, ... radii
gr, gz = np.meshgrid(np.linspace(-1.3 * R, 1.3 * R, 261), np.linspace(-0.9 * R, 0.9 * R, 181))
flux = gr**2 * np.asarray(u(mesh(np.abs(gr.ravel()), gz.ravel()))).reshape(gr.shape)
levels = sorted(r**2 * u(mesh(r, 0)) for r in R * np.arange(0.1, 0.9, 0.15))
ax.contour(gr, gz, flux, levels=levels, colors=INK, linewidths=0.7, zorder=4)
for r in (-R, R):
    for z in (-R / 2, R / 2):
        px.QuickSketch().RectangleC((r, z), c, c).Face().Plot(ax=ax, face_col=INK, line_col=INK, legend=False)
ax.set_aspect("equal", adjustable="box")
ax.set_xlim(-1.35 * R, 1.35 * R)
ax.set_ylim(-0.9 * R, 0.9 * R)
save_figure(fig, "helmholtz-coil", bare=True, bars={1: "|B| / B0"})
plt.close(fig)
```