---
jupyter:
  description: 'Sharp or blunt lightning rods with planex2d and NGSolve: the electric field around rods shaped as
    exact half ellipses, the tip field against E0 (1 + a/b), the field above the tips, where a blunt rod overtakes
    a sharp one, and a rod on the ridge of a house.'
---

# Sharp or blunt lightning rods

Why are lightning rods pointed? At the surface of a conductor, the electric field is strongest where the surface is
curved most sharply, and Benjamin Franklin's rods ended in sharp points. For a rod shaped as a half ellipse of height $a$ and
half width $b$ on flat ground, the field at the tip has a closed form. With its mirror image below the ground, the rod
is a conducting elliptic cylinder in a uniform field $E_0$, and in elliptic coordinates $\xi, \eta$ its potential is
$\varphi = E_0 c \cos\eta \,(\cosh\xi - \cosh\xi_0\, e^{\xi_0 - \xi})$, with $c = \sqrt{a^2 - b^2}$ and the rod's
surface at $\xi_0$. At the tip, this gives

$$E_{tip} = E_0 \left(1 + \frac{a}{b}\right),$$

so the sharper the rod, the stronger the field at its tip, without a limit. Field tests tell a different story. In
seven years on a mountain, Moore, Rison, Mathis and Aulich (2000) saw twelve blunt rods struck and none of the sharp
rods beside them, and their calculations showed why: a few centimeters above the tip, a blunt rod's field is
stronger than a sharper rod's of the same height.

Here, NGSolve solves the electric potential around four rods of the same height, from a half circle to a sharp half
ellipse. The tip fields are checked against $1 + a/b$, the fields above the tips are compared, and a rod goes onto the
ridge of a house.

![The electric field around a blunt and a sharp lightning rod of the same height](../_static/gallery/lightning-rod.png){.result}

```python
import math

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

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


def sampled(E, box, n=241):
    """the field strength E for Row and PlotField, at points of the air inside the view box ((h0, h1), (v0, v1)), n
    along h"""
    def field(face, plot=True, **_):
        (h0, h1), (v0, v1) = box
        gh, gv = np.meshgrid(np.linspace(h0, h1, n), np.linspace(v0, v1, int(n * (v1 - v0) / (h1 - h0))))
        uv = np.c_[gh.ravel(), gv.ravel()]
        uv = uv[face.Distance(uv, signed=True) < -0.01]
        return {"points": [tuple(p) for p in uv], "values": np.asarray(E(uv[:, 0], uv[:, 1])).ravel(),
                "label": "field in units of E0"}
    return field
```

## The rods

Every rod is 10 high, and the half widths run from 10, a half circle, down to 1.25. `px.Ellipse` draws each rod
exactly, standing upright with `angle=90`, and the ground cuts it in half. The air is a box 600 high and 1200 wide,
large against the rods. The rod's edges and the box's ground and top are named before the Boolean operation, which
keeps the names. The rod's edges get a mesh size of a quarter of the radius of curvature at the tip, $b^2/a$:

```python
a, L = 10.0, 600.0


def air(b):
    rod = px.Face([px.Ellipse((0, 0), a, b, angle=90)])
    rod.edges.name = "rod"
    rod.edges.maxh = b**2 / a / 4
    box = px.QuickSketch((-L, 0)).Rectangle(2 * L, L).Face()
    box.edges(near=(0, 0)).name = "ground"
    box.edges(near=(0, L)).name = "top"
    return (box - rod).Face()
```

## The field

The rod and the ground are at the potential 0, and the top of the box at $E_0 L$ with $E_0 = 1$, so the field far from
the rod is uniform. NGSolve solves Laplace's equation for the potential with elements of order 3 on a mesh curved to
the ellipse, and `field` returns the strength of the field at given points:

```python
def field(b, around=air):
    mesh = ngs.Mesh(px.OCCGeometry(around(b)).GenerateMesh(maxh=L / 10, grading=0.25))
    mesh.Curve(3)
    V = ngs.H1(mesh, order=3, dirichlet="rod|ground|top")
    u, v = V.TnT()
    phi = ngs.GridFunction(V)
    phi.Set(L, definedon=mesh.Boundaries("top"))
    K = ngs.BilinearForm(ngs.grad(u) * ngs.grad(v) * ngs.dx).Assemble()
    phi.vec.data += K.mat.Inverse(V.FreeDofs(), inverse="umfpack") * (-K.mat * phi.vec)
    strength = ngs.Norm(ngs.grad(phi))
    return lambda h, v: strength(mesh(h, v))


fields = {b: field(b) for b in (10, 5, 2.5, 1.25)}
for b, E in fields.items():
    print("b", b, "| tip field", E(0.0, a + 1e-6), "| exact", 1 + a / b)
```

Around the blunt and the sharp rod, the field is strongest at the tip (Figure 1). The color scale stops at 4, below
the sharp rod's 9:

```python nbsphinx="hide-input"
box = ((-20, 20), (0, 25))
Row(*[(air(b), f"Half width {b}", {"field": sampled(fields[b], box), "cmap": MAP.reversed(), "clim": (0, 4),
                                    "view": box}) for b in (10, 1.25)],
    caption="The field around a blunt and a sharp rod of the same height")
```

## Above the tip

Along the axis above the rod, the exact solution gives the field as well, with $y = c\cosh\xi$:
$E/E_0 = 1 + \cosh\xi_0\, e^{\xi_0 - \xi} / \sinh\xi$, and for the half circle $1 + a^2/y^2$. The sharp rod's field
is the highest at its tip and falls fastest. A third of the rod's height above the tip, the half circle's field has
overtaken it:

```python
def exact(y, b):
    if b == a:
        return 1 + a**2 / y**2
    c = math.sqrt(a**2 - b**2)
    xi0, xi = math.acosh(a / c), math.acosh(y / c)
    return 1 + math.cosh(xi0) * math.exp(xi0 - xi) / math.sinh(xi)


above = np.linspace(a + 0.01, 3 * a, 2000)
blunt, sharp = (fields[b](0 * above, above).ravel() for b in (10, 1.25))
crossing = above[np.argmax(blunt > sharp)] - a
exact_crossing = above[np.argmax([exact(y, 10) > exact(y, 1.25) for y in above])] - a
print("the half circle overtakes the sharp rod", round(crossing, 2), "above the tip, exact", round(exact_crossing, 2))
```

The computed fields of the two rods follow the exact solution, and the crossing is marked (Figure 2):

```python nbsphinx="hide-input"
few = above[::100]
Chart("height above the tip", "field in units of E0",
      dots={f"finite elements, half width {b}": (few - a, fields[b](0 * few, few).ravel()) for b in (10, 1.25)},
      lines={f"exact, half width {b}": (above - a, [exact(y, b) for y in above]) for b in (10, 1.25)},
      marks={"the crossing": (crossing, blunt[np.argmax(blunt > sharp)])},
      title="The field on the axis above a blunt and a sharp rod")
```

The sharp tip concentrates its field into a small region, and the blunt rod spreads its field farther into the air
above it. This is the effect that Moore and his colleagues found in their calculations for round rods. The model here
is planar, a ridge rather than a round rod, but the comparison comes out the same way.

## On a roof

On a building, the rod stands on the ridge. A house 10 wide, with walls 6 high and a roof pitched at 35°, gets a rod 2
high and 0.2 in half width on its ridge. The house and the rod are one conductor at the potential 0, named `rod`
before they are joined, and the same `field` solves the air around them:

```python
ridge = 6 + 5 * math.tan(math.radians(35))


def building(b, height=2.0):
    walls = px.QuickSketch((-5, 0)).Line((5, 0)).Line((5, 6)).Line((0, ridge)).Line((-5, 6)).Close().Face()
    walls.edges.name = "rod"
    rod = px.Face([px.Ellipse((0, ridge), height, b, angle=90)])
    rod.edges.name = "rod"
    rod.edges.maxh = b**2 / height / 4
    return (walls + rod).Face()


def house(b):
    box = px.QuickSketch((-L, 0)).Rectangle(2 * L, L).Face()
    box.edges(near=(0, 0)).name = "ground"
    box.edges(near=(0, L)).name = "top"
    return (box - building(b)).Face()


on_roof = field(0.2, house)
print("tip field on the ridge", on_roof(0.0, ridge + 2 + 1e-6), "| the same rod on flat ground", 1 + 2 / 0.2)
```

The house raises the field around its ridge, so the rod's tip sees about twice the field of the same rod on flat
ground (Figure 3). Away from the house, the field returns to $E_0$:

```python nbsphinx="hide-input"
view = ((-12, 12), (0, 16))
Row((house(0.2), "A rod on the ridge", {"field": sampled(on_roof, view), "cmap": MAP.reversed(), "clim": (0, 4),
                                        "view": view}),
    caption="The field around a house with a lightning rod on its ridge")
```

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
from gallery import save_figure
from docs_figures import FACE, INK
frame = px.QuickSketch((-12, 0)).Rectangle(24, 16).Face()
fig, ax = plt.subplots(figsize=(6, 4))
px.PlotField(sampled(on_roof, view, 481), house(0.2) * frame, ax=ax, colormap=MAP.reversed(), min=0, max=4)
px.Plot(building(0.2), ax=ax, face_col=FACE, line_col=INK, arc_col=INK, spline_col=INK)
ax.set_aspect("equal", adjustable="box")
ax.set_xlim(-12, 12)
ax.set_ylim(0, 16)
save_figure(fig, "lightning-rod", bare=True, bars={1: "field in units of E0"})
```