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%.

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