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_, \

46 _reverange, _xinstanceof

47 from pygeodesy.constants import INF, _K0_UTM, PI, PI_2, _0_0s, _0_0, \

48 _1_0, _90_0, _copysignINF

49 from pygeodesy.datums import Datum, _spherical_datum, _WGS84, _EWGS84

50# from pygeodesy.ellipsoids import _EWGS84 # from .datums

51 from pygeodesy.errors import _ValueError, _xkwds_pop2, _Xorder

52 from pygeodesy.fmath import hypot, hypot1

53 from pygeodesy.fsums import fsum1f_

54 from pygeodesy.interns import _COMMASPACE_, _singular_

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

56 _polynomial, _unsigned2

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

58 from pygeodesy.named import _NamedBase, pairs, _ALL_LAZY

59 from pygeodesy.namedTuples import Forward4Tuple, Reverse4Tuple

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

61 _update_all

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

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

64 from pygeodesy.utily import atan1d, _loneg, sincos2, sincos2d_

65

66 from cmath import polar

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

68

69__all__ = _ALL_LAZY.ktm

70__version__ = '24.06.11'

71

72

73 class KTMError(_ValueError):

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

75 '''

76 pass

77

78

79 class 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

429 def _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 but stil mentioned in Note 4 of the comments before

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

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

440 '''

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

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

443 return complex(r, j)

444

445

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

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

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

449 '''

450 x = complex(xi_, eta_)

451 if E.f: # isEllipsoidal

452 s, c = sincos2( xi_ * 2)

453 sh, ch = _sinhcosh2(eta_ * 2)

454 n = -s

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

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

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

458

459 y0 = y1 = \

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

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

462 if isodd(n):

463 Cn = C[n]

464 y0 = complex(Cn) # +0j

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

466 n -= 1

467 _c = _cma

468 while n > 0:

469 Cn = C[n]

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

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

472 n -= 1

473 Cn = C[n]

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

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

476 n -= 1

477 # assert n == 0

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

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

480

481 # Gauss-Schreiber to Gauss-Krueger TM

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

483 k, g = polar(c)

484 g = degrees(g)

485 else: # E.isSpherical

486 g, k = _0_0, _1_0

487

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

489

490

491 def _sinhcosh2(x):

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

493 '''

494 return sinh(x), cosh(x)

495

496

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

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

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

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

501 '''

502 Cs = _Coeffs[m]

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

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

505 n = n_ = E.n

506 if n: # isEllipsoidal

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

508 # constant, it cancels when evaluating a definite

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

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

511 _X, _p = X.append, _polynomial

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

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

514 j = i + r + 1

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

516 i = j + 1

517 n_ *= n

518 X = tuple(X)

519 else: # isSpherical

520 X = _0_0s(m + 1)

521 return X

522

523

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

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

526 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4

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

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

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

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

531 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5

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

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

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

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

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

537 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6

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

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

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

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

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

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

544 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7

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

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

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

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

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

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

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

552 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8

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

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

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

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

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

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

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

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

561}

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

563 2: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2

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

565 3: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3

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

567 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4

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

569}

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

571 4: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4

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

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

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

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

576 5: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5

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

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

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

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

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

582 6: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6

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

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

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

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

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

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

589 7: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7

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

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

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

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

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

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

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

597 8: ( # GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8

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

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

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

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

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

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

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

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

606}

607

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

609

610 if __name__ == '__main__':

611

612 from pygeodesy.internals import _usage

613 from sys import argv, exit as _exit

614

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

616

617# **) MIT License

618#

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

620#

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

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

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

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

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

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

627#

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

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

630#

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

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

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

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

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

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

637# OTHER DEALINGS IN THE SOFTWARE.