docs: def sphdist (lat1, lon1, lat2, lon2):
      """Calculate the distance 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.
  
      Returns the separation 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 distance is computed with the "specialized Vincenty formula". Faster
      but more error-prone formulae are possible; see Wikipedia on Great-circle
      Distance.
  
      """
      cd = np.cos (lon2 - lon1)
      sd = np.sin (lon2 - lon1)
      c2 = np.cos (lat2)
      c1 = np.cos (lat1)
      s2 = np.sin (lat2)
      s1 = np.sin (lat1)
      a = np.sqrt ((c2 * sd)**2 + (c1 * s2 - s1 * c2 * cd)**2)
      b = s1 * s2 + c1 * c2 * cd
      return np.arctan2 (a, b)