Eclipses

Chapter 5

Besselian elements

So far we have studied the occurrence of eclipses by computing angles between the bodies as seen from the centre of the Earth. But nobody lives at the centre of the Earth. Next we have to work out how the rays of light meet the surface of the Earth.

Let us start by asking, at some instant, where a solar eclipse is at maximum. That is, where the axis of the Moon’s shadow meets the surface of the Earth. The tools we have at hand are enough to fix the shadow axis: it is the line between the centresActually it is not. In time we will learn to worry about the fact that the roughly spherical shape of the Moon does not sit quite exactly around its centre of mass. of the Sun and the Moon.

Let us look at how this would appear from space.

Figure 5.1:

The situation on 2027-08-02 at 10:05 UTC, when the shadow axis meets the surface of the Earth at 25.6° N, 32.7° E, in Egypt. Apart from the exaggerated sizes and distances of the shadow and Moon, everything else in the figure is from exact computation.

We can, in principle, compute the intersection of a given line and the surface of an ellipsoid in space. Once the direction of the shadow axis (the dashed line in the figure) is known, this would give us its intersection in geographic coordinates. By repeating a large pile of such calculations, we could plot many points as a function of time, and the central line of the eclipse on a map.

Going straight from the axis vector to geographic coordinates leads to fairly complicated expressions, though, and the method turns out more inconvenient than it needs to be, especially given the further calculations we will soon have the need for. So we should simplify the problem geometrically somehow.

5.1 The fundamental plane

This method was introduced in 1824 by the mathematician Friedrich Wilhelm Bessel (1784–1846), director of the Königsberg observatory. The shadow is not followed on the surface of a rotating, oblate Earth. It is instead taken to a plane that is fixed at the centre of the Earth and always turns to stay perpendicular to the shadow axis. On that plane the shadow is a circle, and every quantity describing it changes slowly and fairly evenly. Instead of metres or kilometres, we can use more convenient units of length.

Computing the details of an eclipse then goes roughly as follows. First a set of Besselian elements is drawn up to describe the course of the eclipse on this fundamental plane, the coordinates of the observer are transformed into the same coordinate frame, the equations and conditions for the various phases of the eclipse are written down, and the local points and curves solved from them are projected back into geographic coordinates.

Figure 5.2:

The Besselian fundamental plane is fixed at the centre of the Earth and always turns to stay perpendicular to the axis of the Moon’s shadow. Note that the coordinate lines drawn on the Earth here do not correspond to real parallels and meridians. The Earth may be in any orientation.

5.2 The elements of the shadow

The elements and give the position of the centre of the shadows, and and their radii. The angles of the shadow cones with respect to the shadow axis are and . The signs are chosen so that , and are always positive. is negative for a total eclipse and positive for an annular one. A negative sign means that the vertex of the umbral cone is behind the fundamental plane.

Figure 5.3:

The Besselian elements of the shadows. The centre of the Earth is at the origin and its radius is 1. The centres of the shadows are at . The radius of the penumbra is and that of the umbra . The angles of the shadow cones are and .

5.3 The elements of the Earth’s orientation

In order to express our results in geographic coordinates, we need two more elements: and , which describe the orientation of the Earth in the system. They express where the z-axis of the system points on the celestial sphere, or, putting it a little too simply, at which geographic coordinates the z-axis meets the surface of the Earth. That simplification would hold only for an Earth without flattening, that is, a spherical one.

The declination of the axis, , varies either side of zero with the seasons. It is nearly the same as the geocentric declination of the Sun, but not exactly, because the Sun does not lie quite on the z-axis of the system. Near the equinoxes .

The right ascension of the axis would describe its direction on the celestial sphere. Because the elements should be tied to the orientation of the Earth, we use the hour angle instead. The hour angle is the direction of the axis in the sky at Greenwich: when the prime meridian points along the z-axis.

Figure 5.4:

The Besselian elements and for the orientation of the Earth. Here they are drawn on an Earth assumed to be spherical. The Earth is in fact oblate, and these elements should be given as a direction on the celestial sphere, that is, on a star chart. The plane and Earth can be turned separately.

5.4 Computing the elements

We start with a couple of definitions, and derive the expressions for our elements from them.

Define a geocentric equatorial frame with its origin at the centre of the Earth and the unit vector pointing at the north pole. The vectors and lie in the plane of the equator, with pointing at the vernal equinox. Only is needed for the projection that follows.

Figure 5.5:

The geocentric equatorial frame. points at the north pole. and lie in the plane of the equator and turn with the direction of the vernal equinox, but their orientation does not matter to us. The plane drawn with a dashed line is the Besselian fundamental plane, into which the frame is transformed next.

One thing to keep in mind about the equatorial frame of figure 5.5: it does not follow the rotation of the Earth but is fixedPrecession moves this frame too with respect to the stars, but the timescale of that is counted in tens of thousands of years. Even over shorter intervals, however, it quickly becomes important to compute the ephemeris at the epoch of the observation rather than, say, that of the year 2000. to the stars. Because Besselian elements exist for computing what happens on the surface of the Earth, a little later we will have to use hour angles instead of right ascension.

Let the positions of the Moon and the Sun in the geocentric frame be and . The vector from the Moon to the Sun is , so the z-axis of the Besselian frame, which lies along the shadow axis, is given by the unit vector

The unit vectors of the fundamental plane are now, to the east and to the northThe job of the north pole vector is mainly to fix the directions of the Besselian fundamental plane so that “up”, the direction of the y-axis, is defined towards the north. If , that is, if the axis of the Moon’s shadow pointed straight at the north pole along the axis of the Earth, their cross product would be zero and there would be no way to know which way “up” has to be. This cannot happen on Earth, since the axis tilts only and nowhere near , but computing Besselian elements for Uranus, say, would call for some alternative frame. Helpfully, Uranus does not seem to have moons orbiting around its pole either.

Figure 5.6:

Forming the basis , , of the fundamental plane from the known shadow axis and the north pole of the Earth . Note that the vectors and (not drawn) could turn freely around without changing the result , , .

The position of the Moon in Besselian coordinates now comes simply from projecting

Figure 5.7:

The position vectors of the Moon and the Sun, the shadow axis, and the unit vectors of the Besselian system.

Another way to describe the same thing is to define the matrix

which is nothing but the rotation matrix from one frame to the other:

Nothing more is needed to compute the position of the Moon, and with it the centre of the shadows; Skyfield gives these vectors as they are:

Code for the position of the Moon in Besselian coordinates at the moment of the 2027 eclipse
from numpy import cross
from numpy.linalg import norm
from skyfield.api import load
from skyfield.framelib import true_equator_and_equinox_of_date as of_date

R_EARTH = 6378.1366           # km, IAU equatorial radius

ts = load.timescale()
eph = load('de440s.bsp')
earth, sun, moon = eph['earth'], eph['sun'], eph['moon']

t = ts.utc(2027, 8, 2, 10, 5, 19)
here = earth.at(t)

def apparent(body):
    """Geocentric apparent position in Earth radii, equator and equinox of date."""
    return here.observe(body).apparent().frame_xyz(of_date).km / R_EARTH

S, M = apparent(sun), apparent(moon)

k = (S - M) / norm(S - M)          # the shadow axis, towards the Sun
i = cross([0, 0, 1], k)            # where the plane cuts the equator, east positive
i /= norm(i)
j = cross(k, i)                    # completes the triad, north positive

x, y, z = M @ i, M @ j, M @ k      # the Moon in the fundamental frame
# alternative formulation of the same thing:
# B = column_stack([i, j, k])        # the fundamental basis in equatorial coordinates
# x, y, z = B.T @ M                  # the Moon in the fundamental frame
print(f'x = {x:+.5f}   y = {y:+.5f}   z = {z:.3f}')
x = +0.03896   y = +0.13730   z = 56.032

Note that x and y are well below 1, so the umbra falls somewhere on the Earth. The z-coordinate of the Moon is not one of the Besselian elements, but it is needed in the calculations that follow.

Alternative way to compute x, y and z, from right ascension, declination and distance

If rectangular position vectors are not at hand, only geocentric right ascensions, declinations and distances, the elements follow from the same inputs by trigonometry. This is the form the expressions take in the literature, so they are worth recognising. The same calculations are in chapter 11.3.2.2 of the Explanatory Supplement, although the presentation there is no clearer than this one.

First the direction of the shadow axis, from the known directions and distances of the Moon and the Sun. The position vectors are

Here , and are the geocentric right ascension, declination and distance of the Sun, and , and the same quantities for the Moon. The angles have to be computed at the epoch of the observation, and they have to be apparent, that is, the speed of light has to be accounted for. Distances are given in Earth radii.

Writing the shadow axis as , its right ascension and declination are

Once and are known, the position of the Moon on the fundamental plane follows directly from its own right ascension and declination:

where

The latter form is how the distance of the Moon is given in older sources: the tabulated quantity is the horizontal parallax , and the distance in Earth radii is the reciprocal of its sine.

The expressions are the same three projections as above. They follow from substituting the components of the basis vectors

into the dot products , and , and combining the terms in and using

So both angles appear only in their difference, and the choice of the zero point of right ascension does not affect the result. The result can be checked against , which holds because the projection only rotates the frame.

Next the elements for the orientation of the Earth. Having already computed the shadow axis in the equatorial rectangular frame of the Earth, the declination on the celestial sphere of the direction it points to is

Figure 5.8:

With pointing at the celestial pole in the equatorial frame, the declination of the shadow axis is .

After the declination, the sideways direction: the right ascension and the hour angle of the shadow axis. The right ascension comes from projecting the shadow axis onto the equator and measuring the angle to the direction of the vernal equinox, the vector .

Once a dayOnce a sidereal day, that is. the Greenwich meridian points at the vernal equinox, and Greenwich apparent sidereal time (GAST) is zero degrees, or zero hours, as sidereal time is usually given. At other times GAST is measured as an angle counterclockwise from the x-axis.

Figure 5.9:

The right ascension of the shadow axis on the celestial sphere is the angle . To find the hour angle (not in the figure) we still need the apparent sidereal time at Greenwich (G), GAST.

Finally, the Greenwich hour angle of the shadow axis (also GHA) is

Figure 5.10:

The equatorial plane of the previous figure, seen from the north. Given the sidereal time at Greenwich, GAST, and the right ascension of the shadow axis, the hour angle of the axis is their difference.

Code for the elements d and μ at the same instant
# continues from the previous snippet: t and k are the ones computed there
from numpy import arcsin, arctan2, degrees

d = degrees(arcsin(k[2]))               # the axis's declination
a = degrees(arctan2(k[1], k[0]))        # and its right ascension
gast = t.gast * 15                      # sidereal time, hours to degrees
mu = (gast - a + 180) % 360 - 180       # the hour angle at Greenwich, -180..180

print(f'd  = {d:+.5f} deg')
print(f'mu = {mu:+.5f} deg')
d  = +17.76138 deg
mu = -30.24843 deg

Now for the shadow cones.

Every quantity the shadow calculation needs from the ephemeris has already been computed above: the distance between the Sun and the Moon , and the perpendicular distance of the Moon from the fundamental plane .

Two constants are still needed: the radii of the Sun and the Moon as multiples of the equatorial radius of the Earth. The radius of the Sun and that of the Earth give . For the Moon the adopted radius is , or . We will come back to this number later, when it turns out that the choice needs care in computing the course of totality.

The vertex of the penumbral cone is between the Sun and the Moon. The vertex of the umbral cone is behind the Moon.

Figure 5.11:

The angles of the shadow cones.

For the vertex of either cone

For the penumbra , and for the umbra . Writing these out gives the sines of the angles for both shadows

The same equation gives the distance of the vertex from the centre of the Moon, . The penumbral vertex lies from the Moon towards the Sun and the umbral one in the opposite direction.

The only elements still missing are the radii of these cones on the fundamental plane. First the distances of the vertices from the plane.

The distance of the umbral vertex is negative in a total eclipse, when the vertex is behind the plane.

Figure 5.12:

The sizes of the shadows on the fundamental plane. They are computed from the vertices and and the angles and . Moving the fundamental plane shows how turns positive for an annular eclipse.

The radius of a cone grows away from the vertex in proportion to , so the radii of the shadows on the fundamental plane are

We now know the position and size of both shadows on the fundamental plane. An observer, however, is generally not on that planeOnly near sunrise and sunset., but some way above it. For this reason the Besselian elements also carry the angles of the cones.

Let the observer be at a height (zeta) above the plane. In the parallel plane through the observing site, the local radius of either shadow is then

The radii on the fundamental plane are published as elements in their own right. Once they are known, the local radius is the only calculation the angle is needed for. That is why the elements give the angles directly as the tangents and rather than as angles themselves:

Now the numerical values of the remaining elements, at the instant of the earlier examples.

Code for the shadow radii and the cone angles at the same instant
# continues from the previous snippets: S, M, z and R_EARTH carry over
from numpy import sqrt

k_s = 696000 / R_EARTH             # the Sun's radius, in Earth radii
k_m = 0.2725076                    # the Moon's, adopted

g = norm(S - M)                    # the distance between the centres
sin_f1 = (k_s + k_m) / g
sin_f2 = (k_s - k_m) / g
tan_f1 = sin_f1 / sqrt(1 - sin_f1**2)
tan_f2 = sin_f2 / sqrt(1 - sin_f2**2)

l1 = (z + k_m / sin_f1) * tan_f1   # c1 tan f1
l2 = (z - k_m / sin_f2) * tan_f2   # c2 tan f2

print(f'l1 = {l1:+.5f}   l2 = {l2:+.5f}')
print(f'tan f1 = {tan_f1:.7f}   tan f2 = {tan_f2:.7f}')
l1 = +0.53062   l2 = -0.01569
tan f1 = 0.0046065   tan f2 = 0.0045835

That gives all of the elements at 2027-08-02 10:05:19 UTC:

+0.03896
+0.13730
+17.76138°
−30.24843°
+0.53062
−0.01569
0.0046065
0.0045835

5.5 The elements as a function of time

We can now compute a full set of elements for any instant we choose. The same calculation can be repeated across the whole eclipse at regular intervals, at as fine a time step as we like. The table below, for the day of the 2027 eclipse, is computed at half-hour intervals, but we could just as well compute the elements every second.

time
08:00−1.0988+0.577817.783−61.5820.53054−0.015760.004610.00458
08:30−0.8264+0.472517.778−54.0810.53057−0.015740.004610.00458
09:00−0.5540+0.367017.772−46.5800.53059−0.015710.004610.00458
09:30−0.2817+0.261617.767−39.0790.53061−0.015700.004610.00458
10:00−0.0093+0.156017.762−31.5780.53062−0.015690.004610.00458
10:30+0.2630+0.050417.757−24.0770.53062−0.015680.004610.00458
11:00+0.5354−0.055317.752−16.5760.53062−0.015690.004610.00458
11:30+0.8076−0.161017.747−9.0750.53061−0.015700.004610.00458
12:00+1.0799−0.266817.742−1.5740.53060−0.015710.004610.00458
12:30+1.3521−0.372617.737+5.9270.53057−0.015730.004610.00458

This is enough to animate the course of the eclipse on the fundamental plane.

Figure 5.13:

The shadow crossing the fundamental plane and the Earth, to scale, on 2027-08-02. The rotation of the Earth during the eclipse is shown by projecting the geographic coordinate lines onto the fundamental plane. Note that if the shadows were really projected onto a map, they would grow a little from the sizes shown here, which are computed on the plane.

5.6 The elements as polynomials

The Besselian method above, with the elements on a plane, made the calculation much easier. Although a computer does it quickly, grinding through the ephemerides of the bodies for every second may still be wasted work. In the era of hand computation, and of slow computers, the usual way to speed it up was to compute the ephemeris at a few points only and fit a low-degree polynomial to the results.

The method gives accurate results and simplifies the calculation today as well. The animations on this page, for instance, are computed on the front end, in your own browser, where evaluating a few precomputed polynomials is far more convenient than chewing through the ephemerides of the bodies.

Let us look more closely at how the elements behave, one or two at a time. First and , which move thousands of kilometres during the eclipse, and evenly by the look of it. But how evenly? Fit a straight line to the computed points and see what residuals such a linear approximation leaves.

Figure 5.14:

On the left, the elements and at ten-minute intervals, and the straight line fitted through each set of points. On the right, the residuals of the elements against the line fitted through them.

The residual plot shows that the values depart from the linear estimate by several kilometres at worst. That is not accurate enough. Let us see what the residuals would be for a second and a third degree polynomial fit.

Figure 5.15:

The residuals of the second and third degree fits for the same points.

The second degree fit still leaves several hundred metres of residual, which is not accurate enough either when we want the events of a total solar eclipse to the second. At the third degree the difference is down to centimetres, which is enough for us. There is no sense in going higher than that; the ephemerides themselves will certainly contribute a larger error.

Next and . The declination changes very little during the eclipse, staying around . The hour angle changes much more, very evenly with the rotation of the Earth, at very close to per hour. Fitting a first and a second degree polynomial to both shows that at the first degree the residuals are of the order of metres measured on the surface of the Earth, and at the second degree centimetres.

Figure 5.16:

The first and second degree residuals for the elements and , in metres on the surface of the Earth.

Recall that the hour angle was computed as . Of the two, GAST describes the rotation of the Earth and is linear to very high accuracy. The residuals that remain come from the slow change in the direction of the shadow axis, in its right ascension in the case of .

Meeus (1989) uses the second degree for the declination and the first for the hour angle. For our own purposes the second degree does for both. A third degree term would be nothing but noise.

The shadow radii and behave in an interesting way. The two shadows are of very different sizes, but there is practically no difference in how they change. Over the eclipse both change by about a kilometre, growing first and then shrinking.

A first degree fit does not capture this, but the second degree fits to within centimetres.

Figure 5.17:

The first and second degree residuals for the shadow radii and , in metres. The points are practically identical, so is drawn as crosses.

In the earlier table the tangents of the angles and do not change at all at the precision shown there. Subtracting from each series its own mean brings the change into view. It has a linear or a second degree shape, but its amplitude is hundredths of an arcsecond. These angles can be treated as constants for any sensible purpose.

Figure 5.18:

The size of the change in the cone angles and , in arcseconds. is drawn as crosses, because the points are practically identical.

We now know which degree to use for each element in the polynomial fit. One more shortcut, as in the old days. Above we fitted the polynomials to a dense set of points at 10 minute intervals. Since our fits are of the third degree at most, the number of points used can be cut to five. The traditional way is to look at the moment of greatest eclipse, take the whole hour nearest to itMore usually the epoch chosen is that whole hour in terrestrial time (TT). This is in practice the same thing as UTC without leap seconds. By the time of writing the difference between them has grown to a little over a minute. More on this peculiarity of timekeeping in later chapters. The choice of instant makes no difference to the result. The chosen epoch does of course have to be given when the elements are published, so that they can be used. , and compute the fitting points one and two hours before and after that hour.

Figure 5.19:

The five whole hours chosen for the polynomial fit. The example is the position of the umbra on the fundamental plane.

The results can now be given as a table of polynomials. These are our own Besselian elements for the total solar eclipse of August 2027, valid for 2027-08-02 10:00 UTC .

elementconstant
−0.0000092−0.00004460.5447104−0.0093024
0.0000038−0.0001217−0.21116260.1560073
−0.0000039−0.010181217.7622869
0.000002015.0020965−31.5777830
−0.00001280.00001360.5306203
−0.00001280.0000135−0.0156866
0.0046065
0.0045835

Next we will see that a set of element numbers like this is useful for a wide range of calculations, and that the numbers work as they are in many existing tools.