Python geomath.Math类代码示例

 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

 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)

 def C4coeff(self):
   """Private: return coefficients for C4"""
   coeff = [
     97, 15015,
     1088, 156, 45045,
     -224, -4784, 1573, 45045,
     -10656, 14144, -4576, -858, 45045,
     64, 624, -4576, 6864, -3003, 15015,
     100, 208, 572, 3432, -12012, 30030, 45045,
     1, 9009,
     -2944, 468, 135135,
     5792, 1040, -1287, 135135,
     5952, -11648, 9152, -2574, 135135,
     -64, -624, 4576, -6864, 3003, 135135,
     8, 10725,
     1856, -936, 225225,
     -8448, 4992, -1144, 225225,
     -1440, 4160, -4576, 1716, 225225,
     -136, 63063,
     1024, -208, 105105,
     3584, -3328, 1144, 315315,
     -128, 135135,
     -2560, 832, 405405,
     128, 99099,
   ]
   o = 0; k = 0
   for l in range(Geodesic.nC4_): # l is index of C4[l]
     for j in range(Geodesic.nC4_ - 1, l - 1, -1): # coeff of eps^j
       m = Geodesic.nC4_ - j - 1 # order of polynomial in n
       self._C4x[k] = Math.polyval(m, coeff, o, self._n) / coeff[o + m + 1]
       k += 1
       o += m + 2

 def C3coeff(self):
   """Private: return coefficients for C3"""
   coeff = [
     3, 128,
     2, 5, 128,
     -1, 3, 3, 64,
     -1, 0, 1, 8,
     -1, 1, 4,
     5, 256,
     1, 3, 128,
     -3, -2, 3, 64,
     1, -3, 2, 32,
     7, 512,
     -10, 9, 384,
     5, -9, 5, 192,
     7, 512,
     -14, 7, 512,
     21, 2560,
   ]
   o = 0; k = 0
   for l in range(1, Geodesic.nC3_): # l is index of C3[l]
     for j in range(Geodesic.nC3_ - 1, l - 1, -1): # coeff of eps^j
       m = min(Geodesic.nC3_ - j - 1, j) # order of polynomial in n
       self._C3x[k] = Math.polyval(m, coeff, o, self._n) / coeff[o + m + 1]
       k += 1
       o += m + 2

 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

 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

 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

 def C4f(self, eps, c):
   """Private: return C4"""
   # Evaluate C4 coeffs by Horner's method
   # Elements c[0] thru c[nC4_ - 1] are set
   mult = 1
   o = 0
   for l in range(Geodesic.nC4_): # l is index of C4[l]
     m = Geodesic.nC4_ - l - 1    # order of polynomial in eps
     c[l] = mult * Math.polyval(m, self._C4x, o, eps)
     o += m + 1
     mult *= eps

 def C3f(self, eps, c):
   """Private: return C3"""
   # Evaluate C3
   # Elements c[1] thru c[nC3_ - 1] are set
   mult = 1
   o = 0
   for l in range(1, Geodesic.nC3_): # l is index of C3[l]
     m = Geodesic.nC3_ - l - 1       # order of polynomial in eps
     mult *= eps
     c[l] = mult * Math.polyval(m, self._C3x, o, eps)
     o += m + 1

  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()

 def Add(self, y):
   """Add a value"""
   # Here's Shewchuk's solution...
   # hold exact sum as [s, t, u]
   y, u = Math.sum(y, self._t)             # Accumulate starting at
   self._s, self._t = Math.sum(y, self._s) # least significant end
   # Start is _s, _t decreasing and non-adjacent.  Sum is now (s + t + u)
   # exactly with s, t, u non-adjacent and in decreasing order (except
   # for possible zeros).  The following code tries to normalize the
   # result.  Ideally, we want _s = round(s+t+u) and _u = round(s+t+u -
   # _s).  The follow does an approximate job (and maintains the
   # decreasing non-adjacent property).  Here are two "failures" using
   # 3-bit floats:
   #
   # Case 1: _s is not equal to round(s+t+u) -- off by 1 ulp
   # [12, -1] - 8 -> [4, 0, -1] -> [4, -1] = 3 should be [3, 0] = 3
   #
   # Case 2: _s+_t is not as close to s+t+u as it shold be
   # [64, 5] + 4 -> [64, 8, 1] -> [64,  8] = 72 (off by 1)
   #                    should be [80, -7] = 73 (exact)
   #
   # "Fixing" these problems is probably not worth the expense.  The
   # representation inevitably leads to small errors in the accumulated
   # values.  The additional errors illustrated here amount to 1 ulp of
   # the less significant word during each addition to the Accumulator
   # and an additional possible error of 1 ulp in the reported sum.
   #
   # Incidentally, the "ideal" representation described above is not
   # canonical, because _s = round(_s + _t) may not be true.  For
   # example, with 3-bit floats:
   #
   # [128, 16] + 1 -> [160, -16] -- 160 = round(145).
   # But [160, 0] - 16 -> [128, 16] -- 128 = round(144).
   #
   if self._s == 0:            # This implies t == 0,
     self._s = u               # so result is u
   else:
     self._t += u              # otherwise just accumulate u to t.

 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

 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

 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

 def A3coeff(self):
   """Private: return coefficients for A3"""
   coeff = [
     -3, 128,
     -2, -3, 64,
     -1, -3, -1, 16,
     3, -1, -2, 8,
     1, -1, 2,
     1, 1,
   ]
   o = 0; k = 0
   for j in range(Geodesic.nA3_ - 1, -1, -1): # coeff of eps^j
     m = min(Geodesic.nA3_ - j - 1, j) # order of polynomial in n
     self._A3x[k] = Math.polyval(m, coeff, o, self._n) / coeff[o + m + 1]
     k += 1
     o += m + 2