Coverage for pygeodesy/ktm.py: 98%

210 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-07 07:28 -0400

1 

2# -*- coding: utf-8 -*- 

3 

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>}. 

8 

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}. 

12 

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>}. 

19 

20This is a straight transcription of the formulas in this paper with the 

21following exceptions: 

22 

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%, 

26 

27 - Use Newton's method instead of plain iteration to solve for latitude 

28 in terms of isometric latitude in the Reverse method, 

29 

30 - Use of Horner's representation for evaluating polynomials and Clenshaw's 

31 method for summing trigonometric series, 

32 

33 - Several modifications of the formulas to improve the numerical accuracy, 

34 

35 - Evaluating the convergence and scale using the expression for the 

36 projection or its inverse. 

37 

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 

44 

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_ 

63 

64from cmath import phase 

65from math import atan2, asinh, cos, cosh, degrees, fabs, sin, sinh, sqrt, tanh 

66 

67__all__ = _ALL_LAZY.ktm 

68__version__ = '23.07.01' 

69 

70 

71class KTMError(_ValueError): 

72 '''Error raised for L{KTransverseMercator} and L{KTransverseMercator.forward} issues. 

73 ''' 

74 pass 

75 

76 

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. 

81 

82 Transverse Mercator projection based on Krüger's method which evaluates the 

83 projection and its inverse in terms of a series. 

84 

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. 

91 

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. 

95 

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. 

103 

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 

114 

115 def __init__(self, a_earth=_EWGS84, f=None, lon0=0, k0=_K0_UTM, name=NN, 

116 raiser=False, **TMorder): 

117 '''New L{KTransverseMercator}. 

118 

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}. 

130 

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) 

145 

146 @Property_RO 

147 def _Alp(self): 

148 return _Xs(_AlpCoeffs, self.TMorder, self.ellipsoid) 

149 

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 

161 

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 

166 

167 @Property 

168 def ellipsoid(self): 

169 '''Get the ellipsoid (L{Ellipsoid}). 

170 ''' 

171 return self._E 

172 

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 

182 

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 

188 

189 a = equatoradius 

190 

191 @Property_RO 

192 def flattening(self): 

193 '''Get the C{ellipsoid}'s flattening (C{float}). 

194 ''' 

195 return self.ellipsoid.f 

196 

197 f = flattening 

198 

199 def forward(self, lat, lon, lon0=None, name=NN): 

200 '''Forward projection, from geographic to transverse Mercator. 

201 

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}). 

206 

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. 

212 

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 

224 

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 

247 

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) 

256 

257 return Forward4Tuple(x, y, _norm180(g), k, name=name or self.name) 

258 

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 

264 

265 @k0.setter # PYCHOK setter! 

266 def k0(self, k0): 

267 '''Set the central scale factor (C{float}), aka I{C{scale0}}. 

268 

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 

276 

277 @Property_RO 

278 def _k0_a1(self): 

279 '''(INTERNAL) Cache C{k0 * _b1 * equatoradius}. 

280 ''' 

281 return self._k0_b1 * self.equatoradius 

282 

283 @Property_RO 

284 def _k0_b1(self): 

285 '''(INTERNAL) Cache C{k0 * _b1}. 

286 ''' 

287 return self.k0 * self._b1 

288 

289 @property_doc_(''' the central meridian (C{degrees180}).''') 

290 def lon0(self): 

291 '''Get the central meridian (C{degrees180}). 

292 ''' 

293 return self._lon0 

294 

295 @lon0.setter # PYCHOK setter! 

296 def lon0(self, lon0): 

297 '''Set the central meridian (C{degrees180}). 

298 

299 @raise KTMError: Invalid B{C{lon0}}. 

300 ''' 

301 self._lon0 = _norm180(Degrees(lon0=lon0, Error=KTMError)) 

302 

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 

308 

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) 

315 

316 def reverse(self, x, y, lon0=None, name=NN): 

317 '''Reverse projection, from transverse Mercator to geographic. 

318 

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}). 

322 

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 

331 

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 

352 

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) 

359 

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) 

363 

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 

369 

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 

378 

379 def toStr(self, **kwds): 

380 '''Return a C{str} representation. 

381 

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

388 

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) 

395 

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) 

403 

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) 

426 

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 

432 

433 return x.real, x.imag, g, k 

434 

435 

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}. 

439 

440 @see: CPython function U{_Py_c_prod<https://GitHub.com/python/ 

441 cpython/blob/main/Objects/complexobject.c>}. 

442 

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) 

451 

452 

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 

478 

479 

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} 

563 

564assert set(_AlpCoeffs.keys()) == set(_BetCoeffs.keys()) 

565 

566if __name__ == '__main__': 

567 

568 from pygeodesy.interns import _usage 

569 from sys import argv, exit as _exit 

570 

571 _exit(_usage(*argv).replace('.ktm', '.etm -series')) 

572 

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.