docs: def sphofs (lat1, lon1, r, pa, tol=1e-2, rmax=None):
      """Offset from one location on the sphere to another.
  
      This function is given a start location, expressed as a latitude and
      longitude, a distance to offset, and a direction to offset (expressed as a
      bearing, AKA position angle). It uses these to compute a final location.
      This function mirrors :func:`sphdist` and :func:`sphbear` such that::
  
        # If:
        r = sphdist (lat1, lon1, lat2a, lon2a)
        pa = sphbear (lat1, lon1, lat2a, lon2a)
        lat2b, lon2b = sphofs (lat1, lon1, r, pa)
        # Then lat2b = lat2a and lon2b = lon2a
  
      Arguments are:
  
      lat1
        The latitude of the start location.
      lon1
        The longitude of the start location.
      r
        The distance to offset by.
      pa
        The position angle (“PA” or bearing) to offset towards.
      tol
        The tolerance for the accuracy of the calculation.
      rmax
        The maximum allowed offset distance.
  
      Returns a pair ``(lat2, lon2)``. All arguments and the return values are
      measured in radians. The arguments may be vectors. The PA sign convention
      is astronomical, measuring orientation east from north.
  
      Note that the ordering of the arguments and return values 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 offset is computed naively as::
  
        lat2 = lat1 + r * cos (pa)
        lon2 = lon1 + r * sin (pa) / cos (lat2)
  
      This will fail for large offsets. Error checking can be done in two ways.
      If *tol* is not None, :func:`sphdist` is used to calculate the actual
      distance between the two locations, and if the magnitude of the fractional
      difference between that and *r* is larger than *tol*, :exc:`ValueError` is
      raised. This will add an overhead to the computation that may be
      significant if you're going to be calling this function a lot.
  
      Additionally, if *rmax* is not None, magnitudes of *r* greater than *rmax*
      are rejected. For reference, an *r* of 0.2 (~11 deg) gives a maximum
      fractional distance error of ~3%.
  
      """
      if rmax is not None and np.abs (r) > rmax:
          raise ValueError ('sphofs radius value %f is too big for '
                            'our approximation' % r)
  
      lat2 = lat1 + r * np.cos (pa)
      lon2 = lon1 + r * np.sin (pa) / np.cos (lat2)
  
      if tol is not None:
          s = sphdist (lat1, lon1, lat2, lon2)
          if np.any (np.abs ((s - r) / s) > tol):
              raise ValueError ('sphofs approximation broke down '
                                '(%s %s %s %s %s %s %s)' % (lat1, lon1,
                                                            lat2, lon2,
                                                            r, s, pa))
  
      return lat2, lon2