---
jupyter:
  description: 'Gabriel''s horn with planex2d and OpenCASCADE: y = 1/x as one exact rational B-spline, the region under
    it revolved; the volume tends to pi while the surface grows like 2 pi ln A without bound, both checked against
    their integrals.'
---

# Gabriel's horn

Can a solid with an infinite surface have a finite volume? The hyperbola $y = 1/x$, for $x$ from 1 on, turned about the
$x$ axis, bounds such a solid, Gabriel's horn. Torricelli showed that its volume is $\pi$, though it is infinitely long, and
its surface is infinite. The horn could be filled with $\pi$ of paint, but no amount of paint would cover its surface,
which is the painter's paradox.

The example cuts the horn off at $x = A$, for $A$ from 10 to 1000. The hyperbola is a conic, so one rational quadratic
B-spline traces it exactly. The region under it must have the area $\ln A$. Revolved, its volume must be
$\pi(1 - 1/A)$, its surface the integral $2\pi \int \tfrac{1}{x} \sqrt{1 + 1/x^4}\,dx$, and the surface must grow
without bound, like $2\pi \ln A$.

![Gabriel's horn, cut off at x = 10](../_static/gallery/gabriels-horn.png){.result}

```python
import math
import netgen.occ as occ
from scipy.integrate import quad
from scipy.optimize import brentq
import planex as px
```

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

## The hyperbola

A rational quadratic B-spline from $P_0$ to $P_2$, its middle pole $P_1$ where the tangents at the two ends cross,
traces a conic, and the middle weight $w$ says which one. On the line from the chord's middle $M$ to $P_1$, the curve
passes through the point $M + \tfrac{w}{1 + w}\,(P_1 - M)$, so the weight follows from where that line meets the
hyperbola, $s = w/(1 + w)$:

```python
def hyperbola(x0, x1):
    """y = 1/x from x0 to x1 as one rational quadratic B-spline"""
    p0, p2 = (x0, 1 / x0), (x1, 1 / x1)
    p1 = (2 * x0 * x1 / (x0 + x1), 2 / (x0 + x1))  # where the tangents at both ends cross
    m = ((p0[0] + p2[0]) / 2, (p0[1] + p2[1]) / 2)
    s = brentq(lambda s: (m[0] + s * (p1[0] - m[0])) * (m[1] + s * (p1[1] - m[1])) - 1, 0, 1)
    return px.BSpline([p0, p1, p2], [0, 0, 0, 1, 1, 1], 2, weights=[1, s / (1 - s), 1])


def region(A):
    """the region under the hyperbola from 1 to A"""
    return px.Shape([px.Line((1, 0), (A, 0)), px.Line((A, 0), (A, 1 / A)), hyperbola(1, A), px.Line((1, 1), (1, 0))]).Face()


for A in (10, 100, 1000):
    print("A", A, "| area under the hyperbola", region(A).mass, "ln A", math.log(A))
```

The region from 1 to 10 is turned about the $x$ axis into the horn (Figure 1):

```python
horn = region(10).ToOCCFace().Revolve(occ.Axis((0, 0, 0), occ.X), 360)
```

```python nbsphinx="hide-input"
Row((region(10), "The region under y = 1/x", {}),
    (horn, "Turned about the x axis", {"views": [(horn, "shade"), (horn, "ink")], "angles": (-30, 0, 20)}),
    caption="The region under the hyperbola from 1 to 10 and the horn")
```

## Volume and surface

The horn's volume is $\pi \int (1/x)^2\,dx = \pi(1 - 1/A)$. Its surface, without the two end discs, is the face whose area is
largest. Both are compared with their integrals, and the surface also with $2\pi \ln A$, its lower bound:

```python
for A in (10, 100, 1000):
    solid = region(A).ToOCCFace().Revolve(occ.Axis((0, 0, 0), occ.X), 360)
    side = max(f.mass for f in solid.faces)
    exact = 2 * math.pi * quad(lambda x: math.sqrt(1 + x**-4) / x, 1, A, limit=200)[0]
    print("A", A, "| volume", solid.mass, "exact", math.pi * (1 - 1 / A))
    print("     surface", side, "exact", exact, "| less 2 pi ln A", side - 2 * math.pi * math.log(A))
```

The volumes agree with $\pi(1 - 1/A)$ to about six digits, the surfaces with their integral to all digits. The volume
settles toward $\pi$, while the surface keeps the same lead over $2\pi \ln A$, about 0.71, and grows with $\ln A$
without bound. Each factor of 10 in length adds $2\pi \ln 10$, about 14.5, to the surface and almost nothing to the
volume.

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
from docs_figures import ViewFigure
from gallery import save_figure
fig, aspect = ViewFigure([(horn, "shade"), (horn, "ink")], (-30, 0, 20))
save_figure(fig, "gabriels-horn", dpi=160)
```