my $bearing = atan2($lon2 - $lon1, aths($lat2) - aths($lat1));
sub aths { # atanh(sin(x))
   my ($x) = @_;
   my $y = sin($x);
   return 0.5*log((1 + $y) / (1 - $y));
}