Now that is an interesting topic to discuss! I know what you mean BUT I had the following in place and ran it against the naive/fast atan2 approach and for where I am I found the latter produced a result that if I looked in the sky it gave a better answer! Maybe magnetic deviation? Maybe because we’re looking at very short ranges here so curvature of earth actually is not the right way to go???
Would love to understand! At ranges of 20nm or less I was seeing deviations between the two of 10 or so degrees but agreement at compass “quarters” N E S W.
$dLon = ($lon1-$lon_home);
$y = sin($dLon) * cos($lat1);
$x = cos($lat_home)*sin($lat1) - sin($lat_home)*cos($lat1)*cos($dLon);
$brng=(atan2($y, $x)) / $p; // i.e. * 180 / PI;
$nearest_bearing1=round(360 - (($brng + 360) % 360));
echo("Bearing 1: ".$nearest_bearing1."°<br/>");