Calculate Future Solar Eclipses

According to Herodotus, who wrote The History around 430 BCE, a total solar eclipse interrupted a battle between the Lydians and the Medes:

[J]ust as the battle was growing warm, day was on a sudden changed into night. This event had been foretold by Thales, the Milesian, who forewarned the Ionians of it, fixing for it the very year in which it actually took place. The Medes and Lydians, when they observed the change, ceased fighting, and were alike anxious to have terms of peace agreed on.

(Source: Internet Classic Archives)

This “battle of the eclipse” took place in what is now Turkey, and the eclipse itself could have been the one that occurred on 28 May 585 BCE. Unfortunately, if Thales wrote anything it didn’t survive the centuries. So, we don’t know if he really did predict the eclipse or by what means he did it. All we have is Herodotus to go on. In any case, it’s the first mention in history of anyone predicting an eclipse so it’s interesting on that basis alone.

Wouldn’t it be cool if we could figure out when future solar eclipses occurred for us in the next decade or so? That’s what we’re going to do in this blog post so read on.

What is an Eclipse?

I highly recommend reading the article “What Are Solar Eclipses?” on timeanddate.com if you’re just a little fuzzy about it all. Having said that, let me try my hand at a brief explanation. If you want to dive right into the algorithm, then skip down to the next section.

Have you ever had someone shine a bright flashlight in your face and you’ve instinctively raised a hand to shield your eyes? If so, then you’ve just modeled a solar eclipse. The light is the Sun, your eyes the Earth, and your hand is like the Moon in between them.

Now imagine you are facing a round light bulb that is several meters away. The bulb has an actual diameter of several centimeters. Let’s say 10 cm. But because you are several meters away it will look smaller. This is the bulb’s angular diameter or it’s apparent size from your perspective. Let’s say it’s angular diameter from your perspective is 2 cm. Suppose you now hold up your thumb to block the light. Your thumb might be big enough to cover the entire 2 cm apparent size of the bulb. If so, this is like a total solar eclipse. Suppose your thumb is smaller than 2 cm and light escapes around the edges in all directions. This is like an annular solar eclipse. Or your thumb could be smaller or larger than 2 cm but off center so that it covers only one side of the bulb. That’s a partial solar eclipse which looks like a cookie with a bite out of it. There’s a fourth type of solar eclipse called a hybrid which is not common.

The ancients believed that celestial objects moved in perfect circles and were fixed in their orbits for all time. If so, that would have made these calculations trivial! As it turns out bodies move in elliptical orbits, the moon has its own eccentric orbit of about 0.05 off the plane of the ecliptic, and the gravity of the objects in our solar system pull and tug at each other so that predicting their movements from one decade to another is very difficult.

But if you know where the Sun, Moon, and Earth are at a certain point in time, and you know how they move relative to each other in space, then you can calculate their movements over a period of months and years. But just like in a foot race where you need a starting point and a stopwatch, to calculate these movements you need a starting point and a stopwatch to act as a baseline. This baseline is called the “standard” equinox or the J2000.0 epoch. The starting point for J2000.0 is defined as the moment of 12:00 noon UTC on 1 January 2000.

A solar eclipse can only occur during a new moon phase when the moon is said to be in inferior conjunction. So, our task is to evaluate every new moon over the next several years and calculate the positions of the Earth, Moon, and Sun in relation to each other on each of those dates to see if an eclipse event occurs.

The Algorithm

The approach I’ll use is from the book Astronomical Algorithms (2nd Ed) by Jean Meeus. He describes a method for doing the calculations based on the Bureau des Longitudes’ ELP-2000/82 theory developed primarily by Jean Chapront and Michelle Chapront-Touzé. Here are the steps we need to take:

1. For each new moon, calculate the positions of the Sun, Moon, and Earth on that date. These values are captured in the variables: k, T, JDE, E, Sm, Mm, F, O, Fp, and Ap.

2. Calculate the least distance from the axis of the Moon’s shadow to the center of the Earth. Variables: W and g.

3. Calculate the Moon’s umbral cone in the fundamental plane. Variable: u.

4. If the Moon’s position is such that some type of eclipse is visible from the Earth’s surface then determine the eclipse type: Partial, Total, Annular, or Hybrid (sometimes called Annular-Total). If it’s partial, then calculate the magnitude of the partial eclipse. Finally, calculate the time of greatest eclipse and convert to UTC.

The program I wrote for this post (checked into my GitHub repo here) will follow these steps and display every solar eclipse with some basic metadata.

First thing I did was to set up a helper class that contained every new moon between the years 2022 and 2035 along with a method for enumerating each date:

internal class NewMoonData
{
    // All new moon dates from 2022 thru 2035. Source: timeanddate.com
    private static readonly List<DateTime> _newMoons = new()
    {
        new DateTime(2022, 1, 2),
        new DateTime(2022, 2, 1),
        new DateTime(2022, 3, 2),
        new DateTime(2022, 4, 1),
        new DateTime(2022, 4, 30),
        new DateTime(2022, 5, 30),

        ...

        new DateTime(2035, 10, 1),
        new DateTime(2035, 10, 31),
        new DateTime(2035, 11, 29),
        new DateTime(2035, 12, 29)
    };

    public static IEnumerable<DateTime> EachNewMoon(DateTime from, DateTime thru)
    {
        return _newMoons.Where(m => m.Year >= from.Year && m.Year <= thru.Year);
    }
}

With the NewMoonData helper class in place we can set our date range and enumerate through the list of new moons:

// examine all new moons in date range
DateTime startDate = new(2022, 1, 1);
DateTime endDate = new(2035, 12, 31);

foreach (var dt in NewMoonData.EachNewMoon(startDate, endDate))
{
    /* Algorithm */
}

OK with that bit of plumbing out of the way let’s go step by step. Step 1 is to calculate the positions of the Sun, Earth, and Moon. I mentioned above that our starting point is the standard equinox of J2000.0. The first new moon after that baseline was 6 Jan 2000. The variable k represents that date and so it’s set to k = 0. For each new moon we calculate k and T and then use those two values to get the JDE. This is the time when the moon is in conjunction and an eclipse is possible. I’ll detail the calculations for k, T, and JDE below. If you prefer to just look at them in the finished app go to my repo and clone the C# solution.

Step 1

// for example the new moon on 2 Jan 2022

DateTime d = new(2022, 1, 2);

// get its decimal year
double daysInYear = (DateTime.IsLeapYear(d.Year)) ? 366 : 365.2425;
var fraction = d.DayOfYear / daysInYear;
double dy = d.Year + fraction;

// k
double decimalYear = (dy - 2000) * 12.3685;
int k = (int)Math.Floor(decimalYear);

// Time T
double T = k / 1236.85;

// JDE of moon's conjunction
double JDE = 2451550.09766 +
    (29.530588861 * k) +
    (0.00015437 * (T * T)) -
    (0.000000150 * (T * T * T)) +
    (0.00000000073 * (T * T * T * T));

Next up for Step 1 we need E (Earth’s eccentricity in its orbit), Sm (Sun’s mean anomaly at time JDE), Mm (Moon’s mean anomaly at time JDE), F (Moon’s argument of latitude), O (longitude of Moon’s ascending node), and finally Fp and Ap (angles to obtain time of maximum eclipse). Here are those calculations.

// E
double E = 1 - 0.002516 * T - 0.0000074 * T * T;

// Sm 
// (not shown: reducing the angle and converting to radians)
double Sm = 2.5534 +
    (29.10535670 * k) -
    (0.0000014 * (T * T)) -
    (0.00000011 * (T * T * T));

// Mm
// (not shown: reducing the angle and converting to radians)
double Mm = 201.5643 +
    (385.81693528 * k) +
    (0.0107582 * (T * T)) +
    (0.00001238 * (T * T * T)) -
    (0.000000058 * (T * T * T * T));

// F
double F = 160.7108 +
    (390.67050284 * k) -
    (0.0016118 * (T * T)) -
    (0.00000227 * (T * T * T)) +
    (0.000000011 * (T * T * T * T));

// O (Ω) 
double O = 124.7746 -
    (1.56375588 * k) +
    (0.0020672 * (T * T)) +
    (0.00000215 * (T * T * T));

// Fp
double Fp = F - (0.02665 * Math.Sin(O));

// Ap
double Ap = 299.77 + (0.107408 * k) - (0.009173 * (T * T));

All of those variables in Step 1 give us the positions of the Sun, Earth, and Moon relative to each other on the new moon date in question (2 Jan 2022).

Step 2

This step involves calculating some intermediate variables to get to g (γ or gamma). Meeus writes that gamma (g) “represents the least distance from the axis of the Moon’s shadow to the center of the Earth, in units of the equatorial radius of the Earth.” Yes, that’s a mouthful. It is in reference to the solution that the Prussian mathematician Friedrich Bessel came up with for calculating eclipses. Without getting into too much detail think about how complicated it would have been in Bessel’s day to calculate an oval-shaped shadow on the Earth’s curved surface. What Bessel did to simplify the problem was to insert a two-dimensional plane that ran through the center of the Earth and lay perpendicular to the line created by the Moon’s shadow. This is called the fundamental plane. Now the problem is simple geometry because it’s easier to work with a circle on a Cartesian grid. Here’s a picture from eclipse2024.org and their web site has an excellent explanation of Bessel’s solution.

So that’s what Meeus means when he refers to g being the least distance to the “center of the Earth.” We need three preliminary variables P, Q, and W (which are angular measurements of the moon’s shadow) and then we can calculate g:

// P
double P =
    (0.2070 * E * Math.Sin(Sm)) +
    (0.0024 * E * Math.Sin(Sm * 2)) -
    (0.0392 * Math.Sin(Mm)) +
    (0.0116 * Math.Sin(Mm * 2)) -
    (0.0073 * E * Math.Sin(Mm + Sm)) +
    (0.0067 * E * Math.Sin(Mm - Sm)) +
    (0.0118 * Math.Sin(Fp * 2));

// Q
double Q =
    5.2207 -
    (0.0048 * E * Math.Cos(Sm)) +
    (0.0020 * E * Math.Cos(Sm * 2)) -
    (0.3299 * Math.Cos(Mm)) -
    (0.0060 * E * Math.Cos(Mm + Sm)) +
    (0.0041 * E * Math.Cos(Mm - Sm));

// W
double W = Math.Abs(Math.Cos(Fp));

// g
double g = (P * Math.Cos(Fp) + Q * Math.Sin(Fp)) * (1 - 0.0048 * W);

Gamma g will be positive if the shadow hits the Northern Hemisphere and negative if it hits the Southern Hemisphere.

Step 3

Next, we calculate the Moon’s umbral cone on the fundamental plane as shown in the picture above.

// u
double u =
    0.0059 +
    0.0046 * E * Math.Cos(Sm) -
    0.0182 * Math.Cos(Mm) +
    0.0004 * Math.Cos(Mm * 2) -
    0.0005 * Math.Cos(Sm + Mm);

Step 4

Now we’re in a position to be able to determine if an eclipse occurs on the new moon date in question! Ignoring the sign of g for a moment, we take its absolute value to see if there can even be an eclipse:

If $$| g | > 1.5433 + u$$ then no eclipse is visible. If $$0.9972 <= | g | >= 1.5433 + u$$ then the eclipse is not central but partial. If | g | is outside those bounds and is central then it is either total (u < 0), annular (u > 0.0047), or one of annular or hybrid. Here’s the C# code to wade through these possibilities. I hate the nested if statements, but you get the point:

double absG = Math.Abs(g);
double mag = 0.0;

if (absG <= 1.5433 + u) // then some type of eclipse is visible from the Earth's surface
{
    string eclipseType = "Unknown";
    if (absG >= 0.9972 && absG <= 1.5433 + u)
    {
        // the eclipse is not central but partial
        eclipseType = "Partial";
    }
    else
    {
        // the eclipse is central
        if (u < 0)
        {
            eclipseType = "Total";
        }
        else if (u > 0.0047)
        {
            eclipseType = "Annular";
        }
        else // between 0 and 0.0047 so either annular or hybrid
        {
            var w = 0.00464 * Math.Sqrt(1 - g * g);
            eclipseType = (u < w) ? "Hybrid" : "Annular";
        }
    }
}

From here I can convert to the time of the greatest eclipse and display the same metadata that is published by NASA in the Five Millennium Canon of Solar Eclipses. (See their technical publication here.) If I run the program here’s the displayed output:

Here’s a snapshot of the maps from the NASA Canon for the last four dates from the output above:

Notice how the eclipse of 9 Mar 2035 has g = -0.4388. The negative value tells us the eclipse is in the Southern Hemisphere. And the total eclipse on 2 Sep 2035 has g = 0.3721 so it’s in the Northern Hemisphere.

Improvements and Next Steps

For this project I went out to timeanddate.com and grabbed all the new moon dates and put them in a data structure. That’s just fine for the sake of this blog post. But it would be an improvement to calculate the moon phases instead. Khalid Abuhahmeh did just that, so his solution is worth a look. Another good solution written in PHP is here.

Clear skies!