docs: 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