Law of cosines on a sphere

The previous post looked at the analog of the Pythagorean theorem on a sphere. This post looks at the law of cosines on a sphere.

Yesterday we looked at triangles on a sphere with sides a and b meeting at a right angle and hypotenuse c. Denote the angle opposite a side with the capital version of the same letter, so C is the angle opposite c. We assumed C is a right angle, but now we will remove that assumption.

The spherical analog of the law of cosines says

cos(c) = cos(a) cos(b) + sin(a) sin(b) cos(C).

Note that we have two kinds of angles: the arcs that make up the sides, and the angle formed by the intersection of arcs [1]. So the cos(C) term at the end is different animal from the other terms in the equation.

If C is a right angle, cos(C) = 0 and the second term drops out, leaving us with the spherical counterpart of the Pythagorean theorem. But we do not require that C be a right angle.

Application to air distance

Suppose we want to find how far a plane would travel between Los Angeles (LAX) and Houston (IAH), assuming it takes a great circle path. The lat/long coordinates of the two airports are (33.94°, -118.41°) and (29.98°, -95.34°).

Los Angeles and Houston are approximately at the same latitude, but even if they were at exactly the same latitude we couldn’t just find the flight distance by finding the length of the arc of constant latitude between them because that arc would not be part of a great circle.

To find the distance between LAX and IAH we form a triangle with vertices at LAX, IAH, and the north pole. Call the arc from LAX to the north pole a and the arc from IAH to the north pole b. Since latitude is the angle up from the equator, the angle down from the pole is 90° minus latitude. So

a = 90° – 33.94° = 56.06°

and

b = 90° – 29.98° = 60.02°

The arcs a and b meet at the angle C equal to the differences in the two longitudes. That is,

C = 118.41° – 95.34° = 23.07°.

The law of cosines

cos(c) = cos(a) cos(b) + sin(a) sin(b) cos(C)

reduces to

cos(c) = cos(56.06°) cos(60.02°) + sin(56.06°) sin(60.02°) cos(23.07°)

in our problem.

This tells us

cos(c) = 0.9402

from which we find c = 0.3477 radians or 19.92°. Assume the earth is a perfect sphere of radius r = 3959 miles. The arc length we’re after is r times c in radians, or 1376 miles.

Related posts

[1] In the language of differential geometry, the arcs are geodesics on the sphere and the angles of intersection are measured in the tangent space at the intersection point.

Pythagoras on a sphere

Suppose you drive a distance a, turn right, then drive a distance b. How far are you from where you started?

If a and b are relatively small, then the answer is given by the Pythagorean theorem. Small here means small relative to the radius of the Earth. Any distance you drive is small enough: it matters more that roads are not perfectly straight than that the Earth is not perfectly flat.

Now suppose you fly a distance a, turn right, then fly a distance b. Now how far are you from where you started? Now the curvature of the earth might matter a good deal, and the Pythagorean theorem could give you a poor answer.

Spherical version of the Pythagorean theorem

If a and b are moderately large, we need to look at spherical triangles, not plane triangles. A “straight” line on a sphere is a piece of a great circle, a circle whose center is the center of the sphere. A spherical triangle has three straight sides, straight in the sense of a straight line on a sphere. We need to consider a spherical triangle with sides a and b that meet at right angles, and we want to find the length of the hypotenuse c of the spherical triangle.

The spherical analog of the Pythagorean theorem says

cos(c) = cos(a) cos(b).

This looks weird the first time you see it, because you’re taking the cosine of a side rather than an angle. But on a sphere, the sides are angles. For example, if a side of the triangle is on a meridian, we could measure its length in degrees of latitude. But the theorem above does not require sides to be lined up with any coordinate system. It only requires that the sides a and b meet at a right angle.

Small spherical triangles

Let’s suppose for a moment that a and b are small when measured in radians, small enough that a², b², and ab are negligible compared to a and b. Then we can use the approximation

cos(x) ≈ 1 – x²/2

in the theorem above to show that

c² ≈ a² + b².

This is reassuring since we know the Pythagorean theorem holds for human-scale distances on our planet.

Example with a big triangle

Quito to Nairobi to Jerusalem and back to Quito

Now let’s look at a big triangle. The vertices have longitude and latitude (-80°, 0°), (35°, 0°), and (35°, 30°). These correspond approximately to Quito, Nairobi, and Jerusalem. So a = 115° and b = 30°. We can solve for c and find it equals 111°.

A minute of latitude is a nautical mile; the nautical mile is defined the way to make these kinds of calculations easier. A nautical mile is about 15% longer than a mile.

A degree of a great circle is 60 nautical miles. So our hypotenuse of 111 degrees is 6,660 nautical miles. If we had used the Pythagorean theorem, we would have overestimated this distance as 7,130 nautical miles.

Notice that the triangle with vertices at Quito, Nairobi, and Jerusalem doesn’t look like a triangle when projected flat. Two legs look straight because they follow lines of constant latitude or longitude. But the flight from Jerusalem to Quito is just as straight as the other flights if you look at it on a globe.

Not so great circles

In the example above, I chose Quito and Nairobi because they were near the equator. The equator is a great circle.

Suppose I’d chosen Houston and Wuhan instead. Both cities are at about 30° latitude. I could not find the distance between the two cities by simply subtracting their longitudes because the 30th parallel is not a great circle. Spherical geometry is generally more involved than in the example above.

A great circle path from Houston to Wuhan would not follow the 30th parallel but would actually cross Alaska. The great circle distance is 6,800 nautical miles. Flying along a path of 30° latitude would add an extra 1,000 nautical miles to the flight.

Related posts

A tale of three cities

Big blue marble

Pick three cities and form a spherical triangle by connecting each pair of cities with the shortest arc between them. How might you find the area of this triangle?

For this post, I’ll assume the earth is perfectly spherical. (Taking into account that the earth is slightly oblate makes the problem much more complicated. Maybe a post for another time.)

Spherical excess

If you draw a triangle on the plane, the interior angles add up to 180°. But if you draw a triangle on a sphere, the interior angles add up to more than 180°. The amount by which the sum of the interior angles exceeds 180° is called the spherical excess. It turns out that the area of a spherical triangle is proportional to its spherical excess.

For example, if a spherical triangle is small relative to the radius of the sphere, it’s nearly a Euclidean triangle, and the spherical excess is negligible. But consider a much bigger spherical triangle, one with a vertex at the north pole and vertices at two points 90° apart on the equator. This triangle has three right angles, so the sum of the interior angles is 270° and the spherical excess is 90°.

Relating spherical excess and area

If we measure spherical excess E in radians, then the area of a spherical triangle is

A = ER²

where R is the radius of the sphere. This remarkable theorem was discovered by Thomas Harriot in 1603.

Let’s go back to the example above of a triangle running from the north pole to two points 90° apart on the equator. This triangle takes up 1/8 of the earth’s surface, and area of the earth is 4πR², and so the triangle has area πR²/2. In this case the spherical excess E was π/2, and so we could have come to the same result by multiplying E by R².

Computing area

We can find the area of a spherical triangle by measuring the angle of arcs between pairs of points. We’ll compute area assuming the sphere has radius 1. (Otherwise temporarily assume radius 1 and then multiply by R² when we’re done.)

The previous post explains how to find a parameterization for the great circle connecting points on a sphere. We can take derivatives to find tangent vectors where the great circles intersect, and compute the angle between these tangents to find the interior angles of our triangles.

The previous post showed that if the central angle between two vectors A and B is θ then

cos(t) A + sin(t) (cot θ A – csc θ B)

parameterizes a great circle including A and B. This curve passes through A at time t = 0 with velocity

cot θ A – csc θ B

and so this gives us a tangent vector at A in the direction of B. Repeat the exercise for a great circle between A and C. Then the cosine of the angle of intersection is the two normalized tangent vectors. We can thus obtain all the interior angles, and from there we can compute the spherical excess and thus the area of the spherical triangle.

Related posts

Maidenhead geocode system

The Maidenhead Locator System encodes a pair of longitude and latitude coordinates in a slightly complicated but ingenious way. Amateur radio operators using this geocoding system to describe locations.

The Wikipedia article on the subject describes the what of the system, but I’d like to say more about the why of the system. I’ll also go through an example in great detail. Continue reading

Spherical trig, Research Triangle, and Mathematica

This post will look at the triangle behind North Carolina’s Research Triangle using Mathematica’s geographic functions.

Spherical triangles

A spherical triangle is a triangle drawn on the surface of a sphere. It has three vertices, given by points on the sphere, and three sides. The sides of the triangle are portions of great circles running between two vertices. A great circle is a circle of maximum radius, a circle with the same center as the sphere.

An interesting aspect of spherical geometry is that both the sides and angles of a spherical triangle are angles. Because the sides of a spherical triangle are arcs, they have angular measure, the angle formed by connecting each vertex to the center of the sphere. The arc length of a side is its angular measure times the radius of the sphere.

Denote the three vertices by AB, and C. Denote the side opposite A by a, etc. Denote the angles at AB, and C by α, β, and γ respectively.

Research Triangle

Research Triangle is a (spherical!) triangle with vertices formed by Duke University, North Carolina State University, and University of North Carolina at Chapel Hill.

(That was the origin of the name, though it’s now applied more loosely to the general area around these three universities.)

We’ll take as our vertices

  • A = UNC Chapel Hill (35.9046 N, 79.0468 W)
  • B = Duke University in Durham (36.0011 N, 78.9389 W),
  • C = NCSU in Raleigh (35.7872 N, 78.6705 W)

Research Triangle

Mathematica

We’ll illustrate several features of Mathematica using the spherical triangle corresponding to Research Triangle.

Map

The map above was produced with the following Mathematica code.

    ptA = GeoPosition[{35.9046, -79.0468}]
    ptB = GeoPosition[{36.0011, -78.9389}]
    ptC = GeoPosition[{35.7872, -78.6705}]

    GeoGraphics[{Red, PointSize[Large], 
        Point[ptA], Point[ptB], Point[ptC]}, 
        GeoScaleBar -> "Imperial", 
        GeoRange -> 29000]

Note that longitude is in degrees east, so the longitudes above are negative.

Distance

To find the distance between two locations, you can use the function GeoDistance. For example,

    GeoDistance[ptA, ptB]

tells me that the distance between UNC and Duke is 8.99185 miles. I assume it displays miles by default based on my location in the US, though the GeoRange option above defaults to meters. You can make the system of units explicit. For example

    GeoDistance[ptA, ptB, UnitSystem -> "Metric"]

returns 14.741 km.

If we want to find the length of the side c in degrees, we can use the Earth’s radius.

    r = PlanetData["Earth", "Radius"]
    c = GeoDistance[ptA, ptB] 360 / (2 Pi r)

This shows that c is 0.13014°. Calculating the other sides similarly shows a = 0.30504° and b = 0.32739°.

Angles

Calling GeoDirection[ptA, ptB] returns 42.2432°, which says we need to head at an angle of about 42° to walk from UNC to Duke.

The code

    GeoDirection[ptA, ptB] - GeoDirection[ptA, ptC]

shows that the angle α is 68.6128°. (The code returns the negative of this angle because the angle is clockwise.) Similarly we find β = 87.9808° and γ = 23.4068.

The angles add up to 180° only because our triangle is small compared to the earth’s total area. The actual sum should be slightly more than 180°, but we’ve not retained enough precision to detect the difference. In general the “spherical excess,” i.e. the amount by which the sum of the angles exceed 180°, is proportional to the area of the triangle.

More spherical trig posts

Ellipsoid distance on Earth

globe

To first approximation, Earth is a sphere. But it bulges at the equator, and to second approximation, Earth is an oblate spheroid. Earth is not exactly an oblate spheroid either, but the error in the oblate spheroid model is about 100x smaller than the error in the spherical model.

Finding the distance between two points on a sphere is fairly simple. Here’s a calculator to compute the distance, and here’s a derivation of the formula used in the calculator.

Finding the distance between two points on an ellipsoid is much more complicated. (A spheroid is a kind of ellipsoid.) Wikipedia gives a description of Vincenty’s algorithm for finding the distance between two points on Earth using an oblate spheroid model (specifically WGS-84). I’ll include a Python implementation below.

Comparison with spherical distance

How much difference does it make when you calculate difference on an oblate spheroid rather than a sphere? To address that question I looked at the coordinates of several cities around the world using the CityData function in Mathematica. Latitude is in degrees north of the equator and longitude is in degrees east of the prime meridian.

    |--------------+--------+---------|
    | City         |    Lat |    Long |
    |--------------+--------+---------|
    | Houston      |  29.78 |  -95.39 |
    | Caracas      |  10.54 |  -66.93 |
    | London       |  51.50 |   -0.12 |
    | Tokyo        |  35.67 |  139.77 |
    | Delhi        |  28.67 |   77.21 |
    | Honolulu     |  21.31 | -157.83 |
    | Sao Paulo    | -23.53 |  -46.63 |
    | New York     |  40.66 |  -73.94 |
    | Los Angeles  |  34.02 | -118.41 |
    | Cape Town    | -33.93 |   18.46 |
    | Sydney       | -33.87 |  151.21 |
    | Tromsø       |  69.66 |   18.94 |
    | Singapore    |   1.30 |  103.85 |
    |--------------+--------+---------|

Here are the error extremes.

The spherical model underestimates the distance from London to Tokyo by 12.88 km, and it overestimates the distance from London to Cape Town by 45.40 km.

The relative error is most negative for London to New York (-0.157%) and most positive for Tokyo to Sidney (0.545%).

Update: The comparison above used the same radius for both the spherical and ellipsoidal models. As suggested in the comments, a better comparison would use the mean radius (average of equatorial and polar radii) in the spherical model rather than the equatorial radius.

With that change the most negative absolute error is between Tokyo and New York at -30,038 m. The most negative relative error is between London and New York at -0.324%.

The largest positive error, absolute and relative, is between Tokyo and Sydney. The absolute error is 29,289 m and the relative error is 0.376%.

Python implementation

The code below is a direct implementation of the equations in the Wikipedia article.

Note that longitude and latitude below are assumed to be in radians. You can convert from degrees to radians with SciPy’s deg2rad function.

from scipy import sin, cos, tan, arctan, arctan2, arccos, pi

def spherical_distance(lat1, long1, lat2, long2):
    phi1 = 0.5*pi - lat1
    phi2 = 0.5*pi - lat2
    r = 0.5*(6378137 + 6356752) # mean radius in meters
    t = sin(phi1)*sin(phi2)*cos(long1-long2) + cos(phi1)*cos(phi2)
    return r * arccos(t)

def ellipsoidal_distance(lat1, long1, lat2, long2):

    a = 6378137.0 # equatorial radius in meters 
    f = 1/298.257223563 # ellipsoid flattening 
    b = (1 - f)*a 
    tolerance = 1e-11 # to stop iteration

    phi1, phi2 = lat1, lat2
    U1 = arctan((1-f)*tan(phi1))
    U2 = arctan((1-f)*tan(phi2))
    L1, L2 = long1, long2
    L = L2 - L1

    lambda_old = L + 0

    while True:
    
        t = (cos(U2)*sin(lambda_old))**2
        t += (cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(lambda_old))**2
        sin_sigma = t**0.5
        cos_sigma = sin(U1)*sin(U2) + cos(U1)*cos(U2)*cos(lambda_old)
        sigma = arctan2(sin_sigma, cos_sigma) 
    
        sin_alpha = cos(U1)*cos(U2)*sin(lambda_old) / sin_sigma
        cos_sq_alpha = 1 - sin_alpha**2
        cos_2sigma_m = cos_sigma - 2*sin(U1)*sin(U2)/cos_sq_alpha
        C = f*cos_sq_alpha*(4 + f*(4-3*cos_sq_alpha))/16
    
        t = sigma + C*sin_sigma*(cos_2sigma_m + C*cos_sigma*(-1 + 2*cos_2sigma_m**2))
        lambda_new = L + (1 - C)*f*sin_alpha*t
        if abs(lambda_new - lambda_old) <= tolerance:
            break
        else:
            lambda_old = lambda_new

    u2 = cos_sq_alpha*((a**2 - b**2)/b**2)
    A = 1 + (u2/16384)*(4096 + u2*(-768+u2*(320 - 175*u2)))
    B = (u2/1024)*(256 + u2*(-128 + u2*(74 - 47*u2)))
    t = cos_2sigma_m + 0.25*B*(cos_sigma*(-1 + 2*cos_2sigma_m**2))
    t -= (B/6)*cos_2sigma_m*(-3 + 4*sin_sigma**2)*(-3 + 4*cos_2sigma_m**2)
    delta_sigma = B * sin_sigma * t
    s = b*A*(sigma - delta_sigma)

    return s

Projecting the globe onto regular solids

I was playing around with some geographic features of Mathematica this morning and ran across an interesting example in the documentation for the PolyhedronProjection function given here. Here’s what you get when you project a map of the earth onto each of the five regular (Platonic) solids.

dodecahedron

icosahedron
octahedron
cube
tetrahedron

How the images were made

At first I right-clicked on each image and saved as graphic files. This produced low quality images, even when I saved as SVG. I got better resolution by using the Export command and specifying the ImageSize and ImageResolution options.

The default view of the tetrahedron was face-on, so it looked like a flat triangle. By changing the ViewPoint I was able to get something that’s more clearly three dimensional.

By the way, you can use PolyhedronProjection to project onto more exotic polyhedra than the Platonic solids. For example,

    Export["rhomb.png", 
        PolyhedronProjection["ParagyrateDiminishedRhombicosidodecahedron"], 
        ImageResolution -> 72, 
        ImageSize -> Large]

produces this:

ParagyrateDiminishedRhombicosidodecahedron

More regular solid posts

Equal Earth map projection

There’s no perfectly satisfying way to map the globe on to a flat surface. Every projection has its advantages and disadvantages. The Mercator projection, for example, is much maligned for the way it distorts area, but it has the property that lines of constant bearing correspond to straight lines on the map. Obviously this is convenient if you’re sailing without GPS. But for contemporary use, say in a classroom, minimizing area distortion is often a higher priority than keeping bearing lines straight.

Bojan Šavrič, Tom Patterson, and Bernhard Jenny have developed a new map projection called Equal Earth that nicely balances several competing criteria, including aesthetics.

Equal earth projection

The Equal Earth projection satisfies the following mathematical criteria:

  1. Equal area: The area of a region on the globe is proportional to the area of its projection.
  2. Straight parallels: Lines of equal latitude project on to horizontal lines.
  3. Regularly distributed meridians: At a given latitude, the spacing between lines of longitude are equal.
  4. Bilateral symmetry: The projection is symmetric with respect to the x-axis (equator) and y-axis (prime meridian).

Here are the equations for the Equal Earth projection.

\begin{eqnarray*} \sin\theta &=& \frac{\sqrt{3}}{2} \sin \phi \\ x &=& \frac{2\sqrt{3}}{3} \frac{\lambda \cos \theta}{a_1+ 3a_2 \theta^2 + \theta^6(7a_3 + 9a_4\theta^2)} \\ y &=& \theta (a_1 + a_2\theta^2 + \theta^6(a_3 + a_4\theta^2)) \end{eqnarray*}

Here λ and φ are longitude and latitude respectively. The parametric latitude θ is introduced for convenience and is simply a rescaling of latitude φ.

The parameters ai are given below.

    a_1 =  1.340264 
    a_2 = -0.081106 
    a_3 =  0.000893 
    a_4 =  0.003796 

You can find more in their paper: Bojan Šavrič, Tom Patterson, and Bernhard Jenny. The Equal Earth map projection. International Journal of Geographical Information Science. https://doi.org/10.1080/13658816.2018.1504949.

More geography posts

Misplacing a continent

There are many conventions for describing points on a sphere. For example, does latitude zero refer to the North Pole or the equator? Mathematicians tend to prefer the former and geoscientists the latter. There are also varying conventions for longitude.

Volker Michel describes this clash of conventions colorfully in his book on constructive approximation.

Many mathematicians have faced weird jigsaw puzzles with misplaced continents after using a data set from a geoscientist. If you ever get such figures, too, or if you are, for example, desperately searching South America in a data set but cannot find it, remember the remark you have just read to solve your problem.

More geography posts

Ellipsoid surface area

How much difference does the earth’s equatorial bulge make in its surface area?

To first approximation, the earth is a sphere. The next step in sophistication is to model the earth as an ellipsoid.

The surface area of an ellipsoid with semi-axes abc is

A = 2\pi \left( c^2 + \frac{ab}{\sin\phi} \left( E(\phi, k) \sin^2\phi + F(\phi, k) \cos^2 \phi\right)\right)

where

\cos \phi = \frac{c}{a}

and

m = k^2 = \frac{a^2(b^2 - c^2)}{b^2(a^2 - c^2)}

The functions E and F are incomplete elliptic integrals

 F(\phi, k) = \int_0^\phi \frac{d\theta}{\sqrt{1 - k^2 \sin^2\theta}}

and

E(\phi, k) = \int_0^\phi \sqrt{1 - k^2 \sin^2\theta}\,d\theta

implemented in SciPy as ellipeinc and ellipkinc. Note that the SciPy functions take m as their second argument rather its square root k.

For the earth, a = b and so m = 1.

The following Python code computes the ratio of earth’s surface area as an ellipsoid to its area as a sphere.

from numpy import pi, sin, cos, arccos
from scipy.special import ellipkinc, ellipeinc

# values in meters based on GRS 80
# https://en.wikipedia.org/wiki/GRS_80
equatorial_radius = 6378137
polar_radius = 6356752.314140347

a = b = equatorial_radius
c = polar_radius

phi = arccos(c/a)
# in general, m = (a**2 * (b**2 - c**2)) / (b**2 * (a**2 - c**2))
m = 1

temp = ellipeinc(phi, m)*sin(phi)**2 + ellipkinc(phi, m)*cos(phi)**2
ellipsoid_area = 2*pi*(c**2 + a*b*temp/sin(phi))

# sphere with radius equal to average of polar and equatorial
r = 0.5*(a+c)
sphere_area = 4*pi*r**2

print(ellipsoid_area/sphere_area)

This shows that the ellipsoid model leads to 0.112% more surface area relative to a sphere.

Source: See equation 19.33.2 here.

Update: It was suggested in the comments that it would be better to compare the ellipsoid area to that of a sphere of the same volume. So instead of using the average of the polar and equatorial radii, one would take the geometric mean of the polar radius and two copies of the equatorial radius. Using that radius, the ellipsoid has 0.0002% more area than the sphere.

Update 2: This post gives a simple but accurate approximation for the area of an ellipsoid.