---
jupyter:
  description: 'The Halbach array with planex2d and NGSolve: one wavelength of a linear array of magnet blocks as a
    periodic cell, its magnetostatic field, and the field on the strong and the weak side against
    Br (1 - e^(-kd)) sin(pi/M)/(pi/M) e^(-ky) for four and eight blocks per wavelength.'
---

# The Halbach array

Can a row of magnets have its field on one side only? In a Halbach array, the magnetization turns from block to block,
by a quarter turn in the common version. On one side of the array the fields of the blocks add up, and on the other
side they cancel. An ideal array, whose magnetization turns continuously along it, has no field at all on its weak
side. On its strong side, the field has the constant magnitude $B = B_r\,(1 - e^{-kd})\,e^{-ky}$, for the remanence
$B_r$, the thickness $d$, the wavenumber $k = 2\pi/\lambda$ of the wavelength $\lambda$ and the distance $y$ from the
array. With $M$ blocks per wavelength, the field is lower by the factor $\sin(\pi/M)/(\pi/M)$, and a harmonic of order
$M - 1$ reaches the weak side. Mallinson described this one-sided flux, and Halbach used such arrays for the magnets
of particle accelerators.

The example solves one wavelength of an endless array as a periodic cell, with four and with eight blocks per
wavelength. The field above and below the array must agree with these formulas.

![The field lines and the field strength of a Halbach array with four blocks per wavelength, strong above and weak below](../_static/gallery/halbach-array.png){.result}

```python
import math

import matplotlib.pyplot as plt
import ngsolve as ngs
import planex as px
```

```python nbsphinx="hidden"
from docs_figures import Row
```

## The cell

The array has a wavelength of $\lambda = 40\,\text{mm}$ and is $d = 10\,\text{mm}$ thick, a quarter of the wavelength, as with square blocks
and four blocks per wavelength. Block $j$ is magnetized at the angle $360^\circ j/M$. The cell glues the $M$ blocks to air
$1.5\,\lambda$ deep above and below the array, where the field has faded to almost nothing. Its left and right sides are
named `left` and `right` and identified, so that the mesh repeats from cell to cell, and the far edges above and
below are named `far` (Figure 1):

```python
lam, d, Br = 40.0, 10.0, 1.3  # the wavelength and the thickness in mm, the remanence in T
k = 2 * math.pi / lam
depth = 1.5 * lam


def cell(M):
    """one wavelength of the array with M blocks, the air above and below, its sides identified"""
    width = lam / M
    blocks = [px.QuickSketch((j * width, 0)).Rectangle(width, d).Face() for j in range(M)]
    for j, block in enumerate(blocks):
        block.name = f"block{j}"
    air = [px.QuickSketch((0, d)).Rectangle(lam, depth).Face(), px.QuickSketch((0, -depth)).Rectangle(lam, depth).Face()]
    for part in air:
        part.name = "air"
    shape = px.Glue(blocks + air)
    for edge in shape.edges:
        (x0, y0), (x1, y1) = edge.bounding_box
        if x1 < 1e-9:
            edge.name = "left"
        elif x0 > lam - 1e-9:
            edge.name = "right"
        elif y1 < -depth + 1e-9 or y0 > d + depth - 1e-9:
            edge.name = "far"
    shape.Identify("left", "right")
    return shape


array = cell(4)
print("faces", [face.name for face in array.faces])
```

```python nbsphinx="hide-input"
arrows = ["→", "↑", "←", "↓"]
view = px.QuickSketch((0, -15)).Rectangle(lam, d + 30).Face()
Row((px.Glue([face * view for face in array.faces]), "Four blocks per wavelength",
     {"labels": [(arrows[j], (lam / 8 + j * lam / 4, d / 2)) for j in range(4)]}),
    caption="One wavelength of the array, with the direction of magnetization of each block")
```

## The field

With the field $B = \operatorname{curl} A$ of the vector potential $A$ along $z$, and blocks with the permeability of
air, the magnetostatic equations become
$\int \nabla A \cdot \nabla v = \int (B_{r,x}\,\partial v/\partial y - B_{r,y}\,\partial v/\partial x)$ for every test
function $v$, with the remanence $(B_{r,x}, B_{r,y})$ of each block. NGSolve's `Periodic` space joins the identified
sides, $A$ vanishes on the far edges, and $B = (\partial A/\partial y, -\partial A/\partial x)$:

```python
def field(M, order=4):
    """the mesh, the vector potential A and the flux density B of the array with M blocks per wavelength"""
    mesh = ngs.Mesh(px.OCCGeometry(cell(M)).GenerateMesh(maxh=lam / 16))
    mesh.Curve(order)
    fes = ngs.Periodic(ngs.H1(mesh, order=order, dirichlet="far"))
    u, v = fes.TnT()
    turn = {f"block{j}": 2 * math.pi * j / M for j in range(M)}
    bx = mesh.MaterialCF({name: Br * math.cos(a) for name, a in turn.items()}, default=0)
    by = mesh.MaterialCF({name: Br * math.sin(a) for name, a in turn.items()}, default=0)
    a = ngs.BilinearForm(ngs.grad(u) * ngs.grad(v) * ngs.dx).Assemble()
    f = ngs.LinearForm((bx * ngs.grad(v)[1] - by * ngs.grad(v)[0]) * ngs.dx).Assemble()
    A = ngs.GridFunction(fes)
    A.vec.data = a.mat.Inverse(fes.FreeDofs(), inverse="sparsecholesky") * f.vec
    return mesh, A, ngs.CF((ngs.grad(A)[1], -ngs.grad(A)[0]))


mesh, A, B = field(4)
```

In 2D, the field lines are the lines of constant $A$, and drawn at equal steps of $A$, they lie the closer together the
stronger the field is. They close above the array, and hardly any reach below it:

```python nbsphinx="hide-input"
import numpy as np


def field_lines(ax, mesh, A, levels=28):
    """the lines of constant A over the view, at equal steps"""
    X, Y = np.meshgrid(np.linspace(0, lam, 161), np.linspace(-15, d + 15, 201))
    values = A(mesh(X.flatten(), Y.flatten())).reshape(X.shape)
    ax.contour(X, Y, values, levels=levels, colors="#22303c", linewidths=0.7, linestyles="solid")


fig, ax = plt.subplots(figsize=(6, 6))
px.PlotField(B, mesh, px.Unify(array.faces) * view, ax=ax, min=0, max=1, colormap="default_r")
field_lines(ax, mesh, A)
ax.set_title("|B| in T and the field lines, four blocks per wavelength")
plt.show()
```

## Against the formulas

On the strong side, the harmonic of order $M + 1$ adds a ripple that fades as $e^{-(M+1)ky}$, so a little away from the
array, $|B|$ is the same all along it. On the weak side, only the harmonic of order $n = M - 1$ is left, with the
magnitude $B_r\,(1 - e^{-nkd})\,\lvert\sin(n\pi/M)/(n\pi/M)\rvert\,e^{-nky}$. The smallest and the largest $|B|$ along
the array at the distances $\lambda/8$, $\lambda/4$ and $\lambda/2$ above and below it are compared with the formulas:

```python
for M in (4, 8):
    mesh, A, B = field(M)
    strong = Br * (1 - math.exp(-k * d)) * math.sin(math.pi / M) / (math.pi / M)
    n = M - 1
    weak = Br * (1 - math.exp(-n * k * d)) * abs(math.sin(n * math.pi / M) / (n * math.pi / M))
    for y in (lam / 8, lam / 4, lam / 2):
        above = [ngs.Norm(B)(mesh(lam * i / 40, d + y)) for i in range(40)]
        below = [ngs.Norm(B)(mesh(lam * i / 40, -y)) for i in range(40)]
        print("M", M, "| y", y, "| above", round(min(above), 5), "to", round(max(above), 5), "formula",
              round(strong * math.exp(-k * y), 5), "| below", round(min(below), 6), "to", round(max(below), 6),
              "formula", round(weak * math.exp(-n * k * y), 6))
```

The table shows two things. First, the weak side is weak indeed. A quarter wavelength from the array with four
blocks, the field is 0.0035 T below it and 0.193 T above it, and with eight blocks the field below falls to a few
microtesla. Second, more blocks gain little on the strong side. Eight blocks per wavelength give 0.209 T instead of
0.193 T, 8 % more, since the factor $\sin(\pi/M)/(\pi/M)$ only rises from 0.900 to 0.974. Most arrays therefore use four
blocks per wavelength.

See [Periodic identifications](../integrations/netgen/identifications.md) for `Identify` and periodic meshes.

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
from gallery import save_figure
from matplotlib.patches import FancyArrowPatch
from docs_figures import INK
mesh, A, B = field(4)
fig, ax = plt.subplots(figsize=(6, 6))
px.PlotField(B, mesh, px.Unify(array.faces) * view, ax=ax, min=0, max=1, colormap="default_r")
field_lines(ax, mesh, A)
for j in range(4):  # block j is magnetized at the angle 360° j / M
    cx, cy, turn = lam / 8 + j * lam / 4, d / 2, math.radians(90 * j)
    dx, dy = 0.3 * lam / 4 * math.cos(turn), 0.3 * lam / 4 * math.sin(turn)
    ax.add_patch(FancyArrowPatch((cx - dx, cy - dy), (cx + dx, cy + dy), arrowstyle="-|>", mutation_scale=16,
                                 color=INK, lw=2, zorder=6))
save_figure(fig, "halbach-array", bare=True, panels=[0], bars={1: "|B| (T)"})
plt.close(fig)
```