Coverage for pygeodesy/triaxials/triaxial5.py: 90%

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

4u'''Triaxal ellipsoid classes L{Triaxial} and I{unordered} L{Triaxial_} and Jacobi conformal projections

5L{Conformal} and L{ConformalSphere}, transcoded from I{Karney}'s GeographicLib 2.5.2 C++ class U{JacobiConformal

6<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1JacobiConformal.html#details>} to pure

7Python and miscellaneous classes L{BetaOmega2Tuple}, L{BetaOmega3Tuple} and L{Conformal2Tuple}, I{all kept

8for backward copability}.

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

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

12I{experimental} documentation.

14@see: U{Geodesics on a triaxial ellipsoid<https://WikiPedia.org/wiki/Geodesics_on_an_ellipsoid#

15 Geodesics_on_a_triaxial_ellipsoid>} and U{Triaxial coordinate systems and their geometrical

16 interpretation<https://OLD.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.

18@var Triaxials.Amalthea: Triaxial(name='Amalthea', a=125000, b=73000, c=64000, e2ab=0.658944, e2bc=0.231375493, e2ac=0.737856, volume=2446253479595252, area=93239507787.490341187, R2=86138.05359954)

19@var Triaxials.Ariel: Triaxial(name='Ariel', a=581100, b=577900, c=577700, e2ab=0.01098327, e2bc=0.000692042, e2ac=0.011667711, volume=812633172614203904, area=4211301462766.580078125, R2=578899.578791275)

20@var Triaxials.Earth: Triaxial(name='Earth', a=6378173.435, b=6378103.9, c=6356754.399999999, e2ab=0.000021804, e2bc=0.006683418, e2ac=0.006705077, volume=1083208241574987694080, area=510065911057440.9375, R2=6371008.987886564)

21@var Triaxials.Enceladus: Triaxial(name='Enceladus', a=256600, b=251400, c=248300, e2ab=0.040119337, e2bc=0.024509841, e2ac=0.06364586, volume=67094551514082248, area=798618496278.596679688, R2=252095.300756832)

22@var Triaxials.Europa: Triaxial(name='Europa', a=1564130, b=1561230, c=1560930, e2ab=0.003704694, e2bc=0.000384275, e2ac=0.004087546, volume=15966575194402123776, area=30663773697323.51953125, R2=1562096.533153486)

23@var Triaxials.Io: Triaxial(name='Io', a=1829400, b=1819300, c=1815700, e2ab=0.011011391, e2bc=0.003953651, e2ac=0.014921506, volume=25313121117889765376, area=41691875849096.734375, R2=1821464.812747882)

24@var Triaxials.Mars: Triaxial(name='Mars', a=3394600, b=3393300, c=3376300, e2ab=0.000765776, e2bc=0.009994646, e2ac=0.010752768, volume=162907283585817247744, area=144249140795107.4375, R2=3388064.624110653)

25@var Triaxials.Mimas: Triaxial(name='Mimas', a=207400, b=196800, c=190600, e2ab=0.09960581, e2bc=0.062015624, e2ac=0.155444317, volume=32587072869017956, area=493855762247.691772461, R2=198241.75359411)

26@var Triaxials.Miranda: Triaxial(name='Miranda', a=240400, b=234200, c=232900, e2ab=0.050915557, e2bc=0.011070811, e2ac=0.061422691, volume=54926187094835456, area=698880863325.757080078, R2=235828.692095158)

27@var Triaxials.Moon: Triaxial(name='Moon', a=1735550, b=1735324, c=1734898, e2ab=0.000260419, e2bc=0.000490914, e2ac=0.000751206, volume=21886698675223740416, area=37838824729886.09375, R2=1735257.329122863)

28@var Triaxials.Tethys: Triaxial(name='Tethys', a=535600, b=528200, c=525800, e2ab=0.027441672, e2bc=0.009066821, e2ac=0.036259685, volume=623086233855821440, area=3528073490771.393554688, R2=529863.348254881)

29@var Triaxials.WGS84_3: Triaxial(name='WGS84_3', a=6378171.36, b=6378101.609999999, c=6356751.84, e2ab=0.000021871, e2bc=0.006683505, e2ac=0.00670523, volume=1083207064030173855744, area=510065541435967.5, R2=6371006.679496506)

30@var Triaxials.WGS84_35: Triaxial(name='WGS84_35', a=6378172, b=6378102, c=6356752.314245179, e2ab=0.00002195, e2bc=0.006683478, e2ac=0.006705281, volume=1083207319768789942272, area=510065621722018.25, R2=6371007.180905545)

31@var Triaxials.WGS84_3r: Triaxial(name='WGS84_3r', a=6378172, b=6378102, c=6356752, e2ab=0.00002195, e2bc=0.006683577, e2ac=0.00670538, volume=1083207266220584468480, area=510065604942135.875, R2=6371007.076110449)

32'''

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

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

36from pygeodesy.angles import _SinCos2, Property_RO

37from pygeodesy.basics import _isin, isLatLon

38from pygeodesy.constants import EPS, EPS0, EPS02, _EPS2e4, INT0, \

39 _isfinite, isnear1, _over, _SQRT2_2, \

40 _0_0, _0_5, _1_0, _N_1_0

41from pygeodesy.datums import Datum, _spherical_datum, _WGS84, Fmt

42# from pygeodesy.elliptic import Elliptic # _MODS

43from pygeodesy.errors import _AssertionError, _ValueError, _xkwds_pop2

44from pygeodesy.fmath import Fdot, fdot, hypot, hypot_, fabs, sqrt

45from pygeodesy.fsums import fsumf_, fsum1f_

46from pygeodesy.interns import NN, _beta_, _distant_, _DMAIN_, _finite_, _height_, \

47 _inside_, _near_, _negative_, _not_, _null_, _opposite_, \

48 _outside_, _too_, _x_, _y_

49from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _FOR_DOCS

50from pygeodesy.named import _name__, _Pass

51from pygeodesy.namedTuples import LatLon3Tuple, _NamedTupleTo, Vector2Tuple, \

52 Vector3Tuple, Vector4Tuple

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

54# from pygeodesy.streprs import Fmt # from .datums

55from pygeodesy.triaxials.bases import Conformal5Tuple, _HeightINT0, _hypot2_1, \

56 _not_ordered_, _OrderedTriaxialBase, _over0, \

57 _otherV3d_, _over02, _sqrt0, TriaxialError, \

58 _Triaxial3Base, _TriaxialsBase, _UnOrderedTriaxialBase

59from pygeodesy.units import Degrees, Height_, Lat, Lon, Meter, Radians, Radius_, Scalar_

60from pygeodesy.utily import atan2, atan2d

61from pygeodesy.vector3d import _otherV3d, Vector3d

63# from math import fabs, sqrt # from .fmath

65__all__ = _ALL_LAZY.triaxials_triaxial5

66__version__ = '26.02.15'

68_omega_ = 'omega'

69_TRIPS = 359 # Eberly 1074?

72class _NamedTupleToX(_NamedTupleTo): # in .testNamedTuples

73 '''(INTERNAL) Base class for L{BetaOmega2Tuple},

74 L{BetaOmega3Tuple} and L{Conformal2Tuple}.

75 '''

76 def _toDegrees(self, name, **toDMS_kwds):

77 '''(INTERNAL) Convert C{self[0:2]} to L{Degrees} or C{toDMS}.

78 '''

79 return self._toX3U(_NamedTupleTo._Degrees3, Degrees, name, *self, **toDMS_kwds)

81 def _toRadians(self, name):

82 '''(INTERNAL) Convert C{self[0:2]} to L{Radians}.

83 '''

84 return self._toX3U(_NamedTupleTo._Radians3, Radians, name, *self)

86 def _toX3U(self, _X3, U, name, a, b, *c, **kwds):

87 a, b, s = _X3(self, a, b, **kwds)

88 if s is None or name:

89 n = self._name__(name)

90 s = self.classof(a, b, *c, name=n).reUnit(U, U).toUnits()

91 return s

94class BetaOmega2Tuple(_NamedTupleToX):

95 '''2-Tuple C{(beta, omega)} with I{ellipsoidal} lat- and

96 longitude C{beta} and C{omega} both in L{Radians} (or

97 L{Degrees}).

98 '''

99 _Names_ = (_beta_, _omega_)

100 _Units_ = (_Pass, _Pass)

101

102 def toDegrees(self, name=NN, **toDMS_kwds):

103 '''Convert this L{BetaOmega2Tuple} to L{Degrees} or C{toDMS}.

104

105 @kwarg name: Optional C{B{name}=NN} (C{str}).

106

107 @return: L{BetaOmega2Tuple}C{(beta, omega)} with C{beta} and

108 C{omega} both in L{Degrees} or as L{toDMS} strings

109 provided some B{C{toDMS_kwds}} keyword arguments are

110 specified.

111 '''

112 return self._toDegrees(name, **toDMS_kwds)

113

114 def toRadians(self, **name):

115 '''Convert this L{BetaOmega2Tuple} to L{Radians}.

116

117 @kwarg name: Optional C{B{name}=NN} (C{str}).

118

119 @return: L{BetaOmega2Tuple}C{(beta, omega)} with C{beta} and C{omega}

120 both in L{Radians}.

121 '''

122 return self._toRadians(name)

123

124

125class BetaOmega3Tuple(_NamedTupleToX):

126 '''3-Tuple C{(beta, omega, height)} with I{ellipsoidal} lat- and

127 longitude C{beta} and C{omega} both in L{Radians} (or L{Degrees})

128 and the C{height}, rather the (signed) I{distance} to the triaxial's

129 surface (measured along the radial line to the triaxial's center)

130 in C{meter}, conventionally.

131 '''

132 _Names_ = BetaOmega2Tuple._Names_ + (_height_,)

133 _Units_ = BetaOmega2Tuple._Units_ + ( Meter,)

134

135 def toDegrees(self, name=NN, **toDMS_kwds):

136 '''Convert this L{BetaOmega3Tuple} to L{Degrees} or C{toDMS}.

137

138 @kwarg name: Optional C{B{name}=NN} (C{str}).

139

140 @return: L{BetaOmega3Tuple}C{(beta, omega, height)} with

141 C{beta} and C{omega} both in L{Degrees} or as

142 L{toDMS} strings provided some B{C{toDMS_kwds}}

143 keyword arguments are specified.

144 '''

145 return self._toDegrees(name, **toDMS_kwds)

146

147 def toRadians(self, **name):

148 '''Convert this L{BetaOmega3Tuple} to L{Radians}.

149

150 @kwarg name: Optional C{B{name}=NN} (C{str}), overriding this name.

151

152 @return: L{BetaOmega3Tuple}C{(beta, omega, height)} with C{beta}

153 and C{omega} both in L{Radians}.

154 '''

155 return self._toRadians(name)

156

157 def to2Tuple(self, **name):

158 '''Reduce this L{BetaOmega3Tuple} to a L{BetaOmega2Tuple}.

159

160 @kwarg name: Optional C{B{name}=NN} (C{str}), overriding this name.

161 '''

162 return BetaOmega2Tuple(*self[:2], name=self._name__(name))

163

164

165class Conformal2Tuple(_NamedTupleToX):

166 '''2-Tuple C{(x, y)} with a I{Jacobi Conformal} C{x} and C{y}

167 projection, both in L{Radians} (or L{Degrees}).

168 '''

169 _Names_ = (_x_, _y_)

170 _Units_ = (_Pass, _Pass)

171

172 def toDegrees(self, name=NN, **toDMS_kwds):

173 '''Convert this L{Conformal2Tuple} to L{Degrees} or C{toDMS}.

174

175 @kwarg name: Optional C{B{name}=NN} (C{str}).

176

177 @return: L{Conformal2Tuple}C{(x, y)} with C{x} and C{y} both

178 in L{Degrees} or as L{toDMS} strings provided some

179 B{C{toDMS_kwds}} keyword arguments are specified.

180 '''

181 return self._toDegrees(name, **toDMS_kwds)

182

183 def toRadians(self, **name):

184 '''Convert this L{Conformal2Tuple} to L{Radians}.

185

186 @kwarg name: Optional C{B{name}=NN} (C{str}).

187

188 @return: L{Conformal2Tuple}C{(x, y)} with C{x} and C{y} both in L{Radians}.

189 '''

190 return self._toRadians(name)

191

192 def to5Tuple(self, b_conformal, **z_scale_name):

193 '''Return this L{Conformal2Tuple} as a L{Conformal5Tuple}.

194

195 @arg b_conformal: Middle semi-axis (C{meter}, conventionally)

196 or the original L{Conformal} of this 2-tuple.

197 @kwarg z_scale_name: Optional C{B{z}=0} meter, C{B{scale}=NAN}

198 and C{B{name}=NN} (C{str}).

199

200 @return: A L{Conformal5Tuple}.

201 '''

202 if isinstance(b_conformal, Conformal):

203 b = b_conformal.b

204 else:

205 b = Radius_(b=b_conformal)

206 x, y = self.toRadians()

207 x, y = _over(x, b), _over(y, b)

208 return Conformal5Tuple(x, y, **z_scale_name)

209

210

211class Triaxial_(_UnOrderedTriaxialBase):

212 '''I{Unordered} triaxial ellipsoid.

213

214 Triaxial ellipsoids with right-handed semi-axes C{a}, C{b} and C{c}, oriented

215 such that the large principal ellipse C{ab} is the equator I{Z}=0, I{beta}=0,

216 while the small principal ellipse C{ac} is the prime meridian, plane I{Y}=0,

217 I{omega}=0.

218

219 The four umbilic points, C{abs}(I{omega}) = C{abs}(I{beta}) = C{PI/2}, lie on

220 the middle principal ellipse C{bc} in plane I{X}=0, I{omega}=C{PI/2}.

221

222 @note: I{Geodetic} C{lat}- and C{lon}gitudes are in C{degrees}, I{geodetic}

223 C{phi} and C{lam}bda are in C{radians}, but I{ellipsoidal} lat- and

224 longitude C{beta} and C{omega} are in L{Radians} by default (or in

225 L{Degrees} if converted).

226 '''

227 if _FOR_DOCS:

228 __init__ = _UnOrderedTriaxialBase.__init__

229

230

231class Triaxial(_OrderedTriaxialBase):

232 '''I{Ordered} triaxial ellipsoid.

233

234 @see: L{Triaxial_} for more information.

235 '''

236 if _FOR_DOCS:

237 __init__ = _OrderedTriaxialBase.__init__

238

239 def forwardBetaOmega(self, beta, omega, height=0, **unit_name):

240 '''Convert I{ellipsoidal} lat- C{beta}, longitude C{omega} and C{height}

241 to cartesian.

242

243 @arg beta: Ellipsoidal latitude (C{Ang} or B{C{unit}}).

244 @arg omega: Ellipsoidal longitude (C{Ang} or B{C{unit}}).

245 @kwarg height: Height above or below the triaxial's surface (C{meter},

246 same units as this triaxial's semi-axes.

247 @kwarg unit_name: Optional C{B{name}=NN} (C{str}) and scalar

248 C{B{unit}=}L{Radians} (or L{Degrees}).

249

250 @return: A L{Vector3Tuple}C{(x, y, z)}.

251

252 @see: Method L{Triaxial.reverseBetaOmega} and U{equations (23-25)<https://

253 OLD.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.

254 '''

255 unit, name = _xkwds_pop2(unit_name, unit=Radians)

256 if height:

257 z = self._Height(height) + self.c

258 if z > 0:

259 z2 = z**2

260 x = z * _sqrt0(_1_0 + self._a2c2 / z2)

261 y = z * _sqrt0(_1_0 + self._b2c2 / z2)

262 else:

263 x = y = z = _0_0

264 else:

265 x, y, z = self._abc3

266 if z: # and x and y:

267 sa, ca = _SinCos2(beta, unit)

268 sb, cb = _SinCos2(omega, unit)

269

270 r = self._a2b2_a2c2

271 x *= cb * (_sqrt0(ca**2 + sa**2 * r) if r else fabs(ca))

272 y *= ca * sb

273 z *= sa * (_sqrt0(_1_0 - cb**2 * r) if r else _1_0)

274 return Vector3Tuple(x, y, z, **name)

275

276 def forwardBetaOmega_(self, sbeta, cbeta, somega, comega, **name):

277 '''Convert I{ellipsoidal} lat- and longitude C{beta} and C{omega}

278 to cartesian coordinates I{on this triaxial's surface}.

279

280 @arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}).

281 @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}).

282 @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}).

283 @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}).

284 @kwarg name: Optional C{B{name}=NN} (C{str}).

285

286 @return: A L{Vector3Tuple}C{(x, y, z)} on the surface.

287

288 @raise TriaxialError: This triaxial is near-spherical.

289

290 @see: Method L{Triaxial.reverseBetaOmega}, U{Triaxial ellipsoid coordinate

291 system<https://WikiPedia.org/wiki/Geodesics_on_an_ellipsoid#

292 Triaxial_ellipsoid_coordinate_system>} and U{equations (23-25)<https://

293 OLD.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.

294 '''

295 t = self._radialTo3(sbeta, cbeta, somega, comega)

296 return Vector3Tuple(*t, **name)

297

298 def forwardCartesian(self, x_xyz, y=None, z=None, normal=True, **eps_name):

299 '''Project any cartesian to a cartesian I{on this triaxial's surface}.

300

301 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},

302 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).

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

304 ignored otherwise.

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

306 @kwarg normal: If C{True}, the projection is C{perpendicular} to the surface,

307 otherwise C{radial} to the center of this triaxial (C{bool}).

308 @kwarg eps_name: Root finder tolerance C{B{eps}=EPS} and optional

309 C{B{name}="height4"} (C{str}).

310

311 @return: A L{Vector4Tuple}C{(x, y, z, h)}.

312

313 @see: Method L{Triaxial.reverseCartesian} to reverse the projection and

314 function L{height4<triaxials.triaxial5.height4>} for more details.

315 '''

316 return self.height4(x_xyz, y, z, normal=normal, **eps_name)

317

318 def forwardLatLon(self, lat, lon, height=0, **unit_name):

319 '''Convert I{geodetic} lat-, longitude and height to cartesian.

320

321 @arg lat: Geodetic latitude (C{Ang} or B{C{unit}}).

322 @arg lon: Geodetic longitude (C{Ang} or B{C{unit}}).

323 @arg height: Height above the ellipsoid (C{meter}, same units

324 as this triaxial's semi-axes).

325 @kwarg unit_name: Optional scalar C{B{unit}=}L{Degrees} and

326 C{B{name}=NN} (C{str}).

327

328 @return: A L{Vector3Tuple}C{(x, y, z)}.

329

330 @see: Method L{Triaxial.reverseLatLon} and U{equations (9-11)<https://

331 OLD.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.

332 '''

333 unit, name = _xkwds_pop2(unit_name, unit=Degrees)

334 return self._forwardLatLon3(height, name, *(_SinCos2(lat, unit) +

335 _SinCos2(lon, unit)))

336

337 def forwardLatLon_(self, slat, clat, slon, clon, height=0, **name):

338 '''Convert I{geodetic} lat-, longitude and height to cartesian.

339

340 @arg slat: Geodetic latitude C{sin(lat)} (C{scalar}).

341 @arg clat: Geodetic latitude C{cos(lat)} (C{scalar}).

342 @arg slon: Geodetic longitude C{sin(lon)} (C{scalar}).

343 @arg clon: Geodetic longitude C{cos(lon)} (C{scalar}).

344 @arg height: Height above the ellipsoid (C{meter}, same units

345 as this triaxial's semi-axes).

346 @kwarg name: Optional C{B{name}=NN} (C{str}).

347

348 @return: A L{Vector3Tuple}C{(x, y, z)}.

349

350 @see: Method L{Triaxial.reverseLatLon} and U{equations (9-11)<https://

351 OLD.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.

352 '''

353 sa, ca = self._norm2(slat, clat)

354 sb, cb = self._norm2(slon, clon)

355 return self._forwardLatLon3(height, name, sa, ca, sb, cb)

356

357 def _forwardLatLon3(self, height, name, sa, ca, sb, cb): # name always **name

358 '''(INTERNAL) Helper for C{.forwardLatLon} and C{.forwardLatLon_}.

359 '''

360 h = self._Height(height)

361 x = ca * cb

362 y = ca * sb

363 z = sa

364 # 1 - (1 - (c/a)**2) * sa**2 - (1 - (b/a)**2) * ca**2 * sb**2

365 t = fsumf_(_1_0, -self.e2ac * z**2, -self.e2ab * y**2)

366 n = self.a / _sqrt0(t) # prime vertical

367 x *= h + n

368 y *= h + n * self._b2_a2

369 z *= h + n * self._c2_a2

370 return Vector3Tuple(x, y, z, **name)

371

372 def _Height(self, height):

373 '''(INTERNAL) Validate a C{height}.

374 '''

375 return Height_(height=height, low=-self.c, Error=TriaxialError)

376

377 def reverseBetaOmega(self, x_xyz, y=None, z=None, **name):

378 '''Convert cartesian to I{ellipsoidal} lat- and longitude, C{beta}, C{omega}

379 and height.

380

381 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},

382 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).

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

384 ignored otherwise.

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

386 @kwarg name: Optional C{B{name}=NN} (C{str}).

387

388 @return: A L{BetaOmega3Tuple}C{(beta, omega, height)} with C{beta} and

389 C{omega} in L{Radians} and (radial) C{height} in C{meter}, same

390 units as this triaxial's semi-axes.

391

392 @see: Methods L{Triaxial.forwardBetaOmega} and L{Triaxial.forwardBetaOmega_}

393 and U{equations (21-22)<https://OLD.Topo.Auth.GR/wp-content/uploads/

394 sites/111/2021/12/09_Panou.pdf>}.

395 '''

396 v = _otherV3d_(x_xyz, y, z)

397 a, b, h = _reverseLatLon3(v, atan2, v, self.forwardBetaOmega_)

398 return BetaOmega3Tuple(Radians(beta=a), Radians(omega=b), h, **name)

399

400 def reverseCartesian(self, x_xyz, y=None, z=None, height=0, normal=True, eps=_EPS2e4, **name):

401 '''"Unproject" a cartesian I{off} this triaxial's surface.

402

403 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},

404 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).

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

406 ignored otherwise.

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

408 @kwarg height: Height above or below this triaxial's surface (C{meter}, same

409 units as this triaxial's semi-axes).

410 @kwarg normal: If C{True}, B{C{height}} is C{perpendicular} to the surface,

411 otherwise C{radial} to the center of this triaxial (C{bool}).

412 @kwarg eps: Tolerance for on-surface test (C{scalar}), see method

413 L{Triaxial.sideOf}.

414 @kwarg name: Optional C{B{name}=NN} (C{str}).

415

416 @return: A L{Vector3Tuple}C{(x, y, z)}.

417

418 @raise TrialError: Cartesian B{C{x_xyz}} or C{(x, y, z)} not on this triaxial's

419 surface.

420

421 @see: Methods L{Triaxial.forwardCartesian} and L{Triaxial.height4}.

422 '''

423 h, name = _xkwds_pop2(name, h=height) # h=height for backward compatibility

424 v = _otherV3d_(x_xyz, y, z, **name)

425 _ = self._sideOn(v, eps=eps)

426 h = _HeightINT0(h)

427 if h:

428 if normal:

429 v = v.plus(self.normal3d(*v.xyz, length=h))

430 elif v.length > EPS0:

431 v = v.times(_1_0 + (h / v.length))

432 return v.xyz # Vector3Tuple

433

434 def reverseLatLon(self, x_xyz, y=None, z=None, **name):

435 '''Convert cartesian to I{geodetic} lat-, longitude and height.

436

437 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},

438 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).

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

440 ignored otherwise.

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

442 @kwarg name: Optional C{B{name}=NN} (C{str}).

443

444 @return: A L{LatLon3Tuple}C{(lat, lon, height)} with C{lat} and C{lon}

445 in C{degrees} and (radial) C{height} in C{meter}, same units

446 as this triaxial's semi-axes.

447

448 @see: Methods L{Triaxial.forwardLatLon} and L{Triaxial.forwardLatLon_}

449 and U{equations (4-5)<https://OLD.Topo.Auth.GR/wp-content/uploads/

450 sites/111/2021/12/09_Panou.pdf>}.

451 '''

452 v = _otherV3d_(x_xyz, y, z)

453 s = v.times_(self._c2_a2, # == 1 - e_sub_x**2

454 self._c2_b2, # == 1 - e_sub_y**2

455 _1_0)

456 a, b, h = _reverseLatLon3(s, atan2d, v, self.forwardLatLon_)

457 return LatLon3Tuple(Lat(a), Lon(b), h, **name)

458

459

460class TriaxialB(_Triaxial3Base):

461 '''I{Ordered} triaxial ellipsoid specified by its middle semi-axis and shape.

462

463 @see: L{Triaxial} for details and more information.

464 '''

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

466 '''New L{TriaxialB} triaxial.

467

468 @arg b: Middle semi-axis (C{meter}, conventionally).

469 @kwarg e2: Excentricty I{squared} (C{scalar, e2 = (a**2 - c**2) / B{b}**2}).

470 @kwarg k2: Oblateness I{squared} (C{scalar, k2 = (C{b}**2 - c**2) / (a**2 - c**2)}).

471 @kwarg kp2: Prolateness I{squared} (C{scalar, kp2 = (a**2 - C{b}**2) / (a**2 - c**2)}).

472 @kwarg name: Optional C{B{name}=NN} (C{str}).

473

474 @note: The semi-axes must be ordered as C{B{a} >= B{b} >= B{c} > 0} but

475 may be spherical, C{B{e2} == 0}, i.e. C{B{a} == B{c}}. In that case

476 C{B{k2} == 0} represents a I{prolate sphere}, C{B{k2} == 1} an

477 I{oblate sphere}, otherwise a I{triaxial sphere} with C{0 < B{k2} < 1}.

478

479 @note: The B{C{k2}} and B{C{kp2}} arguments are normalized, C{B{k2} + B{kp2} == 1}.

480

481 @raise TriaxialError: Semi-axes unordered or invalid.

482 '''

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

484

485

486class Conformal(Triaxial):

487 '''This is a I{Jacobi Conformal} projection of a triaxial ellipsoid to a plane where

488 the C{X} and C{Y} grid lines are straight.

489

490 I{Ellipsoidal} coordinates I{beta} and I{omega} are converted to Jacobi Conformal

491 I{y} respectively I{x} separately. Jacobi's coordinates have been multiplied

492 by C{sqrt(B{a}**2 - B{c}**2) / (2 * B{b})} so that the customary results are

493 returned in the case of an ellipsoid of revolution.

494

495 Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2014-2024) and

496 licensed under the MIT/X11 License.

497

498 @note: This constructor can I{not be used to specify a sphere}, see alternate

499 L{ConformalSphere}.

500

501 @see: L{Triaxial}, C++ class U{JacobiConformal<https://GeographicLib.SourceForge.io/

502 C++/doc/classGeographicLib_1_1JacobiConformal.html#details>}, U{Jacobi's conformal

503 projection<https://GeographicLib.SourceForge.io/C++/doc/jacobi.html>} and Jacobi,

504 C. G. J. I{U{Vorlesungen über Dynamik<https://Books.Google.com/books?

505 id=ryEOAAAAQAAJ&pg=PA212>}}, page 212ff.

506 '''

507 if _FOR_DOCS:

508 __init__ = Triaxial.__init__

509

510 @Property_RO

511 def _a2_b(self):

512 return self._a2_b2 * self.b

513

514 @Property_RO

515 def _c2_b(self):

516 return self._c2_b2 * self.b

517

518 def x(self, omega, unit=Radians):

519 '''Compute a I{Jacobi Conformal} C{x} projection.

520

521 @arg omega: Ellipsoidal longitude (C{Ang} or B{C{unit}}).

522 @kwarg unit: Unit of scalar B{C{omega}} (or C{Degrees}).

523

524 @return: The C{x} projection (L{Meter}), same units as

525 this triaxial's semi-axes.

526 '''

527 s, c = _SinCos2(omega, unit)

528 return Meter(x=self._x(s, c, self._a2_b))

529

530 def _x(self, s, c, a2_b_):

531 '''(INTERNAL) Helper for C{.x}, C{.xR_} and C{.xy}.

532 '''

533 s, c = self._norm2(s, c, self.a)

534 return self._xE.fPi(s, c) * a2_b_

535

536 @Property_RO

537 def _xE(self):

538 '''(INTERNAL) Get the x-elliptic function.

539 '''

540 k2, kp2 = self._k2_kp2E

541 # -a2b2 / b2 == (b2 - a2) / b2 == 1 - a2 / b2 == 1 - a2_b2

542 return _MODS.elliptic.Elliptic(k2, _1_0 - self._a2_b2, kp2, self._a2_b2)

543

544 def xR(self, omega, unit=Radians):

545 '''Compute a I{Jacobi Conformal} C{x} projection.

546

547 @arg omega: Ellipsoidal longitude (C{Ang} or B{C{unit}}).

548 @kwarg unit: Unit of scalar B{C{omega}} (or C{Degrees}).

549

550 @return: The C{x} projection (L{Radians}).

551 '''

552 return self.xR_(*_SinCos2(omega, unit))

553

554 def xR_(self, somega, comega):

555 '''Compute a I{Jacobi Conformal} C{x} projection.

556

557 @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}).

558 @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}).

559

560 @return: The C{x} projection (L{Radians}).

561 '''

562 return Radians(x=self._x(somega, comega, self._a2_b2))

563

564 def xy(self, beta, omega, **unit_name):

565 '''Compute a I{Jacobi Conformal} C{x} and C{y} projection.

566

567 @arg beta: Ellipsoidal latitude (C{Ang} or B{C{unit}}).

568 @arg omega: Ellipsoidal longitude (C{Ang} or B{C{unit}}).

569 @kwarg unit_name: Optional scalar C{B{unit}=}L{Radians} and

570 name (C{str}), overriding C{B{name}="xyR2"}.

571

572 @return: A (L{Vector2Tuple}C{(x, y)}), both in C{Meter},

573 same units as this triaxial's semi-axes..

574 '''

575 unit, name = _xkwds_pop2(unit_name, unit=Radians)

576 return Vector2Tuple(self.x(omega, unit=unit),

577 self.y(beta, unit=unit),

578 name=_name__(name, name__=self.xy))

579

580 @Property_RO

581 def xyQ2(self):

582 '''Get the I{Jacobi Conformal} quadrant size in C{meter}

583 (L{Vector2Tuple}C{(x, y)}).

584 '''

585 return Vector2Tuple(Meter(x=self._a2_b * self._xE.cPi),

586 Meter(y=self._c2_b * self._yE.cPi),

587 name=Conformal.xyQ2.name)

588

589 @Property_RO

590 def xyQR2(self):

591 '''Get the I{Jacobi Conformal} quadrant size in C{Radians}

592 (L{Conformal2Tuple}C{(x, y)}).

593 '''

594 return Conformal2Tuple(Radians(x=self._a2_b2 * self._xE.cPi),

595 Radians(y=self._c2_b2 * self._yE.cPi),

596 name=Conformal.xyQR2.name)

597

598 def xyR2(self, beta, omega, **unit_name):

599 '''Compute a I{Jacobi Conformal} C{x} and C{y} projection.

600

601 @arg beta: Ellipsoidal latitude (C{Ang} or B{C{unit}}).

602 @arg omega: Ellipsoidal longitude (C{Ang} or B{C{unit}}).

603 @kwarg unit_name: Optional scalar C{B{unit}=}L{Radians} and

604 name (C{str}), overriding C{B{name}="xyR2"}.

605

606 @return: A L{Conformal2Tuple}C{(x, y)}, both in C{Radians}.

607 '''

608 unit, name = _xkwds_pop2(unit_name, unit=Radians)

609 sb_cb_so_co = _SinCos2(beta, unit) + _SinCos2(omega, unit)

610 return self.xyR2_(*sb_cb_so_co, name=_name__(name, name__=self.xyR2))

611

612 def xyR2_(self, sbeta, cbeta, somega, comega, **name):

613 '''Compute a I{Jacobi Conformal} C{x} and C{y} projection.

614

615 @arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}).

616 @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}).

617 @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}).

618 @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}).

619 @kwarg name: Optional name (C{str}), overriding C{B{name}="xyR2_"}.

620

621 @return: A L{Conformal2Tuple}C{(x, y)}.

622 '''

623 return Conformal2Tuple(self.xR_(somega, comega),

624 self.yR_(sbeta, cbeta),

625 name=_name__(name, name__=self.xyR2_))

626

627 def y(self, beta, unit=Radians):

628 '''Compute a I{Jacobi Conformal} C{y} projection.

629

630 @arg beta: Ellipsoidal latitude (C{Ang} or B{C{unit}}).

631 @kwarg unit: Unit of scalar B{C{beta}} (or C{Degrees}).

632

633 @return: The C{y} projection (L{Meter}), same units as

634 this triaxial's semi-axes.

635 '''

636 s, c = _SinCos2(beta, unit)

637 return Meter(y=self._y(s, c, self._c2_b))

638

639 def _y(self, s, c, c2_b_):

640 '''(INTERNAL) Helper for C{.y}, C{.yR_} and C{.xy}.

641 '''

642 s, c = self._norm2(s, c, self.c)

643 return self._yE.fPi(s, c) * c2_b_

644

645 @Property_RO

646 def _yE(self):

647 '''(INTERNAL) Get the y-elliptic function.

648 '''

649 k2, kp2 = self._k2_kp2E

650 # b2c2 / b2 == (b2 - c2) / b2 == 1 - c2 / b2 == e2bc

651 return _MODS.elliptic.Elliptic(kp2, self.e2bc, k2, self._c2_b2)

652

653 def yR(self, beta, unit=Radians):

654 '''Compute a I{Jacobi Conformal} C{y} projection.

655

656 @arg beta: Ellipsoidal latitude (C{Ang} or B{C{unit}}).

657 @kwarg unit: Unit of scalar B{C{beta}} (or C{Degrees}).

658

659 @return: The C{y} projection (L{Radians}).

660 '''

661 return self.yR_(*_SinCos2(beta, unit))

662

663 def yR_(self, sbeta, cbeta):

664 '''Compute a I{Jacobi Conformal} C{y} projection.

665

666 @arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}).

667 @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}).

668

669 @return: The C{y} projection (L{Radians}).

670 '''

671 return Radians(y=self._y(sbeta, cbeta, self._c2_b2))

672

673

674class ConformalSphere(Conformal):

675 '''Alternate, I{Jacobi Conformal projection} on a I{spherical} triaxial.

676

677 @see: L{Conformal<triaxials.triaxial5.Conformal>} for more information.

678 '''

679 _ab = _bc = 0

680

681 def __init__(self, radius_conformal, ab=0, bc=0, **name):

682 '''New L{ConformalSphere}.

683

684 @arg radius_conformal: Radius (C{scalar}, conventionally in C{meter})

685 or an other L{ConformalSphere} or L{Conformal}.

686 @kwarg ab: Relative magnitude of C{B{a} - B{b}} (C{meter}, same units

687 as C{scalar B{radius}}.

688 @kwarg bc: Relative magnitude of C{B{b} - B{c}} (C{meter}, same units

689 as C{scalar B{radius}}.

690 @kwarg name: Optional C{B{name}=NN} (C{str}).

691

692 @raise TriaxialError: Invalid B{C{radius_conformal}}, negative B{C{ab}},

693 negative B{C{bc}} or C{(B{ab} + B{bc})} not positive.

694

695 @note: If B{C{radius_conformal}} is a L{ConformalSphere} and if B{C{ab}}

696 and B{C{bc}} are both zero or C{None}, the B{C{radius_conformal}}'s

697 C{ab}, C{bc}, C{a}, C{b} and C{c} are copied.

698 '''

699 try:

700 r = radius_conformal

701 if isinstance(r, Triaxial): # ordered only

702 t = r._abc3

703 j = isinstance(r, ConformalSphere) and not bool(ab or bc)

704 else:

705 t = (Radius_(radius=r),) * 3

706 j = False

707 self._ab = r.ab if j else Scalar_(ab=ab) # low=0

708 self._bc = r.bc if j else Scalar_(bc=bc) # low=0

709 if (self.ab + self.bc) <= 0:

710 raise ValueError(_negative_)

711 a, _, c = self._abc3 = t

712 if not (a >= c and _isfinite(self._a2b2)

713 and _isfinite(self._a2c2)):

714 raise ValueError(_not_(_finite_))

715 except (TypeError, ValueError) as x:

716 raise TriaxialError(radius=r, ab=ab, bc=bc, cause=x)

717 if name:

718 self.name = name

719

720 @Property_RO

721 def ab(self):

722 '''Get relative magnitude C{a - b} (C{meter}, same units as B{C{a}}).

723 '''

724 return self._ab

725

726 @Property_RO

727 def _a2b2(self):

728 '''(INTERNAL) Get C{a**2 - b**2} == ab * (a + b).

729 '''

730 a, b, _ = self._abc3

731 return self.ab * (a + b)

732

733 @Property_RO

734 def _a2c2(self):

735 '''(INTERNAL) Get C{a**2 - c**2} == a2b2 + b2c2.

736 '''

737 return self._a2b2 + self._b2c2

738

739 @Property_RO

740 def bc(self):

741 '''Get relative magnitude C{b - c} (C{meter}, same units as B{C{a}}).

742 '''

743 return self._bc

744

745 @Property_RO

746 def _b2c2(self):

747 '''(INTERNAL) Get C{b**2 - c**2} == bc * (b + c).

748 '''

749 _, b, c = self._abc3

750 return self.bc * (b + c)

751

752 @Property_RO

753 def radius(self):

754 '''Get radius (C{meter}, conventionally).

755 '''

756 return self.a

757

758

759def _hartzell3(pov, los, Tun): # in .Ellipsoid.hartzell4, .formy.hartzell

760 '''(INTERNAL) Hartzell's "Satellite Line-of-Sight Intersection ...",

761 formula from a Point-Of-View to an I{un-/ordered} Triaxial.

762 '''

763 def _toUvwV3d(los, pov):

764 try: # pov must be LatLon or Cartesian if los is a Los

765 v = los.toUvw(pov)

766 except (AttributeError, TypeError):

767 v = _otherV3d(los=los)

768 return v

769

770 p3 = _otherV3d(pov=pov.toCartesian() if isLatLon(pov) else pov)

771 if los is True: # normal

772 a, b, c, d, i = _plumbTo5(p3.x, p3.y, p3.z, Tun)

773 return type(p3)(a, b, c), d, i

774

775 u3 = p3.negate() if los is False or los is None else _toUvwV3d(los, pov)

776

777 a, b, c, T = Tun._ordered4

778

779 a2 = T.a2 # largest, factored out

780 b2, p2 = (T.b2, T._b2_a2) if b != a else (a2, _1_0)

781 c2, q2 = (T.c2, T._c2_a2) if c != a else (a2, _1_0)

782

783 p3 = T._order3d(p3)

784 u3 = T._order3d(u3).unit() # unit vector, opposing signs

785

786 x2, y2, z2 = p3.x2y2z23 # p3.times_(p3).xyz3

787 ux, vy, wz = u3.times_(p3).xyz3

788 u2, v2, w2 = u3.x2y2z23 # u3.times_(u3).xyz3

789

790 t = (p2 * c2), c2, b2

791 m = fdot(t, u2, v2, w2) # a2 factored out

792 if m < EPS0: # zero or near-null LOS vector

793 raise _ValueError(_near_(_null_))

794

795 r = fsumf_(b2 * w2, c2 * v2, -v2 * z2, vy * wz * 2,

796 -w2 * y2, -u2 * y2 * q2, -u2 * z2 * p2, ux * wz * 2 * p2,

797 -w2 * x2 * p2, b2 * u2 * q2, -v2 * x2 * q2, ux * vy * 2 * q2)

798 if r > 0: # a2 factored out

799 r = sqrt(r) * b * c # == a * a * b * c / a2

800 elif r < 0: # LOS pointing away from or missing the triaxial

801 raise _ValueError(_opposite_ if max(ux, vy, wz) > 0 else _outside_)

802

803 d = Fdot(t, ux, vy, wz).fadd_(r).fover(m) # -r for antipode, a2 factored out

804 if d > 0: # POV inside or LOS outside or missing the triaxial

805 s = fsumf_(_N_1_0, _over(x2, a2), _over(y2, b2), _over(z2, c2)) # like _sideOf

806 raise _ValueError(_outside_ if s > 0 else _inside_)

807 elif fsum1f_(x2, y2, z2, -d*d) < 0: # d past triaxial's center

808 raise _ValueError(_too_(_distant_))

809

810 v = p3.minus(u3.times(d)) # cartesian type(pov) or Vector3d

811 h = p3.minus(v).length # distance to pov == -d

812 return T._order3d(v, reverse=True), h, None

813

814

815def hartzell4(pov, los=False, tri_biax=_WGS84, **name):

816 '''Compute the intersection of a tri-/biaxial ellipsoid and a Line-Of-Sight from

817 a Point-Of-View outside.

818

819 @arg pov: Point-Of-View outside the tri-/biaxial (C{Cartesian}, L{Ecef9Tuple},

820 C{LatLon} or L{Vector3d}).

821 @kwarg los: Line-Of-Sight, I{direction} to the tri-/biaxial (L{Los}, L{Vector3d})

822 or C{True} for the I{normal, perpendicular, plumb} to the surface of

823 the tri-/biaxial or C{False} or C{None} to point to its center.

824 @kwarg tri_biax: A triaxial (L{Triaxial}, L{Triaxial_}, L{Conformal} or

825 L{ConformalSphere}) or biaxial ellipsoid (L{Datum},

826 L{Ellipsoid}, L{Ellipsoid2}, L{a_f2Tuple}) or spherical earth

827 radius (C{scalar}, conventionally in C{meter}).

828 @kwarg name: Optional name (C{str}), overriding C{B{name}="hartzell4"}.

829

830 @return: L{Vector4Tuple}C{(x, y, z, h)} on the tri-/biaxial's surface, with C{h}

831 the distance from B{C{pov}} to C{(x, y, z)} I{along the} B{C{los}}, all

832 in C{meter}, conventionally.

833

834 @raise TriaxialError: Invalid B{C{pov}} or B{C{pov}} inside the tri-/biaxial or

835 invalid B{C{los}} or B{C{los}} points outside or away from

836 the tri-/biaxial.

837

838 @raise TypeError: Invalid B{C{tri_biax}}, C{ellipsoid} or C{datum}.

839

840 @see: Class L{pygeodesy3.Los}, functions L{pygeodesy.tyr3d} and L{pygeodesy.hartzell}

841 and U{lookAtSpheroid<https://PyPI.org/project/pymap3d>} and U{"Satellite

842 Line-of-Sight Intersection with Earth"<https://StephenHartzell.Medium.com/

843 satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6>}.

844 '''

845 T = _tri_biaxial(tri_biax, hartzell4)

846 try:

847 v, h, i = _hartzell3(pov, los, T)

848 except Exception as x:

849 raise TriaxialError(pov=pov, los=los, tri_biax=tri_biax, cause=x)

850 n = _name__(name, name__=hartzell4) # typename

851 return Vector4Tuple(v.x, v.y, v.z, h, iteration=i, name=n)

852

853

854def height4(x_xyz, y=None, z=None, tri_biax=_WGS84, normal=True, eps=EPS, **name):

855 '''Compute the projection on and the height above or below a tri-/biaxial ellipsoid's surface.

856

857 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian}, L{Ecef9Tuple},

858 L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).

859 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}, ignored

860 otherwise.

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

862 @kwarg normal: If C{True}, the projection is the I{perpendicular, plumb} to the

863 tri-/biaxial's surface, otherwise the C{radial} line to the center

864 of the tri-/biaxial (C{bool}).

865 @kwarg eps: Tolerance for root finding and validation (C{scalar}), use a negative

866 value to skip validation.

867 @kwarg name: Optional C{B{name}="height4"} (C{str}).

868

869 @return: L{Vector4Tuple}C{(x, y, z, h)} with the cartesian coordinates C{x}, C{y}

870 and C{z} of the projection on or the intersection with and with the height

871 C{h} above or below the tri-/biaxial's surface in C{meter}, conventionally.

872

873 @raise TriaxialError: Non-cartesian B{C{xyz}}, invalid B{C{eps}}, no convergence in

874 root finding or validation failed.

875

876 @see: Methods L{Triaxial.normal3d} and L{Ellipsoid.height4}, I{Eberly}'s U{Distance from a Point to ...

877 <https://www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>} and I{Bektas}'

878 U{Shortest Distance from a Point to Triaxial Ellipsoid<https://www.ResearchGate.net/publication/

879 272149005_SHORTEST_DISTANCE_FROM_A_POINT_TO_TRIAXIAL_ELLIPSOID>}.

880 '''

881 v = _otherV3d_(x_xyz, y, z)

882 T = _tri_biaxial(tri_biax, height4)

883 r = T.isSpherical

884

885 i, h = None, v.length

886 if h < EPS0: # EPS

887 x = y = z = _0_0

888 h -= min(T._abc3) # nearest

889 elif r: # .isSpherical

890 x, y, z = v.times(r / h).xyz3

891 h -= r

892 else:

893 x, y, z = v.xyz3

894 try:

895 if normal: # plumb to surface

896 x, y, z, h, i = _plumbTo5(x, y, z, T, eps=eps)

897 else: # radial to center

898 x, y, z = T._radialTo3(z, hypot(x, y), y, x)

899 h = v.minus_(x, y, z).length

900 except Exception as e:

901 raise TriaxialError(x=x, y=y, z=z, cause=e)

902 if h > 0 and T.sideOf(v, eps=EPS0) < 0:

903 h = -h # inside

904 n = _name__(name, name__=height4) # typename

905 return Vector4Tuple(x, y, z, h, iteration=i, name=n)

906

907

908def _ordered(a, b, c):

909 '''(INTERNAL) Is C{a >= b >= c > 0}?

910 '''

911 if not (_isfinite(a) and a >= b >= c > 0):

912 raise TriaxialError(a=a, b=b, c=c, txt=_not_ordered_)

913

914

915def _plumbTo3(px, py, E, eps=EPS): # in .ellipsoids.Ellipsoid.height4

916 '''(INTERNAL) Nearest point in 1st quadrant on a 2-D ellipse.

917 '''

918 a, b = E.a, E.b

919 if min(px, py, a, b) < EPS0:

920 raise _AssertionError(px=px, py=py, a=a, b=b, E=E)

921

922 a2 = a - b * E.b_a

923 b2 = b - a * E.a_b

924 tx = ty = _SQRT2_2

925 for i in range(16): # max 5

926 ex = tx**3 * a2

927 ey = ty**3 * b2

928

929 qx = px - ex

930 qy = py - ey

931 q = hypot(qx, qy)

932 if q < EPS0: # PYCHOK no cover

933 break

934 r = hypot(ex - tx * a,

935 ey - ty * b) / q

936

937 sx, tx = tx, min(_1_0, max(0, (ex + qx * r) / a))

938 sy, ty = ty, min(_1_0, max(0, (ey + qy * r) / b))

939 t = hypot(ty, tx)

940 if t < EPS0: # PYCHOK no cover

941 break

942 tx = tx / t # /= chokes PyChecker

943 ty = ty / t

944 if fabs(sx - tx) < eps and \

945 fabs(sy - ty) < eps:

946 break

947

948 tx *= a / px

949 ty *= b / py

950 return tx, ty, i # x and y as fractions

951

952

953def _plumbTo4(x, y, a, b, eps=EPS):

954 '''(INTERNAL) Nearest point on and distance to a 2-D ellipse, I{unordered}.

955

956 @see: Function C{_plumbTo3} and I{Eberly}'s U{Distance from a Point to ... an Ellipse ...

957 <https://www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>}.

958 '''

959 if b > a:

960 b, a, d, i = _plumbTo4(y, x, b, a, eps=eps)

961 return a, b, d, i

962

963 if not (b > 0 and _isfinite(a)):

964 raise _ValueError(a=a, b=b)

965

966 i = None

967 if y:

968 if x:

969 u = fabs(x / a)

970 w = fabs(y / b)

971 g = _hypot2_1(u, w)

972 if fabs(g) > EPS02:

973 r = (b / a)**2

974 t, i = _rootNd(_1_0 / r, 0, u, 0, w, g) # eps

975 a = _over(x, t * r + _1_0)

976 b = _over(y, t + _1_0)

977 d = hypot(x - a, y - b)

978 else: # on the ellipse

979 a, b, d = x, y, _0_0

980 else: # x == 0

981 if y < 0:

982 b = -b

983 a = x # _copysign_0_0

984 d = fabs(y - b)

985

986 elif x: # y == 0

987 d, r = None, _over0(a * x, (a + b) * (a - b))

988 if r:

989 a *= r

990 r = _1_0 - r**2

991 if r > EPS02:

992 b *= sqrt(r)

993 d = hypot(x - a, y - b)

994 elif x < 0:

995 a = -a

996 if d is None:

997 b = y # _copysign_0_0

998 d = fabs(x - a)

999

1000 else: # x == y == 0

1001 a = x # _copysign_0_0

1002 d = b

1003

1004 return a, b, d, i

1005

1006

1007def _plumbTo5(x, y, z, Tun, eps=EPS): # in .testTriaxials

1008 '''(INTERNAL) Nearest point on and distance to an I{un-/ordered} triaxial.

1009

1010 @see: I{Eberly}'s U{Distance from a Point to ... an Ellipsoid ...<https://

1011 www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>}.

1012 '''

1013 a, b, c, T = Tun._ordered4

1014 if Tun is not T: # T is ordered, Tun isn't

1015 t = T._order3(x, y, z) + (T,)

1016 a, b, c, d, i = _plumbTo5(*t, eps=eps)

1017 return T._order3(a, b, c, reverse=True) + (d, i)

1018

1019 if not (c > 0 and _isfinite(a)):

1020 raise _ValueError(a=a, b=b, c=c)

1021

1022 if eps > 0:

1023 val = max(eps * 1e8, EPS)

1024 else: # no validation

1025 val, eps = 0, max(EPS0, -eps)

1026

1027 i = None

1028 if z:

1029 if y:

1030 if x:

1031 u = fabs(x / a)

1032 v = fabs(y / b)

1033 w = fabs(z / c)

1034 g = _hypot2_1(u, v, w)

1035 if fabs(g) > EPS02:

1036 r = T._c2_a2 # (c / a)**2

1037 s = T._c2_b2 # (c / b)**2

1038 t, i = _rootNd(_1_0 / r, _1_0 / s, u, v, w, g) # eps

1039 a = _over(x, t * r + _1_0)

1040 b = _over(y, t * s + _1_0)

1041 c = _over(z, t + _1_0)

1042 d = hypot_(x - a, y - b, z - c)

1043 else: # on the ellipsoid

1044 a, b, c, d = x, y, z, _0_0

1045 else: # x == 0

1046 a = x # = _copysign_0_0(x)

1047 b, c, d, i = _plumbTo4(y, z, b, c, eps=eps)

1048 elif x: # y == 0

1049 b = y # = _copysign_0_0(y)

1050 a, c, d, i = _plumbTo4(x, z, a, c, eps=eps)

1051 else: # x == y == 0

1052 if z < 0:

1053 c = -c

1054 a, b, d = x, y, fabs(z - c)

1055

1056 else: # z == 0

1057 u = _over0(a * x, T._a2c2) # (a + c) * (a - c)

1058 v = _over0(b * y, T._b2c2) # (b + c) * (b - c)

1059 s = _hypot2_1(u, v)

1060 if u and v and s < 0:

1061 a *= u

1062 b *= v

1063 c *= sqrt(-s)

1064 d = hypot_(x - a, y - b, c)

1065 else:

1066 c = z # _copysign_0_0(z)

1067 a, b, d, i = _plumbTo4(x, y, a, b, eps=eps)

1068

1069 if val > 0:

1070 _validate(a, b, c, d, T, x, y, z, val)

1071 return a, b, c, d, i

1072

1073

1074def _reverseLatLon3(s, atan2_, v, forward_):

1075 '''(INTERNAL) Helper for C{.reverseBetaOmega} and C{.reverseLatLon}.

1076 '''

1077 x, y, z = s.xyz3

1078 d = hypot( x, y)

1079 h = v.minus_(*forward_(z, d, y, x)).length

1080 return atan2_(z, d), atan2_(y, x), (h or INT0)

1081

1082

1083def _rootNd(r, s, u, v, w, g, eps=EPS0):

1084 '''(INTERNAL) Robust 2-D or 3-D root finder: 2-D if C{s == v == 0} else 3-D root.

1085

1086 @see: I{Eberly}'s U{Robust Root Finders ... and Listing 4<https://

1087 www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>}.

1088 '''

1089 u *= r

1090 v *= s # 0 for 2-D root

1091 t0 = w - _1_0

1092 t1 = _0_0 if g < 0 else (hypot_(u, w, v) - _1_0)

1093 # assert t0 <= t1

1094 for i in range(1, _TRIPS): # 48..58

1095 t = (t1 + t0) * _0_5

1096 e = t1 - t0

1097 if eps > e > -eps or _isin(t, t0, t1):

1098 break

1099 g = fsumf_(_N_1_0, # ~= _hypot2_1

1100 _over02(u, t + r),

1101 _over02(w, t + _1_0), (

1102 _over02(v, t + s) if v else _0_0))

1103 if g > 0:

1104 t0 = t

1105 elif g < 0:

1106 t1 = t

1107 else:

1108 break

1109 else: # PYCHOK no cover

1110 t = Fmt.no_convergence(e, eps)

1111 raise _ValueError(t, txt__=_rootNd)

1112 return t, i

1113

1114

1115def _tri_biaxial(tri_biax, where):

1116 '''(INTERNAL) Get a triaxail for C{tri_biax}.

1117 '''

1118 if isinstance(tri_biax, _UnOrderedTriaxialBase):

1119 T = tri_biax

1120 else:

1121 D = tri_biax if isinstance(tri_biax, Datum) else \

1122 _spherical_datum(tri_biax, name__=where) # typename

1123 T = D.ellipsoid._triaxial

1124 return T

1125

1126

1127def _validate(a, b, c, d, T, x, y, z, val):

1128 '''(INTERNAL) Validate an C{_plumTo5} result.

1129 '''

1130 e = T.sideOf(a, b, c, eps=val)

1131 if e: # not near the ellipsoid's surface

1132 raise _ValueError(a=a, b=b, c=c, d=d,

1133 sideOf=e, eps=val)

1134 if d: # angle of delta and normal vector

1135 m = Vector3d(x, y, z).minus_(a, b, c)

1136 if m.euclid > val:

1137 m = m.unit()

1138 n = T.normal3d(a, b, c)

1139 e = n.dot(m) # n.negate().dot(m)

1140 if not isnear1(fabs(e), eps1=val):

1141 raise _ValueError(n=n, m=m,

1142 dot=e, eps=val)

1143

1144

1145class Triaxials(_TriaxialsBase):

1146 _Triaxial = Triaxial

1147

1148Triaxials = Triaxials(Triaxial, Triaxial_) # PYCHOK singleton

1149'''Some pre-defined L{Triaxial}s, all I{lazily} instantiated.'''

1150Triaxials._assert()

1151

1152if __name__ == _DMAIN_:

1153

1154 from pygeodesy.internals import _pregistry, printf

1155

1156 T = Triaxial_(6378388.0, 6378318.0, 6356911.9461)

1157 t = T.height4(3909863.9271, 3909778.123, 3170932.5016)

1158 printf('# Bektas: %r', t)

1159 # __doc__ of this file, force all into registry

1160 _pregistry(Triaxials)

1161

1162# % python3 -m pygeodesy.triaxials.triaxial5

1163#

1164# Bektas: height4(x=3909251.554667, y=3909165.750567, z=3170432.501602, h=999.999996)

1165

1166# **) MIT License

1167#

1169#

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

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

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

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

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

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

1176#

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

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

1179#

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

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

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

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

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

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

1186# OTHER DEALINGS IN THE SOFTWARE.