Month: April 2020

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).