Coverage for pygeodesy/ktm.py: 98%
210 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-07 07:28 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-07 07:28 -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, isodd, neg, neg_, _reverange
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 _ValueError, _xkwds_get, _Xorder
50from pygeodesy.fmath import fsum1_, hypot, hypot1
51# from pygeodesy.fsums import fsum1_ # from .fmath
52from pygeodesy.interns import NN, _COMMASPACE_, _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.07.01'
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 = \
405 z0 = z1 = complex(0) # 0+j0
406 n = self.TMorder # == len(C) - 1
407 if isodd(n):
408 Cn = C[n]
409 y0 = complex(Cn) # +j0
410 z0 = complex(Cn * (n * 2))
411 n -= 1
412 a = c * 2 # cos(zeta * 2) * 2
413 _c = _C
414 while n > 0:
415 Cn = C[n]
416 y1 = _c(a, y0, y1, Cn)
417 z1 = _c(a, z0, z1, Cn * (n * 2))
418 n -= 1
419 Cn = C[n]
420 y0 = _c(a, y1, y0, Cn)
421 z0 = _c(a, z1, z0, Cn * (n * 2))
422 n -= 1
423 # assert n == 0
424 x = _c(s, y0, -x, _0_0)
425 c = _c(c, z0, z1, _1_0)
427 # Gauss-Schreiber to Gauss-Krueger TM
428 # C{cmath.phase} handles INF, NAN, etc.
429 g, k = degrees(phase(c)), abs(c)
430 else: # isSpherical
431 g, k = _0_0, _1_0
433 return x.real, x.imag, g, k
436def _C(a, b0, b1, Cn):
437 '''(INTERNAL) Accurately compute complex M{a * b0 - b1 + Cn}
438 with complex args C{a}, C{b0} and C{b1} and scalar C{Cn}.
440 @see: CPython function U{_Py_c_prod<https://GitHub.com/python/
441 cpython/blob/main/Objects/complexobject.c>}.
443 @note: Python function C{cmath.fsum} is no longer available,
444 but stil mentioned in Note 4 of the comments before
445 CPython function U{math_fsum<https://GitHub.com/python/
446 cpython/blob/main/Modules/mathmodule.c>}
447 '''
448 r = fsum1_(a.real * b0.real, -a.imag * b0.imag, -b1.real, Cn, floats=True)
449 j = fsum1_(a.real * b0.imag, a.imag * b0.real, -b1.imag, floats=True)
450 return complex(r, j)
453def _Xs(_Coeffs, m, E, RA=False): # in .rhumbx
454 '''(INTERNAL) Compute the C{A}, C{B} or C{RA} terms of order
455 B{C{m}} for I{Krüger} series and I{rhumbx._sincosSeries},
456 return a tuple with C{B{m} + 1} terms C{X}, C{X[0]==0}.
457 '''
458 Cs = _Coeffs[m]
459 assert len(Cs) == (((m + 1) * (m + 4)) if RA else
460 ((m + 3) * m)) // 2
461 n = n_ = E.n
462 if n: # isEllipsoidal
463 X = [0] # X[0] never used, it's just an integration
464 # constant, it cancels when evaluating a definite
465 # integral. Don't bother computing it, it is not
466 # used in C{KTransverseMercator._yxgk4} above nor
467 # in C{rhumbx._sincosSeries}.
468 i = (m + 2) if RA else 0
469 for r in _reverange(m): # [m-1 ... 0]
470 j = i + r + 1
471 X.append(_polynomial(n, Cs, i, j) * n_ / Cs[j])
472 i = j + 1
473 n_ *= n
474 X = tuple(X)
475 else: # isSpherical
476 X = _0_0s(m + 1)
477 return X
480# _Alp- and _BetCoeffs in .rhumbx
481_AlpCoeffs = { # Generated by Maxima on 2015-05-14 22:55:13-04:00
482 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4
483 164, 225, -480, 360, 720, # Alp[1]/n^1, polynomial(n), order 3
484 557, -864, 390, 1440, # Alp[2]/n^2, polynomial(n), order 2
485 -1236, 427, 1680, # PYCHOK Alp[3]/n^3, polynomial(n), order 1
486 49561, 161280), # Alp[4]/n^4, polynomial(n), order 0, count = 14
487 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5
488 -635, 328, 450, -960, 720, 1440, # Alp[1]/n^1, polynomial(n), order 4
489 4496, 3899, -6048, 2730, 10080, # PYCHOK Alp[2]/n^2, polynomial(n), order 3
490 15061, -19776, 6832, 26880, # PYCHOK Alp[3]/n^3, polynomial(n), order 2
491 -171840, 49561, 161280, # Alp[4]/n^4, polynomial(n), order 1
492 34729, 80640), # PYCHOK Alp[5]/n^5, polynomial(n), order 0, count = 20
493 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6
494 31564, -66675, 34440, 47250, -100800, 75600, 151200, # Alp[1]/n^1, polynomial(n), order 5
495 -1983433, 863232, 748608, -1161216, 524160, 1935360, # PYCHOK Alp[2]/n^2, polynomial(n), order 4
496 670412, 406647, -533952, 184464, 725760, # Alp[3]/n^3, polynomial(n), order 3
497 6601661, -7732800, 2230245, 7257600, # Alp[4]/n^4, polynomial(n), order 2
498 -13675556, 3438171, 7983360, # PYCHOK Alp[5]/n^5, polynomial(n), order 1
499 212378941, 319334400), # Alp[6]/n^6, polynomial(n), order 0, count = 27
500 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7
501 1804025, 2020096, -4267200, 2204160, 3024000, -6451200, 4838400, 9676800, # Alp[1]/n^1, polynomial(n), order 6
502 4626384, -9917165, 4316160, 3743040, -5806080, 2620800, 9676800, # Alp[2]/n^2, polynomial(n), order 5
503 -67102379, 26816480, 16265880, -21358080, 7378560, 29030400, # PYCHOK Alp[3]/n^3, polynomial(n), order 4
504 155912000, 72618271, -85060800, 24532695, 79833600, # Alp[4]/n^4, polynomial(n), order 3
505 102508609, -109404448, 27505368, 63866880, # Alp[5]/n^5, polynomial(n), order 2
506 -12282192400, 2760926233, 4151347200, # PYCHOK Alp[6]/n^6, polynomial(n), order 1
507 1522256789, 1383782400), # Alp[7]/n^7, polynomial(n), order 0, count = 35
508 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8
509 -75900428, 37884525, 42422016, -89611200, 46287360, 63504000, -135475200, 101606400, 203212800, # Alp[1]/n^1, polynomial(n), order 7
510 148003883, 83274912, -178508970, 77690880, 67374720, -104509440, 47174400, 174182400, # PYCHOK Alp[2]/n^2, polynomial(n), order 6
511 318729724, -738126169, 294981280, 178924680, -234938880, 81164160, 319334400, # PYCHOK Alp[3]/n^3, polynomial(n), order 5
512 -40176129013, 14967552000, 6971354016, -8165836800, 2355138720, 7664025600, # Alp[4]/n^4, polynomial(n), order 4
513 10421654396, 3997835751, -4266773472, 1072709352, 2490808320, # PYCHOK Alp[5]/n^5, polynomial(n), order 3
514 175214326799, -171950693600, 38652967262, 58118860800, # PYCHOK Alp[6]/n^6, polynomial(n), order 2
515 -67039739596, 13700311101, 12454041600, # PYCHOK Alp[7]/n^7, polynomial(n), order 1
516 1424729850961, 743921418240) # PYCHOK Alp[8]/n^8, polynomial(n), order 0, count = 44
517}
518_B1Coeffs = { # Generated by Maxima on 2015-05-14 22:55:13-04:00
519 2: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2
520 1, 16, 64, 64), # b1 * (n + 1), polynomial(n2), order 2
521 3: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3
522 1, 4, 64, 256, 256), # b1 * (n + 1), polynomial(n2), order 3
523 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4
524 25, 64, 256, 4096, 16384, 16384) # PYCHOK b1 * (n + 1), polynomial(n2), order 4
525}
526_BetCoeffs = { # Generated by Maxima on 2015-05-14 22:55:13-04:00
527 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4
528 -4, 555, -960, 720, 1440, # Bet[1]/n^1, polynomial(n), order 3
529 -437, 96, 30, 1440, # Bet[2]/n^2, polynomial(n), order 2
530 -148, 119, 3360, # Bet[3]/n^3, polynomial(n), order 1
531 4397, 161280), # PYCHOK Bet[4]/n^4, polynomial(n), order 0, count = 14
532 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5
533 -3645, -64, 8880, -15360, 11520, 23040, # Bet[1]/n^1, polynomial(n), order 4
534 4416, -3059, 672, 210, 10080, # PYCHOK Bet[2]/n^2, polynomial(n), order 3
535 -627, -592, 476, 13440, # Bet[3]/n^3, polynomial(n), order 2
536 -3520, 4397, 161280, # Bet[4]/n^4, polynomial(n), order 1
537 4583, 161280), # PYCHOK Bet[5]/n^5, polynomial(n), order 0, count = 20
538 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6
539 384796, -382725, -6720, 932400, -1612800, 1209600, 2419200, # Bet[1]/n^1, polynomial(n), order 5
540 -1118711, 1695744, -1174656, 258048, 80640, 3870720, # PYCHOK Bet[2]/n^2, polynomial(n), order 4
541 22276, -16929, -15984, 12852, 362880, # Bet[3]/n^3, polynomial(n), order 3
542 -830251, -158400, 197865, 7257600, # PYCHOK Bet[4]/n^4, polynomial(n), order 2
543 -435388, 453717, 15966720, # PYCHOK Bet[5]/n^5, polynomial(n), order 1
544 20648693, 638668800), # Bet[6]/n^6, polynomial(n), order 0, count = 27
545 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7
546 -5406467, 6156736, -6123600, -107520, 14918400, -25804800, 19353600, 38707200, # Bet[1]/n^1, polynomial(n), order 6
547 829456, -5593555, 8478720, -5873280, 1290240, 403200, 19353600, # PYCHOK Bet[2]/n^2, polynomial(n), order 5
548 9261899, 3564160, -2708640, -2557440, 2056320, 58060800, # PYCHOK Bet[3]/n^3, polynomial(n), order 4
549 14928352, -9132761, -1742400, 2176515, 79833600, # PYCHOK Bet[4]/n^4, polynomial(n), order 3
550 -8005831, -1741552, 1814868, 63866880, # Bet[5]/n^5, polynomial(n), order 2
551 -261810608, 268433009, 8302694400, # Bet[6]/n^6, polynomial(n), order 1
552 219941297, 5535129600), # PYCHOK Bet[7]/n^7, polynomial(n), order 0, count = 35
553 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8
554 31777436, -37845269, 43097152, -42865200, -752640, 104428800, -180633600, 135475200, 270950400, # Bet[1]/n^1, polynomial(n), order 7
555 24749483, 14930208, -100683990, 152616960, -105719040, 23224320, 7257600, 348364800, # Bet[2]/n^2, polynomial(n), order 6
556 -232468668, 101880889, 39205760, -29795040, -28131840, 22619520, 638668800, # PYCHOK Bet[3]/n^3, polynomial(n), order 5
557 324154477, 1433121792, -876745056, -167270400, 208945440, 7664025600, # Bet[4]/n^4, polynomial(n), order 4
558 457888660, -312227409, -67920528, 70779852, 2490808320, # Bet[5]/n^5, polynomial(n), order 3
559 -19841813847, -3665348512, 3758062126, 116237721600, # PYCHOK Bet[6]/n^6, polynomial(n), order 2
560 -1989295244, 1979471673, 49816166400, # PYCHOK Bet[7]/n^7, polynomial(n), order 1
561 191773887257, 3719607091200) # Bet[8]/n^8, polynomial(n), order 0, count = 44
562}
564assert set(_AlpCoeffs.keys()) == set(_BetCoeffs.keys())
566if __name__ == '__main__':
568 from pygeodesy.interns import _usage
569 from sys import argv, exit as _exit
571 _exit(_usage(*argv).replace('.ktm', '.etm -series'))
573# **) MIT License
574#
575# Copyright (C) 2022-2023 -- mrJean1 at Gmail -- All Rights Reserved.
576#
577# Permission is hereby granted, free of charge, to any person obtaining a
578# copy of this software and associated documentation files (the "Software"),
579# to deal in the Software without restriction, including without limitation
580# the rights to use, copy, modify, merge, publish, distribute, sublicense,
581# and/or sell copies of the Software, and to permit persons to whom the
582# Software is furnished to do so, subject to the following conditions:
583#
584# The above copyright notice and this permission notice shall be included
585# in all copies or substantial portions of the Software.
586#
587# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
588# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
589# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
590# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
591# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
592# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
593# OTHER DEALINGS IN THE SOFTWARE.