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:Karney@Alum.MIT.edu>} (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, isodd, neg, neg_, _reverange

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 # _MODS

49# from pygeodesy.ellipsoids import _EWGS84 # from .karney

50 from pygeodesy.errors import _ValueError, _xkwds_get, _Xorder

51 from pygeodesy.fmath import fsum1_, hypot, hypot1

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

53 from pygeodesy.interns import NN, _COMMASPACE_, _singular_

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

55 _norm180, _polynomial, _unsigned2, _EWGS84

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

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

58 from pygeodesy.namedTuples import Forward4Tuple, Reverse4Tuple

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

60 _update_all

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

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

63 from pygeodesy.utily import atand, sincos2, sincos2d_

64

65 from cmath import phase

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

67

68__all__ = _ALL_LAZY.ktm

69__version__ = '23.08.20'

70

71

72 class KTMError(_ValueError):

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

74 '''

75 pass

76

77

78 class 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 _E = _EWGS84

111 _k0 = _K0_UTM # central scale factor

112 _lon0 = _0_0 # central meridian

113 _mTM = 6

114 _raiser = False # throw Error

115

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

117 raiser=False, **TMorder):

118 '''New L{KTransverseMercator}.

119

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

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

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

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

124 a C{scalar}, ignored otherwise.

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

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

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

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

129 singularities (C{bool}).

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

131

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

133 '''

134 if f is not None:

135 self.ellipsoid = a_earth, f

136 elif a_earth not in (_EWGS84, None):

137 self.ellipsoid = a_earth

138 self.lon0 = lon0

139 self.k0 = k0

140 if name: # PYCHOK no cover

141 self.name = name

142 if raiser:

143 self.raiser = True

144 if TMorder:

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

146

147 @Property_RO

148 def _Alp(self):

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

150

151 @Property_RO

152 def _b1(self):

153 n = self.ellipsoid.n

154 if n: # isEllipsoidal

155 m = self.TMorder // 2

156 B1 = _B1Coeffs[m]

157 m += 1

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

159 else: # isSpherical

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

161 return b1

162

163 @Property_RO

164 def _Bet(self):

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

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

167

168 @Property

169 def ellipsoid(self):

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

171 '''

172 return self._E

173

174 @ellipsoid.setter # PYCHOK setter!

175 def ellipsoid(self, a_earth_f):

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

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

178 '''

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

180 if self._E != E:

181 _update_all(self)

182 self._E = E

183

184 @Property_RO

185 def equatoradius(self):

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

187 '''

188 return self.ellipsoid.a

189

190 a = equatoradius

191

192 @Property_RO

193 def flattening(self):

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

195 '''

196 return self.ellipsoid.f

197

198 f = flattening

199

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

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

202

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

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

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

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

207

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

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

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

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

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

213

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

215 C{True}.

216 '''

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

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

219 lon, _lon = _unsigned2(lon)

220 backside = lon > 90

221 if backside: # PYCHOK no cover

222 lon = _180_0 - lon

223 if lat == 0:

224 _lat = True

225

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

227 E = self.ellipsoid

228 if cphi and lat != 90:

229 t = sphi / cphi

230 tp = E.es_taupf(t)

231 h = hypot(tp, clam)

232 if h:

233 xip = atan2(tp, clam)

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

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

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

237 elif self.raiser:

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

239 else: # PYCHOK no cover

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

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

242 else: # PYCHOK no cover

243 xip, etap = PI_2, _0_0

244 g, k = lon, E.es_c

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

246 g -= t

247 k *= z * self._k0_b1

248

249 if backside: # PYCHOK no cover

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

251 y *= self._k0_a1

252 x *= self._k0_a1

253 if _lat:

254 y, g = neg_(y, g)

255 if _lon:

256 x, g = neg_(x, g)

257

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

259

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

261 def k0(self):

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

263 '''

264 return self._k0 # aka scale0

265

266 @k0.setter # PYCHOK setter!

267 def k0(self, k0):

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

269

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

271 '''

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

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

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

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

276 self._k0 = k0

277

278 @Property_RO

279 def _k0_a1(self):

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

281 '''

282 return self._k0_b1 * self.equatoradius

283

284 @Property_RO

285 def _k0_b1(self):

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

287 '''

288 return self.k0 * self._b1

289

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

291 def lon0(self):

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

293 '''

294 return self._lon0

295

296 @lon0.setter # PYCHOK setter!

297 def lon0(self, lon0):

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

299

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

301 '''

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

303

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

305 def raiser(self):

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

307 '''

308 return self._raiser

309

310 @raiser.setter # PYCHOK setter!

311 def raiser(self, raiser):

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

313 for C{forward} singularities.

314 '''

315 self._raiser = bool(raiser)

316

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

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

319

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

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

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

323

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

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

326 '''

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

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

329 backside = xi > PI_2

330 if backside: # PYCHOK no cover

331 xi = PI - xi

332

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

334 t = self._k0_b1

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

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

337 if c > 0:

338 r = hypot(h, c)

339 else: # PYCHOK no cover

340 r = fabs(h)

341 c = _0_0

342 E = self.ellipsoid

343 if r:

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

345 s = sin(xip) # Newton for tau

346 t = E.es_tauf(s / r)

347 lat = atand(t)

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

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

350 else: # PYCHOK no cover

351 lat, lon = _90_0, _0_0

352 k *= E.es_c

353

354 if backside: # PYCHOK no cover

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

356 if _lat:

357 lat, g = neg_(lat, g)

358 if _lon:

359 lon, g = neg_(lon, g)

360

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

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

363 name=name or self.name)

364

365 @Property

366 def TMorder(self):

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

368 '''

369 return self._mTM

370

371 @TMorder.setter # PYCHOK setter!

372 def TMorder(self, order):

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

374 '''

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

376 if self._mTM != m:

377 _update_all(self)

378 self._mTM = m

379

380 def toStr(self, **kwds):

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

382

383 @arg kwds: Optional, overriding keyword arguments.

384 '''

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

386 if self.name: # PYCHOK no cover

387 d.update(name=self.name)

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

389

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

391 '''(INTERNAL) Complex Clenshaw summation with

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

393 '''

394 def _sinhcosh2(x):

395 return sinh(x), cosh(x)

396

397 x = complex(xi_, eta_)

398 if self.f: # isEllipsoidal

399 s, c = sincos2( xi_ * 2)

400 sh, ch = _sinhcosh2(eta_ * 2)

401 n = -s

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

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

404

405 y0 = y1 = \

406 z0 = z1 = complex(0) # 0+j0

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

408 if isodd(n):

409 Cn = C[n]

410 y0 = complex(Cn) # +j0

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

412 n -= 1

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

414 _c = _C

415 while n > 0:

416 Cn = C[n]

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

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

419 n -= 1

420 Cn = C[n]

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

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

423 n -= 1

424 # assert n == 0

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

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

427

428 # Gauss-Schreiber to Gauss-Krueger TM

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

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

431 else: # isSpherical

432 g, k = _0_0, _1_0

433

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

435

436

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

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

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

440

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

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

443

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

445 but stil mentioned in Note 4 of the comments before

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

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

448 '''

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

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

451 return complex(r, j)

452

453

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

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

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

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

458 '''

459 Cs = _Coeffs[m]

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

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

462 n = n_ = E.n

463 if n: # isEllipsoidal

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

465 # constant, it cancels when evaluating a definite

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

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

468 # in C{rhumbx._sincosSeries}.

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

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

471 j = i + r + 1

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

473 i = j + 1

474 n_ *= n

475 X = tuple(X)

476 else: # isSpherical

477 X = _0_0s(m + 1)

478 return X

479

480

481# _Alp- and _BetCoeffs in .rhumbx

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

483 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4

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

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

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

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

488 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5

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

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

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

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

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

494 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6

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

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

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

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

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

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

501 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7

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

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

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

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

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

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

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

509 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8

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

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

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

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

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

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

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

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

518}

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

520 2: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2

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

522 3: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3

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

524 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4

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

526}

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

528 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4

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

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

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

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

533 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5

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

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

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

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

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

539 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6

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

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

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

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

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

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

546 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7

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

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

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

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

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

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

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

554 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8

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

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

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

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

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

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

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

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

563}

564

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

566

567 if __name__ == '__main__':

568

569 from pygeodesy.interns import _usage

570 from sys import argv, exit as _exit

571

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

573

574# **) MIT License

575#

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

577#

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

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

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

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

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

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

584#

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

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

587#

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

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

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

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

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

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

594# OTHER DEALINGS IN THE SOFTWARE.