Angular spearation between two points on a unit sphere.

    This will be an angle between [0, π] radians.

    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
    -----
    The great cicle angular separation of the second point from the
    first is returned as an angle in radians. the return value is
    always in the range [0, π].

    Results agree with those from SLALIB routine sla_dsep. See
    test_sep_against_slalib_dsep() in test_angles.py.


    Examples
    --------
    >>> r2d(sep(0, 0, 0, d2r(90.0)))
    90.0
    >>> r2d(sep(0, d2r(45.0), 0, d2r(90.0)))
    45.00000000000001
    >>> r2d(sep(0, d2r(-45.0), 0, d2r(90.0)))
    135.0

    >>> r2d(sep(0, d2r(-90.0), 0, d2r(90.0)))
    180.0
    >>> r2d(sep(d2r(45.0), d2r(-90.0), d2r(45.0), d2r(90.0)))
    180.0
    >>> r2d(sep(0, 0, d2r(90.0), 0))
    90.0

    >>> r2d(sep(0, d2r(45.0), d2r(90.0), d2r(45.0)))
    60.00000000000001
    >>> import math
    >>> 90.0 * math.cos(d2r(45.0))  # Distance along latitude circle.
    63.63961030678928

def sep(a1, b1, a2, b2):
    """Angular spearation between two points on a unit sphere.

    This will be an angle between [0, π] radians.

    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
    -----
    The great cicle angular separation of the second point from the
    first is returned as an angle in radians. the return value is
    always in the range [0, π].

    Results agree with those from SLALIB routine sla_dsep. See
    test_sep_against_slalib_dsep() in test_angles.py.


    Examples
    --------
    >>> r2d(sep(0, 0, 0, d2r(90.0)))
    90.0
    >>> r2d(sep(0, d2r(45.0), 0, d2r(90.0)))
    45.00000000000001
    >>> r2d(sep(0, d2r(-45.0), 0, d2r(90.0)))
    135.0

    >>> r2d(sep(0, d2r(-90.0), 0, d2r(90.0)))
    180.0
    >>> r2d(sep(d2r(45.0), d2r(-90.0), d2r(45.0), d2r(90.0)))
    180.0
    >>> r2d(sep(0, 0, d2r(90.0), 0))
    90.0

    >>> r2d(sep(0, d2r(45.0), d2r(90.0), d2r(45.0)))
    60.00000000000001
    >>> import math
    >>> 90.0 * math.cos(d2r(45.0))  # Distance along latitude circle.
    63.63961030678928
    """
    # Tolerance to decide if the calculated separation is zero.
    tol = 1e-15

    v = CartesianVector.from_spherical(1.0, a1, b1)
    v2 = CartesianVector.from_spherical(1.0, a2, b2)
    d = v.dot(v2)
    c = v.cross(v2).mod

    res = math.atan2(c, d)

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