Coverage for pygeodesy/ktm.py: 96%

232 statements  

« prev     ^ index     » next       coverage.py v7.6.1, created at 2024-11-12 16:17 -0500

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_pop2, _Xorder 

52from pygeodesy.fmath import fdot_, hypot, hypot1 

53from pygeodesy.interns import _COMMASPACE_, _singular_ 

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

55 _polynomial, _unsigned2 

56# from pygeodesy.lazily import _ALL_LAZY # from .named 

57from pygeodesy.named import _NamedBase, pairs, _ALL_LAZY 

58from pygeodesy.namedTuples import Forward4Tuple, Reverse4Tuple 

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

60 _update_all 

61# from pygeodesy.streprs import pairs # from .named 

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

63from pygeodesy.utily import atan1d, _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__ = '24.11.11' 

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 _datum = _WGS84 

111 _k0 = _K0_UTM # central scale factor 

112 _lat0 = _0_0 # central parallel 

113 _lon0 = _0_0 # central meridian 

114 _mTM = 6 

115 _raiser = False # throw Error 

116 

117 def __init__(self, a_earth=_EWGS84, f=None, lon0=0, k0=_K0_UTM, 

118 raiser=False, **TMorder_name): 

119 '''New L{KTransverseMercator}. 

120 

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

122 L{Datum}, 2-tuple (C{a, f})) or the equatorial radius (C{meter}). 

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

124 is C{meter}, ignored otherwise. 

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

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

127 @kwarg raiser: If C{True}, throw a L{KTMError} for C{forward} singularities (C{bool}). 

128 @kwarg TMorder_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument 

129 C{B{TMorder}=6} for the order of this L{KTransverseMercator}, see 

130 property C{TMorder}. 

131 

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

133 ''' 

134 if TMorder_name: 

135 M = self._mTM 

136 m, name = _xkwds_pop2(TMorder_name, TMorder=M) 

137 if m != M: 

138 self.TMorder = m 

139 if name: 

140 self.name = name 

141 

142 if f is not None: 

143 self.ellipsoid = a_earth, f 

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

145 pass 

146 elif isinstance(a_earth, Datum): 

147 self.datum = a_earth 

148 else: 

149 self.ellipsoid = a_earth 

150 

151 self.lon0 = lon0 

152 self.k0 = k0 

153 if raiser: 

154 self.raiser = True 

155 

156 @Property_RO 

157 def _Alp(self): 

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

159 

160 @Property_RO 

161 def _b1(self): 

162 n = self.ellipsoid.n 

163 if n: # isEllipsoidal 

164 m = self.TMorder // 2 

165 B1 = _B1Coeffs[m] 

166 m += 1 

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

168 else: # isSpherical 

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

170 return b1 

171 

172 @Property_RO 

173 def _Bet(self): 

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

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

176 

177 @property 

178 def datum(self): 

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

180 ''' 

181 return self._datum 

182 

183 @datum.setter # PYCHOK setter! 

184 def datum(self, datum): 

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

186 ''' 

187 _xinstanceof(Datum, datum=datum) 

188 if self._datum != datum: 

189 _update_all(self) 

190 self._datum = datum 

191 

192 @Property 

193 def ellipsoid(self): 

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

195 ''' 

196 return self.datum.ellipsoid 

197 

198 @ellipsoid.setter # PYCHOK setter! 

199 def ellipsoid(self, a_earth_f): 

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

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

202 ''' 

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

204 

205 @Property_RO 

206 def equatoradius(self): 

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

208 ''' 

209 return self.ellipsoid.a 

210 

211 a = equatoradius 

212 

213 @Property_RO 

214 def flattening(self): 

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

216 ''' 

217 return self.ellipsoid.f 

218 

219 f = flattening 

220 

221 def forward(self, lat, lon, lon0=None, **name): 

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

223 

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

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

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

227 @kwarg name: Optional C{B{name}=NN} (C{str}). 

228 

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

230 C{easting} and C{northing} in C{meter}, unfalsed, the 

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

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

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

234 

235 @raise KTMError: For singularities, iff property C{raiser} is C{True}. 

236 ''' 

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

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

239 lon, _lon = _unsigned2(lon) 

240 backside = lon > 90 

241 if backside: # PYCHOK no cover 

242 lon = _loneg(lon) 

243 if lat == 0: 

244 _lat = True 

245 

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

247 E = self.ellipsoid 

248 if cphi and lat != 90: 

249 t = sphi / cphi 

250 tp = E.es_taupf(t) 

251 h = hypot(tp, clam) 

252 if h: 

253 xip = atan2(tp, clam) 

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

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

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

257 elif self.raiser: 

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

259 else: # PYCHOK no cover 

260 xip, etap = _0_0, _copysignINF(slam) 

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

262 else: # PYCHOK no cover 

263 xip, etap = PI_2, _0_0 

264 g, k = lon, E.es_c 

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

266 g -= d 

267 k *= t * self._k0_b1 

268 

269 if backside: # PYCHOK no cover 

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

271 y *= self._k0_a1 

272 x *= self._k0_a1 

273 if _lat: 

274 y, g = neg_(y, g) 

275 if _lon: 

276 x, g = neg_(x, g) 

277 return Forward4Tuple(x, y, _norm180(g), k, name=self._name__(name)) 

278 

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

280 def k0(self): 

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

282 ''' 

283 return self._k0 # aka scale0 

284 

285 @k0.setter # PYCHOK setter! 

286 def k0(self, k0): 

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

288 

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

290 ''' 

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

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

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

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

295 self._k0 = k0 

296 

297 @Property_RO 

298 def _k0_a1(self): 

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

300 ''' 

301 return self._k0_b1 * self.equatoradius 

302 

303 @Property_RO 

304 def _k0_b1(self): 

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

306 ''' 

307 return self.k0 * self._b1 

308 

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

310 def lon0(self): 

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

312 ''' 

313 return self._lon0 

314 

315 @lon0.setter # PYCHOK setter! 

316 def lon0(self, lon0): 

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

318 

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

320 ''' 

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

322 

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

324 def raiser(self): 

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

326 ''' 

327 return self._raiser 

328 

329 @raiser.setter # PYCHOK setter! 

330 def raiser(self, raiser): 

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

332 for C{forward} singularities. 

333 ''' 

334 self._raiser = bool(raiser) 

335 

336 def reset(self, lat0, lon0): 

337 '''Set the central parallel and meridian. 

338 

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

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

341 

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

343 parallel and meridian. 

344 

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

346 ''' 

347 t = self._lat0, self.lon0 

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

349 self. lon0 = lon0 

350 return t 

351 

352 def reverse(self, x, y, lon0=None, **name): 

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

354 

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

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

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

358 @kwarg name: Optional C{B{name}=NN} (C{str}). 

359 

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

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

362 ''' 

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

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

365 backside = xi > PI_2 

366 if backside: # PYCHOK no cover 

367 xi = PI - xi 

368 

369 E = self.ellipsoid 

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

371 t = self._k0_b1 

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

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

374 if c > 0: 

375 r = hypot(h, c) 

376 else: # PYCHOK no cover 

377 r = fabs(h) 

378 c = _0_0 

379 if r: 

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

381 s = sin(xip) # Newton for tau 

382 t = E.es_tauf(s / r) 

383 lat = atan1d(t) 

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

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

386 else: # PYCHOK no cover 

387 lat = _90_0 

388 lon = _0_0 

389 k *= E.es_c 

390 

391 if backside: # PYCHOK no cover 

392 lon, g = _loneg(lon), _loneg(g) 

393 if _lat: 

394 lat, g = neg_(lat, g) 

395 if _lon: 

396 lon, g = neg_(lon, g) 

397 lat += self._lat0 

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

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

400 name=self._name__(name)) 

401 

402 @Property 

403 def TMorder(self): 

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

405 ''' 

406 return self._mTM 

407 

408 @TMorder.setter # PYCHOK setter! 

409 def TMorder(self, order): 

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

411 ''' 

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

413 if self._mTM != m: 

414 _update_all(self) 

415 self._mTM = m 

416 

417 def toStr(self, **kwds): 

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

419 

420 @arg kwds: Optional, overriding keyword arguments. 

421 ''' 

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

423 if self.name: # PYCHOK no cover 

424 d.update(name=self.name) 

425 return _COMMASPACE_.join(pairs(d, **kwds)) 

426 

427 

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

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

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

431 

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

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

434 

435 @note: Python function C{cmath.fsum} no longer exists. 

436 ''' 

437 r = fdot_(a.real, b0.real, -a.imag, b0.imag, -b1.real, _1_0, start=Cn) 

438 j = fdot_(a.real, b0.imag, a.imag, b0.real, -b1.imag, _1_0) 

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 while n > 0: 

464 Cn = C[n] 

465 y1 = _cma(a, y0, y1, Cn) 

466 z1 = _cma(a, z0, z1, Cn * (n * 2)) 

467 n -= 1 

468 Cn = C[n] 

469 y0 = _cma(a, y1, y0, Cn) 

470 z0 = _cma(a, z1, z0, Cn * (n * 2)) 

471 n -= 1 

472 # assert n == 0 

473 x = _cma(s, y0, -x, _0_0) 

474 c = _cma(c, z0, z1, _1_0) 

475 

476 # Gauss-Schreiber to Gauss-Krueger TM 

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

478 k, g = polar(c) 

479 g = degrees(g) 

480 else: # E.isSpherical 

481 g, k = _0_0, _1_0 

482 

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

484 

485 

486def _sinhcosh2(x): 

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

488 ''' 

489 return sinh(x), cosh(x) 

490 

491 

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

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

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

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

496 ''' 

497 Cs = _Coeffs[m] 

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

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

500 n = n_ = E.n 

501 if n: # isEllipsoidal 

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

503 # constant, it cancels when evaluating a definite 

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

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

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

507 _p = _polynomial 

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

509 j = i + r + 1 

510 X.append(_p(n, Cs, i, j) * n_ / Cs[j]) 

511 i = j + 1 

512 n_ *= n 

513 X = tuple(X) 

514 else: # isSpherical 

515 X = _0_0s(m + 1) 

516 return X 

517 

518 

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

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

521 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 

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

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

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

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

526 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 

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

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

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

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

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

532 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 

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

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

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

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

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

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

539 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 

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

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

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

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

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

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

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

547 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 

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

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

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

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

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

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

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

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

556} 

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

558 2: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2 

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

560 3: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3 

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

562 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4 

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

564} 

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

566 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 

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

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

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

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

571 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 

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

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

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

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

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

577 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 

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

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

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

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

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

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

584 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 

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

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

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

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

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

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

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

592 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 

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

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

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

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

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

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

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

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

601} 

602 

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

604 

605if __name__ == '__main__': 

606 

607 from pygeodesy.internals import _usage 

608 from sys import argv, exit as _exit 

609 

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

611 

612# **) MIT License 

613# 

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

615# 

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

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

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

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

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

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

622# 

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

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

625# 

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

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

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

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

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

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

632# OTHER DEALINGS IN THE SOFTWARE.