Find bearing/position angle between two points on a unit sphere.

    Parameters
    ----------
    a1, b1 : float
        Longitude-like and latitude-like angles defining the first
        point. Both are in radians.

    a2, b2 : float
        Longitude-like and latitude-like angles defining the second
        point. Both are in radians.

    Notes
    -----
    Position angle of the second point with respect to the first
    is returned in radians. Position angle is calculated clockwise
    and counter-clockwise from the direction towards the North
    pole. It is between [0 and π] if the second point is in the
    eastern hemisphere w.r.t the first, and between (0, -π) if
    the second point is in the western hemisphere w.r.t the first.

    .. warning::

        If the first point is at the pole then bearing is undefined and
        0 is returned.

    Results agree with those from SLALIB rountine sla_dbear. See
    test_bear_against_slalib_dbear() in test_angles.py.

    Examples
    --------
    >>> from angles import bear, r2d, d2r
    >>> bear(0, 0, 0, -d2r(90.0))
    3.141592653589793
    >>> bear(0, -d2r(90.0), 0, 0)
    0.0
    >>> bear(0, -d2r(45.0), 0, 0)
    0.0
    >>> bear(0, -d2r(89.678), 0, 0)
    0.0

    >>> r2d(bear(d2r(45.0), d2r(45.0), d2r(46.0), d2r(45.0)))
    89.64644212193384

    >>> r2d(bear(d2r(45.0), d2r(45.0), d2r(44.0), d2r(45.0)))
    -89.64644212193421

def bear(a1, b1, a2, b2):
    """Find bearing/position angle between two points on a unit sphere.

    Parameters
    ----------
    a1, b1 : float
        Longitude-like and latitude-like angles defining the first
        point. Both are in radians.

    a2, b2 : float
        Longitude-like and latitude-like angles defining the second
        point. Both are in radians.

    Notes
    -----
    Position angle of the second point with respect to the first
    is returned in radians. Position angle is calculated clockwise
    and counter-clockwise from the direction towards the North
    pole. It is between [0 and π] if the second point is in the
    eastern hemisphere w.r.t the first, and between (0, -π) if
    the second point is in the western hemisphere w.r.t the first.

    .. warning::

        If the first point is at the pole then bearing is undefined and
        0 is returned.

    Results agree with those from SLALIB rountine sla_dbear. See
    test_bear_against_slalib_dbear() in test_angles.py.

    Examples
    --------
    >>> from angles import bear, r2d, d2r
    >>> bear(0, 0, 0, -d2r(90.0))
    3.141592653589793
    >>> bear(0, -d2r(90.0), 0, 0)
    0.0
    >>> bear(0, -d2r(45.0), 0, 0)
    0.0
    >>> bear(0, -d2r(89.678), 0, 0)
    0.0

    >>> r2d(bear(d2r(45.0), d2r(45.0), d2r(46.0), d2r(45.0)))
    89.64644212193384

    >>> r2d(bear(d2r(45.0), d2r(45.0), d2r(44.0), d2r(45.0)))
    -89.64644212193421

    """
    # Find perpendicular to the plane containing the base and
    # z-axis. Then find the perpendicular to the plane containing
    # the base and the target. The angle between these two is the
    # position angle or bearing of the target w.r.t the base. Check
    # sign of the z component of the latter vector to determine
    # quadrant: 1st and 2nd quadrants are +ve while 3rd and 4th are
    # negative.
    #
    # Tolerance to decide if first is on the pole and also to decide if
    # the calculated bearing is zero.
    tol = 1e-15

    v1 = CartesianVector.from_spherical(1.0, a1, b1)
    v2 = CartesianVector.from_spherical(1.0, a2, b2)

    # Z-axis
    v0 = CartesianVector.from_spherical(r=1.0, alpha=0.0, delta=d2r(90.0))

    if abs(v1.cross(v0).mod) < tol:
        # The first point is on the pole. Bearing is undefined.
        warnings.warn(
            "First point is on the pole. Bearing undefined.")
        return 0.0

    # Vector perpendicular to great circle containing two points.
    v12 = v1.cross(v2)

    # Vector perpendicular to great circle containing base and
    # Z-axis.
    v10 = v1.cross(v0)

    # Find angle between these two vectors.
    dot = v12.dot(v10)
    cross = v12.cross(v10).mod
    x = math.atan2(cross, dot)

    # If z is negative then we are in the 3rd or 4th quadrant.
    if v12.z < 0:
        x = -x

    if abs(x) < tol:
        return 0.0
    else:
        return x