Coverage for pygeodesy/ktm.py: 97%

230 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-05-02 14:35 -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_, \ 

46 _reverange, _xinstanceof 

47from pygeodesy.constants import INF, _K0_UTM, PI, PI_2, _0_0s, _0_0, \ 

48 _1_0, _90_0, _copysignINF 

49from pygeodesy.datums import Datum, _spherical_datum, _WGS84, _EWGS84 

50# from pygeodesy.ellipsoids import _EWGS84 # from .datums 

51from pygeodesy.errors import _ValueError, _xkwds_get, _Xorder 

52from pygeodesy.fmath import hypot, hypot1 

53from pygeodesy.fsums import fsum1f_ 

54from pygeodesy.interns import NN, _COMMASPACE_, _singular_ 

55from pygeodesy.karney import _atan2d, _diff182, _fix90, _norm180, \ 

56 _polynomial, _unsigned2, _NamedBase 

57from pygeodesy.lazily import _ALL_LAZY, _pairs 

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

59from pygeodesy.namedTuples import Forward4Tuple, Reverse4Tuple 

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

61 _update_all 

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

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

64from pygeodesy.utily import atan1d, _loneg, sincos2, sincos2d_ 

65 

66from cmath import polar 

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

68 

69__all__ = _ALL_LAZY.ktm 

70__version__ = '24.01.02' 

71 

72 

73class KTMError(_ValueError): 

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

75 ''' 

76 pass 

77 

78 

79class KTransverseMercator(_NamedBase): 

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

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

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

83 

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

85 projection and its inverse in terms of a series. 

86 

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

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

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

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

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

92 degrees for the WGS84 ellipsoid. 

93 

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

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

96 over the entire ellipsoid. 

97 

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

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

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

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

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

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

104 are can be simply included by the calling function. 

105 

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

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

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

109 from true North. 

110 ''' 

111 _datum = _WGS84 

112 _k0 = _K0_UTM # central scale factor 

113 _lat0 = _0_0 # central parallel 

114 _lon0 = _0_0 # central meridian 

115 _mTM = 6 

116 _raiser = False # throw Error 

117 

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

119 raiser=False, **TMorder): 

120 '''New L{KTransverseMercator}. 

121 

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

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

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

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

126 a C{scalar}, ignored otherwise. 

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

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

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

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

131 singularities (C{bool}). 

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

133 

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

135 ''' 

136 if f is not None: 

137 self.ellipsoid = a_earth, f 

138 elif a_earth in (_EWGS84, _WGS84, None): 

139 pass 

140 elif isinstance(a_earth, Datum): 

141 self.datum = a_earth 

142 else: 

143 self.ellipsoid = a_earth 

144 self.lon0 = lon0 

145 self.k0 = k0 

146 if name: # PYCHOK no cover 

147 self.name = name 

148 if raiser: 

149 self.raiser = True 

150 if TMorder: 

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

152 

153 @Property_RO 

154 def _Alp(self): 

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

156 

157 @Property_RO 

158 def _b1(self): 

159 n = self.ellipsoid.n 

160 if n: # isEllipsoidal 

161 m = self.TMorder // 2 

162 B1 = _B1Coeffs[m] 

163 m += 1 

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

165 else: # isSpherical 

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

167 return b1 

168 

169 @Property_RO 

170 def _Bet(self): 

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

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

173 

174 @property 

175 def datum(self): 

176 '''Get this rhumb's datum (L{Datum}). 

177 ''' 

178 return self._datum 

179 

180 @datum.setter # PYCHOK setter! 

181 def datum(self, datum): 

182 '''Set this rhumb's datum (L{Datum}). 

183 ''' 

184 _xinstanceof(Datum, datum=datum) 

185 if self._datum != datum: 

186 _update_all(self) 

187 self._datum = datum 

188 

189 @Property 

190 def ellipsoid(self): 

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

192 ''' 

193 return self.datum.ellipsoid 

194 

195 @ellipsoid.setter # PYCHOK setter! 

196 def ellipsoid(self, a_earth_f): 

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

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

199 ''' 

200 self.datum = _spherical_datum(a_earth_f, Error=KTMError) 

201 

202 @Property_RO 

203 def equatoradius(self): 

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

205 ''' 

206 return self.ellipsoid.a 

207 

208 a = equatoradius 

209 

210 @Property_RO 

211 def flattening(self): 

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

213 ''' 

214 return self.ellipsoid.f 

215 

216 f = flattening 

217 

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

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

220 

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

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

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

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

225 

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

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

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

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

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

231 

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

233 C{True}. 

234 ''' 

235 lat, _lat = _unsigned2(_fix90(lat - self._lat0)) 

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

237 lon, _lon = _unsigned2(lon) 

238 backside = lon > 90 

239 if backside: # PYCHOK no cover 

240 lon = _loneg(lon) 

241 if lat == 0: 

242 _lat = True 

243 

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

245 E = self.ellipsoid 

246 if cphi and lat != 90: 

247 t = sphi / cphi 

248 tp = E.es_taupf(t) 

249 h = hypot(tp, clam) 

250 if h: 

251 xip = atan2(tp, clam) 

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

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

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

255 elif self.raiser: 

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

257 else: # PYCHOK no cover 

258 xip, etap = _0_0, _copysignINF(slam) 

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

260 else: # PYCHOK no cover 

261 xip, etap = PI_2, _0_0 

262 g, k = lon, E.es_c 

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

264 g -= d 

265 k *= t * self._k0_b1 

266 

267 if backside: # PYCHOK no cover 

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

269 y *= self._k0_a1 

270 x *= self._k0_a1 

271 if _lat: 

272 y, g = neg_(y, g) 

273 if _lon: 

274 x, g = neg_(x, g) 

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

276 

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

278 def k0(self): 

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

280 ''' 

281 return self._k0 # aka scale0 

282 

283 @k0.setter # PYCHOK setter! 

284 def k0(self, k0): 

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

286 

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

288 ''' 

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

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

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

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

293 self._k0 = k0 

294 

295 @Property_RO 

296 def _k0_a1(self): 

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

298 ''' 

299 return self._k0_b1 * self.equatoradius 

300 

301 @Property_RO 

302 def _k0_b1(self): 

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

304 ''' 

305 return self.k0 * self._b1 

306 

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

308 def lon0(self): 

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

310 ''' 

311 return self._lon0 

312 

313 @lon0.setter # PYCHOK setter! 

314 def lon0(self, lon0): 

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

316 

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

318 ''' 

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

320 

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

322 def raiser(self): 

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

324 ''' 

325 return self._raiser 

326 

327 @raiser.setter # PYCHOK setter! 

328 def raiser(self, raiser): 

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

330 for C{forward} singularities. 

331 ''' 

332 self._raiser = bool(raiser) 

333 

334 def reset(self, lat0, lon0): 

335 '''Set the central parallel and meridian. 

336 

337 @arg lat0: Latitude of the central parallel (C{degrees90}). 

338 @arg lon0: Longitude of the central parallel (C{degrees180}). 

339 

340 @return: 2-Tuple C{(lat0, lon0)} of the previous central 

341 parallel and meridian. 

342 

343 @raise KTMError: Invalid B{C{lat0}} or B{C{lon0}}. 

344 ''' 

345 t = self._lat0, self.lon0 

346 self._lat0 = _fix90(Degrees(lat0=lat0, Error=KTMError)) 

347 self. lon0 = lon0 

348 return t 

349 

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

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

352 

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

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

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

356 

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

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

359 ''' 

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

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

362 backside = xi > PI_2 

363 if backside: # PYCHOK no cover 

364 xi = PI - xi 

365 

366 E = self.ellipsoid 

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

368 t = self._k0_b1 

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

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

371 if c > 0: 

372 r = hypot(h, c) 

373 else: # PYCHOK no cover 

374 r = fabs(h) 

375 c = _0_0 

376 if r: 

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

378 s = sin(xip) # Newton for tau 

379 t = E.es_tauf(s / r) 

380 lat = atan1d(t) 

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

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

383 else: # PYCHOK no cover 

384 lat = _90_0 

385 lon = _0_0 

386 k *= E.es_c 

387 

388 if backside: # PYCHOK no cover 

389 lon, g = _loneg(lon), _loneg(g) 

390 if _lat: 

391 lat, g = neg_(lat, g) 

392 if _lon: 

393 lon, g = neg_(lon, g) 

394 lat += self._lat0 

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

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

397 name=name or self.name) 

398 

399 @Property 

400 def TMorder(self): 

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

402 ''' 

403 return self._mTM 

404 

405 @TMorder.setter # PYCHOK setter! 

406 def TMorder(self, order): 

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

408 ''' 

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

410 if self._mTM != m: 

411 _update_all(self) 

412 self._mTM = m 

413 

414 def toStr(self, **kwds): 

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

416 

417 @arg kwds: Optional, overriding keyword arguments. 

418 ''' 

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

420 if self.name: # PYCHOK no cover 

421 d.update(name=self.name) 

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

423 

424 

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

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

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

428 

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

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

431 

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

433 but stil mentioned in Note 4 of the comments before 

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

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

436 ''' 

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

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

439 return complex(r, j) 

440 

441 

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

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

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

445 ''' 

446 x = complex(xi_, eta_) 

447 if E.f: # isEllipsoidal 

448 s, c = sincos2( xi_ * 2) 

449 sh, ch = _sinhcosh2(eta_ * 2) 

450 n = -s 

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

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

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

454 

455 y0 = y1 = \ 

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

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

458 if isodd(n): 

459 Cn = C[n] 

460 y0 = complex(Cn) # +0j 

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

462 n -= 1 

463 _c = _cma 

464 while n > 0: 

465 Cn = C[n] 

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

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

468 n -= 1 

469 Cn = C[n] 

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

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

472 n -= 1 

473 # assert n == 0 

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

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

476 

477 # Gauss-Schreiber to Gauss-Krueger TM 

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

479 k, g = polar(c) 

480 g = degrees(g) 

481 else: # E.isSpherical 

482 g, k = _0_0, _1_0 

483 

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

485 

486 

487def _sinhcosh2(x): 

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

489 ''' 

490 return sinh(x), cosh(x) 

491 

492 

493def _Xs(_Coeffs, m, E, RA=False): # in .rhumb.ekx 

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

495 B{C{m}} for I{Krüger} series and I{rhumb.ekx._sincosSeries}, 

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

497 ''' 

498 Cs = _Coeffs[m] 

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

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

501 n = n_ = E.n 

502 if n: # isEllipsoidal 

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

504 # constant, it cancels when evaluating a definite 

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

506 # in C{_Cyxgk4} above and C{rhumb.ekx._sincosSeries}. 

507 _X, _p = X.append, _polynomial 

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

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

510 j = i + r + 1 

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

512 i = j + 1 

513 n_ *= n 

514 X = tuple(X) 

515 else: # isSpherical 

516 X = _0_0s(m + 1) 

517 return X 

518 

519 

520# _Alp- and _BetCoeffs in .rhumb.ekx, .rhumb.bases 

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

522 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 

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

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

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

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

527 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 

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

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

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

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

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

533 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 

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

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

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

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

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

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

540 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 

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

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

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

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

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

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

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

548 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 

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

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

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

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

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

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

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

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

557} 

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

559 2: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2 

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

561 3: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3 

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

563 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4 

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

565} 

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

567 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 

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

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

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

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

572 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 

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

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

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

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

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

578 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 

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

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

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

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

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

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

585 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 

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

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

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

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

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

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

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

593 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 

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

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

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

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

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

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

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

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

602} 

603 

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

605 

606if __name__ == '__main__': 

607 

608 from pygeodesy.interns import _usage 

609 from sys import argv, exit as _exit 

610 

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

612 

613# **) MIT License 

614# 

615# Copyright (C) 2022-2024 -- mrJean1 at Gmail -- All Rights Reserved. 

616# 

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

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

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

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

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

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

623# 

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

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

626# 

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

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

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

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

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

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

633# OTHER DEALINGS IN THE SOFTWARE.