Coverage for pygeodesy/ktm.py: 96%

234 statements  

« prev     ^ index     » next       coverage.py v7.6.0, created at 2024-08-02 18:24 -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.07.16' 

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}, L{a_f2Tuple}, 

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

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

125 is C{meter}, ignored otherwise. 

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

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

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

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

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

131 property C{TMorder}. 

132 

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

134 ''' 

135 if TMorder_name: 

136 M = self._mTM 

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

138 if m != M: 

139 self.TMorder = m 

140 if name: 

141 self.name = name 

142 

143 if f is not None: 

144 self.ellipsoid = a_earth, f 

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

146 pass 

147 elif isinstance(a_earth, Datum): 

148 self.datum = a_earth 

149 else: 

150 self.ellipsoid = a_earth 

151 

152 self.lon0 = lon0 

153 self.k0 = k0 

154 if raiser: 

155 self.raiser = True 

156 

157 @Property_RO 

158 def _Alp(self): 

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

160 

161 @Property_RO 

162 def _b1(self): 

163 n = self.ellipsoid.n 

164 if n: # isEllipsoidal 

165 m = self.TMorder // 2 

166 B1 = _B1Coeffs[m] 

167 m += 1 

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

169 else: # isSpherical 

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

171 return b1 

172 

173 @Property_RO 

174 def _Bet(self): 

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

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

177 

178 @property 

179 def datum(self): 

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

181 ''' 

182 return self._datum 

183 

184 @datum.setter # PYCHOK setter! 

185 def datum(self, datum): 

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

187 ''' 

188 _xinstanceof(Datum, datum=datum) 

189 if self._datum != datum: 

190 _update_all(self) 

191 self._datum = datum 

192 

193 @Property 

194 def ellipsoid(self): 

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

196 ''' 

197 return self.datum.ellipsoid 

198 

199 @ellipsoid.setter # PYCHOK setter! 

200 def ellipsoid(self, a_earth_f): 

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

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

203 ''' 

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

205 

206 @Property_RO 

207 def equatoradius(self): 

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

209 ''' 

210 return self.ellipsoid.a 

211 

212 a = equatoradius 

213 

214 @Property_RO 

215 def flattening(self): 

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

217 ''' 

218 return self.ellipsoid.f 

219 

220 f = flattening 

221 

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

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

224 

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

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

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

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

229 

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

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

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

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

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

235 

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

237 ''' 

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

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

240 lon, _lon = _unsigned2(lon) 

241 backside = lon > 90 

242 if backside: # PYCHOK no cover 

243 lon = _loneg(lon) 

244 if lat == 0: 

245 _lat = True 

246 

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

248 E = self.ellipsoid 

249 if cphi and lat != 90: 

250 t = sphi / cphi 

251 tp = E.es_taupf(t) 

252 h = hypot(tp, clam) 

253 if h: 

254 xip = atan2(tp, clam) 

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

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

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

258 elif self.raiser: 

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

260 else: # PYCHOK no cover 

261 xip, etap = _0_0, _copysignINF(slam) 

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

263 else: # PYCHOK no cover 

264 xip, etap = PI_2, _0_0 

265 g, k = lon, E.es_c 

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

267 g -= d 

268 k *= t * self._k0_b1 

269 

270 if backside: # PYCHOK no cover 

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

272 y *= self._k0_a1 

273 x *= self._k0_a1 

274 if _lat: 

275 y, g = neg_(y, g) 

276 if _lon: 

277 x, g = neg_(x, g) 

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

279 

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

281 def k0(self): 

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

283 ''' 

284 return self._k0 # aka scale0 

285 

286 @k0.setter # PYCHOK setter! 

287 def k0(self, k0): 

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

289 

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

291 ''' 

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

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

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

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

296 self._k0 = k0 

297 

298 @Property_RO 

299 def _k0_a1(self): 

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

301 ''' 

302 return self._k0_b1 * self.equatoradius 

303 

304 @Property_RO 

305 def _k0_b1(self): 

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

307 ''' 

308 return self.k0 * self._b1 

309 

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

311 def lon0(self): 

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

313 ''' 

314 return self._lon0 

315 

316 @lon0.setter # PYCHOK setter! 

317 def lon0(self, lon0): 

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

319 

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

321 ''' 

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

323 

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

325 def raiser(self): 

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

327 ''' 

328 return self._raiser 

329 

330 @raiser.setter # PYCHOK setter! 

331 def raiser(self, raiser): 

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

333 for C{forward} singularities. 

334 ''' 

335 self._raiser = bool(raiser) 

336 

337 def reset(self, lat0, lon0): 

338 '''Set the central parallel and meridian. 

339 

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

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

342 

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

344 parallel and meridian. 

345 

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

347 ''' 

348 t = self._lat0, self.lon0 

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

350 self. lon0 = lon0 

351 return t 

352 

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

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

355 

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

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

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

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

360 

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

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

363 ''' 

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

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

366 backside = xi > PI_2 

367 if backside: # PYCHOK no cover 

368 xi = PI - xi 

369 

370 E = self.ellipsoid 

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

372 t = self._k0_b1 

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

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

375 if c > 0: 

376 r = hypot(h, c) 

377 else: # PYCHOK no cover 

378 r = fabs(h) 

379 c = _0_0 

380 if r: 

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

382 s = sin(xip) # Newton for tau 

383 t = E.es_tauf(s / r) 

384 lat = atan1d(t) 

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

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

387 else: # PYCHOK no cover 

388 lat = _90_0 

389 lon = _0_0 

390 k *= E.es_c 

391 

392 if backside: # PYCHOK no cover 

393 lon, g = _loneg(lon), _loneg(g) 

394 if _lat: 

395 lat, g = neg_(lat, g) 

396 if _lon: 

397 lon, g = neg_(lon, g) 

398 lat += self._lat0 

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

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

401 name=self._name__(name)) 

402 

403 @Property 

404 def TMorder(self): 

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

406 ''' 

407 return self._mTM 

408 

409 @TMorder.setter # PYCHOK setter! 

410 def TMorder(self, order): 

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

412 ''' 

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

414 if self._mTM != m: 

415 _update_all(self) 

416 self._mTM = m 

417 

418 def toStr(self, **kwds): 

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

420 

421 @arg kwds: Optional, overriding keyword arguments. 

422 ''' 

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

424 if self.name: # PYCHOK no cover 

425 d.update(name=self.name) 

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

427 

428 

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

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

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

432 

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

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

435 

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

437 ''' 

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

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

440 return complex(r, j) 

441 

442 

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

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

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

446 ''' 

447 x = complex(xi_, eta_) 

448 if E.f: # isEllipsoidal 

449 s, c = sincos2( xi_ * 2) 

450 sh, ch = _sinhcosh2(eta_ * 2) 

451 n = -s 

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

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

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

455 

456 y0 = y1 = \ 

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

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

459 if isodd(n): 

460 Cn = C[n] 

461 y0 = complex(Cn) # +0j 

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

463 n -= 1 

464 _c = _cma 

465 while n > 0: 

466 Cn = C[n] 

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

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

469 n -= 1 

470 Cn = C[n] 

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

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

473 n -= 1 

474 # assert n == 0 

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

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

477 

478 # Gauss-Schreiber to Gauss-Krueger TM 

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

480 k, g = polar(c) 

481 g = degrees(g) 

482 else: # E.isSpherical 

483 g, k = _0_0, _1_0 

484 

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

486 

487 

488def _sinhcosh2(x): 

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

490 ''' 

491 return sinh(x), cosh(x) 

492 

493 

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

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

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

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

498 ''' 

499 Cs = _Coeffs[m] 

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

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

502 n = n_ = E.n 

503 if n: # isEllipsoidal 

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

505 # constant, it cancels when evaluating a definite 

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

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

508 _X, _p = X.append, _polynomial 

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

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

511 j = i + r + 1 

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

513 i = j + 1 

514 n_ *= n 

515 X = tuple(X) 

516 else: # isSpherical 

517 X = _0_0s(m + 1) 

518 return X 

519 

520 

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

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

523 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 

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

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

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

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

528 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 

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

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

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

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

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

534 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 

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

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

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

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

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

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

541 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 

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

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

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

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

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

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

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

549 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 

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

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

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

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

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

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

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

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

558} 

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

560 2: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2 

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

562 3: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3 

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

564 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4 

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

566} 

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

568 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 

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

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

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

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

573 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 

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

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

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

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

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

579 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 

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

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

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

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

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

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

586 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 

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

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

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

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

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

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

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

594 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 

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

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

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

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

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

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

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

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

603} 

604 

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

606 

607if __name__ == '__main__': 

608 

609 from pygeodesy.internals import _usage 

610 from sys import argv, exit as _exit 

611 

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

613 

614# **) MIT License 

615# 

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

617# 

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

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

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

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

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

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

624# 

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

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

627# 

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

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

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

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

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

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

634# OTHER DEALINGS IN THE SOFTWARE.