Coverage for pygeodesy/triaxials/conformal3.py: 88%

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

4u'''I{Jacobi Conformal projection} classes L{Conformal3}, L{Conformal3B} and L{Conformal3Sphere} on

5triaxial ellipsoids and spheres using L{Ang}, L{Deg}, L{Rad} lat-, longitude, heading and meridian

6(convergence) angles.

8Transcoded to pure Python from I{Karney}'s GeographicLib 2.7 C++ class U{Conformal3<https://

9GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Triaxial_1_1Conformal3.html>}.

11Copyright (C) U{Charles Karney<malto:Karney@Alum.MIT.edu>} (2024-2025) and licensed under the MIT/X11

12License. For more information, see the U{GeographicLib 2.7<https://GeographicLib.SourceForge.io/>}

13documentation.

14'''

15# make sure int/int division yields float quotient, see .basics

16from __future__ import division as _; del _ # noqa: E702 ;

18from pygeodesy.angles import Ang, _Ang3Tuple, isAng, Property_RO

19from pygeodesy.basics import _copysign, signBit

20from pygeodesy.constants import EPS, INF, NAN, PI_2, \

21 _copysign_1_0, _over, _1_over, remainder, \

22 _0_0, _0_01, _0_5, _0_75, _1_0, _2_0, _4_0

23from pygeodesy.constants import _16_0 # PYCHOK used!

24# from pygeodesy.elliptic import Elliptic # _MODS

25from pygeodesy.fmath import hypot1, _ALL_LAZY, Fmt, _MODS

26from pygeodesy.interns import _DMAIN_, _scale_

27# from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS # from .fmath

28# from pygeodesy.named import _NamedTuple, _Pass # from .namedTuples

29from pygeodesy.namedTuples import Vector2Tuple, _Pass

30# from pygeodesy.props import Property_RO # from .triaxials.angles

31# from pygeodesy.streprs import Fmt # from .fmath

32from pygeodesy.triaxials.bases import _bet_, Conformal5Tuple, LLK, _llk_, \

33 _MAXIT, _omg_, TriaxialError, \

34 _Triaxial3Base, Vector3d

35from pygeodesy.triaxials.triaxial3 import Triaxial3B

36from pygeodesy.units import Easting, Northing, Radians, Radius_, Scalar

37from pygeodesy.utily import sincos2

38# from pygeodesy.vector3d import Vector3d # from .triaxials.bases

40from math import atan, atanh, exp, fabs, log, sinh, sqrt

42__all__ = _ALL_LAZY.triaxials_conformal3

43__version__ = '26.02.15'

45_gam_ = 'gam'

46_LOG2MIN = log(EPS) * _2_0

47# _N = (log(_4_0) - log(PI)) / (PI_2 - log(_4_0))

48_NLOG2 = -log(_2_0)

49# _B = exp(_N * PI_2) - pow(_4_0, _N)

52class BetOmgGam5Tuple(_Ang3Tuple):

53 '''5-Tuple C{(bet, omg, gam, scale, llk)} with I{ellipsoidal} lat-

54 C{bet}, longitude C{omg} and meridian convergence C{gam} all

55 L{Ang}les, C{scale} I{and kind} C{llk} I{set to} C{LLK.CONFORMAL}.

56 '''

57 _Names_ = (_bet_, _omg_, _gam_, _scale_, _llk_)

58 _Units_ = ( Ang, Ang, _Pass, Scalar, _Pass)

60 def __new__(cls, bet, omg, gam, scale=NAN, llk=None, unit=None, **kwds): # **iteration_name

61 args = bet, omg, gam, scale, (llk or LLK.CONFORMAL)

62 self = _Ang3Tuple.__new__(cls, args, **kwds)

63 return self.toUnit(unit) if unit else self

66class Conformal3(_Triaxial3Base):

67 '''I{Jacobi Conformal} projection of triaxial ellipsoid using class L{Ang}

68 lat- and longitudes.

70 @see: L{Triaxial<triaxials.triaxial5.Triaxial>} for details.

71 '''

72 @Property_RO

73 def _a_b(self):

74 return self.a / self.b

76 @Property_RO

77 def _a2_b(self):

78 return self.a * self._a_b # == self.b * self._a2_b2

80 @Property_RO

81 def _b_a(self):

82 return self.b / self.a # = sqrt(self._b2_a2)

84 @Property_RO

85 def _b_c(self):

86 return self.b / self.c

88 @Property_RO

89 def _c_b(self):

90 return self.c / self.b # = sqrt(self._c2_b2)

92 @Property_RO

93 def _c2_b(self):

94 return self.c * self._c_b # == self.b * self._c2_b2

96 @Property_RO

97 def _C3S(self):

98 '''(INTERNAL) Cache the I{equivalent Conformal Sphere}.

99 '''

100 return self.equi3Sphere(*self.xyQ2, name=self.name)

101

102 def equi3Sphere(self, x, y, **name):

103 '''Get this projection's I{equivalent Conformal Sphere}.

104

105 @arg x: Quadrant x length, easting (C{meter}).

106 @arg y: Quadrant y length, northing (C{meter}).

107

108 @return: The C{Comformal3Sphere} of this projection.

109

110 @see: Classes L{Conformal3Sphere<triaxials.spheres.Conformal3Sphere>} and

111 L{ConformalSphere<triaxials.triaxial5.ConformalSphere>} and I{Karney}'s

112 GeographicLib 2.7 C++ U{Triaxial::Conformal3<https://GeographicLib.

113 SourceForge.io/C++/doc/classGeographicLib_1_1Triaxial_1_1Conformal3.html>}

114 for more information.

115 '''

116 x, y = self.xyQ2

117 # Find C{b, k2, kp2} s.t. C{b * K(kp2) = x},

118 # C{b * K(k2) = y} and C{x * K(k2) - y * K(kp2) = 0}.

119 _EF = _MODS.elliptic.Elliptic

120 _xy = x < y

121 if _xy:

122 x, y = y, x

123 # Now x >= y, k2 <= 1/2

124 if x != y:

125 s = x + y

126 ny = _over(y, s)

127 if ny:

128 nx = _over(x, s)

129 # assert nx != ny

130 # Find initial guess assume K(k2) = pi/2, so K(kp2) = nx/ny * pi/2.

131 # Invert using approximate k(K) given in https://arxiv.org/abs/2505.17159v4

132 KK = _over(nx, ny) * PI_2

133 # k2 = _16_0 / pow(exp(_N * KK) - _B, _2_0 / _N)

134 # Alternatively using KK = 1/2*log(16/kp) A+S 17.3.26

135 k2 = min(_0_5, exp(-KK * _2_0) * _16_0) # Make sure guess is sane

136 logk2 = log(k2)

137 if logk2 > _LOG2MIN:

138 # Solve for log(k2) to preserve relative accuracy for tiny k2.

139 def _k2s2(logk2):

140 k2 = exp(logk2)

141 kp2 = _1_0 - k2

142 xS = _EF(kp2, 0, k2)

143 yS = _EF(k2, 0, kp2)

144 f = nx * yS.cK - ny * xS.cK

145 fp = (nx * (yS.cE - kp2 * yS.cK) +

146 ny * (xS.cE - k2 * xS.cK)) / (kp2 * _2_0)

147 return f, fp

148

149 logk2, _, _, _ = _root4(_k2s2, 0, logk2, _LOG2MIN, _NLOG2)

150 k2 = exp(logk2)

151 # otherwise accept the asymptotic result

152 kp2 = _1_0 - k2

153 else:

154 k2, kp2 = _0_0, _1_0

155 else:

156 k2 = kp2 = _0_5

157 # b * K(kp2) = x

158 # b * K(k2) = y

159

160 K = _EF(k2, 0, kp2).cK

161 b = _over(y, K)

162 if _xy:

163 k2, kp2 = kp2, k2

164 return Conformal3Sphere(b, k2, kp2, **name)

165

166 def forwardBetOmg(self, bet, omg, M=False, **unit_name):

167 '''Compute the projection to this conformal triaxial.

168

169 @arg bet: Ellipsoidal latitude (C{Ang} or B{C{unit}}).

170 @arg omg: Ellipsoidal longitude (C{Ang} or B{C{unit}})..

171 @kwarg M: If C{True}, compute and include the scale (C{bool}).

172 @kwarg unit_name: Optional C{B{unit}=}L{Radians} and C{B{name}=NN} (C{str}).

173

174 @return: A L{Conformal5Tuple}C{(x, y, z, scale, llk)} with C{z = INT0}

175 I{always} and C{scale} if C{B{M}=True}, otherwise C{scale = NAN}.

176 '''

177 bet, omg, name = _bet_omg_name(bet, omg, **unit_name)

178 m = _1_over(self._invScale(bet, omg)) if M else NAN

179 return Conformal5Tuple(self._x(omg), self._y(bet), scale=m, **name)

180

181 def forwardOther(self, other, bet, omg, M=False, **unit_name):

182 '''Compute the projection to an other conformal triaxial.

183

184 @arg other: A conformal triaxial (C{Conformal3}).

185 @arg bet: Ellipsoidal latitude (C{Ang} or B{C{unit}}).

186 @arg omg: Ellipsoidal longitude (C{Ang} or B{C{unit}}).

187 @kwarg M: If C{True}, compute and include the scale (C{bool}).

188 @kwarg unit_name: Optional C{B{unit}=}L{Radians} and C{B{name}=NN} (C{str}).

189

190 @return: A L{BetOmgGam5Tuple}C{(bet, omg, gam, scale, llk)} with C{scale}

191 if C{B{M}=True}, otherwise C{scale = NAN}.

192 '''

193 if not isinstance(other, Conformal3):

194 raise TriaxialError(other=other)

195 bet, omg, name = _bet_omg_name(bet, omg, **unit_name)

196 m = other._C3S.b / self._C3S.b

197 ct, v, ma = self.forwardSphere3(bet, omg, M=M)

198 ct = Vector3d(ct).times(m)

199 bet, omg, gam, mb, _ = other.reverseSphere(ct, dir3d=v, M=M)

200 if M:

201 m *= _over(ma, mb)

202 else:

203 m = NAN

204 return BetOmgGam5Tuple(bet, omg, gam, m, **name)

205

206 def forwardSphere3(self, bet, omg, M=False, **unit_name):

207 '''Compute the projection to and direction on the equivalent I{Conformal

208 Sphere}.

209

210 @arg bet: Ellipsoidal latitude (C{Ang} or B{C{unit}}).

211 @arg omg: Ellipsoidal longitude (C{Ang} or B{C{unit}}).

212 @kwarg M: If C{True}, compute and include the scale (C{bool}).

213 @kwarg unit_name: Optional C{B{unit}=}L{Radians} and C{B{name}=NN} (C{str}).

214

215 @return: 3-Tuple C{(cartesian, direction, scale)} with a C{cartesian}

216 L{Cartesian5Tuple}C{(x, y, z, h, llk)} on and C{direction} a

217 C{Vector3d} due North and tangent to the I{Conformal Sphere}

218 and C{scale} if C{B{M}=True}, otherwise C{scale = NAN}.

219 '''

220 if M:

221 bet, omg, _ = _bet_omg_name(bet, omg, **unit_name)

222 x, y, _, _, _ = self.forwardBetOmg(bet, omg, M=False, **unit_name)

223 S = self._C3S

224 bets = _invF(S._yE, y / S.b)

225 omgs = _invF(S._xE, x / S.b).shift(-1)

226 ct, dir3d = S.forwardBetOmgAlp2(bets, omgs, Ang.N(), **unit_name)

227 if M:

228 ma = self._invScale(bet, omg)

229 mb = _invScale(S, bets, omgs)

230 m = self._sphScale(ma, mb)

231 else:

232 m = NAN

233 return ct, dir3d, m

234

235 def _invScale(self, bet, omg):

236 # scale helper

237 return _invScale(self, bet, omg.shift(+1))

238

239 def reverseBetOmg(self, x_cf, y=None, M=False, **unit_name):

240 '''Reverse a projection from this conformal triaxial.

241

242 @arg x_cf: Easting (C{scalar}) or a conformal projection (C{Conformal5Tuple}).

243 @kwarg y: Northing (C{scalar}), required if B{C{x_cf}} is C{scalar},

244 ignored otherwise.

245 @kwarg M: If C{True}, compute and include the scale (C{bool}).

246 @kwarg unit_name: Optional C{B{unit}=}L{Radians} and C{B{name}=NN} (C{str}).

247

248 @return: A L{BetOmgGam5Tuple}C{(bet, omg, gam, scale, llk)} with

249 C{gam} set to C{None} and C{scale} only if C{B{M}=True},

250 otherwise C{scale is NAN}.

251 '''

252 x, y = _cf2en(x_cf, y)

253 bet = _invPi(self._yE, y / self._c2_b).mod(self._c_b)

254 omg = _invPi(self._xE, x / self._a2_b).mod(self._a_b).shift(-1)

255 m = _1_over(self._invScale(bet, omg)) if M else NAN

256 return BetOmgGam5Tuple(bet, omg, None, m, **unit_name)

257

258 def reverseOther(self, other, beto, omgo, M=False, **unit_name):

259 '''Reverse a projection from an other conformal triaxial.

260

261 @arg other: A conformal triaxial (C{Conformal3}).

262 @arg beto: Ellipsoidal latitude on the B{C{other}} triaxial

263 (C{Ang} or B{C{unit}}).

264 @arg omgo: Ellipsoidal longitude on the B{C{other}} triaxial

265 (C{Ang} or B{C{unit}}).

266 @kwarg M: If C{True}, compute and include the scale (C{bool}).

267 @kwarg unit_name: Optional C{B{unit}=}L{Radians} and C{B{name}=NN} (C{str}).

268

269 @return: A L{BetOmgGam5Tuple}C{(bet, omg, gam, scale, llk)} with

270 C{scale} if C{B{M}=True}, otherwise C{scale = NAN}.

271 '''

272 if not isinstance(other, Conformal3):

273 raise TriaxialError(other=other)

274 bet, omg, gam, m, _ = other.forwardOther(self, beto, omgo, M=M)

275 m = _1_over(m) if M else NAN

276 return BetOmgGam5Tuple(bet, omg, _Neg(gam), m, **unit_name)

277

278 def reverseSphere(self, x_ct, y=None, z=None, dir3d=None, M=False, **unit_name):

279 '''Reverse a projection from this C{Conformal Sphere} to this triaxial.

280

281 @arg x_ct: X component (C{scalar}) of or a cartesian on the C{Conformal

282 Sphere} (C{Cartesian5Tuple}).

283 @kwarg y: Y component (C{scalar}), required if B{C{x_ct}} is C{scalar},

284 ignored otherwise.

285 @kwarg z: Z component (C{scalar}), like B{C{y}}.

286 @kwarg dir3d: The direction (C{Vector3d} or C{None}), reference.

287 @kwarg M: If C{True}, compute and include the scale (C{bool}).

288 @kwarg unit_name: Optional C{B{unit}=}L{Radians} and C{B{name}=NN} (C{str}).

289

290 @return: A L{BetOmgGam5Tuple}C{(bet, omg, gam, scale, llk)} with

291 C{scale} only if C{B{M}=True}, otherwise C{scale = NAN}.

292 '''

293 S = self._C3S

294 bets, omgs, alp, _, _ = S.reverseBetOmgAlp(x_ct, y, z, dir3d=dir3d)

295 _ = Ang._norm(bets, omgs, alp, alt=bool(S.k2 == 0))

296 x = S.b * _F(S._xE, omgs.shift(+1))

297 y = S.b * _F(S._yE, bets)

298 bet, omg, _, _, _ = self.reverseBetOmg(x, y, M=False)

299 if M:

300 mb = _invScale(S, bets, omgs)

301 ma = _invScale(S, bet, omg)

302 m = self._sphScale(ma, mb)

303 else:

304 m = NAN

305 return BetOmgGam5Tuple(bet, omg, _Neg(alp), m, **unit_name)

306

307 def _sphScale(self, ma, mb):

308 '''(INTERNAL) Compute the C{scale}.

309 '''

310 return _over(mb, ma) if ma and mb else self._sphScale_m

311

312 @Property_RO

313 def _sphScale_m(self):

314 '''(INTERNAL) Cache the constant C{scale}.

315 '''

316 e = sqrt(self.e2)

317

318 def _h1sg(atan_):

319 sg = sinh(e * atan_(e))

320 return hypot1(sg) + sg

321

322 k2, kp2 = self._k2_kp2

323 if not kp2: # oblate pole

324 m = self._c_b * _h1sg(atanh)

325 elif not k2: # prolate pole

326 m = _over(self._a_b, _h1sg(atan))

327 else: # trixial umbilical

328 S = self._C3S

329 s = _over(S.k2 * S.kp2, k2 * kp2)

330 m = (sqrt(s) * (self.b / S.b)) if s > 0 else _0_0

331 return m

332

333 def _x(self, omg): # degrees?

334 '''(INTERNAL) Compute an C{x} projection easting.

335 '''

336 omg = omg.shift(+1)

337 _, c, n = omg.scn3

338 if (n or signBit(c)) and not self.k2:

339 x = NAN

340 else:

341 x = _Pi(self._xE, omg.mod(self._b_a))

342 x *= self._a2_b

343 return x

344

345 @Property_RO

346 def xyQ2(self):

347 '''Get the quadrant length in x and y direction (Vector2Tuple{x, y}).

348 '''

349 return Vector2Tuple(self._a2_b * self._xE.cPi,

350 self._c2_b * self._yE.cPi,

351 name=Conformal3.xyQ2.name)

352

353 def _y(self, bet): # degrees?

354 '''(INTERNAL) Compute a C{y} projection northing.

355 '''

356 _, c, n = bet.scn3

357 if (n or signBit(c)) and not self.kp2:

358 y = NAN

359 else:

360 y = _Pi(self._yE, bet.mod(self._b_c))

361 y *= self._c2_b

362 return y

363

364

365class Conformal3B(Conformal3):

366 '''I{Jacobi Conformal projection} on a triaxial ellipsoid

367 specified by its middle semi-axis and shape.

368

369 @see: L{Conformal3} for details and more information.

370 '''

371 def __init__(self, b, e2=_0_0, k2=_1_0, kp2=_0_0, **name):

372 '''New L{Conformal3B} triaxial.

373

374 @note: Use C{B{b}=radius} and C{B{e2}=0} for a conformal

375 I{spherical} projection.

376

377 @see: L{Conformal<triaxials.triaxial5.Conformal.__init__>}.

378 '''

379 self._init_abc3_e2_k2_kp2(Radius_(b=b), e2, k2, kp2, **name)

380

381

382class Conformal3Sphere(Triaxial3B): # note C{Triaxial3}!

383 '''I{Jacobi Conformal projection} on a I{spherical} triaxial.

384

385 @see: Method L{equiv3Sphere<triaxials.conformal3.Conformal3.equi3Sphere>}.

386 '''

387 def __init__(self, radius, k2=_1_0, kp2=_0_0, **name):

388 '''New, L{Conformal3Sphere} instance.

389

390 @see: L{Triaxial3<triaxials.triaxial3.Triaxial3>} for more information.

391 '''

392 self._init_abc3_e2_k2_kp2(Radius_(radius), 0, k2, kp2, **name)

393 # self._xE = _MODS.elliptic.Elliptic(kp2, 0, k2, **name) # _Triaxial3Base._xE

394 # self._yE = _MODS.elliptic.Elliptic(k2, 0, kp2, **name) # _Triaxial3Base._yE

395

396

397def _bet_omg_name(bet, omg, unit=Radians, **name):

398 '''(INTERNAL) Get C{(bet, omg, name)}.

399 '''

400 return (Ang.fromScalar(bet, unit=unit),

401 Ang.fromScalar(omg, unit=unit), name)

402

403

404def _cf2en(x_cf, y):

405 '''(INTERNAL) Get easting C{x} and notrthing C{y}.

406 '''

407 return x_cf[:2] if isinstance(x_cf, Conformal5Tuple) else (

408 Easting(x_cf), Northing(y))

409

410

411def _F(eF, phi): # -> float

412 '''(INTERNAL) Elliptic function C{F(phi)}.

413 '''

414 s, c, n = phi.scn3

415 if (n or signBit(c)) and not eF.kp2:

416 p = NAN

417 else:

418 p = eF.fF(s, c, eF.fDelta(s, c))

419 if n:

420 p += eF.cK * n * _4_0

421 return p

422

423

424def _invF(eF, x): # -> Ang

425 '''(INTERNAL) Inverse elliptic function C{F(x)}.

426 '''

427 r, y = _invRy2(eF, eF.cK, x)

428 if y: # solve eF.fF(phi) = y for phi

429 def _fF2(phi): # -> pair<real, real>

430 s, c = sincos2(phi)

431 f = eF.fF(s, c, eF.fDelta(s, c))

432 fp = sqrt(eF.kp2 + c**2 * eF.k2)

433 return f, _1_over(fp)

434

435 z = fabs(y)

436 z, _, _, _ = _root4(_fF2, z, z * PI_2 / eF.cK)

437 r += Ang.fromRadians(_copysign(z, y))

438 return r

439

440

441def _invPi(eF, x): # -> Ang

442 '''(INTERNAL) Inverse elliptic function C{Pi(x)}.

443 '''

444 r, y = _invRy2(eF, eF.cPi, x)

445 if y: # solve eF.fPi(phi) = y for phi

446 def _fPi2(phi): # -> pair<real, real>

447 s, c = sincos2(phi)

448 f = eF.fPi(s, c, eF.fDelta(s, c))

449 a2 = eF.alpha2

450 fp = (_1_0 - c**2 * a2) if a2 < 0 else \

451 (eF.alphap2 + s**2 * a2)

452 fp *= sqrt(eF.kp2 + c**2 * eF.k2)

453 return f, _1_over(fp)

454

455 z = fabs(y)

456 z, _, _, _ = _root4(_fPi2, z, z * PI_2 / eF.cPi)

457 r += Ang.fromRadians(_copysign(z, y))

458 return r

459

460

461def _invRy2(eF, eF_c, x):

462 # helper for C{_invF} and C{_invPi}

463 if eF.kp2:

464 n = eF_c * _2_0

465 y = remainder(x, n)

466 n = round((x - y) / n) * _2_0

467 else: # eF_c == N-/INF

468 y, n = x, _0_0

469 if not y:

470 y, n = 0, (n if n else y) # signed 0

471 elif fabs(y) == eF_c:

472 y, n = 0, (_copysign_1_0(y) + n)

473 return Ang.cardinal(n), y

474

475

476def _invScale(triax, bet, omg):

477 # helper for triaxial, sphere scale

478 k2, kp2 = triax._k2_kp2

479 return sqrt(k2 * bet.c**2 + kp2 * omg.c**2)

480

481

482def _Neg(ang):

483 return (-ang) if isAng(ang) else (ang or None)

484

485

486def _Pi(eF, phi): # -> float

487 '''(INTERNAL) Elliptic function C{Pi(phi)}.

488 '''

489 s, c, n = phi.scn3

490 if (n or signBit(c)) and not eF.kp2:

491 p = NAN

492 else:

493 p = eF.fPi(s, c, eF.fDelta(s, c))

494 if n:

495 p += eF.cPi * n * _4_0

496 return p

497

498

499def _root4(_fp2, z, x, xa=_0_0, xb=PI_2, xscale=1, zscale=1, s=1, tol=EPS): # int s

500 '''(INTERNAL) Solve v = _fp2(x) - z = 0.

501 '''

502 k = b = C = 0

503 if xa < xb and xa <= x <= xb:

504 # p = PI_2 * 0 #???

505 # def _fp2z(x):

506 # f, fp = _fp2(x)

507 # f -= z

508 # # "DAT ", x, f, fp

509 # return f

510 # a, _, b = map1(_fp2z, xa, x, xb)

511 # if (a * b) > 0:

512 # raise TriaxalError('"DATBAD")

513 # tol = max(tol, EPS) # tol if tol > 0 else EPS

514 vtol = tol * zscale * _0_01

515 xtol = pow(tol, _0_75) * xscale

516 oldv = oldx = olddx = INF

517 for k in range(1, _MAXIT):

518 # TODO: 20 60 -90 180 127.4974 24.6254 2.4377

519 v, vp = _fp2(x)

520 v -= z

521 va = fabs(v)

522 dx = _over(-v, vp)

523 # "XX ", k, (xa - p), (x - p), (xb - p), dx, (x + dx - p), v, vp

524 if not (va > (0 if k < 2 else vtol)):

525 C = 1 # k, va

526 break

527 elif (v * s) > 0:

528 xb = min(xb, x)

529 else:

530 xa = max(xa, x)

531 x += dx

532 dxa = fabs(dx)

533 if x < xa or x > xb or va > oldv or \

534 (k > 2 and dxa > olddx):

535 b += 1 # k, xa, x, xb

536 x = (xa + xb) * _0_5

537 if x == oldx:

538 C = 3 # k, x, dx

539 break

540 elif not dxa > xtol:

541 C = 2 # k, dx, xtol

542 break

543 # "GAPS ", k, dx, (x - xa), (xb - x), oldx, x, (oldx - x)

544 oldx = x

545 oldv = va

546 olddx = dxa * _0_5

547 else:

548 t = Fmt.no_convergence(dx, xtol)

549 raise TriaxialError(x=x, xa=xa, xb=xb, txt=t)

550 else:

551 x = NAN

552 return x, k, b, C

553

554

555if __name__ == _DMAIN_:

556

557 from pygeodesy import Degrees, printf

558 from pygeodesy.triaxials import Triaxials

559

560 # <https://GeographicLib.SourceForge.io/C++/doc/Cart3Convert.1.html>

561 T = Conformal3(Triaxials.WGS84_3)

562 printf(T)

563 # name='WGS84_3', a=6378171.36, b=6378101.609999999, c=6356751.84, ...

564 t = T.forwardBetOmg(Ang.fromDegrees(33.3), Ang.fromDegrees(44.4), M=True)

565 printf((t.x, t.y, t.scale))

566 # (-5077726.431878843, 3922574.862034525, 1.1970322930902468)

567 # -5077732.396 3922571.859 1.1970343759 C++

568 t = T.reverseBetOmg(*t[:2], M=True)

569 printf((t.bet.degrees0, t.omg.degrees0, t.scale))

570 # (33.300000000000004, 44.39999999999999, 1.1970322930902468)

571 # 33.30000000 44.40000000 1.1970343759 C++

572

573 T = Conformal3(Triaxials.WGS84_3r) # rounded

574 printf(T)

575 # name='WGS84_3r', a=6378172, b=6378102, c=6356752, ...

576 t = T.forwardBetOmg(Degrees(33.3), Degrees(44.4), M=True)

577 printf((t.x, t.y, t.scale))

578 # (-5077732.396020001, 3922571.859124643, 1.197034375926918)

579 # -5077732.396 3922571.859 1.1970343759 C++

580 t = T.reverseBetOmg(*t[:2], M=True)

581 printf((t.bet.degrees0, t.omg.degrees0, t.scale))

582 # (33.3, 44.400000000000006, 1.197034375926918)

583 # 33.30000000 44.40000000 1.1970343759 C++

584

585 c = 6378137 * (1 - 1 / (298257223563 / 1000000000))

586 T = Conformal3(6378172, 6378102, c)

587 printf(T)

588 # name='', a=6378172, b=6378102, c=6356752.314245179, ...

589 t = T.forwardBetOmg(Degrees(33.3), Degrees(44.4), M=True)

590 printf((t.x, t.y, t.scale))

591 # (-5077732.418682504, 3922572.0186951878, 1.197034384522207)

592 # -5077732.396 3922571.859 1.1970343759 C++

593 t = T.reverseBetOmg(*t[:2], M=True)

594 printf((t.bet.degrees0, t.omg.degrees0, t.scale))

595 # (33.300000000000004, 44.4, 1.197034384522207)

596 # 33.30000000 44.40000000 1.1970343759 C++

597

598# % python3 -m pygeodesy.triaxials.conformal3

599

600# name='WGS84_3', a=6378171.36, b=6378101.609999999, c=6356751.84, k2=0.996738165, kp2=0.003261835, volume=1083207064030173855744, area=510065541435967.5, R2=6371006.679496506

601# (-5077726.431878843, 3922574.862034525, 1.1970322930902468)

602# (33.300000000000004, 44.39999999999999, 1.1970322930902468)

603

604# name='WGS84_3r', a=6378172, b=6378102, c=6356752, k2=0.996726547, kp2=0.003273453, volume=1083207266220584468480, area=510065604942135.875, R2=6371007.076110449

605# (-5077732.396020001, 3922571.859124643, 1.197034375926918)

606# (33.3, 44.400000000000006, 1.197034375926918)

607

608# name='', a=6378172, b=6378102, c=6356752.314245179, k2=0.996726499, kp2=0.003273501, volume=1083207319768789942272, area=510065621722018.25, R2=6371007.180905545

609# (-5077732.418682504, 3922572.0186951878, 1.197034384522207)

610# (33.300000000000004, 44.4, 1.197034384522207)

611

612# **) MIT License

613#

615#

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

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

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

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

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

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

622#

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

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

625#

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

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

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

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

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

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

632# OTHER DEALINGS IN THE SOFTWARE.