How to use the geographiclib.geomath.Math function in geographiclib

"""Solve geodesic problems"""

  GEOGRAPHICLIB_GEODESIC_ORDER = 6
  nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nA3x_ = nA3_
  nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nC3x_ = (nC3_ * (nC3_ - 1)) // 2
  nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER
  nC4x_ = (nC4_ * (nC4_ + 1)) // 2
  maxit1_ = 20
  maxit2_ = maxit1_ + Math.digits + 10

  tiny_ = math.sqrt(Math.minval)
  tol0_ = Math.epsilon
  tol1_ = 200 * tol0_
  tol2_ = math.sqrt(tol0_)
  tolb_ = tol0_ * tol2_
  xthresh_ = 1000 * tol2_

  CAP_NONE = GeodesicCapability.CAP_NONE
  CAP_C1   = GeodesicCapability.CAP_C1
  CAP_C1p  = GeodesicCapability.CAP_C1p
  CAP_C2   = GeodesicCapability.CAP_C2
  CAP_C3   = GeodesicCapability.CAP_C3
  CAP_C4   = GeodesicCapability.CAP_C4
  CAP_ALL  = GeodesicCapability.CAP_ALL
  CAP_MASK = GeodesicCapability.CAP_MASK

E = Math.copysign(1, self._salp0)          # East or west going?
      # omg12 = omg2 - omg1
      omg12 = (E * (sig12
                    - (math.atan2(          ssig2,       csig2) -
                       math.atan2(    self._ssig1, self._csig1))
                    + (math.atan2(E *       somg2,       comg2) -
                       math.atan2(E * self._somg1, self._comg1)))
               if outmask & Geodesic.LONG_UNROLL
               else math.atan2(somg2 * self._comg1 - comg2 * self._somg1,
                               comg2 * self._comg1 + somg2 * self._somg1))
      lam12 = omg12 + self._A3c * (
        sig12 + (Geodesic._SinCosSeries(True, ssig2, csig2, self._C3a)
                 - self._B31))
      lon12 = math.degrees(lam12)
      lon2 = (self.lon1 + lon12 if outmask & Geodesic.LONG_UNROLL else
              Math.AngNormalize(Math.AngNormalize(self.lon1) +
                                Math.AngNormalize(lon12)))

    if outmask & Geodesic.LATITUDE:
      lat2 = Math.atan2d(sbet2, self._f1 * cbet2)

    if outmask & Geodesic.AZIMUTH:
      azi2 = Math.atan2d(salp2, calp2)

    if outmask & (Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
      B22 = Geodesic._SinCosSeries(True, ssig2, csig2, self._C2a)
      AB2 = (1 + self._A2m1) * (B22 - self._B21)
      J12 = (self._A1m1 - self._A2m1) * sig12 + (AB1 - AB2)
      if outmask & Geodesic.REDUCEDLENGTH:
        # Add parens around (_csig1 * ssig2) and (_ssig1 * csig2) to ensure
        # accurate cancellation in the case of coincident points.
        m12 = self._b * ((      dn2 * (self._csig1 * ssig2) -

#      1/5   157e6 3.8e9 280e6
        ssig2 = self._ssig1 * csig12 + self._csig1 * ssig12
        csig2 = self._csig1 * csig12 - self._ssig1 * ssig12
        B12 = Geodesic._SinCosSeries(True, ssig2, csig2, self._C1a)
        serr = ((1 + self._A1m1) * (sig12 + (B12 - self._B11)) -
                s12_a12 / self._b)
        sig12 = sig12 - serr / math.sqrt(1 + self._k2 * Math.sq(ssig2))
        ssig12 = math.sin(sig12); csig12 = math.cos(sig12)
        # Update B12 below

    # real omg12, lam12, lon12
    # real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2
    # sig2 = sig1 + sig12
    ssig2 = self._ssig1 * csig12 + self._csig1 * ssig12
    csig2 = self._csig1 * csig12 - self._ssig1 * ssig12
    dn2 = math.sqrt(1 + self._k2 * Math.sq(ssig2))
    if outmask & (
      Geodesic.DISTANCE | Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
      if arcmode or abs(self.f) > 0.01:
        B12 = Geodesic._SinCosSeries(True, ssig2, csig2, self._C1a)
      AB1 = (1 + self._A1m1) * (B12 - self._B11)
    # sin(bet2) = cos(alp0) * sin(sig2)
    sbet2 = self._calp0 * ssig2
    # Alt: cbet2 = hypot(csig2, salp0 * ssig2)
    cbet2 = math.hypot(self._salp0, self._calp0 * csig2)
    if cbet2 == 0:
      # I.e., salp0 = 0, csig2 = 0.  Break the degeneracy in this case
      cbet2 = csig2 = Geodesic.tiny_
    # tan(alp0) = cos(sig2)*tan(alp2)
    salp2 = self._salp0; calp2 = self._calp0 * csig2 # No need to normalize

    if outmask & Geodesic.DISTANCE:

# for s and t by r^3 and r, resp.
      S = p * q / 4            # S = r^3 * s
      r2 = Math.sq(r)
      r3 = r * r2
      # The discriminant of the quadratic equation for T3.  This is zero on
      # the evolute curve p^(1/3)+q^(1/3) = 1
      disc = S * (S + 2 * r3)
      u = r
      if disc >= 0:
        T3 = S + r3
        # Pick the sign on the sqrt to maximize abs(T3).  This minimizes loss
        # of precision due to cancellation.  The result is unchanged because
        # of the way the T is used in definition of u.
        T3 += -math.sqrt(disc) if T3 < 0 else math.sqrt(disc) # T3 = (r * t)^3
        # N.B. cbrt always returns the real root.  cbrt(-8) = -2.
        T = Math.cbrt(T3)       # T = r * t
        # T can be zero; but then r2 / T -> 0.
        u += T + (r2 / T if T != 0 else 0)
      else:
        # T is complex, but the way u is defined the result is real.
        ang = math.atan2(math.sqrt(-disc), -(S + r3))
        # There are three possible cube roots.  We choose the root which
        # avoids cancellation.  Note that disc < 0 implies that r < 0.
        u += 2 * r * math.cos(ang / 3)
      v = math.sqrt(Math.sq(u) + q) # guaranteed positive
      # Avoid loss of accuracy when u < 0.
      uv = q / (v - u) if u < 0 else u + v # u+v, guaranteed positive
      w = (uv - q) / (2 * v)               # positive?
      # Rearrange expression for k to avoid loss of accuracy due to
      # subtraction.  Division by 0 not possible because uv > 0, w >= 0.
      k = uv / (math.sqrt(uv + Math.sq(w)) + w) # guaranteed positive
    else:                                       # q == 0 && r <= 0

# for s and t by r^3 and r, resp.
      S = p * q / 4            # S = r^3 * s
      r2 = Math.sq(r)
      r3 = r * r2
      # The discrimant of the quadratic equation for T3.  This is zero on
      # the evolute curve p^(1/3)+q^(1/3) = 1
      disc = S * (S + 2 * r3)
      u = r
      if (disc >= 0):
        T3 = S + r3
        # Pick the sign on the sqrt to maximize abs(T3).  This minimizes loss
        # of precision due to cancellation.  The result is unchanged because
        # of the way the T is used in definition of u.
        T3 += -math.sqrt(disc) if T3 < 0 else math.sqrt(disc) # T3 = (r * t)^3
        # N.B. cbrt always returns the real root.  cbrt(-8) = -2.
        T = Math.cbrt(T3)       # T = r * t
        # T can be zero; but then r2 / T -> 0.
        u += T + (r2 / T if T != 0 else 0)
      else:
        # T is complex, but the way u is defined the result is real.
        ang = math.atan2(math.sqrt(-disc), -(S + r3))
        # There are three possible cube roots.  We choose the root which
        # avoids cancellation.  Note that disc < 0 implies that r < 0.
        u += 2 * r * math.cos(ang / 3)
      v = math.sqrt(Math.sq(u) + q) # guaranteed positive
      # Avoid loss of accuracy when u < 0.
      uv = q / (v - u) if u < 0 else u + v # u+v, guaranteed positive
      w = (uv - q) / (2 * v)               # positive?
      # Rearrange expression for k to avoid loss of accuracy due to
      # subtraction.  Division by 0 not possible because uv > 0, w >= 0.
      k = uv / (math.sqrt(uv + Math.sq(w)) + w) # guaranteed positive
    else:                                       # q == 0 && r <= 0

def InverseStart(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
                   # Scratch areas of the right size
                   C1a, C2a):
    """Private: Find a starting value for Newton's method."""
    # Return a starting point for Newton's method in salp1 and calp1 (function
    # value is -1).  If Newton's method doesn't need to be used, return also
    # salp2 and calp2 and function value is sig12.
    sig12 = -1; salp2 = calp2 = dnm = Math.nan # Return values
    # bet12 = bet2 - bet1 in [0, pi); bet12a = bet2 + bet1 in (-pi, 0]
    sbet12 = sbet2 * cbet1 - cbet2 * sbet1
    cbet12 = cbet2 * cbet1 + sbet2 * sbet1
    # Volatile declaration needed to fix inverse cases
    # 88.202499451857 0 -88.202499451857 179.981022032992859592
    # 89.262080389218 0 -89.262080389218 179.992207982775375662
    # 89.333123580033 0 -89.333123580032997687 179.99295812360148422
    # which otherwise fail with g++ 4.4.4 x86 -O3
    sbet12a = sbet2 * cbet1
    sbet12a += cbet2 * sbet1

    shortline = cbet12 >= 0 and sbet12 < 0.5 and cbet2 * lam12 < 0.5
    omg12 = lam12
    if shortline:
      sbetm2 = Math.sq(sbet1 + sbet2)
      # sin((bet1+bet2)/2)^2

# for s and t by r^3 and r, resp.
      S = p * q / 4            # S = r^3 * s
      r2 = Math.sq(r)
      r3 = r * r2
      # The discriminant of the quadratic equation for T3.  This is zero on
      # the evolute curve p^(1/3)+q^(1/3) = 1
      disc = S * (S + 2 * r3)
      u = r
      if disc >= 0:
        T3 = S + r3
        # Pick the sign on the sqrt to maximize abs(T3).  This minimizes loss
        # of precision due to cancellation.  The result is unchanged because
        # of the way the T is used in definition of u.
        T3 += -math.sqrt(disc) if T3 < 0 else math.sqrt(disc) # T3 = (r * t)^3
        # N.B. cbrt always returns the real root.  cbrt(-8) = -2.
        T = Math.cbrt(T3)       # T = r * t
        # T can be zero; but then r2 / T -> 0.
        u += T + (r2 / T if T != 0 else 0)
      else:
        # T is complex, but the way u is defined the result is real.
        ang = math.atan2(math.sqrt(-disc), -(S + r3))
        # There are three possible cube roots.  We choose the root which
        # avoids cancellation.  Note that disc < 0 implies that r < 0.
        u += 2 * r * math.cos(ang / 3)
      v = math.sqrt(Math.sq(u) + q) # guaranteed positive
      # Avoid loss of accuracy when u < 0.
      uv = q / (v - u) if u < 0 else u + v # u+v, guaranteed positive
      w = (uv - q) / (2 * v)               # positive?
      # Rearrange expression for k to avoid loss of accuracy due to
      # subtraction.  Division by 0 not possible because uv > 0, w >= 0.
      k = uv / (math.sqrt(uv + Math.sq(w)) + w) # guaranteed positive
    else:                                       # q == 0 && r <= 0