Besselian Elements are generally used to predict the circumstances of solar eclipses because using a full ephemeris, and brute forcing the circumstances, is more computationally involved. Besselian Elements are basically an ephemeris valid only over a short period of time (usually about five hours), and will be reduced to a few polynomials which can be solved quite rapidly for many events observers are usually interested in. In the past, it was most common to precompute and publish the elements for programs to use later, but today it is more common to compute the elements on the fly, though they are still published in publications like the Astronomical Almanac.
Below is an example set of elements for the April 8, 2024 total solar eclipse from the Five Millennium Canon of Solar Eclipses.
n x y d l1 l2 μ
0 -0.3182440 0.2197640 7.5862002 0.5358140 -0.0102720 89.591217
1 0.5117116 0.2709589 0.0148440 0.0000618 0.0000615 15.004080
2 0.0000326 -0.0000595 -0.0000020 -0.0000128 -0.0000127 0.000000
3 -0.0000084 -0.0000047 0.0000000 0.0000000 0.0000000 0.000000
tan f1 = 0.0046683 tan f2 = 0.0046450
"Eclipse Predictions by Fred Espenak, NASA's GSFC"
This article describes how to generate the table above. A different article will show how to take the values from the table above and compute things like the duration of totality for a specific location, or the path of the central line, etc.
This article does not directly follow any specific source, and instead pulls from several sources using equations and derivations which are convenient based on the information assumed to be at hand. The primary source followed is "A Manual of Spherical and Practical Astronomy" by William Chauvenet, but many of the equations used are not numbered in that work, so the code primarily cites the "Explanatory Supplement to the Astronomical Almanac 3rd Ed.", Robin Green's "Spherical Astronomy", and W. M. Smart's "Textbook on Spherical Astronomy". Most of the symbols are consistent across publications, but presented in a slightly different manner. Readers interested in the derivation or a detailed understanding of the equations used here are suggested to consult one of those sources.
Technically speaking, the values \(x, y, d, l_1, l_2, \mu, \tan f_1, \tan f_2 \) computed for a specific instant in time are the Besselian Elements. And the table above is a set of coefficients which will be used to compute those elements. But, more often than not, the table of coefficients is referred to as the Besselian Elements. This terminology is used nearly every place I have seen the coefficent table published (the next most common term is "Polynomial Besselian Elements"). That means there is some ambiguity as to what someone is referring to when they say Besselian Elements. And, for lack of better terminology, this article will not try to eliminate that problem. It should be quite obvious from the context as to whether someone is referring to a table of coefficents, or to the elements for a specific instant in time. And, for the sake of brevity, I will often refer to them as just "elements", or "coefficents".
The equations below, which compute the elements for a given instant in time from the position of the Sun, Earth, and Moon, are referred to as the "Fundamental Equations". Another set of equations will compute the elements using the table of coeffients, and these are usually not referred to by a name, which again produces some ambiguity when referring to "computing the Besselian Elements", but again, should be quite obvious from the context. This article will involve only using the fundamental equations, and use of the coefficents will be relegated to an article on computing the local circumstances.
It is assumed that the reader has an ephemeris that can produce the apparent positions of the Sun, Earth, and Moon to the accuracy desired by the reader. The positions should be geocentric apparent positions in Right Ascention, Declincation, and distance. "Apparent position", meaning the positions should include corrections for precession, nutation, light time, annual and dinural abberation, etc, to the precision desired.
It is also assumed that the reader knows an approximate time an eclipse will occur, this would have been achived through a brute force search by computing the positions of the Sun, Earth, and Moon and determining an eclipse ocurrs at that given time. This time will be referred to as \(T0\), and will be published along with the coefficents table (usually as a Julian Date). T0 does not have to be the exact maxium, or center time of the eclipse, it is usually taken as the nearest whole hour which is closes to the time of maximum. All times are measured in hours relative to T0 as \(t\). For example, one hour before T0, t=-1, and two hours after T0 is t=2, half an hour before T0, t=-.5. Some publications refer to this sligly differently as \(t0\) or \(t_0\), but it should be obvious from the context of the publication or explained in the text. A value of \(t\) can be converted to and from a Julian Date using the formula T = T0 + t/24.
Given a T0 time, we use the fundamental equations below to compute the elements for that given instant in time. We also compute the elements for a few more instants, usually 1 and 2 hours before and after T0, resulting in 5 points. Other publications may use different intervals, but these are commonly used. So, the read would compute the geocentric apparent position of the Sun, and Moon for the times: T0 - 2h, T0 - 1h, T0, T0 + 1h, T0 + 2h, and the fundamental equations below computed for each time.
This results in five values for each element \(x, y, d, l_1, l_2, \mu, \tan f_1, \tan f_2 \). We then interpolate these values to produce polynomials, the coefficients of which will become the table of coefficients that will be published as the Besselian Elements for an eclispe.
The values \( \alpha', \delta', r' \) represent the apparent geocentric right ascension, declination, and distance of the Sun. And the values \( \alpha, \delta, r \) represent the same quantaties, but for the Moon. The distances \( r, r' \) should be in units of Earth radii.
The rectangular coordinates of the Sun: $$ \begin{align*} x_s=&r' \cos \delta' \cos \alpha' \\ y_s=&r' \cos \delta' \sin \alpha' \\ z_s=&r' \sin \delta' \\ \end{align*} $$ The rectangular coordinates of the Moon: $$ \begin{align*} x_m=&r \cos \delta \cos \alpha \\ y_m=&r \cos \delta \sin \alpha \\ z_m=&r \sin \delta \\ \end{align*} $$ The vectors \(\vec S \) and \(\vec M\) represent the geocentric postion of the Sun and Moon and the intermediate vector \(\vec G\) represents a vecotr from the Sun to the Moon. And the scalar value \(G\) is the distance from the Sun to the Moon. $$ \begin{align*} \vec S =& [x_s, y_s ,z_s]\\ \vec M =& [x_m, y_m, z_m]\\ \vec G =& \vec S - \vec M \\ G =& |\vec G | \\ \end{align*} $$ The intermediate values \(a\) and \(d\) represent the right ascension and declincation of the Sun as viewed from the Moon. This is referred to as point Z. $$ \begin{align*} \tan a =& \frac{\vec G_y/G}{\vec G_x/G} \\ \sin d =& \vec G_z/G \\ \end{align*} $$ \(\mu\) is the Greenwich Hour Angle of point Z. $$ \begin{align*} \mu =& \theta - a \\ \end{align*} $$ The coordinates \(x, y\) are the intersection of the center of the shadow with the fundamental plane, and \(z\) along the Z axis is the distance to the Moon. $$ \begin{align*} x =& r \cos \delta \sin (\alpha - a) \\ y =& r \sin \delta \cos d - \cos \delta \sin d \cos(\alpha - a) \\ z =& r \sin \delta \sin d + \cos \delta \cos d \cos(\alpha - a) \\ \end{align*} \\ $$ The constants \(d_s\) and \(k\) are the Sun's radius, and the Moon's radius in units of Earth radii. They are computed from accepted constants defined by the IAU. The values below are from IAU Resolution B3 and IAU Resolution C 10. $$ \begin{align*} d_s =& \frac{6.957e8}{6.3781e6} \\ \\ k =& 0.2725076 \\ \end{align*} $$ The remaining elements \(\tan f_1, \tan f_2, l_1, \) and \( l_2 \) are computed as below. Note the intermediate \(\sin f_x\) values are computed separately only for computational effenciency. $$ \begin{align*} \sin f_1 =& \frac{d_s + k}{G} \\ \sin f_2 =& \frac{d_s - k}{G} \\ \\ c_1 =& z + \frac{k}{\sin f_1} \\ c_2 =& z - \frac{k}{\sin f_2} \\ \\ \tan f_1 =& \tan(\arcsin(\sin f_1)) \\ \tan f_2 =& \tan(\arcsin(\sin f_2)) \\ \\ l_1 =& c_1 \tan f_1 \\ l_2 =& c_2 \tan f_2 \\ \\ \end{align*} $$ The resulting set of Besselian Elements are \(x, y, d, l_1, l_2, \mu, \tan f_1, \tan f_2 \).The value \(x, y, d, l_1, l_2, \mu \) from above must be interpolated to produce the Besselian Element Polynomials, and the values \(\tan f_1, \tan f_2 \) should be averaged. The number of points, and degree of interpolation are up to the programmer, but generally they are computed at one hour intervals starting two hours before T0 eclipse and ending two hours after T0, resulting in five sets of data points. \(x,y\) are usually interpolated as third degree polynomials, \(d, l_1, l_2 \) are usually second degree polynomials, and \(\mu\) a first degree.
For simplicity, here we will interpolate all to the third order, and just disregard the unneeded values. This will simplify the code base by eliminating the number of cases at the cost of some performance. The interpolation is done independently for each element, but is the same for all of them. So I will just use \(x\) as an example, and the remainder of the elements can be interpolated in the same manner.
When we compute the circumstances of an eclipse, we compute the value of \(x\) using: $$ x = a_0 + a_1t + a_2t^2 + a_3t^3 $$ When we compute the elements for an instant, we choose the value of \(t\), so the goal is to interpolate for the cooeffients \(a_i\). There are many methods to do this, but we will use Gauss-Jordan Elimination. Since it is a common process detailed in many other places, I will present only a brief overview. The process is quite simple, we create an augmented matrix with all of our values of \(t\) raised to the powers 0 to 3 on the left, and the values of \(x\) on the right. We then perform a number of row opperations to change the left hand side to reduced row-echelon form, which means it has 1's along the diagonal, and zeroes everywhere else. The permitted row operations are: 1. Swap any two rows, 2. Multiply a row by a constant, 3. Add one row to another multiplied by some constant. The generalized matrix is
$$ \begin{bmatrix} t^0_{-2} & t^1_{-2} & t^2_{-2} & t^3_{-2} &\bigm| & x_{t-2} \\ t^0_{-1} & t^1_{-1} & t^2_{-1} & t^3_{-1} &\bigm| & x_{t-1} \\ t^0_{0} & t^1_{0} & t^2_{0} & t^3_{0} &\bigm| & x_t \\ t^0_{1} & t^1_{1} & t^2_{1} & t^3_{1} &\bigm| & x_{t+1} \\ t^0_{2} & t^1_{2} & t^2_{2} & t^3_{2} &\bigm| & x_{t+2} \\ \end{bmatrix} $$which may appear a bit daunting, but all of the values of \(t\) are known in advance and won't change for each eclipse. So the more specific version of the matrix is
$$ \begin{bmatrix} 1 & -2 & 4 & -8 &\bigm| & x_{t-2} \\ 1 & -1 & 1 & -1 &\bigm| & x_{t-1} \\ 1 & 0 & 0 & 0 &\bigm| & x_{t} \\ 1 & 1 & 1 & 1 &\bigm| & x_{t+1} \\ 1 & 2 & 4 & 8 &\bigm| & x_{t+2} \\ \end{bmatrix} $$Since the entire left side is pre-determined, it would be possible to pre-compute a solution. But this procedure can be used in other problems, so we will implement a generalized solution.
Once we have completed all of the row operations, we will be left with a solution in the form
$$ \begin{bmatrix} 1 & 0 & 0 & 0 &\bigm| & a_0 \\ 0 & 1 & 0 & 0 &\bigm| & a_1 \\ 0 & 0 & 1 & 0 &\bigm| & a_2 \\ 0 & 0 & 0 & 1 &\bigm| & a_3 \\ 0 & 0 & 0 & 0 &\bigm| & a_4 \\ \end{bmatrix} $$Where the \(a_i\) values are the coefficients we were solving for. They will be used to compute \(x\) for an instant in the form \(x = a_0 + a_1t + a_2t^2 + a_3t^3\).
We will walk through an example for the April 8, 2024 total solar eclipse. We will use the T0 = 2460409.25 TBD specified in the Five Millennium Canon of Solar Eclipses. For an ephemeris, we will use an extract of JPL DE 405 which contains only the time period necessary. For the radius of the Sun and Earth, we will use the values from IAU Resolution B3: Sun = 6.957x108m, Earth = 6.3781x106m. And for the diameter of the Moon we will use the value from IAU Resolution C 10 which provides the radius of the Moon in Earth diameters: k = 0.2725076. \( r \) and \( r' \) need to be converted from AU to units of Earth radii, so we multiply by \(\frac{1.496e8AU}{6378.1km}\).
We begin by computing the apparent Geocentric Sun and Moon positions for five hours surrounding T0:
| JD | α' | δ' | r' | α | δ | r |
|---|---|---|---|---|---|---|
| 2460409.1666666665 | 17.815452342642654 | 7.555699209337904 | 23489.964031827803 | 16.444559977344433 | 7.226772473505084 | 56.38611743996088 |
| 2460409.2083333335 | 17.853693662945386 | 7.571209938075441 | 23490.243282970325 | 17.005233846554173 | 7.518443264747429 | 56.397479333523016 |
| 2460409.25 | 17.8919368438214 | 7.586716964262048 | 23490.522447361745 | 17.566432059059977 | 7.809282507962961 | 56.409332542485835 |
| 2460409.2916666665 | 17.9301818901371 | 7.60222028069681 | 23490.80152488961 | 18.12817147239214 | 8.099256674796827 | 56.42167550190302 |
| 2460409.3333333335 | 17.968428806746825 | 7.617719880180428 | 23491.080515449772 | 18.69046940619963 | 8.38833258540859 | 56.4345065543071 |
Next, for each of the rows above, we use the equations above to compute the elements for each corresponding instant in time.
| X | Y | D | L1 | L2 | μ | tan F1 | tan F2 |
|---|---|---|---|---|---|---|---|
| -1.3415037973430761 | -0.3223608681175189 | 7.55648551145187 | 0.5355508655802353 | -0.010781309000520586 | 59.58304622079929 | 0.004666380467014835 | 0.004643121853949435 |
| -0.8299397492022488 | -0.05125078860638919 | 7.571334950213696 | 0.5356511688102383 | -0.01068150558143689 | 74.5871310505214 | 0.004666329257016898 | 0.004643070899208225 |
| -0.31825881990462546 | 0.21976896468344168 | 7.586180926048295 | 0.5357259496713748 | -0.010607097317524097 | 89.59121422033694 | 0.004666276981784016 | 0.004643018884541744 |
| 0.1934881223164183 | 0.4906702756870551 | 7.60102350033966 | 0.5357752188783229 | -0.010558073547510437 | 104.59529576917743 | 0.004666223641832559 | 0.004642965810463789 |
| 0.7052509386759692 | 0.761425358720428 | 7.6158627336629 | 0.5357989904821127 | -0.01053442029040278 | 119.59937573775537 | 0.0046661692378181 | 0.004642911677626662 |
For "tan F1/F2", for the final result, we will just take the value corresponding to T0 (some implementations choose to use an average). For the remaining columns, we need to perform interpolation to produce a third degree polynomial. Using X as an example, we create the augmented matrix $$ \begin{bmatrix} 1 & -2 & 4 & -8 &\bigm| & -1.3415037973430761 \\ 1 & -1 & 1 & -1 &\bigm| & -0.8299397492022488 \\ 1 & 0 & 0 & 0 &\bigm | & -0.31825881990462546 \\ 1 & 1 & 1 & 1 &\bigm | & 0.1934881223164183 \\ 1 & 2 & 4 & 8 &\bigm | & 0.7052509386759692 \\ \end{bmatrix} $$ and after Gauss-Jordan elimination we obtain $$ \begin{bmatrix} 1 & 0 & 0 & 0 &\bigm| & -0.3182588 \\ 0 & 1 & 0 & 0 &\bigm| & 0.5117224 \\ 0 & 0 & 1 & 0 &\bigm| & 0.0000330 \\ 0 & 0 & 0 & 1 &\bigm| & -0.0000085 \\ 0 & 0 & 0 & 0 &\bigm| & 0.0 \\ \end{bmatrix} $$
The last entry will generally not be exactly 0, but most implementations only use the first seven decimal places.
Applying the same interpolation process to the remaining columns, we obtain the final results
| x | y | d | l1 | l2 | μ |
|---|---|---|---|---|---|
| -0.3182588 | 0.2197690 | 7.5861809 | 0.5357259 | -0.0106071 | 89.5912142 |
| 0.5117224 | 0.2709652 | 0.0148443 | 0.0000620 | 0.0000617 | 15.0040824 |
| 0.0000330 | -0.0000592 | -0.0000017 | -0.0000128 | -0.0000127 | -0.0000008 |
| -0.0000085 | -0.0000047 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 |
tan F1 = 0.0046663, tan F2 = 0.0046430
Here is an example implementation in Javascript.