1

2# -*- coding: utf-8 -*-

3

4 u'''A pure Python version of I{Karney}'s C++ class U{TransverseMercator

5<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1TransverseMercator.html>}

6 based on I{Krüger} series. See also I{Karney}'s utility U{TransverseMercatorProj

7<https://GeographicLib.SourceForge.io/C++/doc/TransverseMercatorProj.1.html>}.

8

9 Following and further below is a copy of I{Karney}'s U{TransverseMercator.hpp

10<https://GeographicLib.SourceForge.io/C++/doc/TransverseMercator_8hpp_source.html>}

11 file C{Header}.

12

13 This implementation follows closely JHS 154, ETRS89 - I{järjestelmään liittyvät

14 karttaprojektiot, tasokoordinaatistot ja karttalehtijako} (Map projections, plane

15 coordinates, and map sheet index for ETRS89), published by JUHTA, Finnish Geodetic

16 Institute, and the National Land Survey of Finland (2006). The relevant section

17 is available as the U{2008 PDF file

18<http://Docs.JHS-suositukset.FI/jhs-suositukset/JHS154/JHS154_liite1.pdf>}.

19

20 This is a straight transcription of the formulas in this paper with the

21 following 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

38 Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2008-2023)

39 and licensed under the MIT/X11 License. For more information, see the

40 U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.

41'''

42# make sure int/int division yields float quotient

43 from __future__ import division as _; del _ # PYCHOK semicolon

44

45 from pygeodesy.basics import copysign0, isint, isodd, neg, neg_

46 from pygeodesy.constants import INF, _K0_UTM, NINF, PI, PI_2, _0_0s, \

47 _0_0, _1_0, _90_0, _180_0

48# from pygeodesy.datums import _spherical_datum # in KTransverseMercator.ellipsoid.setter

49 from pygeodesy.errors import _or, _ValueError, _xkwds_get

50 from pygeodesy.fmath import fsum1_, hypot, hypot1

51# from pygeodesy.fsums import fsum1_ # from .fmath

52 from pygeodesy.interns import NN, _COMMASPACE_, _not_, _singular_

53 from pygeodesy.karney import _atan2d, _diff182, _EWGS84, _fix90, \

54 _NamedBase, _norm180, _polynomial, _unsigned2

55 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _pairs

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

57 from pygeodesy.namedTuples import Forward4Tuple, Reverse4Tuple

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

59 _update_all

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

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

62 from pygeodesy.utily import atand, sincos2, sincos2d_

63

64 from cmath import phase

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

66

67__all__ = _ALL_LAZY.ktm

68__version__ = '23.05.15'

69

70

71 class KTMError(_ValueError):

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

73 '''

74 pass

75

76

77 class KTransverseMercator(_NamedBase):

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

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

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

81

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

83 projection and its inverse in terms of a series.

84

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

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

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

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

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

90 degrees for the WGS84 ellipsoid.

91

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

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

94 over the entire ellipsoid.

95

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

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

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

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

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

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

102 are can be simply included by the calling function.

103

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

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

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

107 from true North.

108 '''

109 _E = _EWGS84

110 _k0 = _K0_UTM # central scale factor

111 _lon0 = _0_0 # central meridian

112 _mTM = 6

113 _raiser = False # throw Error

114

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

116 raiser=False, **TMorder):

117 '''New L{KTransverseMercator}.

118

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

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

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

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

123 a C{scalar}, ignored otherwise.

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

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

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

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

128 singularities (C{bool}).

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

130

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

132 '''

133 if f is not None:

134 self.ellipsoid = a_earth, f

135 elif a_earth not in (_EWGS84, None):

136 self.ellipsoid = a_earth

137 self.lon0 = lon0

138 self.k0 = k0

139 if name: # PYCHOK no cover

140 self.name = name

141 if raiser:

142 self.raiser = True

143 if TMorder:

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

145

146 @Property_RO

147 def _Alp(self):

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

149

150 @Property_RO

151 def _b1(self):

152 n = self.ellipsoid.n

153 if n: # isEllipsoidal

154 m = self.TMorder // 2

155 B1 = _B1Coeffs[m]

156 m += 1

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

158 else: # isSpherical

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

160 return b1

161

162 @Property_RO

163 def _Bet(self):

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

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

166

167 @Property

168 def ellipsoid(self):

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

170 '''

171 return self._E

172

173 @ellipsoid.setter # PYCHOK setter!

174 def ellipsoid(self, a_earth_f):

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

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

177 '''

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

179 if self._E != E:

180 _update_all(self)

181 self._E = E

182

183 @Property_RO

184 def equatoradius(self):

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

186 '''

187 return self.ellipsoid.a

188

189 a = equatoradius

190

191 @Property_RO

192 def flattening(self):

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

194 '''

195 return self.ellipsoid.f

196

197 f = flattening

198

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

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

201

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

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

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

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

206

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

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

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

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

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

212

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

214 C{True}.

215 '''

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

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

218 lon, _lon = _unsigned2(lon)

219 backside = lon > 90

220 if backside: # PYCHOK no cover

221 lon = _180_0 - lon

222 if lat == 0:

223 _lat = True

224

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

226 E = self.ellipsoid

227 if cphi and lat != 90:

228 t = sphi / cphi

229 tp = E.es_taupf(t)

230 h = hypot(tp, clam)

231 if h:

232 xip = atan2(tp, clam)

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

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

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

236 elif self.raiser:

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

238 else: # PYCHOK no cover

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

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

241 else: # PYCHOK no cover

242 xip, etap = PI_2, _0_0

243 g, k = lon, E.es_c

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

245 g -= t

246 k *= z * self._k0_b1

247

248 if backside: # PYCHOK no cover

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

250 y *= self._k0_a1

251 x *= self._k0_a1

252 if _lat:

253 y, g = neg_(y, g)

254 if _lon:

255 x, g = neg_(x, g)

256

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

258

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

260 def k0(self):

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

262 '''

263 return self._k0 # aka scale0

264

265 @k0.setter # PYCHOK setter!

266 def k0(self, k0):

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

268

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

270 '''

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

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

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

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

275 self._k0 = k0

276

277 @Property_RO

278 def _k0_a1(self):

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

280 '''

281 return self._k0_b1 * self.equatoradius

282

283 @Property_RO

284 def _k0_b1(self):

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

286 '''

287 return self.k0 * self._b1

288

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

290 def lon0(self):

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

292 '''

293 return self._lon0

294

295 @lon0.setter # PYCHOK setter!

296 def lon0(self, lon0):

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

298

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

300 '''

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

302

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

304 def raiser(self):

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

306 '''

307 return self._raiser

308

309 @raiser.setter # PYCHOK setter!

310 def raiser(self, raiser):

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

312 for C{forward} singularities.

313 '''

314 self._raiser = bool(raiser)

315

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

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

318

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

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

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

322

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

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

325 '''

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

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

328 backside = xi > PI_2

329 if backside: # PYCHOK no cover

330 xi = PI - xi

331

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

333 t = self._k0_b1

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

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

336 if c > 0:

337 r = hypot(h, c)

338 else: # PYCHOK no cover

339 r = fabs(h)

340 c = _0_0

341 E = self.ellipsoid

342 if r:

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

344 s = sin(xip) # Newton for tau

345 t = E.es_tauf(s / r)

346 lat = atand(t)

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

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

349 else: # PYCHOK no cover

350 lat, lon = _90_0, _0_0

351 k *= E.es_c

352

353 if backside: # PYCHOK no cover

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

355 if _lat:

356 lat, g = neg_(lat, g)

357 if _lon:

358 lon, g = neg_(lon, g)

359

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

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

362 name=name or self.name)

363

364 @Property

365 def TMorder(self):

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

367 '''

368 return self._mTM

369

370 @TMorder.setter # PYCHOK setter!

371 def TMorder(self, order):

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

373 '''

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

375 if self._mTM != m:

376 _update_all(self)

377 self._mTM = m

378

379 def toStr(self, **kwds):

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

381

382 @arg kwds: Optional, overriding keyword arguments.

383 '''

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

385 if self.name: # PYCHOK no cover

386 d.update(name=self.name)

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

388

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

390 '''(INTERNAL) Complex Clenshaw summation with

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

392 '''

393 def _sinhcosh2(x):

394 return sinh(x), cosh(x)

395

396 x = complex(xi_, eta_)

397 if self.f: # isEllipsoidal

398 s, c = sincos2( xi_ * 2)

399 sh, ch = _sinhcosh2(eta_ * 2)

400 n = -s

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

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

403

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

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

406 if isodd(n):

407 Cn = C[n]

408 y0 = complex(Cn) # +j0

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

410 n -= 1

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

412 _c = _C

413 while n > 0:

414 Cn = C[n]

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

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

417 n -= 1

418 Cn = C[n]

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

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

421 n -= 1

422 # assert n == 0

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

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

425

426 # Gauss-Schreiber to Gauss-Krueger TM

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

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

429 else: # isSpherical

430 g, k = _0_0, _1_0

431

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

433

434

435 def _C(a, b0, b1, Cn):

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

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

438

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

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

441

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

443 but stil mentioned in Note 4 of the comments before

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

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

446 '''

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

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

449 return complex(r, j)

450

451

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

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

454 '''

455 X, m = Xorder.popitem()

456 if m in _Coeffs and isint(m):

457 return m

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

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

460

461

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

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

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

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

466 '''

467 Cs = _Coeffs[m]

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

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

470 n = n_ = E.n

471 if n: # isEllipsoidal

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

473 # constant, it cancels when evaluating a definite

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

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

476 # in C{rhumbx._sincosSeries}.

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

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

479 j = i + r + 1

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

481 i = j + 1

482 n_ *= n

483 X = tuple(X)

484 else: # isSpherical

485 X = _0_0s(m + 1)

486 return X

487

488

489# _Alp- and _BetCoeffs in .rhumbx

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

491 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4

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

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

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

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

496 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5

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

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

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

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

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

502 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6

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

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

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

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

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

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

509 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7

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

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

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

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

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

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

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

517 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8

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

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

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

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

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

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

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

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

526}

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

528 2: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2

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

530 3: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3

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

532 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4

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

534}

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

536 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4

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

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

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

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

541 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5

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

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

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

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

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

547 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6

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

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

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

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

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

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

554 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7

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

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

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

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

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

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

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

562 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8

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

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

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

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

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

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

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

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

571}

572

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

574

575 if __name__ == '__main__':

576

577 from pygeodesy.interns import _usage

578 from sys import argv, exit as _exit

579

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

581

582# **) MIT License

583#

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

585#

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

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

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

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

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

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

592#

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

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

595#

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

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

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

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

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

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

602# OTHER DEALINGS IN THE SOFTWARE.