Calculate the bearing between two locations on a sphere.

    lat1
      The latitude of the first location.
    lon1
      The longitude of the first location.
    lat2
      The latitude of the second location.
    lon2
      The longitude of the second location.
    tol
      Tolerance for checking proximity to poles and rounding to zero.

    The bearing (AKA the position angle, PA) is the orientation of point 2
    with regards to point 1 relative to the longitudinal axis. Returns the
    bearing in radians. All arguments are in radians as well. The arguments
    may be vectors.

    Note that the ordering of the arguments maps to the nonstandard ordering
    ``(Dec, RA)`` in equatorial coordinates. In a spherical projection it maps
    to ``(Y, X)`` which may also be unexpected.

    The sign convention is astronomical: bearings range from -π to π, with
    negative values if point 2 is in the western hemisphere with regards to
    point 1, positive if it is in the eastern. (That is, “east from north”.)
    If point 1 is very near the pole, the bearing is undefined and the result
    is NaN.

    The *tol* argument is used for checking proximity to the poles and for
    rounding the bearing to precisely zero if it's extremely small.

    Derived from ``bear()`` in `angles.py from Prasanth Nair
    <https://github.com/phn/angles>`_. His version is BSD licensed. This one
    is sufficiently different that I think it counts as a separate
    implementation.

def sphbear (lat1, lon1, lat2, lon2, tol=1e-15):
    """Calculate the bearing between two locations on a sphere.

    lat1
      The latitude of the first location.
    lon1
      The longitude of the first location.
    lat2
      The latitude of the second location.
    lon2
      The longitude of the second location.
    tol
      Tolerance for checking proximity to poles and rounding to zero.

    The bearing (AKA the position angle, PA) is the orientation of point 2
    with regards to point 1 relative to the longitudinal axis. Returns the
    bearing in radians. All arguments are in radians as well. The arguments
    may be vectors.

    Note that the ordering of the arguments maps to the nonstandard ordering
    ``(Dec, RA)`` in equatorial coordinates. In a spherical projection it maps
    to ``(Y, X)`` which may also be unexpected.

    The sign convention is astronomical: bearings range from -π to π, with
    negative values if point 2 is in the western hemisphere with regards to
    point 1, positive if it is in the eastern. (That is, “east from north”.)
    If point 1 is very near the pole, the bearing is undefined and the result
    is NaN.

    The *tol* argument is used for checking proximity to the poles and for
    rounding the bearing to precisely zero if it's extremely small.

    Derived from ``bear()`` in `angles.py from Prasanth Nair
    <https://github.com/phn/angles>`_. His version is BSD licensed. This one
    is sufficiently different that I think it counts as a separate
    implementation.

    """
    # cross product on outer axis:
    ocross = lambda a, b: np.cross (a, b, axisa=0, axisb=0, axisc=0)

    # if args have shape S, this has shape (3, S)
    v1 = np.asarray ([np.cos (lat1) * np.cos (lon1),
                      np.cos (lat1) * np.sin (lon1),
                      np.sin (lat1)])

    v2 = np.asarray ([np.cos (lat2) * np.cos (lon2),
                      np.cos (lat2) * np.sin (lon2),
                      np.sin (lat2)])

    is_bad = (v1[0]**2 + v1[1]**2) < tol

    p12 = ocross (v1, v2) # ~"perpendicular to great circle containing points"
    p1z = np.asarray ([v1[1], -v1[0], np.zeros_like (lat1)]) # ~"perp to base and Z axis"
    cm = np.sqrt ((ocross (p12, p1z)**2).sum (axis=0)) # ~"angle between the vectors"
    bearing = np.arctan2 (cm, np.sum (p12 * p1z, axis=0))
    bearing = np.where (p12[2] < 0, -bearing, bearing) # convert to [-pi/2, pi/2]
    bearing = np.where (np.abs (bearing) < tol, 0, bearing) # clamp
    bearing[np.where (is_bad)] = np.nan
    return bearing