1

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

3

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

5 conformal projections L{JacobiConformal} and L{JacobiConformalSpherical}, transcoded from

6 I{Charles Karney}'s C++ class U{JacobiConformal<https://GeographicLib.SourceForge.io/C++/doc/

7 classGeographicLib_1_1JacobiConformal.html#details>} to pure Python and miscellaneous classes

8 L{BetaOmega2Tuple}, L{BetaOmega3Tuple}, L{Jacobi2Tuple} and L{TriaxialError}.

9

10 Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2008-2023). For more information,

11 see the U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.

12

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

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

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

16

17@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.490371704, area_p=93212299402.670425415)

18@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, area_p=4211301574065.829589844)

19@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=510065911057441.0625, area_p=510065915922713.6875)

20@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, area_p=798619018175.109863281)

21@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, area_p=30663773794562.45703125)

22@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.7421875, area_p=41691877397441.2109375)

23@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, area_p=144249144150662.15625)

24@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.691894531, area_p=493857714107.9375)

25@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.756958008, area_p=698881306767.950317383)

26@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, area_p=37838824733332.2265625)

27@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.394042969, area_p=3528074261832.738769531)

28@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.125, area_p=510065626587483.3125)

29'''

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

31 from __future__ import division as _; del _ # PYCHOK semicolon

32

33 from pygeodesy.basics import isscalar, map1, _zip, _ValueError

34 from pygeodesy.constants import EPS, EPS0, EPS02, EPS4, _EPS2e4, INT0, PI2, PI_3, PI4, \

35 _0_0, _0_5, _1_0, _N_2_0, float0_, isfinite, isnear1, \

36 _4_0 # PYCHOK used!

37 from pygeodesy.datums import Datum, _spherical_datum, _WGS84, Ellipsoid, Fmt

38# from pygeodesy.dms import toDMS # _MODS

39# from pygeodesy.ellipsoids import Ellipsoid # from .datums

40# from pygeodesy.elliptic import Elliptic # _MODS

41# from pygeodesy.errors import _ValueError # from .basics

42 from pygeodesy.fmath import Fdot, fdot, fmean_, hypot, hypot_, norm2

43 from pygeodesy.fsums import Fsum, fsumf_, fsum1f_

44 from pygeodesy.interns import NN, _a_, _b_, _beta_, _c_, _distant_, _finite_, \

45 _height_, _inside_, _near_, _not_, _NOTEQUAL_, _null_, \

46 _opposite_, _outside_, _SPACE_, _spherical_, _too_, \

47 _x_, _y_

48# from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS # from .vector3d

49 from pygeodesy.named import _NamedEnum, _NamedEnumItem, _NamedTuple, _Pass, \

50 _lazyNamedEnumItem as _lazy

51 from pygeodesy.namedTuples import LatLon3Tuple, Vector3Tuple, Vector4Tuple

52 from pygeodesy.props import Property_RO, property_RO

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

54 from pygeodesy.units import Degrees, Float, Height_, Meter, Meter2, Meter3, \

55 Radians, Radius, Scalar_

56 from pygeodesy.utily import asin1, atan2d, km2m, m2km, SinCos2, sincos2d_

57 from pygeodesy.vector3d import _otherV3d, Vector3d, _ALL_LAZY, _MODS

58

59 from math import atan2, fabs, sqrt

60

61__all__ = _ALL_LAZY.triaxials

62__version__ = '23.08.09'

63

64_not_ordered_ = _not_('ordered')

65_omega_ = 'omega'

66_TRIPS = 537 # 52..58, Eberly 1074?

67

68

69 class _NamedTupleTo(_NamedTuple): # in .testNamedTuples

70 '''(INTERNAL) Base for C{-.toDegrees}, C{-.toRadians}.

71 '''

72 def _toDegrees(self, a, b, *c, **toDMS_kwds):

73 if toDMS_kwds:

74 toDMS = _MODS.dms.toDMS

75 a = toDMS(a.toDegrees(), **toDMS_kwds)

76 b = toDMS(b.toDegrees(), **toDMS_kwds)

77 elif isinstance(a, Degrees) and \

78 isinstance(b, Degrees):

79 return self

80 else:

81 a, b = a.toDegrees(), b.toDegrees()

82 return self.classof(a, b, *c, name=self.name)

83

84 def _toRadians(self, a, b, *c):

85 return self if isinstance(a, Radians) and \

86 isinstance(b, Radians) else \

87 self.classof(a.toRadians(), b.toRadians(),

88 *c, name=self.name)

89

90

91 class BetaOmega2Tuple(_NamedTupleTo):

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

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

94 L{Degrees}).

95 '''

96 _Names_ = (_beta_, _omega_)

97 _Units_ = (_Pass, _Pass)

98

99 def toDegrees(self, **toDMS_kwds):

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

101

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

103 C{beta} and C{omega} both in L{Degrees}

104 or as a L{toDMS} string provided some

105 B{C{toDMS_kwds}} keyword arguments are

106 specified.

107 '''

108 return _NamedTupleTo._toDegrees(self, *self, **toDMS_kwds)

109

110 def toRadians(self):

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

112

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

114 C{beta} and C{omega} both in L{Radians}.

115 '''

116 return _NamedTupleTo._toRadians(self, *self)

117

118

119 class BetaOmega3Tuple(_NamedTupleTo):

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

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

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

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

124 in C{meter}, conventionally.

125 '''

126 _Names_ = BetaOmega2Tuple._Names_ + (_height_,)

127 _Units_ = BetaOmega2Tuple._Units_ + ( Meter,)

128

129 def toDegrees(self, **toDMS_kwds):

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

131

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

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

134 L{toDMS} string provided some B{C{toDMS_kwds}}

135 keyword arguments are specified.

136 '''

137 return _NamedTupleTo._toDegrees(self, *self, **toDMS_kwds)

138

139 def toRadians(self):

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

141

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

143 C{beta} and C{omega} both in L{Radians}.

144 '''

145 return _NamedTupleTo._toRadians(self, *self)

146

147 def to2Tuple(self):

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

149 '''

150 return BetaOmega2Tuple(*self[:2])

151

152

153 class Jacobi2Tuple(_NamedTupleTo):

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

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

156 '''

157 _Names_ = (_x_, _y_)

158 _Units_ = (_Pass, _Pass)

159

160 def toDegrees(self, **toDMS_kwds):

161 '''Convert this L{Jacobi2Tuple} to L{Degrees} or C{toDMS}.

162

163 @return: L{Jacobi2Tuple}C{(x, y)} with C{x} and C{y}

164 both in L{Degrees} or as a L{toDMS} string

165 provided some B{C{toDMS_kwds}} keyword

166 arguments are specified.

167 '''

168 return _NamedTupleTo._toDegrees(self, *self, **toDMS_kwds)

169

170 def toRadians(self):

171 '''Convert this L{Jacobi2Tuple} to L{Radians}.

172

173 @return: L{Jacobi2Tuple}C{(x, y)} with C{x}

174 and C{y} both in L{Radians}.

175 '''

176 return _NamedTupleTo._toRadians(self, *self)

177

178

179 class Triaxial_(_NamedEnumItem):

180 '''I{Unordered} triaxial ellipsoid and base class.

181

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

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

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

185 I{omega}=0.

186

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

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

189

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

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

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

193 L{Degrees} if converted).

194 '''

195 _ijk = _kji = None

196 _unordered = True

197

198 def __init__(self, a_triaxial, b=None, c=None, name=NN):

199 '''New I{unordered} L{Triaxial_}.

200

201 @arg a_triaxial: Large, C{X} semi-axis (C{scalar}, conventionally in

202 C{meter}) or an other L{Triaxial} or L{Triaxial_} instance.

203 @kwarg b: Middle, C{Y} semi-axis (C{meter}, same units as B{C{a}}), required

204 if C{B{a_triaxial} is scalar}, ignored otherwise.

205 @kwarg c: Small, C{Z} semi-axis (C{meter}, same units as B{C{a}}), required

206 if C{B{a_triaxial} is scalar}, ignored otherwise.

207 @kwarg name: Optional name (C{str}).

208

209 @raise TriaxialError: Invalid semi-axis or -axes.

210 '''

211 try:

212 a = a_triaxial

213 t = a._abc3 if isinstance(a, Triaxial_) else (

214 Radius(a=a), Radius(b=b), Radius(c=c))

215 except (TypeError, ValueError) as x:

216 raise TriaxialError(a=a, b=b, c=c, cause=x)

217 if name:

218 self.name = name

219

220 a, b, c = self._abc3 = t

221 if self._unordered: # == not isinstance(self, Triaxial)

222 s, _, t = sorted(t)

223 if not (isfinite(t) and s > 0):

224 raise TriaxialError(a=a, b=b, c=c) # txt=_invalid_

225 elif not (isfinite(a) and a >= b >= c > 0):

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

227 elif not (a > c and self._a2c2 > 0 and self.e2ac > 0):

228 raise TriaxialError(a=a, c=c, e2ac=self.e2ac, txt=_spherical_)

229

230 def __str__(self):

231 return self.toStr()

232

233 @Property_RO

234 def a(self):

235 '''Get the (largest) C{x} semi-axis (C{meter}, conventionally).

236 '''

237 a, _, _ = self._abc3

238 return a

239

240 @Property_RO

241 def _a2b2(self):

242 '''(INTERNAL) Get C{a**2 - b**2} == E_sub_e**2.

243 '''

244 a, b, _ = self._abc3

245 return ((a - b) * (a + b)) if a != b else _0_0

246

247 @Property_RO

248 def _a2_b2(self):

249 '''(INTERNAL) Get C{(a/b)**2}.

250 '''

251 a, b, _ = self._abc3

252 return (a / b)**2 if a != b else _1_0

253

254 @Property_RO

255 def _a2c2(self):

256 '''(INTERNAL) Get C{a**2 - c**2} == E_sub_x**2.

257 '''

258 a, _, c = self._abc3

259 return ((a - c) * (a + c)) if a != c else _0_0

260

261 @Property_RO

262 def area(self):

263 '''Get the surface area (C{meter} I{squared}).

264 '''

265 c, b, a = sorted(self._abc3)

266 if a > c:

267 a = Triaxial(a, b, c).area if a > b else \

268 Ellipsoid(a, b=c).areax # a == b

269 else: # a == c == b

270 a = Meter2(area=a**2 * PI4)

271 return a

272

273 def area_p(self, p=1.6075):

274 '''I{Approximate} the surface area (C{meter} I{squared}).

275

276 @kwarg p: Exponent (C{scalar} > 0), 1.6 for near-spherical or 1.5849625007

277 for "near-flat" triaxials.

278

279 @see: U{Surface area<https://WikiPedia.org/wiki/Ellipsoid#Approximate_formula>}.

280 '''

281 a, b, c = self._abc3

282 if a == b == c:

283 a *= a

284 else:

285 _p = pow

286 a = _p(fmean_(_p(a * b, p), _p(a * c, p), _p(b * c, p)), _1_0 / p)

287 return Meter2(area_p=a * PI4)

288

289 @Property_RO

290 def b(self):

291 '''Get the (middle) C{y} semi-axis (C{meter}, same units as B{C{a}}).

292 '''

293 _, b, _ = self._abc3

294 return b

295

296 @Property_RO

297 def _b2c2(self):

298 '''(INTERNAL) Get C{b**2 - c**2} == E_sub_y**2.

299 '''

300 _, b, c = self._abc3

301 return ((b - c) * (b + c)) if b != c else _0_0

302

303 @Property_RO

304 def c(self):

305 '''Get the (smallest) C{z} semi-axis (C{meter}, same units as B{C{a}}).

306 '''

307 _, _, c = self._abc3

308 return c

309

310 @Property_RO

311 def _c2_b2(self):

312 '''(INTERNAL) Get C{(c/b)**2}.

313 '''

314 _, b, c = self._abc3

315 return (c / b)**2 if b != c else _1_0

316

317 @Property_RO

318 def e2ab(self):

319 '''Get the C{ab} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (b/a)**2}.

320 '''

321 return Float(e2ab=(_1_0 - self._1e2ab) or _0_0)

322

323 @Property_RO

324 def _1e2ab(self):

325 '''(INTERNAL) Get C{1 - e2ab} == C{(b/a)**2}.

326 '''

327 a, b, _ = self._abc3

328 return (b / a)**2 if a != b else _1_0

329

330 @Property_RO

331 def e2ac(self):

332 '''Get the C{ac} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (c/a)**2}.

333 '''

334 return Float(e2ac=(_1_0 - self._1e2ac) or _0_0)

335

336 @Property_RO

337 def _1e2ac(self):

338 '''(INTERNAL) Get C{1 - e2ac} == C{(c/a)**2}.

339 '''

340 a, _, c = self._abc3

341 return (c / a)**2 if a != c else _1_0

342

343 @Property_RO

344 def e2bc(self):

345 '''Get the C{bc} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (c/b)**2}.

346 '''

347 return Float(e2bc=(_1_0 - self._1e2bc) or _0_0)

348

349 _1e2bc = _c2_b2 # C{1 - e2bc} == C{(c/b)**2}

350

351 @property_RO

352 def _Elliptic(self):

353 '''(INTERNAL) Get class L{Elliptic}, I{once}.

354 '''

355 Triaxial_._Elliptic = E = _MODS.elliptic.Elliptic # overwrite property_RO

356 return E

357

358 def hartzell4(self, pov, los=None, name=NN):

359 '''Compute the intersection of this triaxial's surface with a Line-Of-Sight

360 from a Point-Of-View in space.

361

362 @see: Function L{pygeodesy.hartzell4} for further details.

363 '''

364 return hartzell4(pov, los=los, tri_biax=self, name=name)

365

366 def height4(self, x_xyz, y=None, z=None, normal=True, eps=EPS):

367 '''Compute the projection on and the height of a cartesian above or below

368 this triaxial's surface.

369

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

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

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

373 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.

374 @kwarg normal: If C{True} the projection is perpendicular to (the nearest

375 point on) this triaxial's surface, otherwise the C{radial}

376 line to this triaxial's center (C{bool}).

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

378 negative value to skip validation.

379

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

381 C{x}, C{y} and C{z} of the projection on or the intersection

382 with and with the height C{h} above or below the triaxial's

383 surface in C{meter}, conventionally.

384

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

386 convergence in root finding or validation failed.

387

388 @see: Method L{Ellipsoid.height4} and I{Eberly}'s U{Distance from a Point

389 to ... an Ellipsoid ...<https://www.GeometricTools.com/Documentation/

390 DistancePointEllipseEllipsoid.pdf>}.

391 '''

392 v, r = _otherV3d_(x_xyz, y, z), self.isSpherical

393

394 i, h = None, v.length

395 if h < EPS0: # EPS

396 x = y = z = _0_0

397 h -= min(self._abc3) # nearest

398 elif r: # .isSpherical

399 x, y, z = v.times(r / h).xyz

400 h -= r

401 else:

402 x, y, z = v.xyz

403 try:

404 if normal: # perpendicular to triaxial

405 x, y, z, h, i = _normalTo5(x, y, z, self, eps=eps)

406 else: # radially to triaxial's center

407 x, y, z = self._radialTo3(z, hypot(x, y), y, x)

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

409 except Exception as e:

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

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

412 h = -h # below the surface

413 return Vector4Tuple(x, y, z, h, iteration=i, name=self.height4.__name__)

414

415 @Property_RO

416 def isOrdered(self):

417 '''Is this triaxial I{ordered} and I{not spherical} (C{bool})?

418 '''

419 a, b, c = self._abc3

420 return bool(a >= b > c) # b > c!

421

422 @Property_RO

423 def isSpherical(self):

424 '''Is this triaxial I{spherical} (C{Radius} or INT0)?

425 '''

426 a, b, c = self._abc3

427 return a if a == b == c else INT0

428

429 def normal3d(self, x_xyz, y=None, z=None, length=_1_0):

430 '''Get a 3-D vector perpendicular to at a cartesian on this triaxial's surface.

431

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

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

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

435 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.

436 @kwarg length: Optional length and in-/outward direction (C{scalar}).

437

438 @return: A C{Vector3d(x_, y_, z_)} normalized to B{C{length}}, pointing

439 in- or outward for neg- respectively positive B{C{length}}.

440

441 @note: Cartesian location C{(B{x}, B{y}, B{z})} must be on this triaxial's

442 surface, use method L{Triaxial.sideOf} to validate.

443 '''

444 # n = 2 * (x / a2, y / b2, z / c2)

445 # == 2 * (x, y * a2 / b2, z * a2 / c2) / a2 # iff ordered

446 # == 2 * (x, y / _1e2ab, z / _1e2ac) / a2

447 # == unit(x, y / _1e2ab, z / _1e2ac).times(length)

448 n = self._normal3d.times_(*_otherV3d_(x_xyz, y, z).xyz)

449 if n.length < EPS0:

450 raise TriaxialError(x=x_xyz, y=y, z=z, txt=_null_)

451 return n.times(length / n.length)

452

453 @Property_RO

454 def _normal3d(self):

455 '''(INTERNAL) Get M{Vector3d((d/a)**2, (d/b)**2, (d/c)**2)}, M{d = max(a, b, c)}.

456 '''

457 d = max(self._abc3)

458 t = tuple(((d / x)**2 if x != d else _1_0) for x in self._abc3)

459 return Vector3d(*t, name=self.normal3d.__name__)

460

461 def _norm2(self, s, c, *a):

462 '''(INTERNAL) Normalize C{s} and C{c} iff not already.

463 '''

464 if fabs(_hypot21(s, c)) > EPS02:

465 s, c = norm2(s, c)

466 if a:

467 s, c = norm2(s * self.b, c * a[0])

468 return float0_(s, c)

469

470 def _order3(self, *abc, **reverse): # reverse=False

471 '''(INTERNAL) Un-/Order C{a}, C{b} and C{c}.

472

473 @return: 3-Tuple C{(a, b, c)} ordered by or un-ordered

474 (reverse-ordered) C{ijk} if C{B{reverse}=True}.

475 '''

476 ijk = self._order_ijk(**reverse)

477 return _getitems(abc, *ijk) if ijk else abc

478

479 def _order3d(self, v, **reverse): # reverse=False

480 '''(INTERNAL) Un-/Order a C{Vector3d}.

481

482 @return: Vector3d(x, y, z) un-/ordered.

483 '''

484 ijk = self._order_ijk(**reverse)

485 return v.classof(*_getitems(v.xyz, *ijk)) if ijk else v

486

487 @Property_RO

488 def _ordered4(self):

489 '''(INTERNAL) Helper for C{_hartzell3d2} and C{_normalTo5}.

490 '''

491 def _order2(reverse, a, b, c):

492 '''(INTERNAL) Un-Order C{a}, C{b} and C{c}.

493

494 @return: 2-Tuple C{((a, b, c), ijk)} with C{a} >= C{b} >= C{c}

495 and C{ijk} a 3-tuple with the initial indices.

496 '''

497 i, j, k = 0, 1, 2 # range(3)

498 if a < b:

499 a, b, i, j = b, a, j, i

500 if a < c:

501 a, c, i, k = c, a, k, i

502 if b < c:

503 b, c, j, k = c, b, k, j

504 # reverse (k, j, i) since (a, b, c) is reversed-sorted

505 ijk = (k, j, i) if reverse else (None if i < j < k else (i, j, k))

506 return (a, b, c), ijk

507

508 abc, T = self._abc3, self

509 if not self.isOrdered:

510 abc, ijk = _order2(False, *abc)

511 if ijk:

512 _, kji = _order2(True, *ijk)

513 T = Triaxial_(*abc)

514 T._ijk, T._kji = ijk, kji

515 return abc + (T,)

516

517 def _order_ijk(self, reverse=False):

518 '''(INTERNAL) Get the un-/order indices.

519 '''

520 return self._kji if reverse else self._ijk

521

522 def _radialTo3(self, sbeta, cbeta, somega, comega):

523 '''(INTERNAL) I{Unordered} helper for C{.height4}.

524 '''

525 def _rphi(a, b, sphi, cphi):

526 # <https://WikiPedia.org/wiki/Ellipse#Polar_form_relative_to_focus>

527 # polar form: radius(phi) = a * b / hypot(a * sphi, b * cphi)

528 return (b / hypot(sphi, b / a * cphi)) if a > b else (

529 (a / hypot(cphi, a / b * sphi)) if a < b else a)

530

531 sa, ca = self._norm2(sbeta, cbeta)

532 sb, cb = self._norm2(somega, comega)

533

534 a, b, c = self._abc3

535 if a != b:

536 a = _rphi(a, b, sb, cb)

537 if a != c:

538 c = _rphi(a, c, sa, ca)

539 z, r = c * sa, c * ca

540 x, y = r * cb, r * sb

541 return x, y, z

542

543 def sideOf(self, x_xyz, y=None, z=None, eps=EPS4):

544 '''Is a cartesian above, below or on the surface of this triaxial?

545

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

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

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

549 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.

550 @kwarg eps: Near surface tolerance(C{scalar}).

551

552 @return: C{INT0} if C{(B{x}, B{y}, B{z})} is near this triaxial's surface

553 within tolerance B{C{eps}}, otherwise a neg- or positive C{float}

554 if in- respectively outside this triaxial.

555

556 @see: Methods L{Triaxial.height4} and L{Triaxial.normal3d}.

557 '''

558 return _sideOf(_otherV3d_(x_xyz, y, z).xyz, self._abc3, eps=eps)

559

560 def _sqrt(self, x):

561 '''(INTERNAL) Helper, see L{pygeodesy.sqrt0}.

562 '''

563 if x < 0:

564 raise TriaxialError(Fmt.PAREN(sqrt=x))

565 return _0_0 if x < EPS02 else sqrt(x)

566

567 def toEllipsoid(self, name=NN):

568 '''Convert this triaxial to an L{Ellipsoid}, provided 2 axes match.

569

570 @return: An L{Ellipsoid} with north along this C{Z} axis if C{a == b},

571 this C{Y} axis if C{a == c} or this C{X} axis if C{b == c}.

572

573 @raise TriaxialError: This C{a != b}, C{b != c} and C{c != a}.

574

575 @see: Method L{Ellipsoid.toTriaxial}.

576 '''

577 a, b, c = self._abc3

578 if a == b:

579 b = c # N = c-Z

580 elif b == c: # N = a-X

581 a, b = b, a

582 elif a != c: # N = b-Y

583 t = _SPACE_(_a_, _NOTEQUAL_, _b_, _NOTEQUAL_, _c_)

584 raise TriaxialError(a=a, b=b, c=c, txt=t)

585 return Ellipsoid(a, b=b, name=name or self.name)

586

587 def toStr(self, prec=9, name=NN, **unused): # PYCHOK signature

588 '''Return this C{Triaxial} as a string.

589

590 @kwarg prec: Precision, number of decimal digits (0..9).

591 @kwarg name: Override name (C{str}) or C{None} to exclude

592 this triaxial's name.

593

594 @return: This C{Triaxial}'s attributes (C{str}).

595 '''

596 T = Triaxial_

597 t = T.a,

598 J = JacobiConformalSpherical

599 t += (J.ab, J.bc) if isinstance(self, J) else (T.b, T.c)

600 t += T.e2ab, T.e2bc, T.e2ac

601 J = JacobiConformal

602 if isinstance(self, J):

603 t += J.xyQ2,

604 t += T.volume, T.area

605 return self._instr(name, prec, props=t, area_p=self.area_p())

606

607 @Property_RO

608 def volume(self):

609 '''Get the volume (C{meter**3}), M{4 / 3 * PI * a * b * c}.

610 '''

611 a, b, c = self._abc3

612 return Meter3(volume=a * b * c * PI_3 * _4_0)

613

614

615 class Triaxial(Triaxial_):

616 '''I{Ordered} triaxial ellipsoid.

617

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

619 '''

620 _unordered = False

621

622 def __init__(self, a_triaxial, b=None, c=None, name=NN):

623 '''New I{ordered} L{Triaxial}.

624

625 @arg a_triaxial: Largest semi-axis (C{scalar}, conventionally in C{meter})

626 or an other L{Triaxial} or L{Triaxial_} instance.

627 @kwarg b: Middle semi-axis (C{meter}, same units as B{C{a}}), required

628 if C{B{a_triaxial} is scalar}, ignored otherwise.

629 @kwarg c: Smallest semi-axis (C{meter}, same units as B{C{a}}), required

630 if C{B{a_triaxial} is scalar}, ignored otherwise.

631 @kwarg name: Optional name (C{str}).

632

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

634 must be ellipsoidal, C{B{a} > B{c}}.

635

636 @raise TriaxialError: Semi-axes not ordered, spherical or invalid.

637 '''

638 Triaxial_.__init__(self, a_triaxial, b=b, c=c, name=name)

639

640 @Property_RO

641 def _a2b2_a2c2(self):

642 '''@see: Methods C{.forwardBetaOmega} and C{._k2_kp2}.

643 '''

644 return self._a2b2 / self._a2c2

645

646 @Property_RO

647 def area(self):

648 '''Get the surface area (C{meter} I{squared}).

649

650 @see: U{Surface area<https://WikiPedia.org/wiki/Ellipsoid#Surface_area>}.

651 '''

652 a, b, c = self._abc3

653 if a != b:

654 kp2, k2 = self._k2_kp2 # swapped!

655 aE = self._Elliptic(k2, _0_0, kp2, _1_0)

656 c2 = self._1e2ac # cos(phi)**2 = (c/a)**2

657 s = sqrt(self.e2ac) # sin(phi)**2 = 1 - c2

658 r = asin1(s) # phi = atan2(sqrt(c2), s)

659 b *= fsum1f_(aE.fE(r) * s, c / a * c / b,

660 aE.fF(r) * c2 / s)

661 a = Meter2(area=a * b * PI2)

662 else: # a == b > c

663 a = Ellipsoid(a, b=c).areax

664 return a

665

666 def _exyz3(self, u):

667 '''(INTERNAL) Helper for C{.forwardBetOmg}.

668 '''

669 if u > 0:

670 u2 = u**2

671 x = u * self._sqrt(_1_0 + self._a2c2 / u2)

672 y = u * self._sqrt(_1_0 + self._b2c2 / u2)

673 else:

674 x = y = u = _0_0

675 return x, y, u

676

677 def forwardBetaOmega(self, beta, omega, height=0, name=NN):

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

679 and height to cartesian.

680

681 @arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}).

682 @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}).

683 @arg height: Height above or below the ellipsoid's surface (C{meter}, same

684 units as this triaxial's C{a}, C{b} and C{c} semi-axes).

685 @kwarg name: Optional name (C{str}).

686

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

688

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

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

691 '''

692 if height:

693 h = Height_(height=height, low=-self.c, Error=TriaxialError)

694 x, y, z = self._exyz3(h + self.c)

695 else:

696 x, y, z = self._abc3 # == self._exyz3(self.c)

697 if z: # and x and y:

698 sa, ca = SinCos2(beta)

699 sb, cb = SinCos2(omega)

700

701 r = self._a2b2_a2c2

702 x *= cb * self._sqrt(ca**2 + r * sa**2)

703 y *= ca * sb

704 z *= sa * self._sqrt(_1_0 - r * cb**2)

705 return Vector3Tuple(x, y, z, name=name)

706

707 def forwardBetaOmega_(self, sbeta, cbeta, somega, comega, name=NN):

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

709 to cartesian coordinates I{on the triaxial's surface}.

710

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

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

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

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

715 @kwarg name: Optional name (C{str}).

716

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

718

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

720

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

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

723 Triaxial_ellipsoid_coordinate_system>} and U{expressions (23-25)<https://

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

725 '''

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

727 return Vector3Tuple(*t, name=name)

728

729 def forwardCartesian(self, x_xyz, y=None, z=None, name=NN, **normal_eps):

730 '''Project a cartesian on this triaxial.

731

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

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

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

735 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.

736 @kwarg name: Optional name (C{str}).

737 @kwarg normal_eps: Optional keyword arguments C{B{normal}=True} and

738 C{B{eps}=EPS}, see method L{Triaxial.height4}.

739

740 @see: Method L{Triaxial.height4} for further information and method

741 L{Triaxial.reverseCartesian} to reverse the projection.

742 '''

743 t = self.height4(x_xyz, y, z, **normal_eps)

744 _ = t.rename(name)

745 return t

746

747 def forwardLatLon(self, lat, lon, height=0, name=NN):

748 '''Convert I{geodetic} lat-, longitude and heigth to cartesian.

749

750 @arg lat: Geodetic latitude (C{degrees}).

751 @arg lon: Geodetic longitude (C{degrees}).

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

753 as this triaxial's C{a}, C{b} and C{c} axes).

754 @kwarg name: Optional name (C{str}).

755

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

757

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

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

760 '''

761 return self._forwardLatLon3(height, name, *sincos2d_(lat, lon))

762

763 def forwardLatLon_(self, slat, clat, slon, clon, height=0, name=NN):

764 '''Convert I{geodetic} lat-, longitude and heigth to cartesian.

765

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

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

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

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

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

771 as this triaxial's axes C{a}, C{b} and C{c}).

772 @kwarg name: Optional name (C{str}).

773

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

775

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

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

778 '''

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

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

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

782

783 def _forwardLatLon3(self, h, name, sa, ca, sb, cb):

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

785 '''

786 ca_x_sb = ca * sb

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

788 t = fsumf_(_1_0, -self.e2ac * sa**2, -self.e2ab * ca_x_sb**2)

789 n = self.a / self._sqrt(t) # prime vertical

790 x = (h + n) * ca * cb

791 y = (h + n * self._1e2ab) * ca_x_sb

792 z = (h + n * self._1e2ac) * sa

793 return Vector3Tuple(x, y, z, name=name)

794

795 @Property_RO

796 def _k2_kp2(self):

797 '''(INTERNAL) Get C{k2} and C{kp2} for C{._xE}, C{._yE} and C{.area}.

798 '''

799 # k2 = a2b2 / a2c2 * c2_b2

800 # kp2 = b2c2 / a2c2 * a2_b2

801 # b2 = b**2

802 # xE = Elliptic(k2, -a2b2 / b2, kp2, a2_b2)

803 # yE = Elliptic(kp2, +b2c2 / b2, k2, c2_b2)

804 # aE = Elliptic(kp2, 0, k2, 1)

805 return (self._a2b2_a2c2 * self._c2_b2,

806 self._b2c2 / self._a2c2 * self._a2_b2)

807

808 def _radialTo3(self, sbeta, cbeta, somega, comega):

809 '''(INTERNAL) Convert I{ellipsoidal} lat- and longitude C{beta} and

810 C{omega} to cartesian coordinates I{on the triaxial's surface},

811 also I{ordered} helper for C{.height4}.

812 '''

813 sa, ca = self._norm2(sbeta, cbeta)

814 sb, cb = self._norm2(somega, comega)

815

816 b2_a2 = self._1e2ab # == (b/a)**2

817 c2_a2 = -self._1e2ac # == -(c/a)**2

818 a2c2_a2 = self. e2ac # (a**2 - c**2) / a**2 == 1 - (c/a)**2

819

820 x2 = Fsum(_1_0, -b2_a2 * sa**2, c2_a2 * ca**2).fover(a2c2_a2)

821 z2 = Fsum(c2_a2, sb**2, b2_a2 * cb**2).fover(a2c2_a2)

822

823 x, y, z = self._abc3

824 x *= cb * self._sqrt(x2)

825 y *= ca * sb

826 z *= sa * self._sqrt(z2)

827 return x, y, z

828

829 def reverseBetaOmega(self, x_xyz, y=None, z=None, name=NN):

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

831 and height.

832

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

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

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

836 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.

837 @kwarg name: Optional name (C{str}).

838

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

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

841 units as this triaxial's axes.

842

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

844 and U{Expressions (21-22)<https://www.Topo.Auth.GR/wp-content/uploads/

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

846 '''

847 v = _otherV3d_(x_xyz, y, z)

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

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

850

851 def reverseCartesian(self, x_xyz, y=None, z=None, h=0, normal=True, eps=_EPS2e4, name=NN):

852 '''"Unproject" a cartesian on to a cartesion I{off} this triaxial's surface.

853

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

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

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

857 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.

858 @arg h: Height above or below this triaxial's surface (C{meter}, same units

859 as the axes).

860 @kwarg normal: If C{True} the height is C{normal} to the surface, otherwise

861 C{radially} to the center of this triaxial (C{bool}).

862 @kwarg eps: Tolerance for surface test (C{scalar}).

863 @kwarg name: Optional name (C{str}).

864

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

866

867 @raise TrialError: Cartesian C{(x, y, z)} not on this triaxial's surface.

868

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

870 '''

871 v = _otherV3d_(x_xyz, y, z, name=name)

872 s = _sideOf(v.xyz, self._abc3, eps=eps)

873 if s: # PYCHOK no cover

874 t = _SPACE_((_inside_ if s < 0 else _outside_), self.toRepr())

875 raise TriaxialError(eps=eps, sideOf=s, x=v.x, y=v.y, z=v.z, txt=t)

876

877 if h:

878 if normal:

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

880 elif v.length > EPS0:

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

882 return v.xyz # Vector3Tuple

883

884 def reverseLatLon(self, x_xyz, y=None, z=None, name=NN):

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

886

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

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

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

890 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.

891 @kwarg name: Optional name (C{str}).

892

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

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

895 as this triaxial's axes.

896

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

898 and U{Expressions (4-5)<https://www.Topo.Auth.GR/wp-content/uploads/

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

900 '''

901 v = _otherV3d_(x_xyz, y, z)

902 s = v.times_(self._1e2ac, # == 1 - e_sub_x**2

903 self._1e2bc, # == 1 - e_sub_y**2

904 _1_0)

905 t = self._reverseLatLon3(s, atan2d, v, self.forwardLatLon_)

906 return LatLon3Tuple(*t, name=name)

907

908 def _reverseLatLon3(self, s, atan2_, v, forward_):

909 '''(INTERNAL) Helper for C{.reverseBetOmg} and C{.reverseLatLon}.

910 '''

911 x, y, z = s.xyz

912 d = hypot( x, y)

913 a = atan2_(z, d)

914 b = atan2_(y, x)

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

916 return a, b, h

917

918

919 class JacobiConformal(Triaxial):

920 '''This is a conformal projection of a triaxial ellipsoid to a plane in which the

921 C{X} and C{Y} grid lines are straight.

922

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

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

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

926 returned in the case of an ellipsoid of revolution.

927

928 Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2014-2023) and

929 licensed under the MIT/X11 License.

930

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

932 L{JacobiConformalSpherical}.

933

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

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

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

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

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

939 '''

940

941 @Property_RO

942 def _xE(self):

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

944 '''

945 k2, kp2 = self._k2_kp2

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

947 return self._Elliptic(k2, _1_0 - self._a2_b2, kp2, self._a2_b2)

948

949 def xR(self, omega):

950 '''Compute a Jacobi Conformal C{x} projection.

951

952 @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}).

953

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

955 '''

956 return self.xR_(*SinCos2(omega))

957

958 def xR_(self, somega, comega):

959 '''Compute a Jacobi Conformal C{x} projection.

960

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

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

963

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

965 '''

966 s, c = self._norm2(somega, comega, self.a)

967 return Radians(x=self._xE.fPi(s, c) * self._a2_b2)

968

969 @Property_RO

970 def xyQ2(self):

971 '''Get the Jacobi Conformal quadrant size (L{Jacobi2Tuple}C{(x, y)}).

972 '''

973 return Jacobi2Tuple(Radians(x=self._a2_b2 * self._xE.cPi),

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

975 name=JacobiConformal.xyQ2.name)

976

977 def xyR2(self, beta, omega, name=NN):

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

979

980 @arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}).

981 @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}).

982 @kwarg name: Optional name (C{str}).

983

984 @return: A L{Jacobi2Tuple}C{(x, y)}.

985 '''

986 return self.xyR2_(*(SinCos2(beta) + SinCos2(omega)),

987 name=name or self.xyR2.__name__)

988

989 def xyR2_(self, sbeta, cbeta, somega, comega, name=NN):

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

991

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

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

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

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

996 @kwarg name: Optional name (C{str}).

997

998 @return: A L{Jacobi2Tuple}C{(x, y)}.

999 '''

1000 return Jacobi2Tuple(self.xR_(somega, comega),

1001 self.yR_(sbeta, cbeta),

1002 name=name or self.xyR2_.__name__)

1003

1004 @Property_RO

1005 def _yE(self):

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

1007 '''

1008 kp2, k2 = self._k2_kp2 # swapped!

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

1010 return self._Elliptic(k2, self.e2bc, kp2, self._c2_b2)

1011

1012 def yR(self, beta):

1013 '''Compute a Jacobi Conformal C{y} projection.

1014

1015 @arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}).

1016

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

1018 '''

1019 return self.yR_(*SinCos2(beta))

1020

1021 def yR_(self, sbeta, cbeta):

1022 '''Compute a Jacobi Conformal C{y} projection.

1023

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

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

1026

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

1028 '''

1029 s, c = self._norm2(sbeta, cbeta, self.c)

1030 return Radians(y=self._yE.fPi(s, c) * self._c2_b2)

1031

1032

1033 class JacobiConformalSpherical(JacobiConformal):

1034 '''An alternate, I{spherical} L{JacobiConformal} projection.

1035

1036 @see: L{JacobiConformal} for other and more details.

1037 '''

1038 _ab = _bc = 0

1039

1040 def __init__(self, radius_triaxial, ab=0, bc=0, name=NN):

1041 '''New L{JacobiConformalSpherical}.

1042

1043 @arg radius_triaxial: Radius (C{scalar}, conventionally in

1044 C{meter}) or an other L{JacobiConformalSpherical},

1045 L{JacobiConformal} or ordered L{Triaxial}.

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

1047 same units as C{scalar B{radius}}.

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

1049 same units as C{scalar B{radius}}.

1050 @kwarg name: Optional name (C{str}).

1051

1052 @raise TriaxialError: Invalid B{C{radius_triaxial}}, negative

1053 B{C{ab}}, negative B{C{bc}} or C{(B{ab}

1054 + B{bc})} not positive.

1055

1056 @note: If B{C{radius_triaxial}} is a L{JacobiConformalSpherical}

1057 and if B{C{ab}} and B{C{bc}} are both zero or C{None},

1058 the B{C{radius_triaxial}}'s C{ab}, C{bc}, C{a}, C{b}

1059 and C{c} are copied.

1060 '''

1061 try:

1062 r, j = radius_triaxial, False

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

1064 if (not (ab or bc)) and isinstance(r, JacobiConformalSpherical):

1065 j = True

1066 t = r._abc3

1067 else:

1068 t = (Radius(radius=r),) * 3

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

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

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

1072 raise ValueError('(ab + bc)')

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

1074 if not (a >= c and isfinite(self._a2b2)

1075 and isfinite(self._a2c2)):

1076 raise ValueError(_not_(_finite_))

1077 except (TypeError, ValueError) as x:

1078 raise TriaxialError(radius_triaxial=r, ab=ab, bc=bc, cause=x)

1079 if name:

1080 self.name = name

1081

1082 @Property_RO

1083 def ab(self):

1084 '''Get relative magnitude C{ab} (C{meter}, same units as B{C{a}}).

1085 '''

1086 return self._ab

1087

1088 @Property_RO

1089 def _a2b2(self):

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

1091 '''

1092 a, b, _ = self._abc3

1093 return self.ab * (a + b)

1094

1095 @Property_RO

1096 def _a2c2(self):

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

1098 '''

1099 return self._a2b2 + self._b2c2

1100

1101 @Property_RO

1102 def bc(self):

1103 '''Get relative magnitude C{bc} (C{meter}, same units as B{C{a}}).

1104 '''

1105 return self._bc

1106

1107 @Property_RO

1108 def _b2c2(self):

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

1110 '''

1111 _, b, c = self._abc3

1112 return self.bc * (b + c)

1113

1114 @Property_RO

1115 def radius(self):

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

1117 '''

1118 return self.a

1119

1120

1121 class TriaxialError(_ValueError):

1122 '''Raised for L{Triaxial} issues.

1123 '''

1124 pass # ...

1125

1126

1127 class Triaxials(_NamedEnum):

1128 '''(INTERNAL) L{Triaxial} registry, I{must} be a sub-class

1129 to accommodate the L{_LazyNamedEnumItem} properties.

1130 '''

1131 def _Lazy(self, *abc, **name):

1132 '''(INTERNAL) Instantiate the C{Triaxial}.

1133 '''

1134 a, b, c = map(km2m, abc)

1135 return Triaxial(a, b, c, **name)

1136

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

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

1139# <https://ArxIV.org/pdf/1909.06452.pdf> Table 1 Semi-axes in Km

1140# <https://www.JPS.NASA.gov/education/images/pdf/ss-moons.pdf>

1141# <https://link.Springer.com/article/10.1007/s00190-022-01650-9>

1142_E = _WGS84.ellipsoid

1143 Triaxials._assert( # a (Km) b (Km) c (Km) planet

1144 Amalthea = _lazy('Amalthea', 125.0, 73.0, 64), # Jupiter

1145 Ariel = _lazy('Ariel', 581.1, 577.9, 577.7), # Uranus

1146 Earth = _lazy('Earth', 6378.173435, 6378.1039, 6356.7544),

1147 Enceladus = _lazy('Enceladus', 256.6, 251.4, 248.3), # Saturn

1148 Europa = _lazy('Europa', 1564.13, 1561.23, 1560.93), # Jupiter

1149 Io = _lazy('Io', 1829.4, 1819.3, 1815.7), # Jupiter

1150 Mars = _lazy('Mars', 3394.6, 3393.3, 3376.3),

1151 Mimas = _lazy('Mimas', 207.4, 196.8, 190.6), # Saturn

1152 Miranda = _lazy('Miranda', 240.4, 234.2, 232.9), # Uranus

1153 Moon = _lazy('Moon', 1735.55, 1735.324, 1734.898), # Earth

1154 Tethys = _lazy('Tethys', 535.6, 528.2, 525.8), # Saturn

1155 WGS84_35 = _lazy('WGS84_35', *map1(m2km, _E.a + 35, _E.a - 35, _E.b)))

1156 del _E

1157

1158

1159 def _getitems(items, *indices):

1160 '''(INTERNAL) Get the C{items} at the given I{indices}.

1161

1162 @return: C{Type(items[i] for i in indices)} with

1163 C{Type = type(items)}, any C{type} having

1164 the special method C{__getitem__}.

1165 '''

1166 return type(items)(map(items.__getitem__, indices))

1167

1168

1169 def _hartzell3d2(pov, los, Tun): # MCCABE 13 in .ellipsoidal.hartzell4, .formy.hartzell

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

1171 formula for I{un-/ordered} triaxials.

1172 '''

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

1174

1175 a2 = a**2 # largest, factored out

1176 b2, p2 = (b**2, T._1e2ab) if b != a else (a2, _1_0)

1177 c2, q2 = (c**2, T._1e2ac) if c != a else (a2, _1_0)

1178

1179 p3 = T._order3d(_otherV3d(pov=pov))

1180 u3 = T._order3d(_otherV3d(los=los)) if los else p3.negate()

1181 u3 = u3.unit() # unit vector, opposing signs

1182

1183 x2, y2, z2 = p3.x2y2z2 # p3.times_(p3).xyz

1184 ux, vy, wz = u3.times_(p3).xyz

1185 u2, v2, w2 = u3.x2y2z2 # u3.times_(u3).xyz

1186

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

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

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

1190 raise _ValueError(_near_(_null_))

1191

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

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

1194 -w2 * x2 * p2, -u2 * y2 * q2, -v2 * x2 * q2, ux * vy * 2 * q2)

1195 if r > 0: # a2 factored out

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

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

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

1199

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

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

1202 s = fsumf_(_1_0, x2 / a2, y2 / b2, z2 / c2, _N_2_0) # like _sideOf

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

1204 elif fsum1f_(x2, y2, z2) < d**2: # d past triaxial's center

1205 raise _ValueError(_too_(_distant_))

1206

1207 v = p3.minus(u3.times(d)) # Vector3d

1208 h = p3.minus(v).length # distance to triaxial

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

1210

1211

1212 def hartzell4(pov, los=None, tri_biax=_WGS84, name=NN):

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

1214 from a Point-Of-View outside.

1215

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

1217 or L{Vector3d}).

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

1219 C{None} to point to the tri-/biaxial's center.

1220 @kwarg tri_biax: A triaxial (L{Triaxial}, L{Triaxial_}, L{JacobiConformal} or

1221 L{JacobiConformalSpherical}) or biaxial ellipsoid (L{Datum},

1222 L{Ellipsoid}, L{Ellipsoid2}, L{a_f2Tuple} or C{scalar} radius,

1223 conventionally in C{meter}).

1224 @kwarg name: Optional name (C{str}).

1225

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

1227 C{h} the distance from B{C{pov}} to C{(x, y, z)} along the B{C{los}},

1228 all in C{meter}, conventionally.

1229

1230 @raise TriaxialError: Null B{C{pov}} or B{C{los}}, or B{C{pov}} is inside the

1231 tri-/biaxial or B{C{los}} points outside the tri-/biaxial

1232 or points in an opposite direction.

1233

1234 @raise TypeError: Invalid B{C{pov}} or B{C{los}}.

1235

1236 @see: Function L{pygeodesy.hartzell}, L{pygeodesy.tyr3d} for B{C{los}} and

1237 U{I{Satellite Line-of-Sight Intersection with Earth}<https://StephenHartzell.

1238 Medium.com/satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6>}.

1239 '''

1240 if isinstance(tri_biax, Triaxial_):

1241 T = tri_biax

1242 else:

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

1244 _spherical_datum(tri_biax, name=hartzell4.__name__)

1245 T = D.ellipsoid._triaxial

1246

1247 try:

1248 v, h = _hartzell3d2(pov, los, T)

1249 except Exception as x:

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

1251 return Vector4Tuple(v.x, v.y, v.z, h, name=name or hartzell4.__name__)

1252

1253

1254 def _hypot21(x, y, z=0):

1255 '''(INTERNAL) Compute M{x**2 + y**2 + z**2 - 1} with C{max(fabs(x),

1256 fabs(y), fabs(z))} rarely greater than 1.0.

1257 '''

1258 return fsumf_(_1_0, x**2, y**2, z**2, _N_2_0) if z else \

1259 fsumf_(_1_0, x**2, y**2, _N_2_0)

1260

1261

1262 def _normalTo4(x, y, a, b, eps=EPS):

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

1264

1265 @see: Function C{pygeodesy.ellipsoids._normalTo3} and I{Eberly}'s U{Distance

1266 from a Point to ... an Ellipsoid ...<https://www.GeometricTools.com/

1267 Documentation/DistancePointEllipseEllipsoid.pdf>}.

1268 '''

1269 if a < b:

1270 b, a, d, i = _normalTo4(y, x, b, a, eps=eps)

1271 return a, b, d, i

1272

1273 if not (b > 0 and isfinite(a)):

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

1275

1276 i = None

1277 if y:

1278 if x:

1279 u = fabs(x / a)

1280 v = fabs(y / b)

1281 g = _hypot21(u, v)

1282 if g:

1283 r = (a / b)**2

1284 t, i = _rootXd(r, 0, u, 0, v, g, eps)

1285 a = x / (t / r + _1_0)

1286 b = y / (t + _1_0)

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

1288 else: # on the ellipse

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

1290 else: # x == 0

1291 if y < 0:

1292 b = -b

1293 a, d = x, fabs(y - b)

1294

1295 else: # y == 0

1296 n = a * x

1297 d = (a + b) * (a - b)

1298 if d > fabs(n): # PYCHOK no cover

1299 r = n / d

1300 a *= r

1301 b *= sqrt(_1_0 - r**2)

1302 d = hypot(x - a, b)

1303 else:

1304 if x < 0:

1305 a = -a

1306 b, d = y, fabs(x - a)

1307 return a, b, d, i

1308

1309

1310 def _normalTo5(x, y, z, Tun, eps=EPS): # MCCABE 19

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

1312

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

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

1315 '''

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

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

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

1319 a, b, c, d, i = _normalTo5(*t, eps=eps)

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

1321

1322 if not (isfinite(a) and c > 0):

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

1324

1325 if eps > 0:

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

1327 else: # no validation

1328 val, eps = 0, -eps

1329

1330 i = None

1331 if z:

1332 if y:

1333 if x:

1334 u = fabs(x / a)

1335 v = fabs(y / b)

1336 w = fabs(z / c)

1337 g = _hypot21(u, v, w)

1338 if g:

1339 r = T._1e2ac # (c / a)**2

1340 s = T._1e2bc # (c / b)**2

1341 t, i = _rootXd(_1_0 / r, _1_0 / s, u, v, w, g, eps)

1342 a = x / (t * r + _1_0)

1343 b = y / (t * s + _1_0)

1344 c = z / (t + _1_0)

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

1346 else: # on the ellipsoid

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

1348 else: # x == 0

1349 a = x # 0

1350 b, c, d, i = _normalTo4(y, z, b, c, eps=eps)

1351 elif x: # y == 0

1352 b = y # 0

1353 a, c, d, i = _normalTo4(x, z, a, c, eps=eps)

1354 else: # x == y == 0

1355 if z < 0:

1356 c = -c

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

1358

1359 else: # z == 0

1360 t = False

1361 n = a * x

1362 d = T._a2c2 # (a + c) * (a - c)

1363 if d > fabs(n):

1364 u = n / d

1365 n = b * y

1366 d = T._b2c2 # (b + c) * (b - c)

1367 if d > fabs(n):

1368 v = n / d

1369 n = _hypot21(u, v)

1370 if n < 0:

1371 a *= u

1372 b *= v

1373 c *= sqrt(-n)

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

1375 t = True

1376 if not t:

1377 c = z # 0

1378 a, b, d, i = _normalTo4(x, y, a, b, eps=eps)

1379

1380 if val > 0: # validate

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

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

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

1384 sideOf=e, eps=val)

1385 if d: # angle of delta and normal vector

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

1387 if m.euclid > val:

1388 m = m.unit()

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

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

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

1392 raise _ValueError(n=n, m=m,

1393 dot=e, eps=val)

1394 return a, b, c, d, i

1395

1396

1397 def _otherV3d_(x_xyz, y, z, **name):

1398 '''(INTERNAL) Get a Vector3d from C{x_xyz}, C{y} and C{z}.

1399 '''

1400 return Vector3d(x_xyz, y, z, **name) if isscalar(x_xyz) else \

1401 _otherV3d(x_xyz=x_xyz)

1402

1403

1404 def _rootXd(r, s, u, v, w, g, eps):

1405 '''(INTERNAL) Robust 2d- or 3d-root finder:

1406 2d- if C{s == v == 0} otherwise 3d-root.

1407 '''

1408 _1, __2 = _1_0, _0_5

1409 _a, _h2 = fabs, _hypot21

1410

1411 u *= r

1412 v *= s # 0 for 2d-root

1413 t0 = w - _1

1414 t1 = _0_0 if g < 0 else _h2(u, w, v)

1415 for i in range(1, _TRIPS):

1416 e = _a(t0 - t1)

1417 if e < eps:

1418 break

1419 t = (t0 + t1) * __2

1420 if t in (t0, t1):

1421 break

1422 g = _h2(u / (t + r), w / (t + _1),

1423 (v / (t + s)) if v else 0)

1424 if g > 0:

1425 t0 = t

1426 elif g < 0:

1427 t1 = t

1428 else:

1429 break

1430 else: # PYCHOK no cover

1431 t = Fmt.no_convergence(e, eps)

1432 raise _ValueError(t, txt=_rootXd.__name__)

1433 return t, i

1434

1435

1436 def _sideOf(xyz, abc, eps=EPS): # in .formy

1437 '''(INTERNAL) Helper for C{_hartzell3d2}, M{.sideOf} and M{.reverseCartesian}.

1438

1439 @return: M{sum((x / a)**2 for x, a in zip(xyz, abc)) - 1} or C{INT0},

1440 '''

1441 s = _hypot21(*((x / a) for x, a in _zip(xyz, abc) if a)) # strict=True

1442 return s if fabs(s) > eps else INT0

1443

1444

1445 if __name__ == '__main__':

1446

1447 from pygeodesy import printf

1448 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

1449

1450 # __doc__ of this file, force all into registery

1451 t = [NN] + Triaxials.toRepr(all=True, asorted=True).split(_NL_)

1452 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))

1453

1454# **) MIT License

1455#

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

1457#

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

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

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

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

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

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

1464#

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

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

1467#

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

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

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

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

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

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

1474# OTHER DEALINGS IN THE SOFTWARE.