Derivation of initial bearing calculation

Question

NOTE: I was recommended by another user that this would perhaps fit better in this stack, initial post link.

I needed to make calculations of the initial bearing between two coordinates described in latitude/longitude. I have found one formulation of this calculation as

Note, I tried formatting the equation with LaTeX but I couldn't see how to do that on stack, using $$ didn't seem to do anything.

This calculation can be found on various sites like link 1 2 3 but I have not been able to find a derivation or a source for the calculation, as such I am looking for a derivation or an explanation of the derivation.

Answer 1

I was surprised that I didn't find this ready-made question and answer on this site.

This formula is sometimes taken from http://www.movable-type.co.uk/scripts/latlong.html, which in turn is taken from http://www.edwilliams.org/avform.htm, but in neither of the two places is this derivation really explained.

Where it is done is in Link, which I limit myself to copying because it is self explanatory:

Let's define three unit vectors, each in the direction of the line from the center of the earth to a point on the surface: N in the direction of the north pole, A in the direction of the initial point, and B in the direction of the final point on the course. Then the bearing we seek is the angle between the plane containing N and A, and the plane containing A and B. Thus it equals the angle between vectors perpendicular to these planes, namely, NxA and BxA.

Let point A have latitude lat1 and longitude 0 (we can rotate our coordinate system so this is true), and let point B have latitude lat2 and longitude dlon (the difference between the actual longitudes of A and B). Then we can calculate NxA and BxA:

 N = (0, 0, 1)
 A = (cos(lat1), 0, sin(lat1))
 B = (cos(lat2)*cos(dlon), cos(lat2)*sin(dlon), sin(lat2))

 NxA = (0, cos(lat1), 0)
 BxA = (sin(lat1)*cos(lat2)*sin(dlon),
        cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(dlon),
        -cos(lat1)*cos(lat2)*sin(dlon))

The usual way to find the angle would be to take the dot product of these vectors, divide by the product of the magnitudes of the vectors, and take the arccosine. But that would be quite a mess. Instead we can take advantage of the fact that NxA is parallel to the y axis. The tangent of the angle between BxA and the y axis is the component of BxA in the x-z plane (the square root of the sum of the squares of the x and z components) divided by the y component:

 tan(theta) = sqrt((sin(lat1)*cos(lat2)*sin(dlon))^2 + 
                   (-cos(lat1)*cos(lat2)*sin(dlon))^2) / 
              (cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(dlon))

The numerator simplifies to

 = cos(lat2)*sin(dlon)

Thus

 tan(theta) = cos(lat2)*sin(dlon) / 
              (cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(dlon))

That's the formula. If we take the arctan of both sides, we get a value in the range from -pi/2 to pi/2 radians. The problem is, this doesn't distinguish between opposite directions, NE vs. SW for instance. But the function atan2(y,x) returns the arctan of y/x with adjustments for the signs of x and y so that the angle returned is the angle of the cartesian point (x,y) in polar coordinates -- just what we need here.

The only problem is that it returns a value in the range (-pi, pi]. The function mod(..., 2pi) moves negative angles up to the range (pi, 2pi). To get a final answer in degrees (0 to 360), you must multiply by 180/pi.