Coverage for pygeodesy/ktm.py: 98%

213 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-09-20 13:43 -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:Karney@Alum.MIT.edu>} (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, PI, PI_2, _0_0s, _0_0, \ 

47 _1_0, _90_0, _copysignINF 

48# from pygeodesy.datums import _spherical_datum # _MODS 

49# from pygeodesy.ellipsoids import _EWGS84 # from .karney 

50from pygeodesy.errors import _ValueError, _xkwds_get, _Xorder 

51from pygeodesy.fmath import hypot, hypot1 

52from pygeodesy.fsums import fsum1f_ 

53from pygeodesy.interns import NN, _COMMASPACE_, _singular_ 

54from pygeodesy.karney import _atan2d, _diff182, _fix90, _polynomial, \ 

55 _norm180, _unsigned2, _EWGS84, _NamedBase 

56from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _pairs 

57# from pygeodesy.named import _NamedBase # from .karney 

58from pygeodesy.namedTuples import Forward4Tuple, Reverse4Tuple 

59from pygeodesy.props import property_doc_, Property, Property_RO, \ 

60 _update_all 

61# from pygeodesy.streprs import pairs as _pairs # from .lazily 

62from pygeodesy.units import Degrees, Scalar_, _1mm as _TOL_10 # PYCHOK used! 

63from pygeodesy.utily import atand, _loneg, sincos2, sincos2d_ 

64 

65from cmath import polar 

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

67 

68__all__ = _ALL_LAZY.ktm 

69__version__ = '23.09.07' 

70 

71 

72class KTMError(_ValueError): 

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

74 ''' 

75 pass 

76 

77 

78class KTransverseMercator(_NamedBase): 

79 '''I{Karney}'s C++ class U{TransverseMercator<https://GeographicLib.SourceForge.io/ 

80 C++/doc/classGeographicLib_1_1TransverseMercator.html>} transcoded to pure 

81 Python, following is a partial copy of I{Karney}'s documentation. 

82 

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

84 projection and its inverse in terms of a series. 

85 

86 There's a singularity in the projection at I{phi = 0, lam - lam0 = +/- (1 - e) 

87 90}, about +/- 82.6 degrees for WGS84, where I{e} is the eccentricity. Beyond 

88 this point, the series ceases to converge and the results from this method 

89 will be garbage. I{To be on the safe side, don't use this method if the 

90 angular distance from the central meridian exceeds (1 - 2e) x 90}, about 75 

91 degrees for the WGS84 ellipsoid. 

92 

93 Class L{ExactTransverseMercator} is an alternative implementation of the 

94 projection using I{exact} formulas which yield accurate (to 8 nm) results 

95 over the entire ellipsoid. 

96 

97 The ellipsoid parameters and the central scale are set in the constructor. 

98 The central meridian (which is a trivial shift of the longitude) is specified 

99 as the C{lon0} keyword argument of the L{KTransverseMercator.forward} and 

100 L{KTransverseMercator.reverse} methods. The latitude of origin is taken to 

101 be the equator. There is no provision in this class for specifying a false 

102 easting or false northing or a different latitude of origin. However these 

103 are can be simply included by the calling function. 

104 

105 The L{KTransverseMercator.forward} and L{KTransverseMercator.reverse} methods 

106 also return the meridian convergence C{gamma} and scale C{k}. The meridian 

107 convergence is the bearing of grid North, the C{y axis}, measured clockwise 

108 from true North. 

109 ''' 

110 _E = _EWGS84 

111 _k0 = _K0_UTM # central scale factor 

112 _lon0 = _0_0 # central meridian 

113 _mTM = 6 

114 _raiser = False # throw Error 

115 

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

117 raiser=False, **TMorder): 

118 '''New L{KTransverseMercator}. 

119 

120 @kwarg a_earth: This rhumb's earth (L{Ellipsoid}, L{Ellipsoid2}, 

121 L{a_f2Tuple}, L{Datum}, 2-tuple (C{a, f})) or the 

122 equatorial radius (C{scalar}, C{meter}). 

123 @kwarg f: The ellipsoid's flattening (C{scalar}), iff B{C{a_earth}} is 

124 a C{scalar}, ignored otherwise. 

125 @kwarg lon0: The central meridian (C{degrees180}). 

126 @kwarg k0: Central scale factor (C{scalar}). 

127 @kwarg name: Optional name (C{str}). 

128 @kwarg raiser: If C{True}, throw a L{KTMError} for C{forward} 

129 singularities (C{bool}). 

130 @kwarg TMorder: Keyword argument B{C{TMorder}}, see property C{TMorder}. 

131 

132 @raise KTMError: Invalid B{C{a_earth}}, B{C{f}} or B{C{TMorder}}. 

133 ''' 

134 if f is not None: 

135 self.ellipsoid = a_earth, f 

136 elif a_earth not in (_EWGS84, None): 

137 self.ellipsoid = a_earth 

138 self.lon0 = lon0 

139 self.k0 = k0 

140 if name: # PYCHOK no cover 

141 self.name = name 

142 if raiser: 

143 self.raiser = True 

144 if TMorder: 

145 self.TMorder = _xkwds_get(TMorder, TMorder=self._mTM) 

146 

147 @Property_RO 

148 def _Alp(self): 

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

150 

151 @Property_RO 

152 def _b1(self): 

153 n = self.ellipsoid.n 

154 if n: # isEllipsoidal 

155 m = self.TMorder // 2 

156 B1 = _B1Coeffs[m] 

157 m += 1 

158 b1 = _polynomial(n**2, B1, 0, m) / (B1[m] * (n + _1_0)) 

159 else: # isSpherical 

160 b1 = _1_0 # B1[m - 1] / B1[m1] == 1, always 

161 return b1 

162 

163 @Property_RO 

164 def _Bet(self): 

165 C = _Xs(_BetCoeffs, self.TMorder, self.ellipsoid) 

166 return tuple(map(neg, C)) if self.f else C # negated if isEllipsoidal 

167 

168 @Property 

169 def ellipsoid(self): 

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

171 ''' 

172 return self._E 

173 

174 @ellipsoid.setter # PYCHOK setter! 

175 def ellipsoid(self, a_earth_f): 

176 '''Set this rhumb's ellipsoid (L{Ellipsoid}, L{Ellipsoid2}, L{Datum}, 

177 L{a_f2Tuple} or 2-tuple C{(a, f)}). 

178 ''' 

179 E = _MODS.datums._spherical_datum(a_earth_f, Error=KTMError).ellipsoid 

180 if self._E != E: 

181 _update_all(self) 

182 self._E = E 

183 

184 @Property_RO 

185 def equatoradius(self): 

186 '''Get the C{ellipsoid}'s equatorial radius, semi-axis (C{meter}). 

187 ''' 

188 return self.ellipsoid.a 

189 

190 a = equatoradius 

191 

192 @Property_RO 

193 def flattening(self): 

194 '''Get the C{ellipsoid}'s flattening (C{scalar}). 

195 ''' 

196 return self.ellipsoid.f 

197 

198 f = flattening 

199 

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

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

202 

203 @arg lat: Latitude of point (C{degrees90}). 

204 @arg lon: Longitude of point (C{degrees180}). 

205 @arg lon0: Central meridian of the projection (C{degrees180}). 

206 @kwarg name: Optional name (C{str}). 

207 

208 @return: L{Forward4Tuple}C{(easting, northing, gamma, scale)} 

209 with C{easting} and C{northing} in C{meter}, unfalsed, the 

210 meridian convergence C{gamma} at point in C{degrees180} 

211 and the C{scale} of projection at point C{scalar}. Any 

212 value may be C{NAN}, C{NINF} or C{INF} for singularities. 

213 

214 @raise KTMError: For singularities, iff property C{raiser} is 

215 C{True}. 

216 ''' 

217 lat, _lat = _unsigned2(_fix90(lat)) 

218 lon, _ = _diff182((self.lon0 if lon0 is None else lon0), lon) 

219 lon, _lon = _unsigned2(lon) 

220 backside = lon > 90 

221 if backside: # PYCHOK no cover 

222 lon = _loneg(lon) 

223 if lat == 0: 

224 _lat = True 

225 

226 sphi, cphi, slam, clam = sincos2d_(lat, lon) 

227 E = self.ellipsoid 

228 if cphi and lat != 90: 

229 t = sphi / cphi 

230 tp = E.es_taupf(t) 

231 h = hypot(tp, clam) 

232 if h: 

233 xip = atan2(tp, clam) 

234 etap = asinh(slam / h) # atanh(sin(lam) / cosh(psi)) 

235 g = _atan2d(slam * tp, clam * hypot1(tp)) # Krueger p 22 (44) 

236 k = sqrt(cphi**2 * E.e2 + E.e21) * hypot1(t) / h 

237 elif self.raiser: 

238 raise KTMError(lat=lat, lon=lon, lon0=lon0, txt=_singular_) 

239 else: # PYCHOK no cover 

240 xip, etap = _0_0, _copysignINF(slam) 

241 g, k = copysign0(_90_0, slam), INF 

242 else: # PYCHOK no cover 

243 xip, etap = PI_2, _0_0 

244 g, k = lon, E.es_c 

245 y, x, d, t = _Cyxgk4(E, xip, etap, self._Alp) 

246 g -= d 

247 k *= t * self._k0_b1 

248 

249 if backside: # PYCHOK no cover 

250 y, g = (PI - y), _loneg(g) 

251 y *= self._k0_a1 

252 x *= self._k0_a1 

253 if _lat: 

254 y, g = neg_(y, g) 

255 if _lon: 

256 x, g = neg_(x, g) 

257 

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

259 

260 @property_doc_(''' the central scale factor (C{float}).''') 

261 def k0(self): 

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

263 ''' 

264 return self._k0 # aka scale0 

265 

266 @k0.setter # PYCHOK setter! 

267 def k0(self, k0): 

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

269 

270 @raise KTMError: Invalid B{C{k0}}. 

271 ''' 

272 k0 = Scalar_(k0=k0, Error=KTMError, low=_TOL_10, high=_1_0) 

273 if self._k0 != k0: # PYCHOK no cover 

274 KTransverseMercator._k0_a1._update(self) # redo ._k0_a1 

275 KTransverseMercator._k0_b1._update(self) # redo ._k0_b1 

276 self._k0 = k0 

277 

278 @Property_RO 

279 def _k0_a1(self): 

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

281 ''' 

282 return self._k0_b1 * self.equatoradius 

283 

284 @Property_RO 

285 def _k0_b1(self): 

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

287 ''' 

288 return self.k0 * self._b1 

289 

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

291 def lon0(self): 

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

293 ''' 

294 return self._lon0 

295 

296 @lon0.setter # PYCHOK setter! 

297 def lon0(self, lon0): 

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

299 

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

301 ''' 

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

303 

304 @property_doc_(''' raise a L{KTMError} for C{forward} singularities (C{bool}).''') 

305 def raiser(self): 

306 '''Get the error setting (C{bool}). 

307 ''' 

308 return self._raiser 

309 

310 @raiser.setter # PYCHOK setter! 

311 def raiser(self, raiser): 

312 '''Set the error setting (C{bool}), to C{True} to throw a L{KTMError} 

313 for C{forward} singularities. 

314 ''' 

315 self._raiser = bool(raiser) 

316 

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

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

319 

320 @arg x: Easting of point (C{meter}). 

321 @arg y: Northing of point (C{meter}). 

322 @arg lon0: Central meridian of the projection (C{degrees180}). 

323 

324 @return: L{Reverse4Tuple}C{(lat, lon, gamma, scale)} with 

325 C{lat}- and C{lon}gitude in C{degrees}, I{unfalsed}. 

326 ''' 

327 eta, _lon = _unsigned2(x / self._k0_a1) 

328 xi, _lat = _unsigned2(y / self._k0_a1) 

329 backside = xi > PI_2 

330 if backside: # PYCHOK no cover 

331 xi = PI - xi 

332 

333 E = self.ellipsoid 

334 xip, etap, g, k = _Cyxgk4(E, xi, eta, self._Bet) 

335 t = self._k0_b1 

336 k = (t / k) if k else _copysignINF(t) # _over(t, k) 

337 h, c = sinh(etap), cos(xip) 

338 if c > 0: 

339 r = hypot(h, c) 

340 else: # PYCHOK no cover 

341 r = fabs(h) 

342 c = _0_0 

343 if r: 

344 lon = _atan2d(h, c) # Krueger p 17 (25) 

345 s = sin(xip) # Newton for tau 

346 t = E.es_tauf(s / r) 

347 lat = atand(t) 

348 g += _atan2d(s * tanh(etap), c) # Krueger p 19 (31) 

349 k *= sqrt(E.e2 / (t**2 + _1_0) + E.e21) * hypot1(t) * r 

350 else: # PYCHOK no cover 

351 lat, lon = _90_0, _0_0 

352 k *= E.es_c 

353 

354 if backside: # PYCHOK no cover 

355 lon, g = _loneg(lon), _loneg(g) 

356 if _lat: 

357 lat, g = neg_(lat, g) 

358 if _lon: 

359 lon, g = neg_(lon, g) 

360 

361 lon += self.lon0 if lon0 is None else _norm180(lon0) 

362 return Reverse4Tuple(lat, _norm180(lon), _norm180(g), k, 

363 name=name or self.name) 

364 

365 @Property 

366 def TMorder(self): 

367 '''Get the I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8). 

368 ''' 

369 return self._mTM 

370 

371 @TMorder.setter # PYCHOK setter! 

372 def TMorder(self, order): 

373 '''Set the I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8). 

374 ''' 

375 m = _Xorder(_AlpCoeffs, KTMError, TMorder=order) 

376 if self._mTM != m: 

377 _update_all(self) 

378 self._mTM = m 

379 

380 def toStr(self, **kwds): 

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

382 

383 @arg kwds: Optional, overriding keyword arguments. 

384 ''' 

385 d = dict(ellipsoid=self.ellipsoid, k0=self.k0, TMorder=self.TMorder) 

386 if self.name: # PYCHOK no cover 

387 d.update(name=self.name) 

388 return _COMMASPACE_.join(_pairs(d, **kwds)) 

389 

390 

391def _cma(a, b0, b1, Cn): 

392 '''(INTERNAL) Compute complex M{a * b0 - b1 + Cn} with complex 

393 C{a}, C{b0} and C{b1} and scalar C{Cn}. 

394 

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

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

397 

398 @note: Python function C{cmath.fsum} is no longer available, 

399 but stil mentioned in Note 4 of the comments before 

400 CPython function U{math_fsum<https://GitHub.com/python/ 

401 cpython/blob/main/Modules/mathmodule.c>} 

402 ''' 

403 r = fsum1f_(a.real * b0.real, -a.imag * b0.imag, -b1.real, Cn) 

404 j = fsum1f_(a.real * b0.imag, a.imag * b0.real, -b1.imag) 

405 return complex(r, j) 

406 

407 

408def _Cyxgk4(E, xi_, eta_, C): 

409 '''(INTERNAL) Complex Clenshaw summation with C{B{C}=._Alp} 

410 or C{B{C}=-._Bet}. 

411 ''' 

412 x = complex(xi_, eta_) 

413 if E.f: # isEllipsoidal 

414 s, c = sincos2( xi_ * 2) 

415 sh, ch = _sinhcosh2(eta_ * 2) 

416 n = -s 

417 s = complex(s * ch, c * sh) # sin(zeta * 2) 

418 c = complex(c * ch, n * sh) # cos(zeta * 2) 

419 a = c * 2 # cos(zeta * 2) * 2 

420 

421 y0 = y1 = \ 

422 z0 = z1 = complex(0) # 0+0j 

423 n = len(C) - 1 # == .TMorder 

424 if isodd(n): 

425 Cn = C[n] 

426 y0 = complex(Cn) # +0j 

427 z0 = complex(Cn * (n * 2)) 

428 n -= 1 

429 _c = _cma 

430 while n > 0: 

431 Cn = C[n] 

432 y1 = _c(a, y0, y1, Cn) 

433 z1 = _c(a, z0, z1, Cn * (n * 2)) 

434 n -= 1 

435 Cn = C[n] 

436 y0 = _c(a, y1, y0, Cn) 

437 z0 = _c(a, z1, z0, Cn * (n * 2)) 

438 n -= 1 

439 # assert n == 0 

440 x = _c(s, y0, -x, _0_0) 

441 c = _c(c, z0, z1, _1_0) 

442 

443 # Gauss-Schreiber to Gauss-Krueger TM 

444 # C{cmath.polar} handles INF, NAN, etc. 

445 k, g = polar(c) 

446 g = degrees(g) 

447 else: # E.isSpherical 

448 g, k = _0_0, _1_0 

449 

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

451 

452 

453def _sinhcosh2(x): 

454 '''(INTERNAL) Like C{sincos2}. 

455 ''' 

456 return sinh(x), cosh(x) 

457 

458 

459def _Xs(_Coeffs, m, E, RA=False): # in .rhumbx 

460 '''(INTERNAL) Compute the C{A}, C{B} or C{RA} terms of order 

461 B{C{m}} for I{Krüger} series and I{rhumbx._sincosSeries}, 

462 return a tuple with C{B{m} + 1} terms C{X}, C{X[0]==0}. 

463 ''' 

464 Cs = _Coeffs[m] 

465 assert len(Cs) == (((m + 1) * (m + 4)) if RA else 

466 ((m + 3) * m)) // 2 

467 n = n_ = E.n 

468 if n: # isEllipsoidal 

469 X = [0] # X[0] never used, it's just an integration 

470 # constant, it cancels when evaluating a definite 

471 # integral. Don't bother computing it, it is unused 

472 # in C{_Cyxgk4} above and C{rhumbx._sincosSeries}. 

473 _X, _p = X.append, _polynomial 

474 i = (m + 2) if RA else 0 

475 for r in _reverange(m): # [m-1 ... 0] 

476 j = i + r + 1 

477 _X(_p(n, Cs, i, j) * n_ / Cs[j]) 

478 i = j + 1 

479 n_ *= n 

480 X = tuple(X) 

481 else: # isSpherical 

482 X = _0_0s(m + 1) 

483 return X 

484 

485 

486# _Alp- and _BetCoeffs in .rhumbx, .rhumbBase 

487_AlpCoeffs = { # Generated by Maxima on 2015-05-14 22:55:13-04:00 

488 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 

489 164, 225, -480, 360, 720, # Alp[1]/n^1, polynomial(n), order 3 

490 557, -864, 390, 1440, # Alp[2]/n^2, polynomial(n), order 2 

491 -1236, 427, 1680, # PYCHOK Alp[3]/n^3, polynomial(n), order 1 

492 49561, 161280), # Alp[4]/n^4, polynomial(n), order 0, count = 14 

493 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 

494 -635, 328, 450, -960, 720, 1440, # Alp[1]/n^1, polynomial(n), order 4 

495 4496, 3899, -6048, 2730, 10080, # PYCHOK Alp[2]/n^2, polynomial(n), order 3 

496 15061, -19776, 6832, 26880, # PYCHOK Alp[3]/n^3, polynomial(n), order 2 

497 -171840, 49561, 161280, # Alp[4]/n^4, polynomial(n), order 1 

498 34729, 80640), # PYCHOK Alp[5]/n^5, polynomial(n), order 0, count = 20 

499 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 

500 31564, -66675, 34440, 47250, -100800, 75600, 151200, # Alp[1]/n^1, polynomial(n), order 5 

501 -1983433, 863232, 748608, -1161216, 524160, 1935360, # PYCHOK Alp[2]/n^2, polynomial(n), order 4 

502 670412, 406647, -533952, 184464, 725760, # Alp[3]/n^3, polynomial(n), order 3 

503 6601661, -7732800, 2230245, 7257600, # Alp[4]/n^4, polynomial(n), order 2 

504 -13675556, 3438171, 7983360, # PYCHOK Alp[5]/n^5, polynomial(n), order 1 

505 212378941, 319334400), # Alp[6]/n^6, polynomial(n), order 0, count = 27 

506 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 

507 1804025, 2020096, -4267200, 2204160, 3024000, -6451200, 4838400, 9676800, # Alp[1]/n^1, polynomial(n), order 6 

508 4626384, -9917165, 4316160, 3743040, -5806080, 2620800, 9676800, # Alp[2]/n^2, polynomial(n), order 5 

509 -67102379, 26816480, 16265880, -21358080, 7378560, 29030400, # PYCHOK Alp[3]/n^3, polynomial(n), order 4 

510 155912000, 72618271, -85060800, 24532695, 79833600, # Alp[4]/n^4, polynomial(n), order 3 

511 102508609, -109404448, 27505368, 63866880, # Alp[5]/n^5, polynomial(n), order 2 

512 -12282192400, 2760926233, 4151347200, # PYCHOK Alp[6]/n^6, polynomial(n), order 1 

513 1522256789, 1383782400), # Alp[7]/n^7, polynomial(n), order 0, count = 35 

514 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 

515 -75900428, 37884525, 42422016, -89611200, 46287360, 63504000, -135475200, 101606400, 203212800, # Alp[1]/n^1, polynomial(n), order 7 

516 148003883, 83274912, -178508970, 77690880, 67374720, -104509440, 47174400, 174182400, # PYCHOK Alp[2]/n^2, polynomial(n), order 6 

517 318729724, -738126169, 294981280, 178924680, -234938880, 81164160, 319334400, # PYCHOK Alp[3]/n^3, polynomial(n), order 5 

518 -40176129013, 14967552000, 6971354016, -8165836800, 2355138720, 7664025600, # Alp[4]/n^4, polynomial(n), order 4 

519 10421654396, 3997835751, -4266773472, 1072709352, 2490808320, # PYCHOK Alp[5]/n^5, polynomial(n), order 3 

520 175214326799, -171950693600, 38652967262, 58118860800, # PYCHOK Alp[6]/n^6, polynomial(n), order 2 

521 -67039739596, 13700311101, 12454041600, # PYCHOK Alp[7]/n^7, polynomial(n), order 1 

522 1424729850961, 743921418240) # PYCHOK Alp[8]/n^8, polynomial(n), order 0, count = 44 

523} 

524_B1Coeffs = { # Generated by Maxima on 2015-05-14 22:55:13-04:00 

525 2: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2 

526 1, 16, 64, 64), # b1 * (n + 1), polynomial(n2), order 2 

527 3: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3 

528 1, 4, 64, 256, 256), # b1 * (n + 1), polynomial(n2), order 3 

529 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4 

530 25, 64, 256, 4096, 16384, 16384) # PYCHOK b1 * (n + 1), polynomial(n2), order 4 

531} 

532_BetCoeffs = { # Generated by Maxima on 2015-05-14 22:55:13-04:00 

533 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 

534 -4, 555, -960, 720, 1440, # Bet[1]/n^1, polynomial(n), order 3 

535 -437, 96, 30, 1440, # Bet[2]/n^2, polynomial(n), order 2 

536 -148, 119, 3360, # Bet[3]/n^3, polynomial(n), order 1 

537 4397, 161280), # PYCHOK Bet[4]/n^4, polynomial(n), order 0, count = 14 

538 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 

539 -3645, -64, 8880, -15360, 11520, 23040, # Bet[1]/n^1, polynomial(n), order 4 

540 4416, -3059, 672, 210, 10080, # PYCHOK Bet[2]/n^2, polynomial(n), order 3 

541 -627, -592, 476, 13440, # Bet[3]/n^3, polynomial(n), order 2 

542 -3520, 4397, 161280, # Bet[4]/n^4, polynomial(n), order 1 

543 4583, 161280), # PYCHOK Bet[5]/n^5, polynomial(n), order 0, count = 20 

544 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 

545 384796, -382725, -6720, 932400, -1612800, 1209600, 2419200, # Bet[1]/n^1, polynomial(n), order 5 

546 -1118711, 1695744, -1174656, 258048, 80640, 3870720, # PYCHOK Bet[2]/n^2, polynomial(n), order 4 

547 22276, -16929, -15984, 12852, 362880, # Bet[3]/n^3, polynomial(n), order 3 

548 -830251, -158400, 197865, 7257600, # PYCHOK Bet[4]/n^4, polynomial(n), order 2 

549 -435388, 453717, 15966720, # PYCHOK Bet[5]/n^5, polynomial(n), order 1 

550 20648693, 638668800), # Bet[6]/n^6, polynomial(n), order 0, count = 27 

551 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 

552 -5406467, 6156736, -6123600, -107520, 14918400, -25804800, 19353600, 38707200, # Bet[1]/n^1, polynomial(n), order 6 

553 829456, -5593555, 8478720, -5873280, 1290240, 403200, 19353600, # PYCHOK Bet[2]/n^2, polynomial(n), order 5 

554 9261899, 3564160, -2708640, -2557440, 2056320, 58060800, # PYCHOK Bet[3]/n^3, polynomial(n), order 4 

555 14928352, -9132761, -1742400, 2176515, 79833600, # PYCHOK Bet[4]/n^4, polynomial(n), order 3 

556 -8005831, -1741552, 1814868, 63866880, # Bet[5]/n^5, polynomial(n), order 2 

557 -261810608, 268433009, 8302694400, # Bet[6]/n^6, polynomial(n), order 1 

558 219941297, 5535129600), # PYCHOK Bet[7]/n^7, polynomial(n), order 0, count = 35 

559 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 

560 31777436, -37845269, 43097152, -42865200, -752640, 104428800, -180633600, 135475200, 270950400, # Bet[1]/n^1, polynomial(n), order 7 

561 24749483, 14930208, -100683990, 152616960, -105719040, 23224320, 7257600, 348364800, # Bet[2]/n^2, polynomial(n), order 6 

562 -232468668, 101880889, 39205760, -29795040, -28131840, 22619520, 638668800, # PYCHOK Bet[3]/n^3, polynomial(n), order 5 

563 324154477, 1433121792, -876745056, -167270400, 208945440, 7664025600, # Bet[4]/n^4, polynomial(n), order 4 

564 457888660, -312227409, -67920528, 70779852, 2490808320, # Bet[5]/n^5, polynomial(n), order 3 

565 -19841813847, -3665348512, 3758062126, 116237721600, # PYCHOK Bet[6]/n^6, polynomial(n), order 2 

566 -1989295244, 1979471673, 49816166400, # PYCHOK Bet[7]/n^7, polynomial(n), order 1 

567 191773887257, 3719607091200) # Bet[8]/n^8, polynomial(n), order 0, count = 44 

568} 

569 

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

571 

572if __name__ == '__main__': 

573 

574 from pygeodesy.interns import _usage 

575 from sys import argv, exit as _exit 

576 

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

578 

579# **) MIT License 

580# 

581# Copyright (C) 2022-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

582# 

583# Permission is hereby granted, free of charge, to any person obtaining a 

584# copy of this software and associated documentation files (the "Software"), 

585# to deal in the Software without restriction, including without limitation 

586# the rights to use, copy, modify, merge, publish, distribute, sublicense, 

587# and/or sell copies of the Software, and to permit persons to whom the 

588# Software is furnished to do so, subject to the following conditions: 

589# 

590# The above copyright notice and this permission notice shall be included 

591# in all copies or substantial portions of the Software. 

592# 

593# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 

594# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 

595# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 

596# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 

597# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 

598# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 

599# OTHER DEALINGS IN THE SOFTWARE.