How to use the pygeodesy.fmath.fsum_ function in PyGeodesy

n2, d22 = _nd2(point2, distance2, r, _2_, point1)
    n3, d32 = _nd2(point3, distance3, r, '3', point1, point2)

    # the following uses x,y coordinate system with origin at n1, x axis n1->n2
    y = n3.minus(n1)
    x = n2.minus(n1)

    d = x.length  # distance n1->n2
    if d > EPS_2:  # and y.length > EPS_2:
        X = x.unit()  # unit vector in x direction n1->n2
        i = X.dot(y)  # signed magnitude of x component of n1->n3
        Y = y.minus(X.times(i)).unit()  # unit vector in y direction
        j = Y.dot(y)  # signed magnitude of y component of n1->n3
        if abs(j) > EPS_2:
            # courtesy Carlos Freitas 
            x = fsum_(d12, -d22, d**2) / (2 * d)  # n1->intersection x- and ...
            y = fsum_(d12, -d32, i**2, j**2) / (2 * j) - (x * i / j)  # ... y-component

            n = n1.plus(X.times(x)).plus(Y.times(y))  # .plus(Z.times(z))
            if useZ:  # include non-NaN, non-zero Z component
                z = fsum_(d12, -(x**2), -(y**2))
                if z > EPS:
                    Z = X.cross(Y)  # unit vector perpendicular to plane
                    n = n.plus(Z.times(sqrt(z)))

            if height is None:
                h = fidw((point1.height, point2.height, point3.height),
                         map1(fabs, distance1, distance2, distance3))
            else:
                h = Height(height)
            kwds = _xkwds(LatLon_kwds, height=h, LatLon=LatLon)
            return n.toLatLon(**kwds)  # Nvector(n.x, n.y, n.z).toLatLon(...)

return _latlon3(degrees90(a), degrees180(b), h,
                        intersections2, LatLon, **LatLon_kwds)

    r1, r2, f = _rads3(rad1, rad2, radius)
    if f:  # swapped
        c1, c2 = c2, c1  # PYCHOK swap

    a1, b1 = c1.philam
    a2, b2 = c2.philam

    db, _ = unrollPI(b1, b2, wrap=wrap)
    d = vincentys_(a2, a1, db)  # radians
    if d < max(r1 - r2, EPS):
        raise ValueError(_near_concentric_)

    x = fsum_(r1, r2, -d)  # overlap
    if x > EPS:
        sd, cd, sr1, cr1, _, cr2 = sincos2(d, r1, r2)
        x = sd * sr1
        if abs(x) < EPS:
            raise ValueError(_invalid_)
        x = acos1((cr2 - cd * cr1) / x)  # 0 <= x <= PI

    elif x < 0:
        t = (d * radius) if too_d is None else too_d
        raise ValueError(_too_distant_fmt_ % (t,))

    if height is None:  # "radical height"
        f, _ = _radical2(d, r1, r2)  # "radical ratio"
        h = Height(favg(c1.height, c2.height, f=f))
    else:
        h = Height(height)

r = self._v4t
        if r is None or d != self.datum:
            # 
            E = d.ellipsoid
            x, y, z = self.xyz

            # Kenneth Gade eqn 23
            p = hypot2(x, y) * E.a2_
            q = (z**2 * E.e12) * E.a2_
            r = fsum_(p, q, -E.e4) / 6
            s = (p * q * E.e4) / (4 * r**3)
            t = cbrt(fsum_(1, s, sqrt(s * (2 + s))))

            u = r * fsum_(1, t, 1 / t)
            v = sqrt(u**2 + E.e4 * q)
            w = E.e2 * fsum_(u, v, -q) / (2 * v)

            k = sqrt(fsum_(u, v, w**2)) - w
            if abs(k) < EPS:
                raise _ValueError(origin=self)
            e = k / (k + E.e2)
#           d = e * hypot(x, y)

#           tmp = 1 / hypot(d, z) == 1 / hypot(e * hypot(x, y), z)
            t = hypot_(e * x, e * y, z)  # == 1 / tmp
            if t < EPS:
                raise _ValueError(origin=self)
            h = fsum_(k, E.e2, -1) / k * t

            s = e / t  # == e * tmp

def _radicaline(d, r1, r2):
    # x coord [0..d] of the "radical line", perpendicular to
    # the x-axis line between both centers (0, 0) and (d, 0)
    return fsum_(d**2, r1**2, -(r2**2)) / (2 * d)

if min(p, r) > EPS:
            # parametric latitude (Bowring eqn 17, replaced)
            t = (E.b * z) / (E.a * p) * (1 + E.e22 * E.b / r)
            c = 1 / hypot1(t)
            s = t * c

            # geodetic latitude (Bowring eqn 18)
            a = atan2(z + E.e22 * E.b * s**3,
                      p - E.e2  * E.a * c**3)
            b = atan2(y, x)  # ... and longitude

            # height above ellipsoid (Bowring eqn 7)
            sa, ca = sincos2(a)
#           r = E.a / E.e2s(sa)  # length of normal terminated by minor axis
#           h = p * ca + z * sa - (E.a * E.a / r)
            h = fsum_(p * ca, z * sa, -E.a * E.e2s(sa))

            C, lat, lon = 1, degrees90(a), degrees180(b)

        # see 
        elif p > EPS:  # lat arbitrarily zero
            C, lat, lon, h = 2, 0.0, degrees180(atan2(y, x)), p - E.a

        else:  # polar lat, lon arbitrarily zero
            C, lat, lon, h = 3, copysign(90.0, z), 0.0, abs(z) - E.b

        r = Ecef9Tuple(x, y, z, lat, lon, h, C, None, self.datum)
        return self._xnamed(r, name)

a2, b2 = other.philam
        if abs(b2 - b1) > PI:
            b1 += PI2  # crossing anti-meridian

        a3 = favg(a1, a2)
        b3 = favg(b1, b2)

        f1 = tanPI_2_2(a1)
        if abs(f1) > EPS:
            f2 = tanPI_2_2(a2)
            f = f2 / f1
            if abs(f) > EPS:
                f = log(f)
                if abs(f) > EPS:
                    f3 = tanPI_2_2(a3)
                    b3 = fsum_(b1 * log(f2),
                              -b2 * log(f1), (b2 - b1) * log(f3)) / f

        h = self._havg(other) if height is None else Height(height)
        return self.classof(degrees90(a3), degrees180(b3), height=h)

s = E.e2s2(sa)  # r, v = E.roc2_(sa, _F0); r = v / r
    v = E.a * _F0 / sqrt(s)  # nu
    r = s / E.e12  # nu / rho == v / (v * E.e12 / s) == s / E.e12

    x2 = r - 1  # η2
    ta = tan(a)

    ca3, ca5 = fpowers(ca, 5, 3)  # PYCHOK false!
    ta2, ta4 = fpowers(ta, 4, 2)  # PYCHOK false!

    vsa = v * sa
    I4 = (E.b * _F0 * _M(E.Mabcd, a) + _N0,
         (vsa /   2) * ca,
         (vsa /  24) * ca3 * fsum_(5, -ta2, 9 * x2),
         (vsa / 720) * ca5 * fsum_(61, ta4, -58 * ta2))

    V4 = (_E0,
         (v        * ca),
         (v /   6) * ca3 * (r - ta2),
         (v / 120) * ca5 * fdot((-18, 1, 14, -58), ta2, 5 + ta4, x2, ta2 * x2))

    d, d2, d3, d4, d5, d6 = fpowers(b - _B0, 6)  # PYCHOK false!
    n = fdot(I4, 1, d2, d4, d6)
    e = fdot(V4, 1, d,  d3, d5)

    if Osgr is None:
        r = _EasNor2Tuple(e, n)
    else:
        r = Osgr(e, n, datum=_Datums_OSGB36, **Osgr_kwds)
        if lon is None and isinstance(latlon, _LLEB):
            r._latlon = latlon  # XXX weakref(latlon)?

}.
    '''
    _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
                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

if abs(sr) > EPS:
            r2 = r**2

            p = (s1 + s2)**2 / (1 + cr)
            q = (s1 - s2)**2 / (1 - cr)
            x = p + q
            y = p - q
            s = 8 * r2 / sr

            a = 64 * r + 2 * s * cr  # 16 * r2 / tan(r)
            d = 48 * sr + s  # 8 * r2 / tan(r)
            b = -2 * d
            e = 30 * s2r
            c = fsum_(30 * r, e / 2, s * cr)  # 8 * r2 / tan(r)

            d = fsum_(a * x, b * y, -c * x**2, d * x * y, e * y**2) * E.f / 32.0
            d = fsum_(d, -x * r, 3 * y * sr) * E.f / 4.0
            r += d
    return r