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

#
    # longsign, swapp, latsign register the transformation to bring the
    # coordinates to this canonical form.  In all cases, 1 means no change was
    # made.  We make these transformations so that there are few cases to
    # check, e.g., on verifying quadrants in atan2.  In addition, this
    # enforces some symmetries in the results returned.

    # real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x

    sbet1, cbet1 = Math.sincosd(lat1); sbet1 *= self._f1
    # Ensure cbet1 = +epsilon at poles
    sbet1, cbet1 = Math.norm(sbet1, cbet1); cbet1 = max(Geodesic.tiny_, cbet1)

    sbet2, cbet2 = Math.sincosd(lat2); sbet2 *= self._f1
    # Ensure cbet2 = +epsilon at poles
    sbet2, cbet2 = Math.norm(sbet2, cbet2); cbet2 = max(Geodesic.tiny_, cbet2)

    # If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the
    # |bet1| - |bet2|.  Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is
    # a better measure.  This logic is used in assigning calp2 in Lambda12.
    # Sometimes these quantities vanish and in that case we force bet2 = +/-
    # bet1 exactly.  An example where is is necessary is the inverse problem
    # 48.522876735459 0 -48.52287673545898293 179.599720456223079643
    # which failed with Visual Studio 10 (Release and Debug)

    if cbet1 < -sbet1:
      if cbet2 == cbet1:
        sbet2 = sbet1 if sbet2 < 0 else -sbet1
    else:
      if abs(sbet2) == -sbet1:
        cbet2 = cbet1

if v > 0 and (numit > Geodesic.maxit1_ or
                        calp1/salp1 > calp1b/salp1b):
            salp1b = salp1; calp1b = calp1
          elif v < 0 and (numit > Geodesic.maxit1_ or
                          calp1/salp1 < calp1a/salp1a):
            salp1a = salp1; calp1a = calp1

          numit += 1
          if numit < Geodesic.maxit1_ and dv > 0:
            dalp1 = -v/dv
            sdalp1 = math.sin(dalp1); cdalp1 = math.cos(dalp1)
            nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
            if nsalp1 > 0 and abs(dalp1) < math.pi:
              calp1 = calp1 * cdalp1 - salp1 * sdalp1
              salp1 = nsalp1
              salp1, calp1 = Math.norm(salp1, calp1)
              # In some regimes we don't get quadratic convergence because
              # slope -> 0.  So use convergence conditions based on epsilon
              # instead of sqrt(epsilon).
              tripn = abs(v) <= 16 * Geodesic.tol0_
              continue
          # Either dv was not positive 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 = Math.norm(salp1, calp1)

"""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, 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 = Math.norm(ssig1, csig1)
    # Math.norm(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.
    # sin(alp2) * cos(bet2) = sin(alp0)
    salp2 = salp0 / cbet2 if cbet2 != cbet1 else salp1
    # calp2 = sqrt(1 - sq(salp2))
    #       = sqrt(sq(calp0) - sq(sbet2)) / cbet2
    # and subst for calp0 and rearrange to give (choose positive sqrt
    # to give alp2 in [0, pi/2]).
    calp2 = (math.sqrt(Math.sq(calp1 * cbet1) +
                       ((cbet2 - cbet1) * (cbet1 + cbet2) if cbet1 < -sbet1
                        else (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2
             if cbet2 != cbet1 or abs(sbet2) != -sbet1 else abs(calp1))
    # tan(bet2) = tan(sig2) * cos(alp2)

#
    #     0 <= lon12 <= 180
    #     -90 <= lat1 <= 0
    #     lat1 <= lat2 <= -lat1
    #
    # longsign, swapp, latsign register the transformation to bring the
    # coordinates to this canonical form.  In all cases, 1 means no change was
    # made.  We make these transformations so that there are few cases to
    # check, e.g., on verifying quadrants in atan2.  In addition, this
    # enforces some symmetries in the results returned.

    # real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x

    sbet1, cbet1 = Math.sincosd(lat1); sbet1 *= self._f1
    # Ensure cbet1 = +epsilon at poles
    sbet1, cbet1 = Math.norm(sbet1, cbet1); cbet1 = max(Geodesic.tiny_, cbet1)

    sbet2, cbet2 = Math.sincosd(lat2); sbet2 *= self._f1
    # Ensure cbet2 = +epsilon at poles
    sbet2, cbet2 = Math.norm(sbet2, cbet2); cbet2 = max(Geodesic.tiny_, cbet2)

    # If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the
    # |bet1| - |bet2|.  Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is
    # a better measure.  This logic is used in assigning calp2 in Lambda12.
    # Sometimes these quantities vanish and in that case we force bet2 = +/-
    # bet1 exactly.  An example where is is necessary is the inverse problem
    # 48.522876735459 0 -48.52287673545898293 179.599720456223079643
    # which failed with Visual Studio 10 (Release and Debug)

    if cbet1 < -sbet1:
      if cbet2 == cbet1:
        sbet2 = sbet1 if sbet2 < 0 else -sbet1

"""the longitude of the first point in degrees (readonly)"""
    if Math.isnan(salp1) or Math.isnan(calp1):
      self.azi1 = Math.AngNormalize(azi1)
      self.salp1, self.calp1 = Math.sincosd(Math.AngRound(azi1))
    else:
      self.azi1 = azi1
      """the azimuth at the first point in degrees (readonly)"""
      self.salp1 = salp1
      """the sine of the azimuth at the first point (readonly)"""
      self.calp1 = calp1
      """the cosine of the azimuth at the first point (readonly)"""

    # real cbet1, sbet1
    sbet1, cbet1 = Math.sincosd(Math.AngRound(lat1)); sbet1 *= self._f1
    # Ensure cbet1 = +epsilon at poles
    sbet1, cbet1 = Math.norm(sbet1, cbet1); cbet1 = max(Geodesic.tiny_, cbet1)
    self._dn1 = math.sqrt(1 + geod._ep2 * Math.sq(sbet1))

    # Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),
    self._salp0 = self.salp1 * cbet1 # alp0 in [0, pi/2 - |bet1|]
    # Alt: calp0 = hypot(sbet1, calp1 * cbet1).  The following
    # is slightly better (consider the case salp1 = 0).
    self._calp0 = math.hypot(self.calp1, self.salp1 * sbet1)
    # Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).
    # sig = 0 is nearest northward crossing of equator.
    # With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).
    # With bet1 =  pi/2, alp1 = -pi, sig1 =  pi/2
    # With bet1 = -pi/2, alp1 =  0 , sig1 = -pi/2
    # Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).
    # With alp0 in (0, pi/2], quadrants for sig and omg coincide.
    # No atan2(0,0) ambiguity at poles since cbet1 = +epsilon.
    # With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.

#    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)
    # Sanity check on starting guess.  Backwards check allows NaN through.
    if not (salp1 <= 0):
      salp1, calp1 = Math.norm(salp1, calp1)
    else:
      salp1 = 1; calp1 = 0
    return sig12, salp1, calp1, salp2, calp2, dnm

m12 = 0.0 + m12x          # Convert -0 to 0

    if outmask & Geodesic.AREA:
      # From Lambda12: sin(alp1) * cos(bet1) = sin(alp0)
      salp0 = salp1 * cbet1
      calp0 = math.hypot(calp1, salp1 * sbet1) # calp0 > 0
      # real alp12
      if calp0 != 0 and salp0 != 0:
        # From Lambda12: tan(bet) = tan(sig) * cos(alp)
        ssig1 = sbet1; csig1 = calp1 * cbet1
        ssig2 = sbet2; csig2 = calp2 * cbet2
        k2 = Math.sq(calp0) * self._ep2
        eps = k2 / (2 * (1 + math.sqrt(1 + k2)) + k2)
        # Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0).
        A4 = Math.sq(self.a) * calp0 * salp0 * self._e2
        ssig1, csig1 = Math.norm(ssig1, csig1)
        ssig2, csig2 = Math.norm(ssig2, csig2)
        C4a = list(range(Geodesic.nC4_))
        self._C4f(eps, C4a)
        B41 = Geodesic._SinCosSeries(False, ssig1, csig1, C4a)
        B42 = Geodesic._SinCosSeries(False, ssig2, csig2, C4a)
        S12 = A4 * (B42 - B41)
      else:
        # Avoid problems with indeterminate sig1, sig2 on equator
        S12 = 0.0

      if not meridian and somg12 > 1:
        somg12 = math.sin(omg12); comg12 = math.cos(omg12)

      if (not meridian and
          # omg12 < 3/4 * pi
          comg12 > -0.7071 and   # Long difference not too big

#    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)
    # Sanity check on starting guess.  Backwards check allows NaN through.
    if not (salp1 <= 0):
      salp1, calp1 = Math.norm(salp1, calp1)
    else:
      salp1 = 1; calp1 = 0
    return sig12, salp1, calp1, salp2, calp2, dnm

# iteration.
    # sin(alp2) * cos(bet2) = sin(alp0)
    salp2 = salp0 / cbet2 if cbet2 != cbet1 else salp1
    # calp2 = sqrt(1 - sq(salp2))
    #       = sqrt(sq(calp0) - sq(sbet2)) / cbet2
    # and subst for calp0 and rearrange to give (choose positive sqrt
    # to give alp2 in [0, pi/2]).
    calp2 = (math.sqrt(Math.sq(calp1 * cbet1) +
                       ((cbet2 - cbet1) * (cbet1 + cbet2) if cbet1 < -sbet1
                        else (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2
             if cbet2 != cbet1 or abs(sbet2) != -sbet1 else abs(calp1))
    # tan(bet2) = tan(sig2) * cos(alp2)
    # tan(omg2) = sin(alp0) * tan(sig2).
    ssig2 = sbet2; somg2 = salp0 * sbet2
    csig2 = comg2 = calp2 * cbet2
    ssig2, csig2 = Math.norm(ssig2, csig2)
    # Math.norm(somg2, comg2); -- don't need to normalize!

    # sig12 = sig2 - sig1, limit to [0, pi]
    sig12 = math.atan2(max(0.0, csig1 * ssig2 - ssig1 * csig2),
                                csig1 * csig2 + ssig1 * ssig2)

    # omg12 = omg2 - omg1, limit to [0, pi]
    somg12 = max(0.0, comg1 * somg2 - somg1 * comg2)
    comg12 =          comg1 * comg2 + somg1 * somg2
    # eta = omg12 - lam120
    eta = math.atan2(somg12 * clam120 - comg12 * slam120,
                     comg12 * clam120 + somg12 * slam120)

    # real B312
    k2 = Math.sq(calp0) * self._ep2
    eps = k2 / (2 * (1 + math.sqrt(1 + k2)) + k2)