That all assumes a spherical earth, doesn't it?
But if we want more accurate results we assume the sphere is slightly squashed-- likely we'll use the WGS84 spheroid, where the distance between the poles divided by the equatorial diameter is 297.257223563 divided by 298.257223563.
About azimuth: those PDFs I linked give the azimuth of the geodesic-- the shortest path on the surface of the spheroid between Point 1 and Point 2. It eventually occurred to me that you would do as well or better with the normal-section azimuth, which is the azimuth (at Point 1) of the plane that passes through Points 1 and 2 and is vertical at Point 1. In other words, you're standing vertically at Point 1, facing Point 2-- which way are you facing?
That's much simpler. In the following formulas Z1, Z2 etc are just the intermediate results in the calculation; A and B are constants that depend on the spheroid you choose. For WGS84, A rounds to 1.0067394967 and B rounds to 0.0066943800. Point 1 is at latitude Lat1; Point 2 is at latitude Lat2, DLon degrees to the west of Point 1.
Z1 = A + tangent-squared of Lat2
Z2 = A + tangent-squared of Lat1
Z3 = B x square root of (Z1/Z2)
Z4 = (tangent Lat2) divided by (tangent Lat1)
Z5 = Z3 + (Z4/A)
Z6 = [(cosine DLon) - Z5] x sine Lat1
tangent of azimuth = (sine DLon)/Z6
This doesn't work if Point 1 is on the equator (zero in the denominator); in that case, use
tangent of azimuth = (-sine DLon) x A divided by (tangent Lat2)
Last edited by Tim Zukas; 29th February 2008 at 17:30.