True-North Bearing between two GPS coordinates

Question

I am trying to find the correct equation to find the bearing off of true-north between two GPS coordinates. I have seen a lot of online calculators and a lot of formulas on this GIS stack exchange, but many of them give me different results! What is the most correct method?

It seems the most popular one is this method: Calculate bearing between two decimal GPS coordinates (Arduino / C++)

Another method is this one: Calculate bearing between two decimal GPS coordinates , but, I am getting completely different answers when using this method (off by like tens of degrees), and I can't even get a matching result on the example from that post: "I get 250.20613449 as expected."

On all of my attempts, I am consistently using the degrees notation for lat and lon (I'm not using the hours, minutes, second notation). Also, on all of my attempts, I swap the start/stop position of both coordinates and recalculate just to make sure I'm not solving relatively to the wrong start position.

Online calculators seem to give me a variance as well. Sometimes they match the first method, sometimes they don't.

Can someone help nip this in the bud, once and for all?

Here's an example prob we can focus on:

Start: 42.766698 -75.746404 Stop: 42.761405 -75.757653

If I plot these points on a map and guesstimate the bearing by eye:

Answer 1

In the screenshot included in your question, a straight line is drawn on a Mercator projection from a sphere. The angle drwan is the angle returned by the first question/answer referred.

---

About the first answer (Calculate bearing between two decimal GPS coordinates (Arduino / C++)): It is computing the Rhumb Line (a line with constant bearing) azimuth over a sphere.

The formulas are converting geographic coordinates to a Mercator projection over a unitary sphere, see http://mathworld.wolfram.com/MercatorProjection.html.

Since Mercator projection is conformal, the straight line in the projection is a rhumb line over the sphere. But we do know that the shortest path, over the curved suface, between two points is not a rhumb line.

About the second answer (Calculate bearing between two decimal GPS coordinates): It is computing the intial azimuth of the travel over the Great Circle (the intersection of a plane that contains the center of a sphere and its surface) that contains both points, using spherical trigonometry. You can see Link.

Also, we do know that the Earth is not a sphere (also not an ellipsoid of revolution, but closer).

All depends mainly of two factors, the time that you have to compute the bearing, and the accuracy that you need for the result.

If you want to know the WGS84 ellipsoidal initial azimuth in an online calculator, use GeodSolve to solve the Inverse direct problem of geodesy. If you want to have it installed in your machine, install the GeographicLib Utility programs, and if you want to implement it in your code, just see GeographicLib.