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, 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 # in KTransverseMercator.ellipsoid.setter

49 from pygeodesy.errors import _ValueError, _xkwds_get, _Xorder

50 from pygeodesy.fmath import fsum1_, hypot, hypot1

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

52 from pygeodesy.interns import NN, _COMMASPACE_, _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.07.01'

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 = \

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

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

407 if isodd(n):

408 Cn = C[n]

409 y0 = complex(Cn) # +j0

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

411 n -= 1

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

413 _c = _C

414 while n > 0:

415 Cn = C[n]

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

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

418 n -= 1

419 Cn = C[n]

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

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

422 n -= 1

423 # assert n == 0

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

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

426

427 # Gauss-Schreiber to Gauss-Krueger TM

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

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

430 else: # isSpherical

431 g, k = _0_0, _1_0

432

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

434

435

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

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

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

439

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

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

442

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

444 but stil mentioned in Note 4 of the comments before

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

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

447 '''

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

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

450 return complex(r, j)

451

452

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

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

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

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

457 '''

458 Cs = _Coeffs[m]

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

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

461 n = n_ = E.n

462 if n: # isEllipsoidal

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

464 # constant, it cancels when evaluating a definite

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

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

467 # in C{rhumbx._sincosSeries}.

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

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

470 j = i + r + 1

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

472 i = j + 1

473 n_ *= n

474 X = tuple(X)

475 else: # isSpherical

476 X = _0_0s(m + 1)

477 return X

478

479

480# _Alp- and _BetCoeffs in .rhumbx

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

482 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4

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

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

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

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

487 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5

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

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

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

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

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

493 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6

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

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

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

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

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

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

500 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7

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

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

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

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

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

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

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

508 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8

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

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

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

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

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

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

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

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

517}

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

519 2: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2

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

521 3: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3

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

523 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4

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

525}

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

527 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4

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

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

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

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

532 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5

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

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

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

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

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

538 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6

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

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

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

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

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

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

545 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7

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

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

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

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

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

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

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

553 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8

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

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

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

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

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

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

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

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

562}

563

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

565

566 if __name__ == '__main__':

567

568 from pygeodesy.interns import _usage

569 from sys import argv, exit as _exit

570

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

572

573# **) MIT License

574#

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

576#

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

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

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

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

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

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

583#

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

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

586#

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

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

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

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

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

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

593# OTHER DEALINGS IN THE SOFTWARE.