Coverage for pygeodesy/ktm.py: 96%

235 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-06-01 11: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_, \ 

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 hypot, hypot1 

53from pygeodesy.fsums import fsum1f_ 

54from pygeodesy.interns import _COMMASPACE_, _singular_ 

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

56 _polynomial, _unsigned2 

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

58from pygeodesy.named import _NamedBase, pairs, _ALL_LAZY 

59from pygeodesy.namedTuples import Forward4Tuple, Reverse4Tuple 

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

61 _update_all 

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

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

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, 

119 raiser=False, **TMorder_name): 

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

126 is a C{scalar}, ignored otherwise. 

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

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

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

130 singularities (C{bool}). 

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

132 keyword argument C{B{TMorder}=6} for the order of 

133 this L{KTransverseMercator}, see property C{TMorder}. 

134 

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

136 ''' 

137 if TMorder_name: 

138 M = self._mTM 

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

140 if m != M: 

141 self.TMorder = m 

142 if name: 

143 self.name = name 

144 

145 if f is not None: 

146 self.ellipsoid = a_earth, f 

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

148 pass 

149 elif isinstance(a_earth, Datum): 

150 self.datum = a_earth 

151 else: 

152 self.ellipsoid = a_earth 

153 

154 self.lon0 = lon0 

155 self.k0 = k0 

156 if raiser: 

157 self.raiser = True 

158 

159 @Property_RO 

160 def _Alp(self): 

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

162 

163 @Property_RO 

164 def _b1(self): 

165 n = self.ellipsoid.n 

166 if n: # isEllipsoidal 

167 m = self.TMorder // 2 

168 B1 = _B1Coeffs[m] 

169 m += 1 

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

171 else: # isSpherical 

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

173 return b1 

174 

175 @Property_RO 

176 def _Bet(self): 

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

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

179 

180 @property 

181 def datum(self): 

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

183 ''' 

184 return self._datum 

185 

186 @datum.setter # PYCHOK setter! 

187 def datum(self, datum): 

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

189 ''' 

190 _xinstanceof(Datum, datum=datum) 

191 if self._datum != datum: 

192 _update_all(self) 

193 self._datum = datum 

194 

195 @Property 

196 def ellipsoid(self): 

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

198 ''' 

199 return self.datum.ellipsoid 

200 

201 @ellipsoid.setter # PYCHOK setter! 

202 def ellipsoid(self, a_earth_f): 

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

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

205 ''' 

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

207 

208 @Property_RO 

209 def equatoradius(self): 

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

211 ''' 

212 return self.ellipsoid.a 

213 

214 a = equatoradius 

215 

216 @Property_RO 

217 def flattening(self): 

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

219 ''' 

220 return self.ellipsoid.f 

221 

222 f = flattening 

223 

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

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

226 

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

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

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

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

231 

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

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

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

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

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

237 

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

239 C{True}. 

240 ''' 

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

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

243 lon, _lon = _unsigned2(lon) 

244 backside = lon > 90 

245 if backside: # PYCHOK no cover 

246 lon = _loneg(lon) 

247 if lat == 0: 

248 _lat = True 

249 

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

251 E = self.ellipsoid 

252 if cphi and lat != 90: 

253 t = sphi / cphi 

254 tp = E.es_taupf(t) 

255 h = hypot(tp, clam) 

256 if h: 

257 xip = atan2(tp, clam) 

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

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

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

261 elif self.raiser: 

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

263 else: # PYCHOK no cover 

264 xip, etap = _0_0, _copysignINF(slam) 

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

266 else: # PYCHOK no cover 

267 xip, etap = PI_2, _0_0 

268 g, k = lon, E.es_c 

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

270 g -= d 

271 k *= t * self._k0_b1 

272 

273 if backside: # PYCHOK no cover 

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

275 y *= self._k0_a1 

276 x *= self._k0_a1 

277 if _lat: 

278 y, g = neg_(y, g) 

279 if _lon: 

280 x, g = neg_(x, g) 

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

282 

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

284 def k0(self): 

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

286 ''' 

287 return self._k0 # aka scale0 

288 

289 @k0.setter # PYCHOK setter! 

290 def k0(self, k0): 

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

292 

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

294 ''' 

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

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

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

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

299 self._k0 = k0 

300 

301 @Property_RO 

302 def _k0_a1(self): 

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

304 ''' 

305 return self._k0_b1 * self.equatoradius 

306 

307 @Property_RO 

308 def _k0_b1(self): 

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

310 ''' 

311 return self.k0 * self._b1 

312 

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

314 def lon0(self): 

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

316 ''' 

317 return self._lon0 

318 

319 @lon0.setter # PYCHOK setter! 

320 def lon0(self, lon0): 

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

322 

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

324 ''' 

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

326 

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

328 def raiser(self): 

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

330 ''' 

331 return self._raiser 

332 

333 @raiser.setter # PYCHOK setter! 

334 def raiser(self, raiser): 

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

336 for C{forward} singularities. 

337 ''' 

338 self._raiser = bool(raiser) 

339 

340 def reset(self, lat0, lon0): 

341 '''Set the central parallel and meridian. 

342 

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

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

345 

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

347 parallel and meridian. 

348 

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

350 ''' 

351 t = self._lat0, self.lon0 

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

353 self. lon0 = lon0 

354 return t 

355 

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

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

358 

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

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

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

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

363 

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

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

366 ''' 

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

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

369 backside = xi > PI_2 

370 if backside: # PYCHOK no cover 

371 xi = PI - xi 

372 

373 E = self.ellipsoid 

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

375 t = self._k0_b1 

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

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

378 if c > 0: 

379 r = hypot(h, c) 

380 else: # PYCHOK no cover 

381 r = fabs(h) 

382 c = _0_0 

383 if r: 

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

385 s = sin(xip) # Newton for tau 

386 t = E.es_tauf(s / r) 

387 lat = atan1d(t) 

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

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

390 else: # PYCHOK no cover 

391 lat = _90_0 

392 lon = _0_0 

393 k *= E.es_c 

394 

395 if backside: # PYCHOK no cover 

396 lon, g = _loneg(lon), _loneg(g) 

397 if _lat: 

398 lat, g = neg_(lat, g) 

399 if _lon: 

400 lon, g = neg_(lon, g) 

401 lat += self._lat0 

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

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

404 name=self._name__(name)) 

405 

406 @Property 

407 def TMorder(self): 

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

409 ''' 

410 return self._mTM 

411 

412 @TMorder.setter # PYCHOK setter! 

413 def TMorder(self, order): 

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

415 ''' 

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

417 if self._mTM != m: 

418 _update_all(self) 

419 self._mTM = m 

420 

421 def toStr(self, **kwds): 

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

423 

424 @arg kwds: Optional, overriding keyword arguments. 

425 ''' 

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

427 if self.name: # PYCHOK no cover 

428 d.update(name=self.name) 

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

430 

431 

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

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

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

435 

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

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

438 

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

440 but stil mentioned in Note 4 of the comments before 

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

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

443 ''' 

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

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

446 return complex(r, j) 

447 

448 

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

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

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

452 ''' 

453 x = complex(xi_, eta_) 

454 if E.f: # isEllipsoidal 

455 s, c = sincos2( xi_ * 2) 

456 sh, ch = _sinhcosh2(eta_ * 2) 

457 n = -s 

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

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

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

461 

462 y0 = y1 = \ 

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

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

465 if isodd(n): 

466 Cn = C[n] 

467 y0 = complex(Cn) # +0j 

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

469 n -= 1 

470 _c = _cma 

471 while n > 0: 

472 Cn = C[n] 

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

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

475 n -= 1 

476 Cn = C[n] 

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

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

479 n -= 1 

480 # assert n == 0 

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

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

483 

484 # Gauss-Schreiber to Gauss-Krueger TM 

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

486 k, g = polar(c) 

487 g = degrees(g) 

488 else: # E.isSpherical 

489 g, k = _0_0, _1_0 

490 

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

492 

493 

494def _sinhcosh2(x): 

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

496 ''' 

497 return sinh(x), cosh(x) 

498 

499 

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

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

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

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

504 ''' 

505 Cs = _Coeffs[m] 

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

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

508 n = n_ = E.n 

509 if n: # isEllipsoidal 

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

511 # constant, it cancels when evaluating a definite 

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

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

514 _X, _p = X.append, _polynomial 

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

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

517 j = i + r + 1 

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

519 i = j + 1 

520 n_ *= n 

521 X = tuple(X) 

522 else: # isSpherical 

523 X = _0_0s(m + 1) 

524 return X 

525 

526 

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

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

529 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 

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

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

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

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

534 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 

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

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

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

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

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

540 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 

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

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

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

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

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

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

547 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 

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

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

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

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

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

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

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

555 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 

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

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

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

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

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

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

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

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

564} 

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

566 2: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2 

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

568 3: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3 

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

570 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4 

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

572} 

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

574 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 

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

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

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

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

579 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 

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

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

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

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

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

585 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 

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

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

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

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

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

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

592 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 

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

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

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

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

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

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

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

600 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 

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

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

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

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

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

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

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

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

609} 

610 

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

612 

613if __name__ == '__main__': 

614 

615 from pygeodesy.internals import _usage 

616 from sys import argv, exit as _exit 

617 

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

619 

620# **) MIT License 

621# 

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

623# 

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

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

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

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

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

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

630# 

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

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

633# 

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

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

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

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

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

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

640# OTHER DEALINGS IN THE SOFTWARE.