James Still

In my last post (Calculate Future Solar Eclipses) I showed how to calculate the types of solar eclipses that would occur for a given date range. I took a shortcut and used a NewMoonData class that had a list of all known new moon dates that I culled from the timeanddate.com website. That worked just fine for the purposes of that blog post.

But in the real world who wants to hardcode every single new moon date into the far future? That’s not a realistic solution. I’m following up with an implementation in C# .NET, based on Jean Meeus’ algorithm he adopted from the ELP-2000/82 theory for the Moon. The results show the exact instant of the new moon phase (UTC) for any date from -1999 BCE to +3000 CE, accurate to within a few seconds.

The code is checked into my repo at FindNewMoonConsoleApp. Feel free to clone the code base and take it for a test ride. I won’t do a full walkthrough like with the solar eclipse calculations. I’ll just say that the general approach I took was to start with the first new moon in the J2000.0 epoch (6 Jan 2000) and enumerate over the mean new moon phases that follow. Starting with 6 Jan 2000 the program adds the lunar cycle of 29.5306 days to get 4 Feb 2000. Add the lunar cycle again and you get 5 Mar 2000, and so on for as many months and years as you like.

Sounds easy right? The problem is there are forces of gravity pulling on the Moon and Earth. And none of these bodies move at constant rates. So, the times of the mean phases of the moon just gets us in the ballpark of give or take a day or so. From there we have to apply an algorithm that corrects for all of the perturbations in order to get a more accurate date and time for each new moon.

Suppose I want to know the instant of the new moon in Dec 2022. I run the program and get: 12/23/2022 10:16:47 UTC. The Griffith Observatory, which uses the NASA/JPL Horizons system, rounds it up to 12/23/2022 10:17 UTC.

The program will produce results for every new moon that are accurate to within several seconds for any date range you provide.

Clear skies!

In two-dimensional space (a plane), the shortest distance between two points is a line. The shortest distance between two points on a sphere is the great-circle distance. The shorter arc described by the great circle running through the two points is then calculated by this well known formula below. Given the geographical longitude and latitude \(L_1, \phi_1\) of point P and \(L_2, \phi_2\) of point Q:

$$\cos d = \sin \phi_1 \space \sin \phi_2 + \cos \phi_1 \space \cos \phi_2 \space \cos (L_1 – L_2)$$

Where d is expressed in degrees:

$$d = \arccos (\sin \phi_1 \space \sin \phi_2 + \cos \phi_1 \space \cos \phi_2 \space \cos (L_1 – L_2))$$

I once worked at a logistics company primarily as a database developer. I was asked to calculate distances between various ports in China and the U.S. in order to estimate fuel costs. I implemented the above great-circle algorithm in an Oracle PL/SQL package and it worked just fine. However, the Earth isn’t a sphere. Because of centrifugal force generated from the Earth’s rotation it bulges at the equator. It’s better to say that the Earth is a spheroid or an ellipsoid.

So while the great-circle algorithm is good enough what if we want even greater accuracy? To do that we need to take into account the Earth’s true shape as an ellipsoid. There are two variables we need. We need to know the equatorial radius of the Earth (a) and we need to know its flattening (f). Consider this ellipse that represents the Earth with an exaggerated equatorial bulge:

The equatorial radius AC (or BC) is labeled a and flattening f is the ratio of \(a – b \over a\). These values are constants accepted by the International Astronomical Union (IAU) and I can pull them from the wiki:

a = 6378.140 km

f = 0.00335281

Given these two constants we can now calculate the geodesic distance between two points on the Earth’s surface using Andoyer’s Formula. Here I’ve created the program in C# following Meeus’ Method as described in his book Astronomical Algorithms (p. 85). This method will produce an error of +/- 50 meters.

First there are two helper classes. DoubleExtensions has a method ToRadians that will be necessary for the trigonometric functions:

public static class DoubleExtensions
{
    public static double ToRadians(this double d)
    {
        return d * (Math.PI / 180); // convert degrees to radians
    }

    public static double ToDegrees(this double d)
    {
        return d * (180 / Math.PI); // convert radians to degrees
    }
}

And I created a small struct to wrap the geocoordinates of the two points between which we want to calculate the distance:

public struct GeoCoordinate
{
    public GeoCoordinate(double latitude, double longitude)
    {
        Latitude = latitude;
        Longitude = longitude;
        Name = "Unknown";
    }

    public GeoCoordinate(double latitude, double longitude, string name)
    {
        Latitude = latitude;
        Longitude = longitude;
        Name = name;
    }

    public double Latitude { get; }
    public double Longitude { get; }
    public string Name { get; }
}

With those two helpers in place the main program will calculate the distance between the Observatoire de Paris and the U.S. Naval Observatory in Washington, D.C.

class Program
{
    static void Main(string[] args)
    {
        // constants
        const double a = 6378.140;
        const double f = 0.00335281;

        var p = new GeoCoordinate(
            48.83698, -2.33652, "Observatoire de Paris"); 
        var q = new GeoCoordinate(
            38.921261, 77.06652, "U.S. Naval Observatory");

        double F = (p.Latitude + q.Latitude) / 2;
        double G = (p.Latitude - q.Latitude) / 2;
        double L = (p.Longitude - q.Longitude) / 2;

        double S = Math.Pow(Math.Sin(G.ToRadians()), 2) 
            * Math.Pow(Math.Cos(L.ToRadians()), 2) 
            + Math.Pow(Math.Cos(F.ToRadians()), 2) 
            * Math.Pow(Math.Sin(L.ToRadians()), 2);

        double C = Math.Pow(Math.Cos(G.ToRadians()), 2) 
            * Math.Pow(Math.Cos(L.ToRadians()), 2) 
            + Math.Pow(Math.Sin(F.ToRadians()), 2) 
            * Math.Pow(Math.Sin(L.ToRadians()), 2);

        double O = Math.Atan(Math.Sqrt(S / C));
        double R = Math.Sqrt(S * C) / O;
        double D = 2 * O * a;
        double H1 = (3 * R - 1) / (2 * C);
        double H2 = (3 * R + 1) / (2 * S);

        double s = D * (1 + f * H1 
            * Math.Pow(Math.Sin(F.ToRadians()), 2) 
            * Math.Pow(Math.Cos(G.ToRadians()), 2) 
            - f * H2 
            * Math.Pow(Math.Cos(F.ToRadians()), 2) 
            * Math.Pow(Math.Sin(G.ToRadians()), 2)
            );

        string d = Math.Round(s, 2).ToString();

        Console.WriteLine($"The distance between 
            {p.Name} and {q.Name} is {d} km (+/- 50 m).");
        Console.ReadLine();
    }
}

The output produces the correct result:

The distance between Observatoire de Paris and U.S. Naval Observatory is 6181.63 km (+/- 50 m).