Coverage for pygeodesy/ktm.py: 98%
215 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-05-20 11:54 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-05-20 11:54 -0400
2# -*- coding: utf-8 -*-
4u'''A pure Python version of I{Karney}'s C++ class U{TransverseMercator
5<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1TransverseMercator.html>}
6based on I{Krüger} series. See also I{Karney}'s utility U{TransverseMercatorProj
7<https://GeographicLib.SourceForge.io/C++/doc/TransverseMercatorProj.1.html>}.
9Following and further below is a copy of I{Karney}'s U{TransverseMercator.hpp
10<https://GeographicLib.SourceForge.io/C++/doc/TransverseMercator_8hpp_source.html>}
11file C{Header}.
13This implementation follows closely JHS 154, ETRS89 - I{järjestelmään liittyvät
14karttaprojektiot, tasokoordinaatistot ja karttalehtijako} (Map projections, plane
15coordinates, and map sheet index for ETRS89), published by JUHTA, Finnish Geodetic
16Institute, and the National Land Survey of Finland (2006). The relevant section
17is available as the U{2008 PDF file
18<http://Docs.JHS-suositukset.FI/jhs-suositukset/JHS154/JHS154_liite1.pdf>}.
20This is a straight transcription of the formulas in this paper with the
21following exceptions:
23 - Use of 6th order series instead of 4th order series. This reduces the
24 error to about 5 nm for the UTM range of coordinates (instead of 200 nm),
25 with a speed penalty of only 1%,
27 - Use Newton's method instead of plain iteration to solve for latitude
28 in terms of isometric latitude in the Reverse method,
30 - Use of Horner's representation for evaluating polynomials and Clenshaw's
31 method for summing trigonometric series,
33 - Several modifications of the formulas to improve the numerical accuracy,
35 - Evaluating the convergence and scale using the expression for the
36 projection or its inverse.
38Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2008-2023)
39and licensed under the MIT/X11 License. For more information, see the
40U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.
41'''
42# make sure int/int division yields float quotient
43from __future__ import division as _; del _ # PYCHOK semicolon
45from pygeodesy.basics import copysign0, isint, isodd, neg, neg_
46from pygeodesy.constants import INF, _K0_UTM, NINF, PI, PI_2, _0_0s, \
47 _0_0, _1_0, _90_0, _180_0
48# from pygeodesy.datums import _spherical_datum # in KTransverseMercator.ellipsoid.setter
49from pygeodesy.errors import _or, _ValueError, _xkwds_get
50from pygeodesy.fmath import fsum1_, hypot, hypot1
51# from pygeodesy.fsums import fsum1_ # from .fmath
52from pygeodesy.interns import NN, _COMMASPACE_, _not_, _singular_
53from pygeodesy.karney import _atan2d, _diff182, _EWGS84, _fix90, \
54 _NamedBase, _norm180, _polynomial, _unsigned2
55from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _pairs
56# from pygeodesy.named import _NamedBase # from .karney
57from pygeodesy.namedTuples import Forward4Tuple, Reverse4Tuple
58from pygeodesy.props import property_doc_, Property, Property_RO, \
59 _update_all
60# from pygeodesy.streprs import pairs as _pairs # from .lazily
61from pygeodesy.units import Degrees, Scalar_, _1mm as _TOL_10 # PYCHOK used!
62from pygeodesy.utily import atand, sincos2, sincos2d_
64from cmath import phase
65from math import atan2, asinh, cos, cosh, degrees, fabs, sin, sinh, sqrt, tanh
67__all__ = _ALL_LAZY.ktm
68__version__ = '23.05.15'
71class KTMError(_ValueError):
72 '''Error raised for L{KTransverseMercator} and L{KTransverseMercator.forward} issues.
73 '''
74 pass
77class KTransverseMercator(_NamedBase):
78 '''I{Karney}'s C++ class U{TransverseMercator<https://GeographicLib.SourceForge.io/
79 C++/doc/classGeographicLib_1_1TransverseMercator.html>} transcoded to pure
80 Python, following is a partial copy of I{Karney}'s documentation.
82 Transverse Mercator projection based on Krüger's method which evaluates the
83 projection and its inverse in terms of a series.
85 There's a singularity in the projection at I{phi = 0, lam - lam0 = +/- (1 - e)
86 90}, about +/- 82.6 degrees for WGS84, where I{e} is the eccentricity. Beyond
87 this point, the series ceases to converge and the results from this method
88 will be garbage. I{To be on the safe side, don't use this method if the
89 angular distance from the central meridian exceeds (1 - 2e) x 90}, about 75
90 degrees for the WGS84 ellipsoid.
92 Class L{ExactTransverseMercator} is an alternative implementation of the
93 projection using I{exact} formulas which yield accurate (to 8 nm) results
94 over the entire ellipsoid.
96 The ellipsoid parameters and the central scale are set in the constructor.
97 The central meridian (which is a trivial shift of the longitude) is specified
98 as the C{lon0} keyword argument of the L{KTransverseMercator.forward} and
99 L{KTransverseMercator.reverse} methods. The latitude of origin is taken to
100 be the equator. There is no provision in this class for specifying a false
101 easting or false northing or a different latitude of origin. However these
102 are can be simply included by the calling function.
104 The L{KTransverseMercator.forward} and L{KTransverseMercator.reverse} methods
105 also return the meridian convergence C{gamma} and scale C{k}. The meridian
106 convergence is the bearing of grid North, the C{y axis}, measured clockwise
107 from true North.
108 '''
109 _E = _EWGS84
110 _k0 = _K0_UTM # central scale factor
111 _lon0 = _0_0 # central meridian
112 _mTM = 6
113 _raiser = False # throw Error
115 def __init__(self, a_earth=_EWGS84, f=None, lon0=0, k0=_K0_UTM, name=NN,
116 raiser=False, **TMorder):
117 '''New L{KTransverseMercator}.
119 @kwarg a_earth: This rhumb's earth (L{Ellipsoid}, L{Ellipsoid2},
120 L{a_f2Tuple}, L{Datum}, 2-tuple (C{a, f})) or the
121 equatorial radius (C{scalar}, C{meter}).
122 @kwarg f: The ellipsoid's flattening (C{scalar}), iff B{C{a_earth}} is
123 a C{scalar}, ignored otherwise.
124 @kwarg lon0: The central meridian (C{degrees180}).
125 @kwarg k0: Central scale factor (C{scalar}).
126 @kwarg name: Optional name (C{str}).
127 @kwarg raiser: If C{True}, throw a L{KTMError} for C{forward}
128 singularities (C{bool}).
129 @kwarg TMorder: Keyword argument B{C{TMorder}}, see property C{TMorder}.
131 @raise KTMError: Invalid B{C{a_earth}}, B{C{f}} or B{C{TMorder}}.
132 '''
133 if f is not None:
134 self.ellipsoid = a_earth, f
135 elif a_earth not in (_EWGS84, None):
136 self.ellipsoid = a_earth
137 self.lon0 = lon0
138 self.k0 = k0
139 if name: # PYCHOK no cover
140 self.name = name
141 if raiser:
142 self.raiser = True
143 if TMorder:
144 self.TMorder = _xkwds_get(TMorder, TMorder=self._mTM)
146 @Property_RO
147 def _Alp(self):
148 return _Xs(_AlpCoeffs, self.TMorder, self.ellipsoid)
150 @Property_RO
151 def _b1(self):
152 n = self.ellipsoid.n
153 if n: # isEllipsoidal
154 m = self.TMorder // 2
155 B1 = _B1Coeffs[m]
156 m += 1
157 b1 = _polynomial(n**2, B1, 0, m) / (B1[m] * (n + _1_0))
158 else: # isSpherical
159 b1 = _1_0 # B1[m - 1] / B1[m1] == 1, always
160 return b1
162 @Property_RO
163 def _Bet(self):
164 C = _Xs(_BetCoeffs, self.TMorder, self.ellipsoid)
165 return tuple(map(neg, C)) if self.f else C # negated if isEllispoidal
167 @Property
168 def ellipsoid(self):
169 '''Get the ellipsoid (L{Ellipsoid}).
170 '''
171 return self._E
173 @ellipsoid.setter # PYCHOK setter!
174 def ellipsoid(self, a_earth_f):
175 '''Set this rhumb's ellipsoid (L{Ellipsoid}, L{Ellipsoid2}, L{Datum},
176 L{a_f2Tuple} or 2-tuple C{(a, f)}).
177 '''
178 E = _MODS.datums._spherical_datum(a_earth_f, Error=KTMError).ellipsoid
179 if self._E != E:
180 _update_all(self)
181 self._E = E
183 @Property_RO
184 def equatoradius(self):
185 '''Get the C{ellipsoid}'s equatorial radius, semi-axis (C{meter}).
186 '''
187 return self.ellipsoid.a
189 a = equatoradius
191 @Property_RO
192 def flattening(self):
193 '''Get the C{ellipsoid}'s flattening (C{float}).
194 '''
195 return self.ellipsoid.f
197 f = flattening
199 def forward(self, lat, lon, lon0=None, name=NN):
200 '''Forward projection, from geographic to transverse Mercator.
202 @arg lat: Latitude of point (C{degrees90}).
203 @arg lon: Longitude of point (C{degrees180}).
204 @arg lon0: Central meridian of the projection (C{degrees180}).
205 @kwarg name: Optional name (C{str}).
207 @return: L{Forward4Tuple}C{(easting, northing, gamma, scale)}
208 with C{easting} and C{northing} in C{meter}, unfalsed, the
209 meridian convergence C{gamma} at point in C{degrees180}
210 and the C{scale} of projection at point C{scalar}. Any
211 value may be C{NAN}, C{NINF} or C{INF} for singularities.
213 @raise KTMError: For singularities, iff property C{raiser} is
214 C{True}.
215 '''
216 lat, _lat = _unsigned2(_fix90(lat))
217 lon, _ = _diff182((self.lon0 if lon0 is None else lon0), lon)
218 lon, _lon = _unsigned2(lon)
219 backside = lon > 90
220 if backside: # PYCHOK no cover
221 lon = _180_0 - lon
222 if lat == 0:
223 _lat = True
225 sphi, cphi, slam, clam = sincos2d_(lat, lon)
226 E = self.ellipsoid
227 if cphi and lat != 90:
228 t = sphi / cphi
229 tp = E.es_taupf(t)
230 h = hypot(tp, clam)
231 if h:
232 xip = atan2(tp, clam)
233 etap = asinh(slam / h) # atanh(sin(lam) / cosh(psi))
234 g = _atan2d(slam * tp, clam * hypot1(tp)) # Krueger p 22 (44)
235 k = sqrt(E.e21 + E.e2 * cphi**2) * hypot1(t) / h
236 elif self.raiser:
237 raise KTMError(lat=lat, lon=lon, lon0=lon0, txt=_singular_)
238 else: # PYCHOK no cover
239 xip, etap = _0_0, (NINF if slam < 0 else INF)
240 g, k = copysign0(_90_0, slam), INF
241 else: # PYCHOK no cover
242 xip, etap = PI_2, _0_0
243 g, k = lon, E.es_c
244 y, x, t, z = self._yxgk4(xip, etap, self._Alp)
245 g -= t
246 k *= z * self._k0_b1
248 if backside: # PYCHOK no cover
249 y, g = (PI - y), (_180_0 - g)
250 y *= self._k0_a1
251 x *= self._k0_a1
252 if _lat:
253 y, g = neg_(y, g)
254 if _lon:
255 x, g = neg_(x, g)
257 return Forward4Tuple(x, y, _norm180(g), k, name=name or self.name)
259 @property_doc_(''' the central scale factor (C{float}).''')
260 def k0(self):
261 '''Get the central scale factor (C{float}), aka I{C{scale0}}.
262 '''
263 return self._k0 # aka scale0
265 @k0.setter # PYCHOK setter!
266 def k0(self, k0):
267 '''Set the central scale factor (C{float}), aka I{C{scale0}}.
269 @raise KTMError: Invalid B{C{k0}}.
270 '''
271 k0 = Scalar_(k0=k0, Error=KTMError, low=_TOL_10, high=_1_0)
272 if self._k0 != k0: # PYCHOK no cover
273 KTransverseMercator._k0_a1._update(self) # redo ._k0_a1
274 KTransverseMercator._k0_b1._update(self) # redo ._k0_b1
275 self._k0 = k0
277 @Property_RO
278 def _k0_a1(self):
279 '''(INTERNAL) Cache C{k0 * _b1 * equatoradius}.
280 '''
281 return self._k0_b1 * self.equatoradius
283 @Property_RO
284 def _k0_b1(self):
285 '''(INTERNAL) Cache C{k0 * _b1}.
286 '''
287 return self.k0 * self._b1
289 @property_doc_(''' the central meridian (C{degrees180}).''')
290 def lon0(self):
291 '''Get the central meridian (C{degrees180}).
292 '''
293 return self._lon0
295 @lon0.setter # PYCHOK setter!
296 def lon0(self, lon0):
297 '''Set the central meridian (C{degrees180}).
299 @raise KTMError: Invalid B{C{lon0}}.
300 '''
301 self._lon0 = _norm180(Degrees(lon0=lon0, Error=KTMError))
303 @property_doc_(''' raise a L{KTMError} for C{forward} singularities (C{bool}).''')
304 def raiser(self):
305 '''Get the error setting (C{bool}).
306 '''
307 return self._raiser
309 @raiser.setter # PYCHOK setter!
310 def raiser(self, raiser):
311 '''Set the error setting (C{bool}), to C{True} to throw a L{KTMError}
312 for C{forward} singularities.
313 '''
314 self._raiser = bool(raiser)
316 def reverse(self, x, y, lon0=None, name=NN):
317 '''Reverse projection, from transverse Mercator to geographic.
319 @arg x: Easting of point (C{meter}).
320 @arg y: Northing of point (C{meter}).
321 @arg lon0: Central meridian of the projection (C{degrees180}).
323 @return: L{Reverse4Tuple}C{(lat, lon, gamma, scale)} with
324 C{lat}- and C{lon}gitude in C{degrees}, I{unfalsed}.
325 '''
326 eta, _lon = _unsigned2(x / self._k0_a1)
327 xi, _lat = _unsigned2(y / self._k0_a1)
328 backside = xi > PI_2
329 if backside: # PYCHOK no cover
330 xi = PI - xi
332 xip, etap, g, k = self._yxgk4(xi, eta, self._Bet)
333 t = self._k0_b1
334 k = (t / k) if k else (NINF if t < 0 else INF)
335 h, c = sinh(etap), cos(xip)
336 if c > 0:
337 r = hypot(h, c)
338 else: # PYCHOK no cover
339 r = fabs(h)
340 c = _0_0
341 E = self.ellipsoid
342 if r:
343 lon = _atan2d(h, c) # Krueger p 17 (25)
344 s = sin(xip) # Newton for tau
345 t = E.es_tauf(s / r)
346 lat = atand(t)
347 g += _atan2d(s * tanh(etap), c) # Krueger p 19 (31)
348 k *= sqrt(E.e21 + E.e2 / (t**2 + _1_0)) * hypot1(t) * r
349 else: # PYCHOK no cover
350 lat, lon = _90_0, _0_0
351 k *= E.es_c
353 if backside: # PYCHOK no cover
354 lon, g = (_180_0 - lon), (_180_0 - g)
355 if _lat:
356 lat, g = neg_(lat, g)
357 if _lon:
358 lon, g = neg_(lon, g)
360 lon += self.lon0 if lon0 is None else _norm180(lon0)
361 return Reverse4Tuple(lat, _norm180(lon), _norm180(g), k,
362 name=name or self.name)
364 @Property
365 def TMorder(self):
366 '''Get the I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8).
367 '''
368 return self._mTM
370 @TMorder.setter # PYCHOK setter!
371 def TMorder(self, order):
372 '''Set the I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8).
373 '''
374 m = _Xorder(_AlpCoeffs, KTMError, TMorder=order)
375 if self._mTM != m:
376 _update_all(self)
377 self._mTM = m
379 def toStr(self, **kwds):
380 '''Return a C{str} representation.
382 @arg kwds: Optional, overriding keyword arguments.
383 '''
384 d = dict(ellipsoid=self.ellipsoid, k0=self.k0, TMorder=self.TMorder)
385 if self.name: # PYCHOK no cover
386 d.update(name=self.name)
387 return _COMMASPACE_.join(_pairs(d, **kwds))
389 def _yxgk4(self, xi_, eta_, C):
390 '''(INTERNAL) Complex Clenshaw summation with
391 C{B{C}=_Alp} or C{B{C}=-_Bet}, negated!
392 '''
393 def _sinhcosh2(x):
394 return sinh(x), cosh(x)
396 x = complex(xi_, eta_)
397 if self.f: # isEllipsoidal
398 s, c = sincos2( xi_ * 2)
399 sh, ch = _sinhcosh2(eta_ * 2)
400 n = -s
401 s = complex(s * ch, c * sh) # sin(zeta * 2)
402 c = complex(c * ch, n * sh) # cos(zeta * 2)
404 y0 = y1 = z0 = z1 = complex(0) # 0+j0
405 n = self.TMorder # == len(C) - 1
406 if isodd(n):
407 Cn = C[n]
408 y0 = complex(Cn) # +j0
409 z0 = complex(Cn * (n * 2))
410 n -= 1
411 a = c * 2 # cos(zeta * 2) * 2
412 _c = _C
413 while n > 0:
414 Cn = C[n]
415 y1 = _c(a, y0, y1, Cn)
416 z1 = _c(a, z0, z1, Cn * (n * 2))
417 n -= 1
418 Cn = C[n]
419 y0 = _c(a, y1, y0, Cn)
420 z0 = _c(a, z1, z0, Cn * (n * 2))
421 n -= 1
422 # assert n == 0
423 x = _c(s, y0, -x, _0_0)
424 c = _c(c, z0, z1, _1_0)
426 # Gauss-Schreiber to Gauss-Krueger TM
427 # C{cmath.phase} handles INF, NAN, etc.
428 g, k = degrees(phase(c)), abs(c)
429 else: # isSpherical
430 g, k = _0_0, _1_0
432 return x.real, x.imag, g, k
435def _C(a, b0, b1, Cn):
436 '''(INTERNAL) Accurately compute complex M{a * b0 - b1 + Cn}
437 with complex args C{a}, C{b0} and C{b1} and scalar C{Cn}.
439 @see: CPython function U{_Py_c_prod<https://GitHub.com/python/
440 cpython/blob/main/Objects/complexobject.c>}.
442 @note: Python function C{cmath.fsum} is no longer available,
443 but stil mentioned in Note 4 of the comments before
444 CPython function U{math_fsum<https://GitHub.com/python/
445 cpython/blob/main/Modules/mathmodule.c>}
446 '''
447 r = fsum1_(a.real * b0.real, -a.imag * b0.imag, -b1.real, Cn, floats=True)
448 j = fsum1_(a.real * b0.imag, a.imag * b0.real, -b1.imag, floats=True)
449 return complex(r, j)
452def _Xorder(_Coeffs, Error, **Xorder): # in .rhumbx
453 '''(INTERNAL) Validate C{RAorder} or C{TMorder}.
454 '''
455 X, m = Xorder.popitem()
456 if m in _Coeffs and isint(m):
457 return m
458 t = sorted(map(str, _Coeffs.keys()))
459 raise Error(X, m, txt=_not_(_or(*t)))
462def _Xs(_Coeffs, m, E, RA=False): # in .rhumbx
463 '''(INTERNAL) Compute the C{A}, C{B} or C{RA} terms of order
464 B{C{m}} for I{Krüger} series and I{rhumbx._sincosSeries},
465 return a tuple with C{B{m} + 1} terms C{X}, C{X[0]==0}.
466 '''
467 Cs = _Coeffs[m]
468 assert len(Cs) == (((m + 1) * (m + 4)) if RA else
469 ((m + 3) * m)) // 2
470 n = n_ = E.n
471 if n: # isEllipsoidal
472 X = [0] # X[0] never used, it's just an integration
473 # constant, it cancels when evaluating a definite
474 # integral. Don't bother computing it, it is not
475 # used in C{KTransverseMercator._yxgk4} above nor
476 # in C{rhumbx._sincosSeries}.
477 i = (m + 2) if RA else 0
478 for r in range(m - 1, -1, -1): # [m-1 ... 0]
479 j = i + r + 1
480 X.append(_polynomial(n, Cs, i, j) * n_ / Cs[j])
481 i = j + 1
482 n_ *= n
483 X = tuple(X)
484 else: # isSpherical
485 X = _0_0s(m + 1)
486 return X
489# _Alp- and _BetCoeffs in .rhumbx
490_AlpCoeffs = { # Generated by Maxima on 2015-05-14 22:55:13-04:00
491 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4
492 164, 225, -480, 360, 720, # Alp[1]/n^1, polynomial(n), order 3
493 557, -864, 390, 1440, # Alp[2]/n^2, polynomial(n), order 2
494 -1236, 427, 1680, # PYCHOK Alp[3]/n^3, polynomial(n), order 1
495 49561, 161280), # Alp[4]/n^4, polynomial(n), order 0, count = 14
496 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5
497 -635, 328, 450, -960, 720, 1440, # Alp[1]/n^1, polynomial(n), order 4
498 4496, 3899, -6048, 2730, 10080, # PYCHOK Alp[2]/n^2, polynomial(n), order 3
499 15061, -19776, 6832, 26880, # PYCHOK Alp[3]/n^3, polynomial(n), order 2
500 -171840, 49561, 161280, # Alp[4]/n^4, polynomial(n), order 1
501 34729, 80640), # PYCHOK Alp[5]/n^5, polynomial(n), order 0, count = 20
502 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6
503 31564, -66675, 34440, 47250, -100800, 75600, 151200, # Alp[1]/n^1, polynomial(n), order 5
504 -1983433, 863232, 748608, -1161216, 524160, 1935360, # PYCHOK Alp[2]/n^2, polynomial(n), order 4
505 670412, 406647, -533952, 184464, 725760, # Alp[3]/n^3, polynomial(n), order 3
506 6601661, -7732800, 2230245, 7257600, # Alp[4]/n^4, polynomial(n), order 2
507 -13675556, 3438171, 7983360, # PYCHOK Alp[5]/n^5, polynomial(n), order 1
508 212378941, 319334400), # Alp[6]/n^6, polynomial(n), order 0, count = 27
509 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7
510 1804025, 2020096, -4267200, 2204160, 3024000, -6451200, 4838400, 9676800, # Alp[1]/n^1, polynomial(n), order 6
511 4626384, -9917165, 4316160, 3743040, -5806080, 2620800, 9676800, # Alp[2]/n^2, polynomial(n), order 5
512 -67102379, 26816480, 16265880, -21358080, 7378560, 29030400, # PYCHOK Alp[3]/n^3, polynomial(n), order 4
513 155912000, 72618271, -85060800, 24532695, 79833600, # Alp[4]/n^4, polynomial(n), order 3
514 102508609, -109404448, 27505368, 63866880, # Alp[5]/n^5, polynomial(n), order 2
515 -12282192400, 2760926233, 4151347200, # PYCHOK Alp[6]/n^6, polynomial(n), order 1
516 1522256789, 1383782400), # Alp[7]/n^7, polynomial(n), order 0, count = 35
517 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8
518 -75900428, 37884525, 42422016, -89611200, 46287360, 63504000, -135475200, 101606400, 203212800, # Alp[1]/n^1, polynomial(n), order 7
519 148003883, 83274912, -178508970, 77690880, 67374720, -104509440, 47174400, 174182400, # PYCHOK Alp[2]/n^2, polynomial(n), order 6
520 318729724, -738126169, 294981280, 178924680, -234938880, 81164160, 319334400, # PYCHOK Alp[3]/n^3, polynomial(n), order 5
521 -40176129013, 14967552000, 6971354016, -8165836800, 2355138720, 7664025600, # Alp[4]/n^4, polynomial(n), order 4
522 10421654396, 3997835751, -4266773472, 1072709352, 2490808320, # PYCHOK Alp[5]/n^5, polynomial(n), order 3
523 175214326799, -171950693600, 38652967262, 58118860800, # PYCHOK Alp[6]/n^6, polynomial(n), order 2
524 -67039739596, 13700311101, 12454041600, # PYCHOK Alp[7]/n^7, polynomial(n), order 1
525 1424729850961, 743921418240) # PYCHOK Alp[8]/n^8, polynomial(n), order 0, count = 44
526}
527_B1Coeffs = { # Generated by Maxima on 2015-05-14 22:55:13-04:00
528 2: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2
529 1, 16, 64, 64), # b1 * (n + 1), polynomial(n2), order 2
530 3: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3
531 1, 4, 64, 256, 256), # b1 * (n + 1), polynomial(n2), order 3
532 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4
533 25, 64, 256, 4096, 16384, 16384) # PYCHOK b1 * (n + 1), polynomial(n2), order 4
534}
535_BetCoeffs = { # Generated by Maxima on 2015-05-14 22:55:13-04:00
536 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4
537 -4, 555, -960, 720, 1440, # Bet[1]/n^1, polynomial(n), order 3
538 -437, 96, 30, 1440, # Bet[2]/n^2, polynomial(n), order 2
539 -148, 119, 3360, # Bet[3]/n^3, polynomial(n), order 1
540 4397, 161280), # PYCHOK Bet[4]/n^4, polynomial(n), order 0, count = 14
541 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5
542 -3645, -64, 8880, -15360, 11520, 23040, # Bet[1]/n^1, polynomial(n), order 4
543 4416, -3059, 672, 210, 10080, # PYCHOK Bet[2]/n^2, polynomial(n), order 3
544 -627, -592, 476, 13440, # Bet[3]/n^3, polynomial(n), order 2
545 -3520, 4397, 161280, # Bet[4]/n^4, polynomial(n), order 1
546 4583, 161280), # PYCHOK Bet[5]/n^5, polynomial(n), order 0, count = 20
547 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6
548 384796, -382725, -6720, 932400, -1612800, 1209600, 2419200, # Bet[1]/n^1, polynomial(n), order 5
549 -1118711, 1695744, -1174656, 258048, 80640, 3870720, # PYCHOK Bet[2]/n^2, polynomial(n), order 4
550 22276, -16929, -15984, 12852, 362880, # Bet[3]/n^3, polynomial(n), order 3
551 -830251, -158400, 197865, 7257600, # PYCHOK Bet[4]/n^4, polynomial(n), order 2
552 -435388, 453717, 15966720, # PYCHOK Bet[5]/n^5, polynomial(n), order 1
553 20648693, 638668800), # Bet[6]/n^6, polynomial(n), order 0, count = 27
554 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7
555 -5406467, 6156736, -6123600, -107520, 14918400, -25804800, 19353600, 38707200, # Bet[1]/n^1, polynomial(n), order 6
556 829456, -5593555, 8478720, -5873280, 1290240, 403200, 19353600, # PYCHOK Bet[2]/n^2, polynomial(n), order 5
557 9261899, 3564160, -2708640, -2557440, 2056320, 58060800, # PYCHOK Bet[3]/n^3, polynomial(n), order 4
558 14928352, -9132761, -1742400, 2176515, 79833600, # PYCHOK Bet[4]/n^4, polynomial(n), order 3
559 -8005831, -1741552, 1814868, 63866880, # Bet[5]/n^5, polynomial(n), order 2
560 -261810608, 268433009, 8302694400, # Bet[6]/n^6, polynomial(n), order 1
561 219941297, 5535129600), # PYCHOK Bet[7]/n^7, polynomial(n), order 0, count = 35
562 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8
563 31777436, -37845269, 43097152, -42865200, -752640, 104428800, -180633600, 135475200, 270950400, # Bet[1]/n^1, polynomial(n), order 7
564 24749483, 14930208, -100683990, 152616960, -105719040, 23224320, 7257600, 348364800, # Bet[2]/n^2, polynomial(n), order 6
565 -232468668, 101880889, 39205760, -29795040, -28131840, 22619520, 638668800, # PYCHOK Bet[3]/n^3, polynomial(n), order 5
566 324154477, 1433121792, -876745056, -167270400, 208945440, 7664025600, # Bet[4]/n^4, polynomial(n), order 4
567 457888660, -312227409, -67920528, 70779852, 2490808320, # Bet[5]/n^5, polynomial(n), order 3
568 -19841813847, -3665348512, 3758062126, 116237721600, # PYCHOK Bet[6]/n^6, polynomial(n), order 2
569 -1989295244, 1979471673, 49816166400, # PYCHOK Bet[7]/n^7, polynomial(n), order 1
570 191773887257, 3719607091200) # Bet[8]/n^8, polynomial(n), order 0, count = 44
571}
573assert set(_AlpCoeffs.keys()) == set(_BetCoeffs.keys())
575if __name__ == '__main__':
577 from pygeodesy.interns import _usage
578 from sys import argv, exit as _exit
580 _exit(_usage(*argv).replace('.ktm', '.etm -series'))
582# **) MIT License
583#
584# Copyright (C) 2022-2023 -- mrJean1 at Gmail -- All Rights Reserved.
585#
586# Permission is hereby granted, free of charge, to any person obtaining a
587# copy of this software and associated documentation files (the "Software"),
588# to deal in the Software without restriction, including without limitation
589# the rights to use, copy, modify, merge, publish, distribute, sublicense,
590# and/or sell copies of the Software, and to permit persons to whom the
591# Software is furnished to do so, subject to the following conditions:
592#
593# The above copyright notice and this permission notice shall be included
594# in all copies or substantial portions of the Software.
595#
596# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
597# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
598# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
599# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
600# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
601# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
602# OTHER DEALINGS IN THE SOFTWARE.