Coverage for pygeodesy/ktm.py: 98%

214 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-05 15:46 -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 

14liittyvät karttaprojektiot, tasokoordinaatistot ja karttalehtijako} (Map 

15projections, plane coordinates, and map sheet index for ETRS89), published 

16by JUHTA, Finnish Geodetic Institute, and the National Land Survey of Finland 

17(2006). The relevant section is 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-2022) 

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

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.03.19' 

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 '''Transverse Mercator projection based on Krüger's method which evaluates the 

79 projection and its inverse in terms of a series. 

80 

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

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

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

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

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

86 degrees for the WGS84 ellipsoid. 

87 

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

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

90 over the entire ellipsoid. 

91 

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

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

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

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

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

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

98 are can be simply included by the calling function. 

99 

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

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

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

103 from true North. 

104 ''' 

105 _E = _EWGS84 

106 _k0 = _K0_UTM # central scale factor 

107 _lon0 = _0_0 # central meridian 

108 _mTM = 6 

109 _raiser = False # throw Error 

110 

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

112 raiser=False, **TMorder): 

113 '''New L{KTransverseMercator}. 

114 

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

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

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

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

119 a C{scalar}, ignored otherwise. 

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

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

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

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

124 singularities (C{bool}). 

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

126 

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

128 ''' 

129 if f is not None: 

130 self.ellipsoid = a_earth, f 

131 elif a_earth not in (_EWGS84, None): 

132 self.ellipsoid = a_earth 

133 self.lon0 = lon0 

134 self.k0 = k0 

135 if name: # PYCHOK no cover 

136 self.name = name 

137 if raiser: 

138 self.raiser = True 

139 if TMorder: 

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

141 

142 @Property_RO 

143 def _Alp(self): 

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

145 

146 @Property_RO 

147 def _b1(self): 

148 n = self.ellipsoid.n 

149 if n: # isEllipsoidal 

150 m = self.TMorder // 2 

151 B1 = _B1Coeffs[m] 

152 m += 1 

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

154 else: # isSpherical 

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

156 return b1 

157 

158 @Property_RO 

159 def _Bet(self): 

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

161 return tuple(map(neg, C)) if self.f else C # negated if isEllispoidal 

162 

163 @Property 

164 def ellipsoid(self): 

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

166 ''' 

167 return self._E 

168 

169 @ellipsoid.setter # PYCHOK setter! 

170 def ellipsoid(self, a_earth_f): 

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

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

173 ''' 

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

175 if self._E != E: 

176 _update_all(self) 

177 self._E = E 

178 

179 @Property_RO 

180 def equatoradius(self): 

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

182 ''' 

183 return self.ellipsoid.a 

184 

185 a = equatoradius 

186 

187 @Property_RO 

188 def flattening(self): 

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

190 ''' 

191 return self.ellipsoid.f 

192 

193 f = flattening 

194 

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

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

197 

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

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

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

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

202 

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

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

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

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

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

208 

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

210 C{True}. 

211 ''' 

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

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

214 lon, _lon = _unsigned2(lon) 

215 backside = lon > 90 

216 if backside: # PYCHOK no cover 

217 lon = _180_0 - lon 

218 if lat == 0: 

219 _lat = True 

220 

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

222 E = self.ellipsoid 

223 if cphi and lat != 90: 

224 t = sphi / cphi 

225 tp = E.es_taupf(t) 

226 h = hypot(tp, clam) 

227 if h: 

228 xip = atan2(tp, clam) 

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

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

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

232 elif self.raiser: 

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

234 else: # PYCHOK no cover 

235 xip, etap = _0_0, (NINF if slam < 0 else INF) 

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

237 else: # PYCHOK no cover 

238 xip, etap = PI_2, _0_0 

239 g, k = lon, E.es_c 

240 y, x, t, z = self._yxgk4(xip, etap, self._Alp) 

241 g -= t 

242 k *= z * self._k0_b1 

243 

244 if backside: # PYCHOK no cover 

245 y, g = (PI - y), (_180_0 - g) 

246 y *= self._k0_a1 

247 x *= self._k0_a1 

248 if _lat: 

249 y, g = neg_(y, g) 

250 if _lon: 

251 x, g = neg_(x, g) 

252 

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

254 

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

256 def k0(self): 

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

258 ''' 

259 return self._k0 # aka scale0 

260 

261 @k0.setter # PYCHOK setter! 

262 def k0(self, k0): 

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

264 

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

266 ''' 

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

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

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

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

271 self._k0 = k0 

272 

273 @Property_RO 

274 def _k0_a1(self): 

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

276 ''' 

277 return self._k0_b1 * self.equatoradius 

278 

279 @Property_RO 

280 def _k0_b1(self): 

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

282 ''' 

283 return self.k0 * self._b1 

284 

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

286 def lon0(self): 

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

288 ''' 

289 return self._lon0 

290 

291 @lon0.setter # PYCHOK setter! 

292 def lon0(self, lon0): 

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

294 

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

296 ''' 

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

298 

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

300 def raiser(self): 

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

302 ''' 

303 return self._raiser 

304 

305 @raiser.setter # PYCHOK setter! 

306 def raiser(self, raiser): 

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

308 for C{forward} singularities. 

309 ''' 

310 self._raiser = bool(raiser) 

311 

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

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

314 

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

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

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

318 

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

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

321 ''' 

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

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

324 backside = xi > PI_2 

325 if backside: # PYCHOK no cover 

326 xi = PI - xi 

327 

328 xip, etap, g, k = self._yxgk4(xi, eta, self._Bet) 

329 t = self._k0_b1 

330 k = (t / k) if k else (NINF if t < 0 else INF) 

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

332 if c > 0: 

333 r = hypot(h, c) 

334 else: # PYCHOK no cover 

335 r = fabs(h) 

336 c = _0_0 

337 E = self.ellipsoid 

338 if r: 

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

340 s = sin(xip) # Newton for tau 

341 t = E.es_tauf(s / r) 

342 lat = atand(t) 

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

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

345 else: # PYCHOK no cover 

346 lat, lon = _90_0, _0_0 

347 k *= E.es_c 

348 

349 if backside: # PYCHOK no cover 

350 lon, g = (_180_0 - lon), (_180_0 - g) 

351 if _lat: 

352 lat, g = neg_(lat, g) 

353 if _lon: 

354 lon, g = neg_(lon, g) 

355 

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

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

358 name=name or self.name) 

359 

360 @Property 

361 def TMorder(self): 

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

363 ''' 

364 return self._mTM 

365 

366 @TMorder.setter # PYCHOK setter! 

367 def TMorder(self, order): 

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

369 ''' 

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

371 if self._mTM != m: 

372 _update_all(self) 

373 self._mTM = m 

374 

375 def toStr(self, **kwds): 

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

377 

378 @arg kwds: Optional, overriding keyword arguments. 

379 ''' 

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

381 if self.name: # PYCHOK no cover 

382 d.update(name=self.name) 

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

384 

385 def _yxgk4(self, xi_, eta_, C): 

386 '''(INTERNAL) Complex Clenshaw summation with 

387 C{B{C}=_Alp} or C{B{C}=-_Bet}, negated! 

388 ''' 

389 def _sinhcosh2(x): 

390 return sinh(x), cosh(x) 

391 

392 x = complex(xi_, eta_) 

393 if self.f: # isEllipsoidal 

394 s, c = sincos2( xi_ * 2) 

395 sh, ch = _sinhcosh2(eta_ * 2) 

396 n = -s 

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

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

399 

400 y0 = y1 = z0 = z1 = complex(0) # 0+j0 

401 n = self.TMorder # == len(C) - 1 

402 if isodd(n): 

403 Cn = C[n] 

404 y0 = complex(Cn) # +j0 

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

406 n -= 1 

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

408 while n > 0: 

409 Cn = C[n] 

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

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

412 n -= 1 

413 Cn = C[n] 

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

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

416 n -= 1 

417 # assert n == 0 

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

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

420 

421 # Gauss-Schreiber to Gauss-Krueger TM 

422 # C{cmath.phase} handles INF, NAN, etc. 

423 g, k = degrees(phase(c)), abs(c) 

424 else: # isSpherical 

425 g, k = _0_0, _1_0 

426 

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

428 

429 

430def _c(a, b0, b1, Cn): 

431 '''(INTERNAL) Accurately compute complex M{a * b0 - b1 + Cn} 

432 with complex args C{a}, C{b0} and C{b1} and scalar C{Cn}. 

433 

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

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

436 

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

438 but stil mentioned in Note 4 of the comments before 

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

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

441 ''' 

442 r = fsum1_(a.real * b0.real, -a.imag * b0.imag, -b1.real, Cn, floats=True) 

443 j = fsum1_(a.real * b0.imag, a.imag * b0.real, -b1.imag, floats=True) 

444 return complex(r, j) 

445 

446 

447def _Xorder(_Coeffs, Error, **Xorder): # in .rhumbx 

448 '''(INTERNAL) Validate C{RAorder} or C{TMorder}. 

449 ''' 

450 X, m = Xorder.popitem() 

451 if m in _Coeffs and isint(m): 

452 return m 

453 t = sorted(map(str, _Coeffs.keys())) 

454 raise Error(X, m, txt=_not_(_or(*t))) 

455 

456 

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

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

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

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

461 ''' 

462 Cs = _Coeffs[m] 

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

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

465 n = n_ = E.n 

466 if n: # isEllipsoidal 

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

468 # constant, it cancels when evaluating a definite 

469 # integral. Don't bother computing it, it is not 

470 # used in C{KTransverseMercator._yxgk4} above nor 

471 # in C{rhumbx._sincosSeries}. 

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

473 for r in range(m - 1, -1, -1): # [m-1 ... 0] 

474 j = i + r + 1 

475 X.append(_polynomial(n, Cs, i, j) * n_ / Cs[j]) 

476 i = j + 1 

477 n_ *= n 

478 X = tuple(X) 

479 else: # isSpherical 

480 X = _0_0s(m + 1) 

481 return X 

482 

483 

484# _Alp- and _BetCoeffs in .rhumbx 

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

486 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 

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

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

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

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

491 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 

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

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

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

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

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

497 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 

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

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

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

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

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

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

504 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 

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

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

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

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

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

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

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

512 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 

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

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

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

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

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

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

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

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

521} 

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

523 2: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2 

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

525 3: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3 

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

527 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4 

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

529} 

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

531 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 

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

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

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

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

536 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 

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

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

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

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

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

542 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 

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

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

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

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

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

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

549 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 

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

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

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

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

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

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

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

557 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 

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

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

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

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

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

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

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

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

566} 

567 

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

569 

570if __name__ == '__main__': 

571 

572 from pygeodesy.interns import _usage 

573 from sys import argv, exit as _exit 

574 

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

576 

577# **) MIT License 

578# 

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

580# 

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

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

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

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

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

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

587# 

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

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

590# 

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

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

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

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

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

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

597# OTHER DEALINGS IN THE SOFTWARE.