How to use the pygeodesy.utily.sincos2 function in PyGeodesy

s2, c2, s1, c1, _, c21 = sincos2(phi2, phi1, lam21)
    E = datum.ellipsoid
    if E.f and abs(c1) > EPS and abs(c2) > EPS:
        r1 = atan(E.b_a * s1 / c1)
        r2 = atan(E.b_a * s2 / c2)

        j = (r2 + r1) / 2.0
        k = (r2 - r1) / 2.0
        sj, cj, sk, ck, sl_2, _ = sincos2(j, k, lam21 / 2.0)

        h = fsum_(sk**2, (ck * sl_2)**2, -(sj * sl_2)**2)
        if EPS < abs(h) < EPS1:
            u = 1 / (1 - h)
            d = 2 * atan(sqrt(h * u))  # == acos(1 - 2 * h)
            sd, cd = sincos2(d)
            if abs(sd) > EPS:
                u = 2 * (sj * ck)**2 * u
                v = 2 * (sk * cj)**2 / h
                x = u + v
                y = u - v

                t = d / sd
                s = 4 * t**2
                e = 2 * cd
                a = s * e
                b = 2 * d
                c = t - (a - e) / 2.0

                s = fsum_(a * x,  c * x**2, -b * y, -e * y**2, s * x * y) * E.f / 16.0
                s = fsum_(t * x, -y, -s) * E.f / 4.0
                return d - s * sd

limit and/or if this and the B{C{other}}
                                 point are coincident or near-antipodal.
        '''
        E = self.ellipsoids(other)

        c1, s1, _ = _r3(self.lat, E.f)
        c2, s2, _ = _r3(other.lat, E.f)

        c1c2, s1c2 = c1 * c2, s1 * c2
        c1s2, s1s2 = c1 * s2, s1 * s2

        dl, _ = unroll180(self.lon, other.lon, wrap=wrap)
        ll = dl = radians(dl)
        for self._iteration in range(self._iterations):
            ll_ = ll
            sll, cll = sincos2(ll)

            ss = hypot(c2 * sll, c1s2 - s1c2 * cll)
            if ss < EPS:  # coincident or antipodal, ...
                if self.isantipodeTo(other, eps=self._epsilon):
                    t = '%r %sto %r' % (self, _antipodal_, other)
                    raise VincentyError(_ambiguous_, txt=t)
                # return zeros like Karney, but unlike Veness
                return Distance3Tuple(0.0, 0, 0)

            cs = s1s2 + c1c2 * cll
            s = atan2(ss, cs)

            sa = c1c2 * sll / ss
            c2a = 1 - sa**2
            if abs(c2a) < EPS:
                c2a = 0  # equatorial line

i = radians(bearing)  # initial bearing (forward azimuth)
        si, ci = sincos2(i)
        s12 = atan2(t1, ci) * 2

        sa = c1 * si
        c2a = 1 - sa**2
        if c2a < EPS:
            c2a = 0
            A, B = 1, 0
        else:  # e22 == (a / b)**2 - 1
            A, B = _p2(c2a * E.e22)

        s = d = distance / (E.b * A)
        for self._iteration in range(self._iterations):
            ss, cs = sincos2(s)
            c2sm = cos(s12 + s)
            s_, s = s, d + _ds(B, cs, ss, c2sm)
            if abs(s - s_) < self._epsilon:
                break
        else:
            raise VincentyError(_no_convergence_, txt=repr(self))  # self.toRepr()

        t = s1 * ss - c1 * cs * ci
        # final bearing (reverse azimuth +/- 180)
        r = degrees360(atan2(sa, -t))

        if llr:
            # destination latitude in [-90, 90)
            a = degrees90(atan2(s1 * cs + c1 * ss * ci,
                                (1 - E.f) * hypot(sa, t)))
            # destination longitude in [-180, 180)

def _reverse(self, x, y, name, LatLon, LatLon_kwds, _c_t, lea):
        '''(INTERNAL) Azimuthal (spherical) reverse C{x, y} to C{lat, lon}.
        '''
        x = Scalar(x, name=_x_)
        y = Scalar(y, name=_y_)

        r = hypot(x, y)
        c, t = _c_t(r / self.radius)
        if t:
            s0, c0 = self._sc0
            sc, cc = sincos2(c)
            k = c / sc
            z = atan2d(x, y)  # (x, y) for azimuth from true North

            lat = degrees(asin1(s0 * cc + c0 * sc * (y / r)))
            if lea or abs(c0) > EPS:
                lon = atan2(x * sc, c0 * cc * r - s0 * sc * y)
            else:
                lon = atan2(x, (y if s0 < 0 else -y))
            lon = _norm180(self.lon0 + degrees(lon))
        else:
            k, z = 1, 0
            lat, lon = self.latlon0

        t = self._toLatLon(lat, lon, LatLon, LatLon_kwds) if LatLon else \
            Azimuthal7Tuple(x, y, lat, lon, z, k, self.datum)
        return _xnamed(t, name or self.name)

       @arg lam2: End longitude (C{radians}).
       @kwarg final: Return final bearing if C{True}, initial
                     otherwise (C{bool}).
       @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}).

       @return: Initial or final bearing (compass C{radiansPI2}) or
                zero if start and end point coincide.
    '''
    if final:
        phi1, lam1, phi2, lam2 = phi2, lam2, phi1, lam1
        r = PI2 + PI
    else:
        r = PI2

    db, _ = unrollPI(lam1, lam2, wrap=wrap)
    sa1, ca1, sa2, ca2, sdb, cdb = sincos2(phi1, phi2, db)

    # see 
    x = ca1 * sa2 - sa1 * ca2 * cdb
    y = sdb * ca2
    return (atan2(y, x) + r) % PI2  # wrapPI2

       @return: Angular distance (C{radians}).

       @see: Functions L{vincentys}, L{cosineAndoyerLambert_},
             L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},
             L{equirectangular_}, L{euclidean_}, L{flatLocal_}/L{hubeny_},
             L{flatPolar_}, L{haversine_} and L{thomas_}.

       @note: Functions L{vincentys_}, L{haversine_} and L{cosineLaw_}
              produce equivalent results, but L{vincentys_} is suitable
              for antipodal points and slightly more expensive (M{3 cos,
              3 sin, 1 hypot, 1 atan2, 6 mul, 2 add}) than L{haversine_}
              (M{2 cos, 2 sin, 2 sqrt, 1 atan2, 5 mul, 1 add}) and
              L{cosineLaw_} (M{3 cos, 3 sin, 1 acos, 3 mul, 1 add}).
    '''
    sa1, ca1, sa2, ca2, sb21, cb21 = sincos2(phi1, phi2, lam21)

    c = ca2 * cb21
    x = sa1 * sa2 + ca1 * c
    y = ca1 * sa2 - sa1 * c
    return atan2(hypot(ca2 * sb21, y), x)

       @kwarg datum: Ellipsoidal datum to use (L{Datum}).

       @return: Angular distance (C{radians}).

       @raise TypeError: Invalid B{C{datum}}.

       @see: Functions L{thomas}, L{cosineAndoyerLambert_},
             L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},
             L{equirectangular_}, L{euclidean_}, L{flatLocal_}/L{hubeny_},
             L{flatPolar_}, L{haversine_} and L{vincentys_} and U{Geodesy-PHP
             }.
    '''
    _xinstanceof(Datum, datum=datum)

    s2, c2, s1, c1, _, c21 = sincos2(phi2, phi1, lam21)
    E = datum.ellipsoid
    if E.f and abs(c1) > EPS and abs(c2) > EPS:
        r1 = atan(E.b_a * s1 / c1)
        r2 = atan(E.b_a * s2 / c2)

        j = (r2 + r1) / 2.0
        k = (r2 - r1) / 2.0
        sj, cj, sk, ck, sl_2, _ = sincos2(j, k, lam21 / 2.0)

        h = fsum_(sk**2, (ck * sl_2)**2, -(sj * sl_2)**2)
        if EPS < abs(h) < EPS1:
            u = 1 / (1 - h)
            d = 2 * atan(sqrt(h * u))  # == acos(1 - 2 * h)
            sd, cd = sincos2(d)
            if abs(sd) > EPS:
                u = 2 * (sj * ck)**2 * u

def _destination2(a, b, r, t):
    '''(INTERNAL) Destination phi- and longitude in C{radians}.

       @arg a: Latitude (C{radians}).
       @arg b: Longitude (C{radians}).
       @arg r: Angular distance (C{radians}).
       @arg t: Bearing (compass C{radians}).

       @return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}).
    '''
    # see 
    sa, ca, sr, cr, st, ct = sincos2(a, r, t)

    a = asin(ct * sr * ca + cr * sa)
    d = atan2(st * sr * ca, cr - sa * sin(a))
    # note, in EdWilliams.org/avform.htm W is + and E is -
    return a, b + d

@arg lat: Latitude at the crossing (C{degrees}).
           @kwarg wrap: Wrap and unroll longitudes (C{bool}).

           @return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or
                    C{None} if the great circle doesn't reach B{C{lat}}.
        '''
        self.others(other)

        a1, b1 = self.philam
        a2, b2 = other.philam

        a = radians(lat)
        db, b2 = unrollPI(b1, b2, wrap=wrap)

        sa,  ca,  sa1, ca1, \
        sa2, ca2, sdb, cdb = sincos2(a, a1, a2, db)

        x = sa1 * ca2 * ca * sdb
        y = sa1 * ca2 * ca * cdb - ca1 * sa2 * ca
        z = ca1 * ca2 * sa * sdb

        h = hypot(x, y)
        if h < EPS or abs(z) > h:
            return None  # great circle doesn't reach latitude

        m = atan2(-y, x) + b1  # longitude at max latitude
        d = acos1(z / h)  # delta longitude to intersections
        return degrees180(m - d), degrees180(m + d)