Python Math.sq方法代码示例

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
  def __init__(self, a, f):
    """
    Construct a Geodesic object for ellipsoid with major radius a and
    flattening f.
    """

    self._a = float(a)
    self._f = float(f) if f <= 1 else 1.0/f
    self._f1 = 1 - self._f
    self._e2 = self._f * (2 - self._f)
    self._ep2 = self._e2 / Math.sq(self._f1) # e2 / (1 - e2)
    self._n = self._f / ( 2 - self._f)
    self._b = self._a * self._f1
    # 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"
    self._etol2 = 0.01 * Geodesic.tol2_ / max(0.1, math.sqrt(abs(self._e2)))
    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 = range(Geodesic.nA3x_)
    self._C3x = range(Geodesic.nC3x_)
    self._C4x = range(Geodesic.nC4x_)
    self.A3coeff()
    self.C3coeff()
    self.C4coeff()

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
 def A1m1f(eps):
   """Private: return A1-1."""
   coeff = [
     1, 4, 64, 0, 256,
   ]
   m = Geodesic.nA1_//2
   t = Math.polyval(m, coeff, 0, Math.sq(eps)) / coeff[m + 1]
   return (t + eps) / (1 - eps)

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
 def A2m1f(eps):
   """Private: return A2-1"""
   coeff = [
     25, 36, 64, 0, 256,
   ]
   m = Geodesic.nA2_//2
   t = Math.polyval(m, coeff, 0, Math.sq(eps)) / coeff[m + 1]
   return t * (1 - eps) - eps

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
 def Astroid(x, y):
   """Private: solve astroid equation."""
   # Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k.
   # This solution is adapted from Geocentric::Reverse.
   p = Math.sq(x)
   q = Math.sq(y)
   r = (p + q - 1) / 6
   if not(q == 0 and r <= 0):
     # Avoid possible division by zero when r = 0 by multiplying equations
     # 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
     # y = 0 with |x| <= 1.  Handle this case directly.
     # for y small, positive root is k = abs(y)/sqrt(1-x^2)
     k = 0
   return k

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
  def __init__(self, a, f):
    """Construct a Geodesic object for ellipsoid with major radius a and
    flattening f.

    """

    self._a = float(a)
    self._f = float(f) if f <= 1 else 1.0/f
    self._f1 = 1 - self._f
    self._e2 = self._f * (2 - self._f)
    self._ep2 = self._e2 / Math.sq(self._f1) # e2 / (1 - e2)
    self._n = self._f / ( 2 - self._f)
    self._b = self._a * self._f1
    # 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()

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
 def C2f(eps, c):
   eps2 = Math.sq(eps)
   d = eps
   c[1] = d*(eps2*(eps2+2)+16)/32
   d *= eps
   c[2] = d*(eps2*(35*eps2+64)+384)/2048
   d *= eps
   c[3] = d*(15*eps2+80)/768
   d *= eps
   c[4] = d*(7*eps2+35)/512
   d *= eps
   c[5] = 63*d/1280
   d *= eps
   c[6] = 77*d/2048

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
 def C1pf(eps, c):
   eps2 = Math.sq(eps)
   d = eps
   c[1] = d*(eps2*(205*eps2-432)+768)/1536
   d *= eps
   c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288
   d *= eps
   c[3] = d*(116-225*eps2)/384
   d *= eps
   c[4] = d*(2695-7173*eps2)/7680
   d *= eps
   c[5] = 3467*d/7680
   d *= eps
   c[6] = 38081*d/61440

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
 def C1f(eps, c):
   eps2 = Math.sq(eps)
   d = eps
   c[1] = d*((6-eps2)*eps2-16)/32
   d *= eps
   c[2] = d*((64-9*eps2)*eps2-128)/2048
   d *= eps
   c[3] = d*(9*eps2-16)/768
   d *= eps
   c[4] = d*(3*eps2-5)/512
   d *= eps
   c[5] = -7*d/1280
   d *= eps
   c[6] = -7*d/2048

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
 def C1pf(eps, c):
   """Private: return C1'"""
   coeff = [
     205, -432, 768, 1536,
     4005, -4736, 3840, 12288,
     -225, 116, 384,
     -7173, 2695, 7680,
     3467, 7680,
     38081, 61440,
   ]
   eps2 = Math.sq(eps)
   d = eps
   o = 0
   for l in range(1, Geodesic.nC1p_ + 1): # l is index of C1p[l]
     m = (Geodesic.nC1p_ - l) // 2 # order of polynomial in eps^2
     c[l] = d * Math.polyval(m, coeff, o, eps2) / coeff[o + m + 1]
     o += m + 2
     d *= eps

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
 def C1f(eps, c):
   """Private: return C1."""
   coeff = [
     -1, 6, -16, 32,
     -9, 64, -128, 2048,
     9, -16, 768,
     3, -5, 512,
     -7, 1280,
     -7, 2048,
   ]
   eps2 = Math.sq(eps)
   d = eps
   o = 0
   for l in range(1, Geodesic.nC1_ + 1): # l is index of C1p[l]
     m = (Geodesic.nC1_ - l) // 2        # order of polynomial in eps^2
     c[l] = d * Math.polyval(m, coeff, o, eps2) / coeff[o + m + 1]
     o += m + 2
     d *= eps

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
 def C2f(eps, c):
   """Private: return C2"""
   coeff = [
     1, 2, 16, 32,
     35, 64, 384, 2048,
     15, 80, 768,
     7, 35, 512,
     63, 1280,
     77, 2048,
   ]
   eps2 = Math.sq(eps)
   d = eps
   o = 0
   for l in range(1, Geodesic.nC2_ + 1): # l is index of C2[l]
     m = (Geodesic.nC2_ - l) // 2        # order of polynomial in eps^2
     c[l] = d * Math.polyval(m, coeff, o, eps2) / coeff[o + m + 1]
     o += m + 2
     d *= eps

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
  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.
    # 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 = Geodesic.SinCosNorm(ssig2, csig2)
    # SinCosNorm(somg2, comg2); -- don't need to normalize!

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

    # omg12 = omg2 - omg1, limit to [0, pi]
    omg12 = math.atan2(max(comg1 * somg2 - somg1 * comg2, 0.0),
                       comg1 * comg2 + somg1 * somg2)
    # real B312, h0
    k2 = Math.sq(calp0) * self._ep2
    eps = k2 / (2 * (1 + math.sqrt(1 + k2)) + k2)
    self.C3f(eps, C3a)
    B312 = (Geodesic.SinCosSeries(True, ssig2, csig2, C3a, Geodesic.nC3_-1) -
            Geodesic.SinCosSeries(True, ssig1, csig1, C3a, Geodesic.nC3_-1))
    h0 = -self._f * self.A3f(eps)
    domg12 = salp0 * h0 * (sig12 + B312)
    lam12 = omg12 + domg12

    if diffp:
      if calp2 == 0:
        dlam12 = - 2 * self._f1 * dn1 / sbet1
      else:
        dummy, dlam12, dummy, dummy, dummy = self.Lengths(
          eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2,
          False, C1a, C2a)
        dlam12 *= self._f1 / (calp2 * cbet2)
    else:
      dlam12 = Math.nan

    return (lam12, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2, eps,
            domg12, dlam12)

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
  def GenInverse(self, lat1, lon1, lat2, lon2, outmask):
    """Private: General version of the inverse problem"""
    a12 = s12 = azi1 = azi2 = m12 = M12 = M21 = S12 = Math.nan # return vals

    outmask &= Geodesic.OUT_MASK
    # Compute longitude difference (AngDiff does this carefully).  Result is
    # in [-180, 180] but -180 is only for west-going geodesics.  180 is for
    # east-going and meridional geodesics.
    lon12 = Math.AngDiff(Math.AngNormalize(lon1), Math.AngNormalize(lon2))
    # If very close to being on the same half-meridian, then make it so.
    lon12 = Geodesic.AngRound(lon12)
    # Make longitude difference positive.
    lonsign = 1 if lon12 >= 0 else -1
    lon12 *= lonsign
    # If really close to the equator, treat as on equator.
    lat1 = Geodesic.AngRound(lat1)
    lat2 = Geodesic.AngRound(lat2)
    # Swap points so that point with higher (abs) latitude is point 1
    swapp = 1 if abs(lat1) >= abs(lat2) else -1
    if swapp < 0:
      lonsign *= -1
      lat2, lat1 = lat1, lat2
    # Make lat1 <= 0
    latsign = 1 if lat1 < 0 else -1
    lat1 *= latsign
    lat2 *= latsign
    # Now we have
    #
    #     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

    phi = lat1 * Math.degree
    # Ensure cbet1 = +epsilon at poles
    sbet1 = self._f1 * math.sin(phi)
    cbet1 = Geodesic.tiny_ if lat1 == -90 else math.cos(phi)
    sbet1, cbet1 = Geodesic.SinCosNorm(sbet1, cbet1)

    phi = lat2 * Math.degree
    # Ensure cbet2 = +epsilon at poles
    sbet2 = self._f1 * math.sin(phi)
    cbet2 = Geodesic.tiny_ if abs(lat2) == 90 else math.cos(phi)
    sbet2, cbet2 = Geodesic.SinCosNorm(sbet2, 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

    dn1 = math.sqrt(1 + self._ep2 * Math.sq(sbet1))
    dn2 = math.sqrt(1 + self._ep2 * Math.sq(sbet2))

    lam12 = lon12 * Math.degree
    slam12 = 0 if lon12 == 180 else math.sin(lam12)
    clam12 = math.cos(lam12)      # lon12 == 90 isn't interesting

    # real a12, sig12, calp1, salp1, calp2, salp2
    # index zero elements of these arrays are unused
    C1a = list(range(Geodesic.nC1_ + 1))
    C2a = list(range(Geodesic.nC2_ + 1))
    C3a = list(range(Geodesic.nC3_))

    meridian = lat1 == -90 or slam12 == 0

    if meridian:

      # Endpoints are on a single full meridian, so the geodesic might lie on
      # a meridian.

      calp1 = clam12; salp1 = slam12 # Head to the target longitude
      calp2 = 1; salp2 = 0           # At the target we're heading north

      # tan(bet) = tan(sig) * cos(alp)
      ssig1 = sbet1; csig1 = calp1 * cbet1
      ssig2 = sbet2; csig2 = calp2 * cbet2

      # sig12 = sig2 - sig1
      sig12 = math.atan2(max(csig1 * ssig2 - ssig1 * csig2, 0.0),
                         csig1 * csig2 + ssig1 * ssig2)

      s12x, m12x, dummy, M12, M21 = self.Lengths(
        self._n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2,
#.........这里部分代码省略.........

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
 def A2m1f(eps):
   eps2 = Math.sq(eps)
   t = eps2*(eps2*(25*eps2+36)+64)/256
   return t * (1 - eps) - eps

# 需要导入模块: from geographiclib.geomath import Math [as 别名]
# 或者: from geographiclib.geomath.Math import sq [as 别名]
  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
      # =  (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2)
      sbetm2 /= sbetm2 + Math.sq(cbet1 + cbet2)
      dnm = math.sqrt(1 + self._ep2 * sbetm2)
      omg12 /= self._f1 * dnm
    somg12 = math.sin(omg12); comg12 = math.cos(omg12)

    salp1 = cbet2 * somg12
    calp1 = (
      sbet12 + cbet2 * sbet1 * Math.sq(somg12) / (1 + comg12) if comg12 >= 0
      else sbet12a - cbet2 * sbet1 * Math.sq(somg12) / (1 - comg12))

    ssig12 = math.hypot(salp1, calp1)
    csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12

    if shortline and ssig12 < self._etol2:
      # really short lines
      salp2 = cbet1 * somg12
      calp2 = sbet12 - cbet1 * sbet2 * (Math.sq(somg12) / (1 + comg12)
                                        if comg12 >= 0 else 1 - comg12)
      salp2, calp2 = Geodesic.SinCosNorm(salp2, calp2)
      # Set return value
      sig12 = math.atan2(ssig12, csig12)
    elif (abs(self._n) >= 0.1 or # Skip astroid calc if too eccentric
          csig12 >= 0 or
          ssig12 >= 6 * abs(self._n) * math.pi * Math.sq(cbet1)):
      # Nothing to do, zeroth order spherical approximation is OK
      pass
    else:
      # Scale lam12 and bet2 to x, y coordinate system where antipodal point
      # is at origin and singular point is at y = 0, x = -1.
      # real y, lamscale, betscale
      # Volatile declaration needed to fix inverse case
      # 56.320923501171 0 -56.320923501171 179.664747671772880215
      # which otherwise fails with g++ 4.4.4 x86 -O3
      # volatile real x
      if self._f >= 0:            # In fact f == 0 does not get here
        # x = dlong, y = dlat
        k2 = Math.sq(sbet1) * self._ep2
        eps = k2 / (2 * (1 + math.sqrt(1 + k2)) + k2)
        lamscale = self._f * cbet1 * self.A3f(eps) * math.pi
        betscale = lamscale * cbet1
        x = (lam12 - math.pi) / lamscale
        y = sbet12a / betscale
      else:                     # _f < 0
        # x = dlat, y = dlong
        cbet12a = cbet2 * cbet1 - sbet2 * sbet1
        bet12a = math.atan2(sbet12a, cbet12a)
        # real m12b, m0, dummy
        # In the case of lon12 = 180, this repeats a calculation made in
        # Inverse.
        dummy, m12b, m0, dummy, dummy = self.Lengths(
          self._n, math.pi + bet12a, sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
          cbet1, cbet2, False, C1a, C2a)
        x = -1 + m12b / (cbet1 * cbet2 * m0 * math.pi)
        betscale = (sbet12a / x if x < -0.01
                    else -self._f * Math.sq(cbet1) * math.pi)
        lamscale = betscale / cbet1
        y = (lam12 - math.pi) / lamscale

      if y > -Geodesic.tol1_ and x > -1 - Geodesic.xthresh_:
        # strip near cut
        if self._f >= 0:
          salp1 = min(1.0, -x); calp1 = - math.sqrt(1 - Math.sq(salp1))
        else:
          calp1 = max((0.0 if x > -Geodesic.tol1_ else -1.0), x)
          salp1 = math.sqrt(1 - Math.sq(calp1))
      else:
        # Estimate alp1, by solving the astroid problem.
        #
        # Could estimate alpha1 = theta + pi/2, directly, i.e.,
        #   calp1 = y/k; salp1 = -x/(1+k);  for _f >= 0
        #   calp1 = x/(1+k); salp1 = -y/k;  for _f < 0 (need to check)
        #
        # However, it's better to estimate omg12 from astroid and use
        # spherical formula to compute alp1.  This reduces the mean number of
        # Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12
#.........这里部分代码省略.........