[[Eclipses]]
![[Attachments/Pasted image 20240321171904.png]]
https://phys.libretexts.org/Bookshelves/Astronomy__Cosmology/Celestial_Mechanics_(Tatum)/01%3A_Numerical_Methods/1.10%3A_1.10-Besselian_Interpolation
https://www.gregmiller.net/riseandsetgraph/index.html
https://www.gregmiller.net/eclipsesGeneratingBesselianElements/generatingBesselianPolynomials.html
![[Attachments/Pasted image 20240322131409.png]]
### Introduction
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](https://eclipse.gsfc.nasa.gov/SEpubs/5MCSE.html).
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.
### Terminology
Technically speaking, the values �,�,�,�1,�2,�,tan�1,tan�2x,y,d,l1,l2,μ,tanf1,tanf2 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.
### Process Overview and Prerequisites
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 �0T0, 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 �0t0 or �0t0, 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 �,�,�,�1,�2,�,tan�1,tan�2x,y,d,l1,l2,μ,tanf1,tanf2. 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.
### Computing the Besselian Elements for a given instant in time
The values �′,�′,�′α′,δ′,r′ represent the apparent geocentric right ascension, declination, and distance of the Sun. And the values �,�,�α,δ,r represent the same quantaties, but for the Moon.
The rectangular coordinates of the Sun:��=�′cos�′cos�′��=�′cos�′sin�′��=�′sin�′xs=ys=zs=r′cosδ′cosα′r′cosδ′sinα′r′sinδ′The rectangular coordinates of the Moon:��=�cos�cos���=�cos�sin���=�sin�xm=ym=zm=rcosδcosαrcosδsinαrsinδThe vectors �⃗S and �⃗M represent the geocentric postion of the Sun and Moon and the intermediate vector �⃗G represents a vecotr from the Sun to the Moon. And the scalar value �G is the distance from the Sun to the Moon.�⃗=[��,��,��]�⃗=[��,��,��]�⃗=�⃗−�⃗�=∣�⃗∣S=M=G=G=[xs,ys,zs][xm,ym,zm]S−M∣G∣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.tan�=�⃗�/��⃗�/�sin�=�⃗�/�tana=sind=Gx/GGy/GGz/G�μ is the Greenwich Hour Angle of point Z.�=�−�μ=θ−aThe 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.�=�cos�sin(�−�)�=�sin�cos�−cos�sin�cos(�−�)�=�sin�sin�+cos�cos�cos(�−�)x=y=z=rcosδsin(α−a)rsinδcosd−cosδsindcos(α−a)rsinδsind+cosδcosdcos(α−a)The constants ��ds 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.��=6.957�86.3781�6�=0.2725076ds=k=6.3781e66.957e80.2725076The remaining elements tan�1,tan�2,�1,tanf1,tanf2,l1, and �2l2 are computed as below. Note the intermediate sin��sinfx values are computed separately only for computational effenciency.sin�1=��+��sin�2=��−���1=�+�sin�1�2=�−�sin�2tan�1=tan(arcsin(sin�1))tan�2=tan(arcsin(sin�2))�1=�1tan�1�2=�2tan�2sinf1=sinf2=c1=c2=tanf1=tanf2=l1=l2=Gds+kGds−kz+sinf1kz−sinf2ktan(arcsin(sinf1))tan(arcsin(sinf2))c1tanf1c2tanf2The resulting set of Besselian Elements are �,�,�,�1,�2,�,tan�1,tan�2x,y,d,l1,l2,μ,tanf1,tanf2.
### Interpolating Elements
The value �,�,�,�1,�2,�x,y,d,l1,l2,μ from above must be interpolated to produce the Besselian Element Polynomials, and the values tan�1,tan�2tanf1,tanf2 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, �,�1,�2d,l1,l2 are usually second degree polynomials, and �μ 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:�=�0+�1�+�2�2+�3�3x=a0+a1t+a2t2+a3t3When we compute the elements for an instant, we choose the value of �t, so the goal is to interpolate for the cooeffients ��ai. There are many methods to do this, but we will use [Gauss-Jordan Elimination](https://en.wikipedia.org/wiki/Gaussian_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
[�−20�−21�−22�−23∣��−2�−10�−11�−12�−13∣��−1�00�01�02�03∣���10�11�12�13∣��+1�20�21�22�23∣��+2]⎣⎡t−20t−10t00t10t20t−21t−11t01t11t21t−22t−12t02t12t22t−23t−13t03t13t23∣∣∣∣∣∣∣∣∣∣xt−2xt−1xtxt+1xt+2⎦⎤
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
[1−24−8∣��−21−11−1∣��−11000∣��1111∣��+11248∣��+2]⎣⎡11111−2−101241014−8−1018∣∣∣∣∣∣∣∣∣∣xt−2xt−1xtxt+1xt+2⎦⎤
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
[1000∣�00100∣�10010∣�20001∣�30000∣�4]⎣⎡10000010000010000010∣∣∣∣∣∣∣∣∣∣a0a1a2a3a4⎦⎤
Where the ��ai values are the coefficients we were solving for. They will be used to compute �x for an instant in the form �=�0+�1�+�2�2+�3�3x=a0+a1t+a2t2+a3t3.
### Example
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.
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[1−24−8∣−1.34150379734307611−11−1∣−0.82993974920224881000∣−0.318258819904625461111∣0.19348812231641831248∣0.7052509386759692]⎣⎡11111−2−101241014−8−1018∣∣∣∣∣∣∣∣∣∣−1.3415037973430761−0.8299397492022488−0.318258819904625460.19348812231641830.7052509386759692⎦⎤and after Gauss-Jordan elimination we obtain[1000∣−0.31825880100∣0.51172240010∣0.00003300001∣−0.00000850000∣0.0]⎣⎡10000010000010000010∣∣∣∣∣∣∣∣∣∣−0.31825880.51172240.0000330−0.00000850.0⎦⎤
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](https://www.gregmiller.net/eclipsesGeneratingBesselianElements/example01/example.html).
### Computing the Besselian Elements for a given instant in time
The values �′,�′,�′α′,δ′,r′ represent the apparent geocentric right ascension, declination, and distance of the Sun. And the values �,�,�α,δ,r represent the same quantaties, but for the Moon.
The rectangular coordinates of the Sun:��=�′cos�′cos�′��=�′cos�′sin�′��=�′sin�′xs=ys=zs=r′cosδ′cosα′r′cosδ′sinα′r′sinδ′The rectangular coordinates of the Moon:��=�cos�cos���=�cos�sin���=�sin�xm=ym=zm=rcosδcosαrcosδsinαrsinδThe vectors �⃗S and �⃗M represent the geocentric postion of the Sun and Moon and the intermediate vector �⃗G represents a vecotr from the Sun to the Moon. And the scalar value �G is the distance from the Sun to the Moon.�⃗=[��,��,��]�⃗=[��,��,��]�⃗=�⃗−�⃗�=∣�⃗∣S=M=G=G=[xs,ys,zs][xm,ym,zm]S−M∣G∣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.tan�=�⃗�/��⃗�/�sin�=�⃗�/�tana=sind=Gx/GGy/GGz/G�μ is the Greenwich Hour Angle of point Z.�=�−�μ=θ−aThe 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.�=�cos�sin(�−�)�=�sin�cos�−cos�sin�cos(�−�)�=�sin�sin�+cos�cos�cos(�−�)x=y=z=rcosδsin(α−a)rsinδcosd−cosδsindcos(α−a)rsinδsind+cosδcosdcos(α−a)The constants ��ds 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.��=6.957�86.3781�6�=0.2725076ds=k=6.3781e66.957e80.2725076The remaining elements tan�1,tan�2,�1,tanf1,tanf2,l1, and �2l2 are computed as below. Note the intermediate sin��sinfx values are computed separately only for computational effenciency.sin�1=��+��sin�2=��−���1=�+�sin�1�2=�−�sin�2tan�1=tan(arcsin(sin�1))tan�2=tan(arcsin(sin�2))�1=�1tan�1�2=�2tan�2sinf1=sinf2=c1=c2=tanf1=tanf2=l1=l2=Gds+kGds−kz+sinf1kz−sinf2ktan(arcsin(sinf1))tan(arcsin(sinf2))c1tanf1c2tanf2The resulting set of Besselian Elements are �,�,�,�1,�2,�,tan�1,tan�2x,y,d,l1,l2,μ,tanf1,tanf2.
### Interpolating Elements
The value �,�,�,�1,�2,�x,y,d,l1,l2,μ from above must be interpolated to produce the Besselian Element Polynomials, and the values tan�1,tan�2tanf1,tanf2 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, �,�1,�2d,l1,l2 are usually second degree polynomials, and �μ 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:�=�0+�1�+�2�2+�3�3x=a0+a1t+a2t2+a3t3When we compute the elements for an instant, we choose the value of �t, so the goal is to interpolate for the cooeffients ��ai. There are many methods to do this, but we will use [Gauss-Jordan Elimination](https://en.wikipedia.org/wiki/Gaussian_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
[�−20�−21�−22�−23∣��−2�−10�−11�−12�−13∣��−1�00�01�02�03∣���10�11�12�13∣��+1�20�21�22�23∣��+2]⎣⎡t−20t−10t00t10t20t−21t−11t01t11t21t−22t−12t02t12t22t−23t−13t03t13t23∣∣∣∣∣∣∣∣∣∣xt−2xt−1xtxt+1xt+2⎦⎤
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
[1−24−8∣��−21−11−1∣��−11000∣��1111∣��+11248∣��+2]⎣⎡11111−2−101241014−8−1018∣∣∣∣∣∣∣∣∣∣xt−2xt−1xtxt+1xt+2⎦⎤
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
[1000∣�00100∣�10010∣�20001∣�30000∣�4]⎣⎡10000010000010000010∣∣∣∣∣∣∣∣∣∣a0a1a2a3a4⎦⎤
Where the ��ai values are the coefficients we were solving for. They will be used to compute �x for an instant in the form �=�0+�1�+�2�2+�3�3x=a0+a1t+a2t2+a3t3.
### Computing the Besselian Elements for a given instant in time
The values �′,�′,�′α′,δ′,r′ represent the apparent geocentric right ascension, declination, and distance of the Sun. And the values �,�,�α,δ,r represent the same quantaties, but for the Moon.
The rectangular coordinates of the Sun:��=�′cos�′cos�′��=�′cos�′sin�′��=�′sin�′xs=ys=zs=r′cosδ′cosα′r′cosδ′sinα′r′sinδ′The rectangular coordinates of the Moon:��=�cos�cos���=�cos�sin���=�sin�xm=ym=zm=rcosδcosαrcosδsinαrsinδThe vectors �⃗S and �⃗M represent the geocentric postion of the Sun and Moon and the intermediate vector �⃗G represents a vecotr from the Sun to the Moon. And the scalar value �G is the distance from the Sun to the Moon.�⃗=[��,��,��]�⃗=[��,��,��]�⃗=�⃗−�⃗�=∣�⃗∣S=M=G=G=[xs,ys,zs][xm,ym,zm]S−M∣G∣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.tan�=�⃗�/��⃗�/�sin�=�⃗�/�tana=sind=Gx/GGy/GGz/G�μ is the Greenwich Hour Angle of point Z.�=�−�μ=θ−aThe 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.�=�cos�sin(�−�)�=�sin�cos�−cos�sin�cos(�−�)�=�sin�sin�+cos�cos�cos(�−�)x=y=z=rcosδsin(α−a)rsinδcosd−cosδsindcos(α−a)rsinδsind+cosδcosdcos(α−a)The constants ��ds 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.��=6.957�86.3781�6�=0.2725076ds=k=6.3781e66.957e80.2725076The remaining elements tan�1,tan�2,�1,tanf1,tanf2,l1, and �2l2 are computed as below. Note the intermediate sin��sinfx values are computed separately only for computational effenciency.sin�1=��+��sin�2=��−���1=�+�sin�1�2=�−�sin�2tan�1=tan(arcsin(sin�1))tan�2=tan(arcsin(sin�2))�1=�1tan�1�2=�2tan�2sinf1=sinf2=c1=c2=tanf1=tanf2=l1=l2=Gds+kGds−kz+sinf1kz−sinf2ktan(arcsin(sinf1))tan(arcsin(sinf2))c1tanf1c2tanf2The resulting set of Besselian Elements are �,�,�,�1,�2,�,tan�1,tan�2x,y,d,l1,l2,μ,tanf1,tanf2.
### Interpolating Elements
The value �,�,�,�1,�2,�x,y,d,l1,l2,μ from above must be interpolated to produce the Besselian Element Polynomials, and the values tan�1,tan�2tanf1,tanf2 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, �,�1,�2d,l1,l2 are usually second degree polynomials, and �μ 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:�=�0+�1�+�2�2+�3�3x=a0+a1t+a2t2+a3t3When we compute the elements for an instant, we choose the value of �t, so the goal is to interpolate for the cooeffients ��ai. There are many methods to do this, but we will use [Gauss-Jordan Elimination](https://en.wikipedia.org/wiki/Gaussian_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
[�−20�−21�−22�−23∣��−2�−10�−11�−12�−13∣��−1�00�01�02�03∣���10�11�12�13∣��+1�20�21�22�23∣��+2]⎣⎡t−20t−10t00t10t20t−21t−11t01t11t21t−22t−12t02t12t22t−23t−13t03t13t23∣∣∣∣∣∣∣∣∣∣xt−2xt−1xtxt+1xt+2⎦⎤
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
[1−24−8∣��−21−11−1∣��−11000∣��1111∣��+11248∣��+2]⎣⎡11111−2−101241014−8−1018∣∣∣∣∣∣∣∣∣∣xt−2xt−1xtxt+1xt+2⎦⎤
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
[1000∣�00100∣�10010∣�20001∣�30000∣�4]⎣⎡10000010000010000010∣∣∣∣∣∣∣∣∣∣a0a1a2a3a4⎦⎤
Where the ��ai values are the coefficients we were solving for. They will be used to compute �x for an instant in the form �=�0+�1�+�2�2+�3�3x=a0+a1t+a2t2+a3t3.
https://peltiertech.com/excel-interpolation-formulas/
https://exceloffthegrid.com/interpolate-values-using-the-forecast-function/
https://www.exceldemy.com/do-linear-interpolation-in-excel/
BESSELI function
Excel for Microsoft 365 Excel for Microsoft 365 for Mac Excel for the web More...
We're offering you a free trial of Microsoft 365
This article describes the formula syntax and usage of the **BESSELI** function in Microsoft Excel.
Description
Returns the modified Bessel function, which is equivalent to the Bessel function evaluated for purely imaginary arguments.
BESSELI(X, N)
The BESSELI function syntax has the following arguments:
- **X** Required. The value at which to evaluate the function.
- **N** Required. The order of the Bessel function. If n is not an integer, it is truncated.
Lunar Eclipses
This Link has an index for all lunar eclipses from -1999 to +3000, predominantly a statistics page but also has this page that contains how to calculate when lunar eclipses are.
There is more than one formula depending on which time frame you are trying to look in.
This is the formula for eclipses between the year 2005 and 2050:
ΔT=62.92+0.32217∗t+0.005589∗t2
Where:
y=year+(month−0.5)/12
t=y−2000
Solar Eclipses
This Link has an index like above but for all of the solar eclipses from -1999 to +3000.
This link has the formula for calculating solar eclipses. This is the formula for between 2005 and 2050:
ΔT=62.92+0.32217∗t+0.005589∗t2
Where:
y=year+(month−0.5)/12
t=y−2000
In order to predict the general characteristics as well as the local circumstances for a solar eclipse, high accuracy ephemerides for both the Sun and the Moon are required. Conventional ephemerides tabulate the positions and distances of these bodies with respect to Earth's center. However, the solar eclipse calculator is primarily interested in the position, dimensions and velocity of the Moon's shadow on Earth's surface.
In 1824, the Prussian astronomer and mathematician Friedrich Bessel introduced a new method for the prediction of solar eclipses. It was so successful that it remains today as the most powerful technique, even with the application of the digital computer. The key to Bessel's method is the expression of the ephemerides of the Sun and Moon in terms of the Moon's shadow with respect to Earth's center. This change in the frame of reference greatly simplifies the mathematics and geometry without any sacrifice in accuracy.
To define the Besselian elements of an eclipse, a plane is passed through the center of Earth which is fixed perpendicular to the axis of the Moon's shadow. This is called the _fundamental plane_ on which on X-Y rectangular coordinate system is constructed with its origin at the geocenter. The axes of this system are oriented with north in the positive Y direction and east in the positive X direction. The Z axis is perpendicular to the fundamental plane and parallel to the shadow axis.
The X-Y coordinates of the shadow axis can now be expressed in units of the equatorial radius of Earth. The radii of the penumbral and umbral shadows on the fundamental plane are identified as L1 and L2, respectively. The direction of the shadow axis on the celestial sphere is defined by its declination _'d'_ and ephemeris hour angle _'µ'_. Finally, the angles which the penumbral and umbral shadow cones make with the shadow axis are expressed as f1 and f2, respectively. These eight parameters and their variations over time serve as the only input needed to characterize a solar eclipse. They may be tabulated at hourly intervals, or expressed as a series of third order polynomials over the several hour period of the eclipse.
The eight Besselian elements needed to characterize a solar eclipse can be summarized as follows:
- _**x, y**_ - Cartesian coordinates of the lunar shadow axis in the _Fundamental Plane_ (in units of Earth's equatorial radius)
- _**L1, L2**_ - Radii of the Moon's _penumbral_ and _umbral/antumbral_ shadows in the _Fundamental Plane_ (in units of Earth's equatorial radius)
- _**d**_ - Declination of the Moon's shadow axis on the celestial sphere
- _**µ**_ - Hour angle of the Moon's shadow axis on the celestial sphere
- _**f1, f2**_ - Angles of the _penumbral_ and _umbral/antumbral_ shadow cones with respect to the axis of the lunar shadow
The details for actual eclipse calculations using the Besselian elements can be found in the references listed below.
![[Attachments/Pasted image 20240320053524.png]]
![[Attachments/Pasted image 20240319032515.png]]
https://en.wikipedia.org/wiki/Ecliptic_coordinate_system
## Remarks
- If x is nonnumeric, BESSELI returns the #VALUE! error value.
- If n is nonnumeric, BESSELI returns the #VALUE! error value.
- If n < 0, BESSELI returns the #NUM! error value.
- The n-th order modified Bessel function of the variable x is:
## Example
Copy the example data in the following table, and paste it in cell A1 of a new Excel worksheet. For formulas to show results, select them, press F2, and then press Enter. If you need to, you can adjust the column widths to see all the data.
https://phys.libretexts.org/Bookshelves/Astronomy__Cosmology/Celestial_Mechanics_(Tatum)/01%3A_Numerical_Methods/1.10%3A_1.10-Besselian_Interpolation
https://www.gregmiller.net/Sky_Path/
https://engineerexcel.com/bilinear-interpolation-excel/
https://exceloffthegrid.com/interpolate-values-using-the-forecast-function/
In the days before the widespread use of high-speed computers, extensive use was commonly made of printed tables of the common mathematical functions. For example, a table of the Bessel function J0(x)�0(�) would indicate
J0(1.7)=0.397 984 859�0(1.7)=0.397 984 859
J0(1.8)=0.339 986 411�0(1.8)=0.339 986 411
If one wanted the Bessel function for x=1.762�=1.762, one would have to interpolate between the tabulated values.
Today it would be easier simply to calculate the Bessel function for any particular desired value of the argument x�, and there is less need today for printed tables or to know how to interpolate. Indeed, most computing systems today have internal routines that will enable one to calculate the commoner functions such as Bessel functions even if one has only a hazy notion of what a Bessel function is.
The need has not entirely passed, however. For example, in orbital calculations, we often need the geocentric coordinates of the Sun. These are not trivial for the nonspecialist to compute, and it may be easier to look them up in _The Astronomical Almanac_, where it is tabulated for every day of the year, such as, for example, July 14 and July 15. But, if one needs y� for July 14.395, how does one interpolate?
In an ideal world, a tabulated function would be tabulated at sufficiently fine intervals so that linear interpolation between two tabulated values would be adequate to return the function to the same number of significant figures as the tabulated points. The world is not perfect, however, and to achieve such perfection, the tabulation interval would have to change as the function changed more or less rapidly. We need to know, therefore, how to do nonlinear interpolation.
Suppose a function y(x)�(�) is tabulated at x=x1�=�1 and x=x2�=�2, the interval x2−x1�2−�1 being δxδ�. If one wishes to find the value of y� at x+θδx�+�δ�, linear interpolation gives
y(x1+θΔx)=y1+θ(y2−y1)=θy2+(1−θ)y1,(1.10.1)(1.10.1)�(�1+�Δ�)=�1+�(�2−�1)=��2+(1−�)�1,
where y1=y(x1)�1=�(�1) and y2=y(x2)�2=�(�2). Here it is assumed that is a fraction between 00 and 11; if θ� is outside this range (that is negative, or greater than 11) we are extrapolating, not interpolating, and that is always a dangerous thing to do.
Let us now look at the situation where linear interpolation is not good enough. Suppose that a function y(x)�(�) is tabulated for four points x1, x2, x3, x4�1, �2, �3, �4 of the argument x�, the corresponding values of the function being y1, y2, y3, y4�1, �2, �3, �4. We wish to evaluate y� for x=x2+θδx�=�2+�δ�, where δxδ� is the interval x2−x1�2−�1 or x3−x2�3−�2 or x4−x3�4−�3. The situation is illustrated in figure I.6AI.6A.
A possible approach would be to fit a polynomial to the four adjacent points:
y=a+bx+cx2+dx3.(1.10.1)(1.10.1)�=�+��+��2+��3.
We write down this Equation for the four adjacent tabulated points and solve for the coefficients, and hence we can evaluate the function for any value of x� that we like in the interval between x1�1 and x4�4. Unfortunately, this might well involve more computational effort than evaluating the original function itself.
FIGURE I.6AFIGURE I.6A
The problem has been solved in a convenient fashion in terms of finite difference calculus, the logical development of which would involve an additional substantial chapter beyond the intended scope of this book. I therefore just provide the method only, without proof.
The essence of the method is to set up a table of differences as illustrated below. The first two columns are x� and y�. The entries in the remaining columns are the differences between the two entries in the column immediately to the left. Thus, for example, δ4.5=y5−y4δ4.5=�5−�4, δ24=δ4.5−δ3.5δ42=δ4.5−δ3.5, etc.
Let us suppose that we want to find y� for a value of x� that is a fraction θ� of the way from x4�4 to x5�5. Bessel's interpolation formula is then
y(x)=12(y4+y5)+B1δ4.5+B2(δ24+δ25)+B3δ34.5+B4(δ44+δ45)+...(1.10.3)(1.10.3)�(�)=12(�4+�5)+�1δ4.5+�2(δ42+δ52)+�3δ4.53+�4(δ44+δ54)+...
Here the Bn�� are the Besselian interpolation coefficients, and the successive terms in parentheses in the expansion are the sums of the numbers in the boxes in the table.
The Besselian coefficients are
Bn(θ)=12(θ+12n−1n)if n is even,(1.10.4)(1.10.4)��(�)=12(�+12�−1�)if n is even,
and
Bn(θ)=θ−12n(θ+12n−32n−1)if n is odd.(1.10.5)(1.10.5)��(�)=�−12�(�+12�−32�−1)if n is odd.
The notation (mn)(��) means the coefficient of x m in the binomial expansion of (1+x)n(1+�)�.
Explicitly,
B1=θ−12(1.10.6)(1.10.6)�1=�−12
B2=12θ(θ−1)/2!=θ(θ−1)/4(1.10.7)(1.10.7)�2=12�(�−1)/2!=�(�−1)/4
B3=(θ−12)θ(θ−1)/3!=θ(0.5+θ(−1.5+θ))/6(1.10.8)(1.10.8)�3=(�−12)�(�−1)/3!=�(0.5+�(−1.5+�))/6
B4=12(θ+1)θ(θ−1)(θ−2)/4!=θ(2+θ(−1+θ(−2+θ)))/48(1.10.9)(1.10.9)�4=12(�+1)�(�−1)(�−2)/4!=�(2+�(−1+�(−2+�)))/48
B5=(θ−12)(θ+1)θ(θ−1)(θ−2)/5!=θ(−1+θ(2.5+θ2(−2.5+θ)))/120(1.10.10)(1.10.10)�5=(�−12)(�+1)�(�−1)(�−2)/5!=�(−1+�(2.5+�2(−2.5+�)))/120
The reader should convince him- or herself that the interpolation formula taken as far as B1�1 is merely linear interpolation. Addition of successively higher terms effectively fits a curve to more and more points around the desired value and more and more accurately reflects the actual change of y� with x�.
The above table is taken from _The Astronomical Almanac_ for 1997, and it shows the y�-coordinate of the Sun for eight consecutive days in July. The first three difference columns are tabulated, and it is clear that further difference columns are unwarranted.
If we want to find the value of y�, for example, for July 4.746, we have θ=0.746�=0.746 and the first three Bessel coefficients are
B1B2B3===+0.246−0.047 371−0.007 768 844�1=+0.246�2=−0.047 371�3=−0.007 768 844
The reader can verify the following calculations for y� from the sum of the first 2, 3 and 4 terms of the Besselian interpolation series formula. The sum of the first two terms is the result of linear interpolation.
Sum of the first 2 terms, y=0.909 580 299�=0.909 580 299
Sum of the first 3 terms, y=0.909 604 723�=0.909 604 723
Sum of the first 4 terms, y=0.909 604 715�=0.909 604 715
Provided the table is not tabulated at inappropriately coarse intervals, one need rarely go past the third Bessel coefficient. In that case an alternative and equivalent interpolation formula (for t=t4+θΔt)�=�4+�Δ�), which avoids having to construct a difference table, is
y(t4+θΔt)=−16θ[(2−θ(3−θ))y3+(1−θ)y6]+12[(2+θ(−1+θ(−2+θ)))y4+θ(2+θ(1−θ))y5].(1.10.2)(1.10.2)�(�4+�Δ�)=−16�[(2−�(3−�))�3+(1−�)�6]+12[(2+�(−1+�(−2+�)))�4+�(2+�(1−�))�5].
Readers should check that this gives the same answer, at the same time noting that the nested parentheses make the calculation very rapid and they are easy to program on either a calculator or a computer.
##### Exercise 1.10.11.10.1: Bessel Coefficients
From the following table, construct a difference table up to fourth differences. Calculate the first four Bessel coefficients for θ=0.73�=0.73. Hence calculate the value of y� for x=0.273�=0.273.
x0.00.10.20.30.40.5y+0.381300+0.285603+0.190092+0.096327+0.008268−0.067725(1.10.3)(1.10.3)��0.0+0.3813000.1+0.2856030.2+0.1900920.3+0.0963270.4+0.0082680.5−0.067725
**Answers**
- B1=+0.23�1=+0.23
- B2−−0.049275�2−−0.049275
- B3=−7.5555×10−3�3=−7.5555×10−3
- B4=+9.021841875×10−3�4=+9.021841875×10−3
- y=0.121289738�=0.121289738
Note: the table was calculated from a formula, and the interpolated answer is correct to nine significant figures.
##### Exercise 1.10.21.10.2: Linear Interpolation vs. Besselian Interpolation
From the following table of sinxsin�, use linear interpolation and Besselian interpolation to estimate sin51∘sin51∘ to three significant figures.
x∘0306090sinx0.00.53–√/2−0.866031.0(1.10.4)(1.10.4)�∘sin�00.0300.5603/2−0.86603901.0
**Answers**
- By linear interpolation, sin51∘=0.634.sin51∘=0.634.
- By Besselian interpolation, sin51∘=0.776.sin51∘=0.776.
The correct value is 0.777. You should be impressed – but there is more on interpolation to come, in Section 1.11.
---
This page titled [1.10: 1.10- Besselian Interpolation](https://phys.libretexts.org/Bookshelves/Astronomy__Cosmology/Celestial_Mechanics_(Tatum)/01%3A_Numerical_Methods/1.10%3A_1.10-Besselian_Interpolation) is shared under a [CC BY-NC 4.0](https://creativecommons.org/licenses/by-nc/4.0) license and was authored, remixed, and/or curated by [Jeremy Tatum](http://orca.phys.uvic.ca/~tatum/index.php) via [source content](http://orca.phys.uvic.ca/~tatum/celmechs.html) that was edited to the style and standards of the LibreTexts platform; a detailed edit history is available upon request.
# Bessel interpolation formula
A formula which is defined as half the sum of the Gauss formula (cf.[Gauss interpolation formula]( "Gauss interpolation formula")) for forward interpolation on the nodes
x0,x0+h,x0‘hx0+nh,x0‘nh,x0+(n+1)h,0,0+,0‘0+,0‘,0+(+1),
at the pointx=x0+th=0+:
G2n+2(x0+th)=f0+f11/2t+f20t(t‘1)2!++(1)(1)2+2(0+)=0+1/21+02(‘1)2!++
+f2n+11/2t(t2‘1)(t2‘n2)(2n+1)!,+1/22+1(2‘1)(2‘2)(2+1)!,
and the Gauss formula of the same order for backward interpolation with respect to the nodex1=x0+h1=0+, i.e. with respect to the population of nodes
x0+h,x0,x0+2h,x0‘hx0+(n+1)h,x0‘nh:0+,0,0+2,0‘0+(+1),0‘:
G2n+2(x0+th)=f1+f11/2(t‘1)+f21t(t‘1)2!++(2)(2)2+2(0+)=1+1/21(‘1)+12(‘1)2!++
+f2n+11/2t(t2‘1)[t2‘(n‘1)2](t‘n)(t‘n‘1)(2n+1)!.+1/22+1(2‘1)[2‘(‘1)2](‘)(‘‘1)(2+1)!.
Putting
f2k1/2=(f2k0+f2k1)2,1/22=(02+12)2,
Bessel's interpolation formula assumes the form ([[1]]([[2]](
B2n+2(x0+th)=(3)(3)2+2(0+)=
=f1/2+f11/2(t‘12)+f21/2t(t‘1)2!++=1/2+1/21(‘12)+1/22(‘1)2!++
+f2n1/2t(t2‘1)[t2‘(n‘1)2](t‘n)(2n)!++1/22(2‘1)[2‘(‘1)2](‘)(2)!+
+f2n+11/2t(t2‘1)[t2‘(n‘1)2](t‘n)(t‘1/2)(2n+1)!.+1/22+1(2‘1)[2‘(‘1)2](‘)(‘1/2)(2+1)!.
Bessel's interpolation formula has certain advantages over Gauss' formulas (1), (2); in particular, if the interpolation is at the middle of the segment, i.e. att=1/2=1/2, all coefficients at the differences of odd orders vanish. If the last term on the right-hand side of (3) is omitted, the polynomialB2n+1(x0+th)2+1(0+), which is not a proper interpolation polynomial (it coincides withf(x)()only in the2n2nodesx0‘(n‘1)hx0+nh0‘(‘1)0+), represents a better estimate of the residual term (cf.[Interpolation formula]( "Interpolation formula")) than the interpolation polynomial of the same degree. Thus, for instance, ifx=x0+th(x0,x1)=0+(0,1), the estimate of the last term using the polynomial which is most frequently employed
B3(x0+th)=f1/2+f11/2(t‘12)+f21/2t(t‘1)2,3(0+)=1/2+1/21(‘12)+1/22(‘1)2,
written with respect to the nodesx0‘h,x0,x0+h,x0+2h0‘,0,0+,0+2, is almost 8 times better than that of the interpolation polynomial written with respect to the nodesx0‘h,x0,x0+h0‘,0,0+orx0,x0+h,x0+2h0,0+,0+2([[2]](
https://www.gregmiller.net/MeeusEclipseExamples/
Below is an example set of elements for the April 8, 2024 total solar eclipse from the [Five Millennium Canon of Solar Eclipses](https://eclipse.gsfc.nasa.gov/SEpubs/5MCSE.html).
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"
### Terminology
Technically speaking, the values �,�,�,�1,�2,�,tan�1,tan�2x,y,d,l1,l2,μ,tanf1,tanf2 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".
![[Attachments/Pasted image 20240403053228.png]]
![[Attachments/Pasted image 20240403053238.png]]
![[Attachments/Pasted image 20240403053249.png]]
![[Attachments/Pasted image 20240403053300.png]]
#### Twitter Response
Did someone ask how I predicted the apparent shadow for this upcoming total solar eclipse on April 8th? https://twitter.com/AntiDisinfo86/status/1770413357209702404
Wrong, Wrong, and Wrong. Oh and Wong again.
Todays lesson is courtesy of @BradleyMcK67 , who likes to repeat talking points with a cursory understanding of the topic. He chose to insist that I 'stole' something from NASA.... let us correct that misunderstanding. For the record, I would not repeat cheap trash FTFE posts, especially since his sole intent with any targeted post is to send his army of trolls to swarm.
The thing I 'stole' was the definition of Besselian Elements. I attached the offending slide.
1 No one predicts anything, there are cycles only, and no predictions.
2 No one at NASA has any responsibility for any eclipse 'predictions', and it says this on every single webpage they mention the word eclipse. (means you have never seen a single webpage about eclipses from NASA)
3 Besselian Elements are absolutely not the result of transform equations, that data is not released. Go and find me one of these transform equations ...The Besselian Elements are only valid for about 5 hours, so I am VERY curious where NASA posted all these for the public to easily use , in 5 hour increments for all of eternity. No, you ACTUALLY just saw a meme from flatearth.ws that said they can calculate the position of the sun and moon perpetually, but that data actually comes from completing the great astrological year of 651 years, so that you can 'predict' in perpetuity, so kind of a chicken or the egg situation there.
4 Try looking into the actual sources, Fred's publications (i know you wont, I attached them at the bottom)
5 Besselian Elements are the result of complicated calculations using a form of interpolation to predict unseen portions of a trend line. I used excel for this, but here is the process. its a decent overview: https://www.gregmiller.net/eclipsesGeneratingBesselianElements/generatingBesselianPolynomials.html
5 Once you look into what the Besselian Elements ACTUALLY ARE, tell me how compromising it was to copy and paste the definition. Especially since it was to show people like you what I had calculated in the first place, since you clearly have no idea. 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.
6 They don't even come from NASA. Moreover, they dont even come from Fred, They come from Fred's Partner Jean Meeus.
7 Here is a video where I outline the process of calculating the Besselian Elements. You wont watch, but that's ok, this post is no longer about you, your purpose has been served. https://x.com/AntiDisinfo86/status/1771400846078587313?s=20
Sources
Fifty Year Canon of Solar Eclipses: 1986 - 2035 https://eclipse.gsfc.nasa.gov/5MCLE/5MKLE-214173.pdf Fifty Year Canon of Lunar Eclipses: 1986 - 2035 https://eclipse.gsfc.nasa.gov/SEpubs/FiftyYearCanonOfSolarEclipses1986-2035_NASA1178a.pdf
The actual sources I listed and referenced when I calculated all the stuff that you clearly had no idea how to do.
https://publish.obsidian.md/shanesql/Eclipse+Sources
FINAL NOTE*** I had to use the celestial sphere to calculate these, and that is exclusively geocentric.
So.. yea.. there you go.
Eclipse Maps: (all flat of course) https://publish.obsidian.md/shanesql/Eclipse+Maps
Eclipse PDFS:
https://publish.obsidian.md/shanesql/Eclipse+PDFs
How eclipses actually work: https://publish.obsidian.md/shanesql/Eclipse+Maps
Besselian Elements & Interpolation
https://publish.obsidian.md/shanesql/Besselian+Elements+%26+Interpolation
Works cited:
Trowbridge, D. (1874). "The Recurrence of Eclipses." The Analyst, 1(1), 3-5. Retrieved from https://www.jstor.org/stable/2635490.
Rehak, A. (2002). "Regression of Lunar nodes."
Dimbleby, J. B. (1894). "All Past Time: All eclipses and transits from creation."
Hayakawa, H. (2021). "Analyses of a Datable Solar Eclipse Record in Maya Classic Period Monumental Inscriptions." Publications of the Astronomical Society of Japan, D.
HUNGER, H. (1988). "ASTRONOMICAL DIARIES AND RELATED TEXT FROM BABYLONIA."
Fendt, W. (2002). "Basics of Positional Astronomy and Ephemerides."
Costantino Sigismondi, P. D. V. (2024). "Eclipses: A Brief History of Celestial Mechanics, Astrometry and Astrophysics." ICRANet International Center for Relativistic Astrophysics Network.
Espenak, F. (2009). "Five Millennium Catalog of Lunar Eclipses."
Espenak, F. (2009). "Five Millennium Catalog of Solar Eclipses."
Jeffreys, H. (2002). "Historical eclipses and Earth’s rotation."
RAPPENGLÜCK, M. A. (2001). "PALAEOLITHIC TIMEKEEPERS LOOKING AT THE GOLDEN GATE."
Trowbridge, D. (1874). "The Recurrence of Eclipses." The Analyst, 1(1), 3-5.
Brad, B. (1995). "The Saros Cycle - Eclipses come in Families." www.BernadetteBrady.com.
Sigismondi, C. (2006). "Solar radius in 2005 and 2006 eclipses." ICRA, International Center for Relativistic Astrophysics.