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

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, _0_0, _0_01, \

21 _0_5, _0_75, _1_0, _2_0, _4_0, \

22 _copysign_1_0, _over, _1_over

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_

27from pygeodesy.karney import _remainder2, _tauf

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

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

30from pygeodesy.namedTuples import Vector2Tuple, _Pass

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

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

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

34 _MAXIT, _omg_, TriaxialError, \

35 _Triaxial3Base, Vector3d

36from pygeodesy.triaxials.triaxial3 import Triaxial3B

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

38from pygeodesy.utily import sincos2

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

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

43__all__ = _ALL_LAZY.triaxials_conformal3

44__version__ = '26.10.30'

46_gam_ = 'gam'

47_LOG2MIN = log(EPS) * _2_0

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

49_NLOG2 = -log(_2_0)

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

51_useTauf = True # see ._invPI and GeographicLib 2.8 Comformal3.cpp

54class BetOmgGam5Tuple(_Ang3Tuple):

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

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

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

58 '''

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

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

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

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

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

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

68class Conformal3(_Triaxial3Base):

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

70 lat- and longitudes.

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

73 '''

74 @Property_RO

75 def _a_b(self):

76 return self.a / self.b

78 @Property_RO

79 def _a2_b(self):

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

82 @Property_RO

83 def _b_a(self):

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

86 @Property_RO

87 def _b_c(self):

88 return self.b / self.c

90 @Property_RO

91 def _c_b(self):

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

94 @Property_RO

95 def _c2_b(self):

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

98 @Property_RO

99 def _C3S(self):

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

101 '''

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

103

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

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

106

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

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

109

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

111

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

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

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

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

116 for more information.

117 '''

118 x, y = self.xyQ2

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

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

121 _EF = _MODS.elliptic.Elliptic

122 _xy = x < y

123 if _xy:

124 x, y = y, x

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

126 if x != y:

127 s = x + y

128 ny = _over(y, s)

129 if ny:

130 nx = _over(x, s)

131 # assert nx != ny

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

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

134 KK = _over(nx, ny) * PI_2

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

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

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

138 logk2 = log(k2)

139 if logk2 > _LOG2MIN:

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

141 def _k2s2(logk2):

142 k2 = exp(logk2)

143 kp2 = _1_0 - k2

144 xS = _EF(kp2, 0, k2)

145 yS = _EF(k2, 0, kp2)

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

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

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

149 return f, fp

150

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

152 k2 = exp(logk2)

153 # otherwise accept the asymptotic result

154 kp2 = _1_0 - k2

155 else:

156 k2, kp2 = _0_0, _1_0

157 else:

158 k2 = kp2 = _0_5

159 # b * K(kp2) = x

160 # b * K(k2) = y

161

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

163 b = _over(y, K)

164 if _xy:

165 k2, kp2 = kp2, k2

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

167

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

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

170

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

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

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

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

175

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

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

178 '''

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

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

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

182

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

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

185

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

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

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

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

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

191

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

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

194 '''

195 if not isinstance(other, Conformal3):

196 raise TriaxialError(other=other)

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

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

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

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

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

202 if M:

203 m *= _over(ma, mb)

204 else:

205 m = NAN

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

207

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

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

210 Sphere}.

211

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

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

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

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

216

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

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

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

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

221 '''

222 if M:

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

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

225 S = self._C3S

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

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

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

229 if M:

230 ma = self._invScale(bet, omg)

231 mb = _invScale(S, bets, omgs)

232 m = self._sphScale(ma, mb)

233 else:

234 m = NAN

235 return ct, dir3d, m

236

237 def _invScale(self, bet, omg):

238 # scale helper

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

240

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

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

243

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

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

246 ignored otherwise.

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

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

249

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

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

252 otherwise C{scale is NAN}.

253 '''

254 x, y = _cf2en(x_cf, y)

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

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

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

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

259

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

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

262

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

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

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

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

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

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

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

270

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

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

273 '''

274 if not isinstance(other, Conformal3):

275 raise TriaxialError(other=other)

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

277 m = _1_over(m) if M else NAN

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

279

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

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

282

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

284 Sphere} (C{Cartesian5Tuple}).

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

286 ignored otherwise.

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

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

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

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

291

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

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

294 '''

295 S = self._C3S

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

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

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

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

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

301 if M:

302 mb = _invScale(S, bets, omgs)

303 ma = _invScale(S, bet, omg)

304 m = self._sphScale(ma, mb)

305 else:

306 m = NAN

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

308

309 def _sphScale(self, ma, mb):

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

311 '''

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

313

314 @Property_RO

315 def _sphScale_m(self):

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

317 '''

318 e = sqrt(self.e2)

319

320 def _h1sg(atan_):

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

322 return hypot1(sg) + sg

323

324 k2, kp2 = self._k2_kp2

325 if not kp2: # oblate pole

326 m = self._c_b * _h1sg(atanh)

327 elif not k2: # prolate pole

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

329 else: # triaxial umbilical

330 S = self._C3S

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

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

333 return m

334

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

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

337 '''

338 omg = omg.shift(+1)

339 _, c, n = omg.scn3

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

341 x = NAN

342 else:

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

344 x *= self._a2_b

345 return x

346

347 @Property_RO

348 def xyQ2(self):

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

350 '''

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

352 self._c2_b * self._yE.cPi,

353 name=Conformal3.xyQ2.name)

354

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

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

357 '''

358 _, c, n = bet.scn3

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

360 y = NAN

361 else:

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

363 y *= self._c2_b

364 return y

365

366

367class Conformal3B(Conformal3):

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

369 specified by its middle semi-axis and shape.

370

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

372 '''

373 def __init__(self, b, e2=_0_0, k2=_1_0, kp2=_0_0, **name): # PYLINT W0231

374 '''New L{Conformal3B} triaxial.

375

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

377 I{spherical} projection.

378

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

380 '''

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

382

383

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

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

386

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

388 '''

389 def __init__(self, radius, k2=_1_0, kp2=_0_0, **name): # PYLINT W0231

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

391

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

393 '''

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

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

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

397

398

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

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

401 '''

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

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

404

405

406def _cf2en(x_cf, y):

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

408 '''

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

410 Easting(x_cf), Northing(y))

411

412

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

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

415 '''

416 s, c, n = phi.scn3

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

418 p = NAN

419 else:

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

421 if n:

422 p += eF.cK * n * _4_0

423 return p

424

425

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

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

428 '''

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

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

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

432 s, c = sincos2(phi)

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

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

435 return f, _1_over(fp)

436

437 z = fabs(y)

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

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

440 return r

441

442

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

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

445 '''

446 a2, ap2 = eF.alpha2, eF.alphap2

447 if _useTauf and not eF.kp2: # i.e. if eF.k2 == 1

448 # see GeographicLib 2.8 Conformal3.cpp

449 es = _copysign(sqrt(fabs(a2)), a2)

450 r = _tauf(sinh(ap2 * x), es)

451 r = Ang.fromRadians(r)

452 else:

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

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

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

456 s, c = sincos2(phi)

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

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

459 ( ap2 + s**2 * a2)

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

461 return f, _1_over(fp)

462

463 z = fabs(y)

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

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

466 return r

467

468

469def _invRy2(eF, eF_c, x):

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

471 if eF.kp2:

472 y, n = _remainder2(x, eF_c)

473 else: # eF_c == N-/INF

474 y, n = x, _0_0

475 if not y:

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

477 elif fabs(y) == eF_c:

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

479 return Ang.cardinal(n), y

480

481

482def _invScale(triax, bet, omg):

483 # helper for triaxial, sphere scale

484 k2, kp2 = triax._k2_kp2

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

486

487

488def _Neg(ang):

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

490

491

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

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

494 '''

495 s, c, n = phi.scn3

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

497 p = NAN

498 else:

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

500 if n:

501 p += eF.cPi * n * _4_0

502 return p

503

504

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

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

507 '''

508 k = b = C = 0

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

510 # p = PI_2 * 0 #???

511 # def _fp2z(x):

512 # f, fp = _fp2(x)

513 # f -= z

514 # # "DAT ", x, f, fp

515 # return f

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

517 # if (a * b) > 0:

518 # raise TriaxalError('"DATBAD")

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

520 vtol = tol * zscale * _0_01

521 xtol = pow(tol, _0_75) * xscale

522 oldv = oldx = olddx = INF

523 for k in range(1, _MAXIT):

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

525 v, vp = _fp2(x)

526 v -= z

527 va = fabs(v)

528 dx = _over(-v, vp)

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

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

531 C = 1 # k, va

532 break

533 elif (v * s) > 0:

534 xb = min(xb, x)

535 else:

536 xa = max(xa, x)

537 x += dx

538 dxa = fabs(dx)

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

540 (k > 2 and dxa > olddx):

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

542 x = (xa + xb) * _0_5

543 if x == oldx:

544 C = 3 # k, x, dx

545 break

546 elif not dxa > xtol:

547 C = 2 # k, dx, xtol

548 break

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

550 oldx = x

551 oldv = va

552 olddx = dxa * _0_5

553 else:

554 t = Fmt.no_convergence(dx, xtol)

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

556 else:

557 x = NAN

558 return x, k, b, C

559

560

561if __name__ == _DMAIN_:

562

563 from pygeodesy import Degrees, printf

564 from pygeodesy.triaxials import Triaxials

565

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

567 T = Conformal3(Triaxials.WGS84_3)

568 printf(T)

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

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

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

572 # (-5077726.431878843, 3922574.862034525, 1.1970322930902468)

573 # -5077732.396 3922571.859 1.1970343759 C++

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

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

576 # (33.300000000000004, 44.39999999999999, 1.1970322930902468)

577 # 33.30000000 44.40000000 1.1970343759 C++

578

579 T = Conformal3(Triaxials.WGS84_3r) # rounded

580 printf(T)

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

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

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

584 # (-5077732.396020001, 3922571.859124643, 1.197034375926918)

585 # -5077732.396 3922571.859 1.1970343759 C++

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

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

588 # (33.3, 44.400000000000006, 1.197034375926918)

589 # 33.30000000 44.40000000 1.1970343759 C++

590

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

592 T = Conformal3(6378172, 6378102, c)

593 printf(T)

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

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

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

597 # (-5077732.418682504, 3922572.0186951878, 1.197034384522207)

598 # -5077732.396 3922571.859 1.1970343759 C++

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

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

601 # (33.300000000000004, 44.4, 1.197034384522207)

602 # 33.30000000 44.40000000 1.1970343759 C++

603

604# % python3 -m pygeodesy.triaxials.conformal3

605

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

607# (-5077726.431878843, 3922574.862034525, 1.1970322930902468)

608# (33.300000000000004, 44.39999999999999, 1.1970322930902468)

609

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

611# (-5077732.396020001, 3922571.859124643, 1.197034375926918)

612# (33.3, 44.400000000000006, 1.197034375926918)

613

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

615# (-5077732.418682504, 3922572.0186951878, 1.197034384522207)

616# (33.300000000000004, 44.4, 1.197034384522207)

617

618# **) MIT License

619#

621#

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

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

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

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

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

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

628#

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

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

631#

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

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

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

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

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

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

638# OTHER DEALINGS IN THE SOFTWARE.