AzimuthModules
(Determining the direction between two points on Earth's surface)

Coded in Mathematica v. 12.2.0.0 by Vince Cronin
Revised February 28, 2021

Introduction

The purpose of this Mathematica notebook is to explain how to compute the direction (azimuth or bearing) from a reference point to another point on Earth’s surface.  Azimuth is measured from true north in a clockwise direction, so the azimuth of due east is 90°, south is 180°, west is 270°, and north is 0° (or 360°).

Assumptions

We assume that Earth is a sphere of unit radius.

Specify the input data

Input only decimal degrees (i.e., not degrees, minutes, seconds).  The sign convention follows:  north latitude is positive, south latitude is negative, east longitude is positive, west longitude is negative.  
The variables refPtLat and refPtLong are the latitude and longitude of the reference point from which the distance and azimuth to the observed point is to be determined.

refPtLat=35;

refPtLong=-119;

The variables otherPtLat and otherPtLong are the latitude and longitude of the other point.

otherPtLat=37;

otherPtLong=-121;

Operate on the input data

User-defined function

The function geog2Cartesian computes the unit location vector of a given point on an assumed-spherical Earth, given that point's latitude and longitude.  The X coordinate is equal to Cos[latitude]*Cos[longitude].  The Y coordinate is equal to Cos[latitude]*Sin[longitude].  The Z coordinate is equal to Sin[latitude].  The input latitude and longitude values must be in decimal degrees (i.e., not degrees, minutes, seconds).  The forms "lat Degree" and "long Degree" in the square brackets invokes the built-in Mathematica function Degree to convert the input data from degrees to radians, as required for the built-in Mathematica trigonometric functions (e.g., Sin, Cos, Tan, ArcSin, ArcCos, ArcTan).

geog2Cartesian[lat_,long_]:={Cos[lat Degree] Cos[long Degree],Cos[lat Degree] Sin[long Degree],Sin[lat Degree]};

Computation

Utilize the user-defined fuction makeVector to convert the input latitude and longitude data into unit location vectors

refPtVect=geog2Cartesian[refPtLat,refPtLong];

otherPtVect=geog2Cartesian[otherPtLat,otherPtLong];

Define the unit location vector associated with Earth' s north pole (northPole)

northPole={0,0,1};

Find the vector that is normal to the plane that contains the given two vectors, using a vector cross product.  The built-in Mathematica function Cross performs this computation.  If you had two vectors α and β, respectively specified by components {(Subscript[α, x]),  (Subscript[α, y]), (Subscript[α, z])} and {(Subscript[β, x]), (Subscript[β, y]), (Subscript[β, z])}, the result of the vector cross product α × β would be a vector γ with components {(Subscript[γ, x]), (Subscript[γ, y]), (Subscript[γ, z])} = {(((Subscript[α, y]) (Subscript[β, z]))-((Subscript[α, z]) (Subscript[β, y]))), (((Subscript[α, z]) (Subscript[β, x]))-((Subscript[α, x]) (Subscript[β, z]))), (((Subscript[α, x]) (Subscript[β, y]))-((Subscript[α, y]) (Subscript[β, x])))} that is perpendicular to the the plane containing α and β.  Given α × β = γ, vector β is a right-handed rotation around γ starting at α.

vectorA = Cross[northPole,refPtVect];

vectorB = Cross[otherPtVect,refPtVect];

Utilize the built-in Mathematica function VectorAngle to find the angle θ1  (expressed in degrees) between the two vectors vectorA and vectorB, which is also the angular difference between [1] the azimuth from the reference point toward north, and [2] the azimuth from the reference point toward the other point.

θ1=VectorAngle[vectorA,vectorB](180/π);

Find the angular distance (θ2), in degrees, between the location vector to the other point (otherPtVect) and the vectorA.  This angle helps us define the bearing from the reference point to the other point.

θ2=VectorAngle[vectorA,otherPtVect](180/π);

Imagine that Earth is divided into two hemispheres bounded by the plane that contains Earth's spin axis and the reference point.  We will specify that the north pole is "up" in this frame of reference, so one hemisphere would be to the right of the reference point and the other would be to the left.
     If the other point is located along the same circle of longitude as the reference point, then θ2 will equal 90° or 180°, and the bearing from the reference point to the other point will either be 0° or 180°.
     If the other point is located on the hemisphere to the right of the reference point, then θ2 will be less than 90°.  In that case, the azimuth or bearing from the reference point to the other point is equal to θ1.  Bearings are expressed in degrees relative to north, where due north is 0°, due east is 90°, due south is 180°, and due west is 270°.
     If the other point is located on the hemisphere to the left of the reference point, then θ2 will be greater than 90°.  In that case, the bearing from the reference point to the other point is equal to (360°- θ1).

The vector vectorC and the angle θ3 are used when the other point is on the same meridian as the reference point, to differentiate between a bearing of 0° (north) or 180° (south) from the reference point to the other point.  If θ3 is less than 90°,  the other point is north of the reference point.  If θ3 is equal to 90°,  the other point is either coincident with the reference point or 180° from the reference point, in which case the bearing from the reference point to the other point is undefined.  For simplicity, when θ3=90°, we set the bearing to 0°.

vectorC=Cross[refPtVect,vectorA];

θ3=VectorAngle[vectorC,otherPtVect](180/π);

The code fragment "If[((θ2 == 90) ∨ (θ2 == 180))" in the following line  literally means "if θ2 equals 90 or if θ2 equals 180."

azimuth1=If[((θ2==90)∨(θ2==180)),If[(θ3<=90),0,180],If[(θ2>90),(360-θ1),θ1]];

azimuth2=N[azimuth1];

ClearAll[northPole,vectorA,vectorB,θ1,θ2,vectorC,θ3]

Modules

The module findAzimuthVectInput[]) acts on unit location vectors to determine the azimuth from a reference point (location vector refPtVect) to a target point (location vector otherPtVect).

findAzimuthVectInput[refPtVect_,otherPtVect_]:=Module[{northPole,vectorA,vectorB,θ1,θ2,vectorC,θ3,azimuth},northPole={0,0,1};vectorA = Cross[northPole,refPtVect];vectorB = Cross[otherPtVect,refPtVect];θ1=VectorAngle[vectorA,vectorB](180/π);θ2=VectorAngle[vectorA,otherPtVect](180/π);vectorC=Cross[refPtVect,vectorA];θ3=VectorAngle[vectorC,otherPtVect](180/π);azimuth=If[((θ2==90)∨(θ2==180)),If[(θ3<=90),0,180],If[(θ2>90),(360-θ1),θ1]];azimuth ];

The module findAzimuthGeogInput[]) acts on two sets of geographic coordinates to determine the azimuth from a reference point (refPtLat and refPtLong) to a target point (otherPtLat and otherPtLong).

findAzimuthGeogInput[refPtLat_,refPtLong_,otherPtLat_,otherPtLong_]:=Module[{refPtVect,otherPtVect,northPole,vectorA,vectorB,θ1,θ2,vectorC,θ3,azimuth},refPtVect={Cos[refPtLat Degree] Cos[refPtLong Degree],Cos[refPtLat Degree] Sin[refPtLong Degree],Sin[refPtLat Degree]};otherPtVect={Cos[otherPtLat Degree] Cos[otherPtLong Degree],Cos[otherPtLat Degree] Sin[otherPtLong Degree],Sin[otherPtLat Degree]};northPole={0,0,1};vectorA = Cross[northPole,refPtVect];vectorB = Cross[otherPtVect,refPtVect];θ1=VectorAngle[vectorA,vectorB](180/π);θ2=VectorAngle[vectorA,otherPtVect](180/π);vectorC=Cross[refPtVect,vectorA];θ3=VectorAngle[vectorC,otherPtVect](180/π);azimuth=If[((θ2==90)∨(θ2==180)),If[(θ3<=90),0,180],If[(θ2>90),(360-θ1),θ1]];azimuth ];

Computation

azimuth3=findAzimuthVectInput[refPtVect,otherPtVect];

azimuth4=N[azimuth3];

azimuth5=findAzimuthGeogInput[refPtLat,refPtLong,otherPtLat,otherPtLong];

azimuth6=N[azimuth5];

Input Summary

Reference point {lat, long}

{refPtLat,refPtLong}

Other point {lat, long}

{otherPtLat,otherPtLong}

Output

The direction from the reference point to the other point (i.e., the azimuth) is given in degrees relative to north, where due north is 0°, due east is 90 °, due south is 180°, and due west is 270 °.  If the modules faithfully reproduce the “longhand” computational process, all three of the output values shown below should be equal to one another.

azimuth2 is the value computed using a “longhand” process involving several global variables and lines of code.

azimuth2

azimuth4 is the value computed using the module findAzimuthVectInput[], whose input data are two unit location vectors.

azimuth4

azimuth6 is the value computed using the module findAzimuthGeogInput[], whose input data are two pairs of geographic coordinates.

azimuth6

References

Davis, H.F., and Snider, A.D., 1987, Introduction to Vector Analysis [5th edition]:  Boston, Massachusets, Allyn and Bacon, 365 p., ISBN 0-205-10263-8.