How to use the gtsam.3rdparty.GeographicLib.python.geographiclib.geodesic.Geodesic function in gtsam

#    2 93315 102225
        #    3 36189  32341
        #    4  5396      7
        #    5   455      1
        #    6    56      0
        #
        # Because omg12 is near pi, estimate work with omg12a = pi - omg12
        k = Geodesic.Astroid(x, y)
        omg12a = lamscale * ( -x * k/(1 + k) if self._f >= 0
                               else -y * (1 + k)/k )
        somg12 = math.sin(omg12a); comg12 = -math.cos(omg12a)
        # Update spherical estimate of alp1 using omg12 instead of lam12
        salp1 = cbet2 * somg12
        calp1 = sbet12a - cbet2 * sbet1 * Math.sq(somg12) / (1 - comg12)
    if salp1 > 0:               # Sanity check on starting guess
      salp1, calp1 = Geodesic.SinCosNorm(salp1, calp1)
    else:
      salp1 = 1; calp1 = 0
    return sig12, salp1, calp1, salp2, calp2, dnm

def Lambda12(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1, diffp,
               # Scratch areas of the right size
               C1a, C2a, C3a):
    """Private: Solve hybrid problem"""
    if sbet1 == 0 and calp1 == 0:
      # Break degeneracy of equatorial line.  This case has already been
      # handled.
      calp1 = -Geodesic.tiny_

    # sin(alp1) * cos(bet1) = sin(alp0)
    salp0 = salp1 * cbet1
    calp0 = math.hypot(calp1, salp1 * sbet1) # calp0 > 0

    # real somg1, comg1, somg2, comg2, omg12, lam12
    # tan(bet1) = tan(sig1) * cos(alp1)
    # tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1)
    ssig1 = sbet1; somg1 = salp0 * sbet1
    csig1 = comg1 = calp1 * cbet1
    ssig1, csig1 = Geodesic.SinCosNorm(ssig1, csig1)
    # SinCosNorm(somg1, comg1); -- don't need to normalize!

    # Enforce symmetries in the case abs(bet2) = -bet1.  Need to be careful
    # about this case, since this can yield singularities in the Newton
    # iteration.

outmask determines which fields get included and if outmask is
    omitted, then only the basic geodesic fields are computed.  The mask
    is an or'ed combination of the following values

      Geodesic.LATITUDE
      Geodesic.LONGITUDE
      Geodesic.AZIMUTH
      Geodesic.REDUCEDLENGTH
      Geodesic.GEODESICSCALE
      Geodesic.AREA
      Geodesic.ALL
    """

    lon1 = Geodesic.CheckPosition(lat1, lon1)
    azi1 = Geodesic.CheckAzimuth(azi1)
    Geodesic.CheckDistance(s12)

    result = {'lat1': lat1, 'lon1': lon1, 'azi1': azi1, 's12': s12}
    a12, lat2, lon2, azi2, s12, m12, M12, M21, S12 = self.GenDirect(
      lat1, lon1, azi1, False, s12, outmask)
    outmask &= Geodesic.OUT_ALL
    result['a12'] = a12
    if outmask & Geodesic.LATITUDE: result['lat2'] = lat2
    if outmask & Geodesic.LONGITUDE: result['lon2'] = lon2
    if outmask & Geodesic.AZIMUTH: result['azi2'] = azi2
    if outmask & Geodesic.REDUCEDLENGTH: result['m12'] = m12
    if outmask & Geodesic.GEODESICSCALE:
      result['M12'] = M12; result['M21'] = M21
    if outmask & Geodesic.AREA: result['S12'] = S12
    return result

# being attached to 0 correctly.  The following ensures the correct
        # behavior.
        if salp12 == 0 and calp12 < 0:
          salp12 = Geodesic.tiny_ * calp1
          calp12 = -1
        alp12 = math.atan2(salp12, calp12)
      S12 += self._c2 * alp12
      S12 *= swapp * lonsign * latsign
      # Convert -0 to 0
      S12 += 0

    # Convert calp, salp to azimuth accounting for lonsign, swapp, latsign.
    if swapp < 0:
      salp2, salp1 = salp1, salp2
      calp2, calp1 = calp1, calp2
      if outmask & Geodesic.GEODESICSCALE:
        M21, M12 = M12, M21

    salp1 *= swapp * lonsign; calp1 *= swapp * latsign
    salp2 *= swapp * lonsign; calp2 *= swapp * latsign

    if outmask & Geodesic.AZIMUTH:
      # minus signs give range [-180, 180). 0- converts -0 to +0.
      azi1 = 0 - math.atan2(-salp1, calp1) / Math.degree
      azi2 = 0 - math.atan2(-salp2, calp2) / Math.degree

    # Returned value in [0, 180]
    return a12, s12, azi1, azi2, m12, M12, M21, S12

def Lengths(self, eps, sig12,
              ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2, scalep,
              # Scratch areas of the right size
              C1a, C2a):
    """Private: return a bunch of lengths"""
    # Return m12b = (reduced length)/_b; also calculate s12b = distance/_b,
    # and m0 = coefficient of secular term in expression for reduced length.
    Geodesic.C1f(eps, C1a)
    Geodesic.C2f(eps, C2a)
    A1m1 = Geodesic.A1m1f(eps)
    AB1 = (1 + A1m1) * (
      Geodesic.SinCosSeries(True, ssig2, csig2, C1a, Geodesic.nC1_) -
      Geodesic.SinCosSeries(True, ssig1, csig1, C1a, Geodesic.nC1_))
    A2m1 = Geodesic.A2m1f(eps)
    AB2 = (1 + A2m1) * (
      Geodesic.SinCosSeries(True, ssig2, csig2, C2a, Geodesic.nC2_) -
      Geodesic.SinCosSeries(True, ssig1, csig1, C2a, Geodesic.nC2_))
    m0 = A1m1 - A2m1
    J12 = m0 * sig12 + (AB1 - AB2)
    # Missing a factor of _b.
    # Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure accurate
    # cancellation in the case of coincident points.
    m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) - csig1 * csig2 * J12
    # Missing a factor of _b
    s12b = (1 + A1m1) * sig12 + AB1
    if scalep:
      csig12 = csig1 * csig2 + ssig1 * ssig2
      t = self._ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2)
      M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1

def Lengths(self, eps, sig12,
              ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2, scalep,
              # Scratch areas of the right size
              C1a, C2a):
    """Private: return a bunch of lengths"""
    # Return m12b = (reduced length)/_b; also calculate s12b = distance/_b,
    # and m0 = coefficient of secular term in expression for reduced length.
    Geodesic.C1f(eps, C1a)
    Geodesic.C2f(eps, C2a)
    A1m1 = Geodesic.A1m1f(eps)
    AB1 = (1 + A1m1) * (
      Geodesic.SinCosSeries(True, ssig2, csig2, C1a, Geodesic.nC1_) -
      Geodesic.SinCosSeries(True, ssig1, csig1, C1a, Geodesic.nC1_))
    A2m1 = Geodesic.A2m1f(eps)
    AB2 = (1 + A2m1) * (
      Geodesic.SinCosSeries(True, ssig2, csig2, C2a, Geodesic.nC2_) -
      Geodesic.SinCosSeries(True, ssig1, csig1, C2a, Geodesic.nC2_))
    m0 = A1m1 - A2m1
    J12 = m0 * sig12 + (AB1 - AB2)
    # Missing a factor of _b.
    # Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure accurate
    # cancellation in the case of coincident points.
    m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) - csig1 * csig2 * J12
    # Missing a factor of _b
    s12b = (1 + A1m1) * sig12 + AB1
    if scalep:

area the area (counter-clockwise traversal positive)

    There is no need to "close" the polygon.  If polyline is set to
    True, then the points define a polyline instead of a polygon, the
    length is returned as the perimeter, and the area is not returned.
    """

    from geographiclib.polygonarea import PolygonArea
    for p in points:
      Geodesic.CheckPosition(p['lat'], p['lon'])
    num, perimeter, area = PolygonArea.Area(self, points, polyline)
    result = {'number': num, 'perimeter': perimeter}
    if not polyline: result['area'] = area
    return result

Geodesic.WGS84 = Geodesic(Constants.WGS84_a, Constants.WGS84_f)

#
        # The histogram of iterations is (m = number of iterations estimating
        # alp1 directly, n = number of iterations estimating via omg12, total
        # number of trials = 148605):
        #
        #  iter    m      n
        #    0   148    186
        #    1 13046  13845
        #    2 93315 102225
        #    3 36189  32341
        #    4  5396      7
        #    5   455      1
        #    6    56      0
        #
        # Because omg12 is near pi, estimate work with omg12a = pi - omg12
        k = Geodesic.Astroid(x, y)
        omg12a = lamscale * ( -x * k/(1 + k) if self._f >= 0
                               else -y * (1 + k)/k )
        somg12 = math.sin(omg12a); comg12 = -math.cos(omg12a)
        # Update spherical estimate of alp1 using omg12 instead of lam12
        salp1 = cbet2 * somg12
        calp1 = sbet12a - cbet2 * sbet1 * Math.sq(somg12) / (1 - comg12)
    if salp1 > 0:               # Sanity check on starting guess
      salp1, calp1 = Geodesic.SinCosNorm(salp1, calp1)
    else:
      salp1 = 1; calp1 = 0
    return sig12, salp1, calp1, salp2, calp2, dnm

# instead of sqrt(epsilon).
              tripn = abs(v) <= 16 * Geodesic.tol0_
              continue
          # Either dv was not postive or updated value was outside legal range.
          # Use the midpoint of the bracket as the next estimate.  This
          # mechanism is not needed for the WGS84 ellipsoid, but it does catch
          # problems with more eccentric ellipsoids.  Its efficacy is such for
          # the WGS84 test set with the starting guess set to alp1 = 90deg:
          # the WGS84 test set: mean = 5.21, sd = 3.93, max = 24
          # WGS84 and random input: mean = 4.74, sd = 0.99
          salp1 = (salp1a + salp1b)/2
          calp1 = (calp1a + calp1b)/2
          salp1, calp1 = Geodesic.SinCosNorm(salp1, calp1)
          tripn = False
          tripb = (abs(salp1a - salp1) + (calp1a - calp1) < Geodesic.tolb_ or
                   abs(salp1 - salp1b) + (calp1 - calp1b) < Geodesic.tolb_)

        s12x, m12x, dummy, M12, M21 = self.Lengths(
          eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2,
          (outmask & Geodesic.GEODESICSCALE) != 0, C1a, C2a)

        m12x *= self._b
        s12x *= self._b
        a12 = sig12 / Math.degree
        omg12 = lam12 - omg12
    # end elif not meridian

    if outmask & Geodesic.DISTANCE:
      s12 = 0 + s12x           # Convert -0 to 0

    if outmask & Geodesic.REDUCEDLENGTH:
      m12 = 0 + m12x           # Convert -0 to 0

# authalic radius squared
    self._c2 = (Math.sq(self._a) + Math.sq(self._b) *
                (1 if self._e2 == 0 else
                 (Math.atanh(math.sqrt(self._e2)) if self._e2 > 0 else
                  math.atan(math.sqrt(-self._e2))) /
                 math.sqrt(abs(self._e2))))/2
    # The sig12 threshold for "really short".  Using the auxiliary sphere
    # solution with dnm computed at (bet1 + bet2) / 2, the relative error in
    # the azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2.
    # (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000.  For a given
    # f and sig12, the max error occurs for lines near the pole.  If the old
    # rule for computing dnm = (dn1 + dn2)/2 is used, then the error increases
    # by a factor of 2.)  Setting this equal to epsilon gives sig12 = etol2.
    # Here 0.1 is a safety factor (error decreased by 100) and max(0.001,
    # abs(f)) stops etol2 getting too large in the nearly spherical case.
    self._etol2 = 0.1 * Geodesic.tol2_ / math.sqrt( max(0.001, abs(self._f)) *
                                                    min(1.0, 1-self._f/2) / 2 )
    if not(Math.isfinite(self._a) and self._a > 0):
      raise ValueError("Major radius is not positive")
    if not(Math.isfinite(self._b) and self._b > 0):
      raise ValueError("Minor radius is not positive")
    self._A3x = list(range(Geodesic.nA3x_))
    self._C3x = list(range(Geodesic.nC3x_))
    self._C4x = list(range(Geodesic.nC4x_))
    self.A3coeff()
    self.C3coeff()
    self.C4coeff()