Gilbert’s magnetic Earth (1600)

William Gilbert carved a sphere of lodestone, his terrella, and carried a small dip needle over it. The needle lay flat at the terrella’s equator, stood vertical at its poles and dipped steadily in between, just as compass needles did at sea. He concluded that the Earth itself is a great magnet. For an axial dipole the dip obeys tan(I) = 2 tan(latitude) and the field at the poles is twice that at the equator.

import matplotlib.pyplot as plt
import numpy as np

from earthsciencekit.geophysics import dipole_field, dipole_inclination, field_elements
from earthsciencekit.units import convert

latitude = np.linspace(-np.pi / 2, np.pi / 2, 181)
inclination = dipole_inclination(latitude)
field = dipole_field(latitude, 1.0, 1.0)
intensity = field_elements(field[..., 0], 0, field[..., 1])[..., 2]
assert np.isclose(intensity[-1] / intensity[90], 2)
np.testing.assert_allclose(np.tan(inclination[1:-1]), 2 * np.tan(latitude[1:-1]))

fig, (left, right) = plt.subplots(1, 2, figsize=(10, 4.5))
circle = np.linspace(0, 2 * np.pi, 361)
left.fill(np.cos(circle), np.sin(circle), color="0.85")
for lat in np.radians(np.arange(-80, 81, 20)):
    north, down = dipole_field(lat, 1.0, 1.0)
    direction = north * np.array([-np.sin(lat), np.cos(lat)]) + down * np.array([-np.cos(lat), -np.sin(lat)])
    direction /= np.hypot(*direction)
    for side in (1, -1):
        centre = 1.15 * np.array([side * np.cos(lat), np.sin(lat)])
        needle = 0.1 * direction * np.array([side, 1])
        left.annotate("", centre + needle, centre - needle, arrowprops={"arrowstyle": "->", "color": "C3"})
left.set(xlim=(-1.4, 1.4), ylim=(-1.4, 1.4), aspect="equal", title="Dip needles around a terrella")
left.axis("off")
right.plot(convert(latitude, "rad", "degree"), convert(inclination, "rad", "degree"), label="inclination")
right.plot(convert(latitude, "rad", "degree"), 45 * intensity / intensity[90], "--", label="intensity (x 45 degrees)")
right.set(xlabel="magnetic latitude [degrees]", ylabel="degrees", title="tan I = 2 tan(latitude)")
right.legend()
fig.tight_layout()

plt.show()