1

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

3

4 u'''Ellipsoidal and spherical earth models.

5

6 Classes L{a_f2Tuple}, L{Ellipsoid} and L{Ellipsoid2}, an L{Ellipsoids} registry and

7 2 dozen functions to convert I{equatorial} radius, I{polar} radius, I{eccentricities},

8 I{flattenings} and I{inverse flattening}.

9

10 See module L{datums} for L{Datum} and L{Transform} information and other details.

11

12 Following is the list of predefined L{Ellipsoid}s, all instantiated lazily.

13

14@var Ellipsoids.Airy1830: Ellipsoid(name='Airy1830', a=6377563.396, f=0.00334085, f_=299.3249646, b=6356256.90923729)

15@var Ellipsoids.AiryModified: Ellipsoid(name='AiryModified', a=6377340.189, f=0.00334085, f_=299.3249646, b=6356034.44793853)

16@var Ellipsoids.ATS1977: Ellipsoid(name='ATS1977', a=6378135, f=0.00335281, f_=298.257, b=6356750.30492159)

17@var Ellipsoids.Australia1966: Ellipsoid(name='Australia1966', a=6378160, f=0.00335289, f_=298.25, b=6356774.71919531)

18@var Ellipsoids.Bessel1841: Ellipsoid(name='Bessel1841', a=6377397.155, f=0.00334277, f_=299.1528128, b=6356078.962818)

19@var Ellipsoids.BesselModified: Ellipsoid(name='BesselModified', a=6377492.018, f=0.00334277, f_=299.1528128, b=6356173.5087127)

20@var Ellipsoids.CGCS2000: Ellipsoid(name='CGCS2000', a=6378137, f=0.00335281, f_=298.2572221, b=6356752.31414036)

21@var Ellipsoids.Clarke1866: Ellipsoid(name='Clarke1866', a=6378206.4, f=0.00339008, f_=294.97869821, b=6356583.8)

22@var Ellipsoids.Clarke1880: Ellipsoid(name='Clarke1880', a=6378249.145, f=0.00340756, f_=293.465, b=6356514.86954978)

23@var Ellipsoids.Clarke1880IGN: Ellipsoid(name='Clarke1880IGN', a=6378249.2, f=0.00340755, f_=293.46602129, b=6356515)

24@var Ellipsoids.Clarke1880Mod: Ellipsoid(name='Clarke1880Mod', a=6378249.145, f=0.00340755, f_=293.46630766, b=6356514.96639549)

25@var Ellipsoids.CPM1799: Ellipsoid(name='CPM1799', a=6375738.7, f=0.00299052, f_=334.39, b=6356671.92557493)

26@var Ellipsoids.Delambre1810: Ellipsoid(name='Delambre1810', a=6376428, f=0.00321027, f_=311.5, b=6355957.92616372)

27@var Ellipsoids.Engelis1985: Ellipsoid(name='Engelis1985', a=6378136.05, f=0.00335282, f_=298.2566, b=6356751.32272154)

28@var Ellipsoids.Everest1969: Ellipsoid(name='Everest1969', a=6377295.664, f=0.00332445, f_=300.8017, b=6356094.667915)

29@var Ellipsoids.Everest1975: Ellipsoid(name='Everest1975', a=6377299.151, f=0.00332445, f_=300.8017255, b=6356098.14512013)

30@var Ellipsoids.Fisher1968: Ellipsoid(name='Fisher1968', a=6378150, f=0.00335233, f_=298.3, b=6356768.33724438)

31@var Ellipsoids.GEM10C: Ellipsoid(name='GEM10C', a=6378137, f=0.00335281, f_=298.2572236, b=6356752.31424783)

32@var Ellipsoids.GPES: Ellipsoid(name='GPES', a=6378135, f=0, f_=0, b=6378135)

33@var Ellipsoids.GRS67: Ellipsoid(name='GRS67', a=6378160, f=0.00335292, f_=298.24716743, b=6356774.51609071)

34@var Ellipsoids.GRS80: Ellipsoid(name='GRS80', a=6378137, f=0.00335281, f_=298.2572221, b=6356752.31414035)

35@var Ellipsoids.Helmert1906: Ellipsoid(name='Helmert1906', a=6378200, f=0.00335233, f_=298.3, b=6356818.16962789)

36@var Ellipsoids.IAU76: Ellipsoid(name='IAU76', a=6378140, f=0.00335281, f_=298.257, b=6356755.28815753)

37@var Ellipsoids.IERS1989: Ellipsoid(name='IERS1989', a=6378136, f=0.00335281, f_=298.257, b=6356751.30156878)

38@var Ellipsoids.IERS1992TOPEX: Ellipsoid(name='IERS1992TOPEX', a=6378136.3, f=0.00335281, f_=298.25722356, b=6356751.61659215)

39@var Ellipsoids.IERS2003: Ellipsoid(name='IERS2003', a=6378136.6, f=0.00335282, f_=298.25642, b=6356751.85797165)

40@var Ellipsoids.Intl1924: Ellipsoid(name='Intl1924', a=6378388, f=0.003367, f_=297, b=6356911.94612795)

41@var Ellipsoids.Intl1967: Ellipsoid(name='Intl1967', a=6378157.5, f=0.0033529, f_=298.24961539, b=6356772.2)

42@var Ellipsoids.Krassovski1940: Ellipsoid(name='Krassovski1940', a=6378245, f=0.00335233, f_=298.3, b=6356863.01877305)

43@var Ellipsoids.Krassowsky1940: Ellipsoid(name='Krassowsky1940', a=6378245, f=0.00335233, f_=298.3, b=6356863.01877305)

44@var Ellipsoids.Maupertuis1738: Ellipsoid(name='Maupertuis1738', a=6397300, f=0.0052356, f_=191, b=6363806.28272251)

45@var Ellipsoids.Mercury1960: Ellipsoid(name='Mercury1960', a=6378166, f=0.00335233, f_=298.3, b=6356784.28360711)

46@var Ellipsoids.Mercury1968Mod: Ellipsoid(name='Mercury1968Mod', a=6378150, f=0.00335233, f_=298.3, b=6356768.33724438)

47@var Ellipsoids.NWL1965: Ellipsoid(name='NWL1965', a=6378145, f=0.00335289, f_=298.25, b=6356759.76948868)

48@var Ellipsoids.OSU86F: Ellipsoid(name='OSU86F', a=6378136.2, f=0.00335281, f_=298.2572236, b=6356751.51693008)

49@var Ellipsoids.OSU91A: Ellipsoid(name='OSU91A', a=6378136.3, f=0.00335281, f_=298.2572236, b=6356751.6165948)

50@var Ellipsoids.Plessis1817: Ellipsoid(name='Plessis1817', a=6376523, f=0.00324002, f_=308.64, b=6355862.93325557)

51@var Ellipsoids.PZ90: Ellipsoid(name='PZ90', a=6378136, f=0.0033528, f_=298.2578393, b=6356751.36174571)

52@var Ellipsoids.SGS85: Ellipsoid(name='SGS85', a=6378136, f=0.00335281, f_=298.257, b=6356751.30156878)

53@var Ellipsoids.SoAmerican1969: Ellipsoid(name='SoAmerican1969', a=6378160, f=0.00335289, f_=298.25, b=6356774.71919531)

54@var Ellipsoids.Sphere: Ellipsoid(name='Sphere', a=6371008.771415, f=0, f_=0, b=6371008.771415)

55@var Ellipsoids.SphereAuthalic: Ellipsoid(name='SphereAuthalic', a=6371000, f=0, f_=0, b=6371000)

56@var Ellipsoids.SpherePopular: Ellipsoid(name='SpherePopular', a=6378137, f=0, f_=0, b=6378137)

57@var Ellipsoids.Struve1860: Ellipsoid(name='Struve1860', a=6378298.3, f=0.00339294, f_=294.73, b=6356657.14266956)

58@var Ellipsoids.WGS60: Ellipsoid(name='WGS60', a=6378165, f=0.00335233, f_=298.3, b=6356783.28695944)

59@var Ellipsoids.WGS66: Ellipsoid(name='WGS66', a=6378145, f=0.00335289, f_=298.25, b=6356759.76948868)

60@var Ellipsoids.WGS72: Ellipsoid(name='WGS72', a=6378135, f=0.00335278, f_=298.26, b=6356750.52001609)

61@var Ellipsoids.WGS84: Ellipsoid(name='WGS84', a=6378137, f=0.00335281, f_=298.25722356, b=6356752.31424518)

62@var Ellipsoids.WGS84_NGS: Ellipsoid(name='WGS84_NGS', a=6378137, f=0.00335281, f_=298.2572221, b=6356752.31414035)

63'''

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

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

66

67 from pygeodesy.basics import copysign0, isbool, isint

68 from pygeodesy.constants import EPS, EPS0, EPS02, EPS1, INF, NINF, PI4, PI_2, PI_3, R_M, R_MA, R_FM, \

69 _EPSqrt, _EPStol as _TOL, _floatuple as _T, _isfinite, _SQRT2_2, \

70 _0_0s, _0_0, _0_5, _1_0, _1_EPS, _2_0, _4_0, _90_0, \

71 _0_25, _3_0 # PYCHOK used!

72 from pygeodesy.errors import _AssertionError, IntersectionError, _ValueError, _xattr, _xkwds_not

73 from pygeodesy.fmath import cbrt, cbrt2, fdot, Fhorner, fpowers, hypot, hypot_, \

74 hypot1, hypot2, sqrt3, Fsum

75# from pygeodesy.fsums import Fsum # from .fmath

76 from pygeodesy.interns import NN, _a_, _Airy1830_, _AiryModified_, _b_, _Bessel1841_, _beta_, \

77 _Clarke1866_, _Clarke1880IGN_, _DOT_, _f_, _GRS80_, _height_, \

78 _Intl1924_, _incompatible_, _invalid_, _Krassovski1940_, \

79 _Krassowsky1940_, _lat_, _meridional_, _negative_, _not_finite_, \

80 _prime_vertical_, _radius_, _Sphere_, _SPACE_, _vs_, _WGS72_, _WGS84_

81# from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS # from .named

82 from pygeodesy.named import _lazyNamedEnumItem as _lazy, _name__, _NamedEnum, \

83 _NamedEnumItem, _NamedTuple, _Pass, _ALL_LAZY, _MODS

84 from pygeodesy.namedTuples import Distance2Tuple, Vector3Tuple, Vector4Tuple

85 from pygeodesy.props import deprecated_Property_RO, Property_RO, property_doc_, \

86 deprecated_property_RO, property_RO

87 from pygeodesy.streprs import Fmt, fstr, instr, strs, unstr

88# from pygeodesy.triaxials import _hartzell3 # _MODS

89 from pygeodesy.units import Bearing_, Distance, Float, Float_, Height, Lamd, Lat, Meter, \

90 Meter2, Meter3, Phi, Phid, Radius, Radius_, Scalar

91 from pygeodesy.utily import atan1, atan1d, atan2b, degrees90, m2radians, radians2m, sincos2d

92

93 from math import asinh, atan, atanh, cos, degrees, exp, fabs, radians, sin, sinh, sqrt, tan

94

95__all__ = _ALL_LAZY.ellipsoids

96__version__ = '24.06.24'

97

98_f_0_0 = Float(f =_0_0) # zero flattening

99_f__0_0 = Float(f_=_0_0) # zero inverse flattening

100# see U{WGS84_f<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Constants.html>}

101_f__WGS84 = Float(f_=_1_0 / (1000000000 / 298257223563)) # 298.25722356299997 vs 298.257223563

102

103

104 def _aux(lat, inverse, auxLat, clip=90):

105 '''Return a named auxiliary latitude in C{degrees}.

106 '''

107 return Lat(lat, clip=clip, name=_lat_ if inverse else auxLat.__name__)

108

109

110 def _s2_c2(phi):

111 '''(INTERNAL) Return 2-tuple C{(sin(B{phi})**2, cos(B{phi})**2)}.

112 '''

113 if phi:

114 s2 = sin(phi)**2

115 if s2 > EPS:

116 c2 = _1_0 - s2

117 if c2 > EPS:

118 if c2 < EPS1:

119 return s2, c2

120 else:

121 return _1_0, _0_0 # phi == PI_2

122 return _0_0, _1_0 # phi == 0

123

124

125 class a_f2Tuple(_NamedTuple):

126 '''2-Tuple C{(a, f)} specifying an ellipsoid by I{equatorial}

127 radius C{a} in C{meter} and scalar I{flattening} C{f}.

128

129 @see: Class L{Ellipsoid2}.

130 '''

131 _Names_ = (_a_, _f_) # name 'f' not 'f_'

132 _Units_ = (_Pass, _Pass)

133

134 def __new__(cls, a, f, **name):

135 '''New L{a_f2Tuple} ellipsoid specification.

136

137 @arg a: Equatorial radius (C{scalar} > 0).

138 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

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

140

141 @return: An L{a_f2Tuple}C{(a, f)} instance.

142

143 @raise UnitError: Invalid B{C{a}} or B{C{f}}.

144

145 @note: C{abs(B{f}) < EPS} is forced to C{B{f}=0}, I{spherical}.

146 Negative C{B{f}} produces a I{prolate} ellipsoid.

147 '''

148 a = Radius_(a=a) # low=EPS, high=None

149 f = Float_( f=f, low=None, high=EPS1)

150 if fabs(f) < EPS: # force spherical

151 f = _f_0_0

152 return _NamedTuple.__new__(cls, a, f, **name)

153

154 @Property_RO

155 def b(self):

156 '''Get the I{polar} radius (C{meter}), M{a * (1 - f)}.

157 '''

158 return a_f2b(self.a, self.f) # PYCHOK .a and .f

159

160 def ellipsoid(self, **name):

161 '''Return an L{Ellipsoid} for this 2-tuple C{(a, f)}.

162

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

164

165 @raise NameError: A registered C{ellipsoid} with the

166 same B{C{name}} already exists.

167 '''

168 return Ellipsoid(self.a, f=self.f, name=self._name__(name)) # PYCHOK .a and .f

169

170 @Property_RO

171 def f_(self):

172 '''Get the I{inverse} flattening (C{scalar}), M{1 / f} == M{a / (a - b)}.

173 '''

174 return f2f_(self.f) # PYCHOK .f

175

176

177 class Circle4Tuple(_NamedTuple):

178 '''4-Tuple C{(radius, height, lat, beta)} of the C{radius} and C{height},

179 both conventionally in C{meter} of a parallel I{circle of latitude} at

180 (geodetic) latitude C{lat} and the I{parametric (or reduced) auxiliary

181 latitude} C{beta}, both in C{degrees90}.

182

183 The C{height} is the (signed) distance along the z-axis between the

184 parallel and the equator. At near-polar C{lat}s, the C{radius} is C{0},

185 the C{height} is the ellipsoid's (signed) polar radius and C{beta}

186 equals C{lat}.

187 '''

188 _Names_ = (_radius_, _height_, _lat_, _beta_)

189 _Units_ = ( Radius, Height, Lat, Lat)

190

191

192 class Curvature2Tuple(_NamedTuple):

193 '''2-Tuple C{(meridional, prime_vertical)} of radii of curvature, both in

194 C{meter}, conventionally.

195 '''

196 _Names_ = (_meridional_, _prime_vertical_)

197 _Units_ = ( Meter, Meter)

198

199 @property_RO

200 def transverse(self):

201 '''Get this I{prime_vertical}, aka I{transverse} radius of curvature.

202 '''

203 return self.prime_vertical

204

205

206 class Ellipsoid(_NamedEnumItem):

207 '''Ellipsoid with I{equatorial} and I{polar} radii, I{flattening}, I{inverse

208 flattening} and other, often used, I{cached} attributes, supporting

209 I{oblate} and I{prolate} ellipsoidal and I{spherical} earth models.

210 '''

211 _a = 0 # equatorial radius, semi-axis (C{meter})

212 _b = 0 # polar radius, semi-axis (C{meter}): a * (f - 1) / f

213 _f = 0 # (1st) flattening: (a - b) / a

214 _f_ = 0 # inverse flattening: 1 / f = a / (a - b)

215

216 _geodsolve = NN # means, use PYGEODESY_GEODSOLVE

217 _KsOrder = 8 # Krüger series order (4, 6 or 8)

218 _rhumbsolve = NN # means, use PYGEODESY_RHUMBSOLVE

219

220 def __init__(self, a, b=None, f_=None, f=None, **name):

221 '''New L{Ellipsoid} from the I{equatorial} radius I{and} either

222 the I{polar} radius or I{inverse flattening} or I{flattening}.

223

224 @arg a: Equatorial radius, semi-axis (C{meter}).

225 @arg b: Optional polar radius, semi-axis (C{meter}).

226 @arg f_: Inverse flattening: M{a / (a - b)} (C{float} >>> 1.0).

227 @arg f: Flattening: M{(a - b) / a} (C{scalar}, near zero for

228 spherical).

229 @kwarg name: Optional, unique C{B{name}=NN} (C{str}).

230

231 @raise NameError: Ellipsoid with the same B{C{name}} already exists.

232

233 @raise ValueError: Invalid B{C{a}}, B{C{b}}, B{C{f_}} or B{C{f}} or

234 B{C{f_}} and B{C{f}} are incompatible.

235

236 @note: M{abs(f_) > 1 / EPS} or M{abs(1 / f_) < EPS} is forced

237 to M{1 / f_ = 0}, spherical.

238 '''

239 ff_ = f, f_ # assertion below

240 n = _name__(**name) if name else NN

241 try:

242 a = Radius_(a=a) # low=EPS

243 if not _isfinite(a):

244 raise ValueError(_SPACE_(_a_, _not_finite_))

245

246 if b: # not in (_0_0, None)

247 b = Radius_(b=b) # low=EPS

248 f = a_b2f(a, b) if f is None else Float(f=f)

249 f_ = f2f_(f) if f_ is None else Float(f_=f_)

250 elif f is not None:

251 f = Float(f=f)

252 b = a_f2b(a, f)

253 f_ = f2f_(f) if f_ is None else Float(f_=f_)

254 elif f_:

255 f_ = Float(f_=f_)

256 b = a_f_2b(a, f_) # a * (f_ - 1) / f_

257 f = f_2f(f_)

258 else: # only a, spherical

259 f_ = f = 0

260 b = a # superfluous

261

262 if not f < _1_0: # sanity check, see .ecef.Ecef.__init__

263 raise ValueError(_SPACE_(_f_, _invalid_))

264 if not _isfinite(b):

265 raise ValueError(_SPACE_(_b_, _not_finite_))

266

267 if fabs(f) < EPS or a == b or not f_: # spherical

268 b = a

269 f = _f_0_0

270 f_ = _f__0_0

271

272 except (TypeError, ValueError) as x:

273 d = _xkwds_not(None, b=b, f_=f_, f=f)

274 t = instr(self, a=a, name=n, **d)

275 raise _ValueError(t, cause=x)

276

277 self._a = a

278 self._b = b

279 self._f = f

280 self._f_ = f_

281

282 self._register(Ellipsoids, n)

283

284 if f and f_: # see .test/testEllipsoidal.py

285 d = dict(eps=_TOL)

286 if None in ff_: # both f_ and f given

287 d.update(Error=_ValueError, txt=_incompatible_)

288 self._assert(_1_0 / f, f_=f_, **d)

289 self._assert(_1_0 / f_, f =f, **d)

290 self._assert(self.b2_a2, e21=self.e21, eps=EPS)

291

292 def __eq__(self, other):

293 '''Compare this and an other ellipsoid.

294

295 @arg other: The other ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).

296

297 @return: C{True} if equal, C{False} otherwise.

298 '''

299 return self is other or (isinstance(other, Ellipsoid) and

300 self.a == other.a and

301 (self.f == other.f or self.b == other.b))

302

303 def __hash__(self):

304 return self._hash # memoized

305

306 @Property_RO

307 def a(self):

308 '''Get the I{equatorial} radius, semi-axis (C{meter}).

309 '''

310 return self._a

311

312 equatoradius = a # = Requatorial

313

314 @Property_RO

315 def a2(self):

316 '''Get the I{equatorial} radius I{squared} (C{meter} I{squared}), M{a**2}.

317 '''

318 return Meter2(a2=self.a**2)

319

320 @Property_RO

321 def a2_(self):

322 '''Get the inverse of the I{equatorial} radius I{squared} (C{meter} I{squared}), M{1 / a**2}.

323 '''

324 return Float(a2_=_1_0 / self.a2)

325

326 @Property_RO

327 def a_b(self):

328 '''Get the ratio I{equatorial} over I{polar} radius (C{float}), M{a / b} == M{1 / (1 - f)}.

329 '''

330 return Float(a_b=self.a / self.b if self.f else _1_0)

331

332 @Property_RO

333 def a2_b(self):

334 '''Get the I{polar} meridional (or polar) radius of curvature (C{meter}), M{a**2 / b}.

335

336 @see: U{Radii of Curvature

337 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}

338 and U{Moritz, H. (1980), Geodetic Reference System 1980

339 <https://WikiPedia.org/wiki/Earth_radius#cite_note-Moritz-2>}.

340

341 @note: Symbol C{c} is used by IUGG and IERS for the U{polar radius of curvature

342 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}, see L{c2}

343 and L{R2} or L{Rauthalic}.

344 '''

345 return Radius(a2_b=self.a2 / self.b if self.f else self.a) # = rocPolar

346

347 @Property_RO

348 def a2_b2(self):

349 '''Get the ratio I{equatorial} over I{polar} radius I{squared} (C{float}),

350 M{(a / b)**2} == M{1 / (1 - e**2)} == M{1 / (1 - e2)} == M{1 / e21}.

351 '''

352 return Float(a2_b2=self.a_b**2 if self.f else _1_0)

353

354 @Property_RO

355 def a_f(self):

356 '''Get the I{equatorial} radius and I{flattening} (L{a_f2Tuple}), see method C{toEllipsoid2}.

357 '''

358 return a_f2Tuple(self.a, self.f, name=self.name)

359

360 @Property_RO

361 def A(self):

362 '''Get the UTM I{meridional (or rectifying)} radius (C{meter}).

363

364 @see: I{Meridian arc unit} U{Q<https://StudyLib.net/doc/7443565/>}.

365 '''

366 A, n = self.a, self.n

367 if n:

368 d = (n + _1_0) * 1048576 / A

369 if d: # use 6 n**2 terms, half-way between the _KsOrder's 4, 6, 8

370 # <https://GeographicLib.SourceForge.io/C++/doc/tmseries30.html>

371 # <https://GeographicLib.SourceForge.io/C++/doc/transversemercator.html> and

372 # <https://www.MyGeodesy.id.AU/documents/Karney-Krueger%20equations.pdf> (3)

373 # A *= fhorner(n**2, 1048576, 262144, 16384, 4096, 1600, 784, 441) / 1048576) / (1 + n)

374 A = Radius(A=Fhorner(n**2, 1048576, 262144, 16384, 4096, 1600, 784, 441).fover(d))

375 return A

376

377 @Property_RO

378 def _albersCyl(self):

379 '''(INTERNAL) Helper for C{auxAuthalic}.

380 '''

381 return _MODS.albers.AlbersEqualAreaCylindrical(datum=self, name=self.name)

382

383 @Property_RO

384 def AlphaKs(self):

385 '''Get the I{Krüger} U{Alpha series coefficients<https://GeographicLib.SourceForge.io/C++/doc/tmseries30.html>} (C{KsOrder}C{-tuple}).

386 '''

387 return self._Kseries( # XXX int/int quotients may require from __future__ import division as _; del _ # PYCHOK semicolon

388 # n n**2 n**3 n**4 n**5 n**6 n**7 n**8

389 _T(1/2, -2/3, 5/16, 41/180, -127/288, 7891/37800, 72161/387072, -18975107/50803200),

390 _T(13/48, -3/5, 557/1440, 281/630, -1983433/1935360, 13769/28800, 148003883/174182400), # PYCHOK unaligned

391 _T(61/240, -103/140, 15061/26880, 167603/181440, -67102379/29030400, 79682431/79833600), # PYCHOK unaligned

392 _T(49561/161280, -179/168, 6601661/7257600, 97445/49896, -40176129013/7664025600), # PYCHOK unaligned

393 _T(34729/80640, -3418889/1995840, 14644087/9123840, 2605413599/622702080), # PYCHOK unaligned

394 _T(212378941/319334400, -30705481/10378368, 175214326799/58118860800), # PYCHOK unaligned

395 _T(1522256789/1383782400, -16759934899/3113510400), # PYCHOK unaligned

396 _T(1424729850961/743921418240)) # PYCHOK unaligned

397

398 @Property_RO

399 def area(self):

400 '''Get the ellipsoid's surface area (C{meter} I{squared}), M{4 * PI * c2}.

401

402 @see: Properties L{areax}, L{c2} and L{R2} and functions

403 L{ellipsoidalExact.areaOf} and L{ellipsoidalKarney.areaOf}.

404 '''

405 return Meter2(area=self.c2 * PI4)

406

407 @Property_RO

408 def areax(self):

409 '''Get the ellipsoid's surface area (C{meter} I{squared}), M{4 * PI * c2x}, more

410 accurate for very I{oblate} ellipsoids.

411

412 @see: Properties L{area}, L{c2x} and L{R2x}, class L{GeodesicExact} and

413 functions L{ellipsoidalExact.areaOf} and L{ellipsoidalKarney.areaOf}.

414 '''

415 return Meter2(areax=self.c2x * PI4)

416

417 def _assert(self, val, eps=_TOL, f0=_0_0, Error=_AssertionError, txt=NN, **name_value):

418 '''(INTERNAL) Assert a C{name=value} vs C{val}.

419 '''

420 for n, v in name_value.items():

421 if fabs(v - val) > eps: # PYCHOK no cover

422 t = (v, _vs_, val)

423 t = _SPACE_.join(strs(t, prec=12, fmt=Fmt.g))

424 t = Fmt.EQUAL(self._DOT_(n), t)

425 raise Error(t, txt=txt or Fmt.exceeds_eps(eps))

426 return Float(v if self.f else f0, name=n)

427 raise Error(unstr(self._DOT_(self._assert.__name__), val,

428 eps=eps, f0=f0, **name_value))

429

430 def auxAuthalic(self, lat, inverse=False):

431 '''Compute the I{authalic} auxiliary latitude or the I{inverse} thereof.

432

433 @arg lat: The geodetic (or I{authalic}) latitude (C{degrees90}).

434 @kwarg inverse: If C{True}, B{C{lat}} is the I{authalic} and

435 return the geodetic latitude (C{bool}).

436

437 @return: The I{authalic} (or geodetic) latitude in C{degrees90}.

438

439 @see: U{Inverse-/AuthalicLatitude<https://GeographicLib.SourceForge.io/

440 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Authalic latitude

441 <https://WikiPedia.org/wiki/Latitude#Authalic_latitude>}, and

442 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, p 16.

443 '''

444 if self.f:

445 f = self._albersCyl._tanf if inverse else self._albersCyl._txif # PYCHOK attr

446 lat = atan1d(f(tan(Phid(lat)))) # PYCHOK attr

447 return _aux(lat, inverse, Ellipsoid.auxAuthalic)

448

449 def auxConformal(self, lat, inverse=False):

450 '''Compute the I{conformal} auxiliary latitude or the I{inverse} thereof.

451

452 @arg lat: The geodetic (or I{conformal}) latitude (C{degrees90}).

453 @kwarg inverse: If C{True}, B{C{lat}} is the I{conformal} and

454 return the geodetic latitude (C{bool}).

455

456 @return: The I{conformal} (or geodetic) latitude in C{degrees90}.

457

458 @see: U{Inverse-/ConformalLatitude<https://GeographicLib.SourceForge.io/

459 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Conformal latitude

460 <https://WikiPedia.org/wiki/Latitude#Conformal_latitude>}, and

461 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 15-16.

462 '''

463 if self.f:

464 f = self.es_tauf if inverse else self.es_taupf # PYCHOK attr

465 lat = atan1d(f(tan(Phid(lat)))) # PYCHOK attr

466 return _aux(lat, inverse, Ellipsoid.auxConformal)

467

468 def auxGeocentric(self, lat, inverse=False):

469 '''Compute the I{geocentric} auxiliary latitude or the I{inverse} thereof.

470

471 @arg lat: The geodetic (or I{geocentric}) latitude (C{degrees90}).

472 @kwarg inverse: If C{True}, B{C{lat}} is the geocentric and

473 return the I{geocentric} latitude (C{bool}).

474

475 @return: The I{geocentric} (or geodetic) latitude in C{degrees90}.

476

477 @see: U{Inverse-/GeocentricLatitude<https://GeographicLib.SourceForge.io/

478 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Geocentric latitude

479 <https://WikiPedia.org/wiki/Latitude#Geocentric_latitude>}, and

480 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 17-18.

481 '''

482 if self.f:

483 f = self.a2_b2 if inverse else self.b2_a2

484 lat = atan1d(f * tan(Phid(lat)))

485 return _aux(lat, inverse, Ellipsoid.auxGeocentric)

486

487 def auxIsometric(self, lat, inverse=False):

488 '''Compute the I{isometric} auxiliary latitude or the I{inverse} thereof.

489

490 @arg lat: The geodetic (or I{isometric}) latitude (C{degrees}).

491 @kwarg inverse: If C{True}, B{C{lat}} is the I{isometric} and

492 return the geodetic latitude (C{bool}).

493

494 @return: The I{isometric} (or geodetic) latitude in C{degrees}.

495

496 @note: The I{isometric} latitude for geodetic C{+/-90} is far

497 outside the C{[-90..+90]} range but the inverse

498 thereof is the original geodetic latitude.

499

500 @see: U{Inverse-/IsometricLatitude<https://GeographicLib.SourceForge.io/

501 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Isometric latitude

502 <https://WikiPedia.org/wiki/Latitude#Isometric_latitude>}, and

503 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 15-16.

504 '''

505 if self.f:

506 r = Phid(lat, clip=0)

507 lat = degrees(atan1(self.es_tauf(sinh(r))) if inverse else

508 asinh(self.es_taupf(tan(r))))

509 # clip=0, since auxIsometric(+/-90) is far outside [-90..+90]

510 return _aux(lat, inverse, Ellipsoid.auxIsometric, clip=0)

511

512 def auxParametric(self, lat, inverse=False):

513 '''Compute the I{parametric} auxiliary latitude or the I{inverse} thereof.

514

515 @arg lat: The geodetic (or I{parametric}) latitude (C{degrees90}).

516 @kwarg inverse: If C{True}, B{C{lat}} is the I{parametric} and

517 return the geodetic latitude (C{bool}).

518

519 @return: The I{parametric} (or geodetic) latitude in C{degrees90}.

520

521 @see: U{Inverse-/ParametricLatitude<https://GeographicLib.SourceForge.io/

522 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Parametric latitude

523 <https://WikiPedia.org/wiki/Latitude#Parametric_(or_reduced)_latitude>},

524 and U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, p 18.

525 '''

526 if self.f:

527 lat = self._beta(Lat(lat), inverse=inverse)

528 return _aux(lat, inverse, Ellipsoid.auxParametric)

529

530 auxReduced = auxParametric # synonymous

531

532 def auxRectifying(self, lat, inverse=False):

533 '''Compute the I{rectifying} auxiliary latitude or the I{inverse} thereof.

534

535 @arg lat: The geodetic (or I{rectifying}) latitude (C{degrees90}).

536 @kwarg inverse: If C{True}, B{C{lat}} is the I{rectifying} and

537 return the geodetic latitude (C{bool}).

538

539 @return: The I{rectifying} (or geodetic) latitude in C{degrees90}.

540

541 @see: U{Inverse-/RectifyingLatitude<https://GeographicLib.SourceForge.io/

542 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Rectifying latitude

543 <https://WikiPedia.org/wiki/Latitude#Rectifying_latitude>}, and

544 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 16-17.

545 '''

546 if self.f:

547 lat = Lat(lat)

548 if 0 < fabs(lat) < _90_0:

549 if inverse:

550 e = self._elliptic_e22

551 d = degrees90(e.fEinv(e.cE * lat / _90_0))

552 lat = self.auxParametric(d, inverse=True)

553 else:

554 lat = _90_0 * self.Llat(lat) / self.L

555 return _aux(lat, inverse, Ellipsoid.auxRectifying)

556

557 @Property_RO

558 def b(self):

559 '''Get the I{polar} radius, semi-axis (C{meter}).

560 '''

561 return self._b

562

563 polaradius = b # = Rpolar

564

565 @Property_RO

566 def b_a(self):

567 '''Get the ratio I{polar} over I{equatorial} radius (C{float}), M{b / a == f1 == 1 - f}.

568

569 @see: Property L{f1}.

570 '''

571 return self._assert(self.b / self.a, b_a=self.f1, f0=_1_0)

572

573 @Property_RO

574 def b2(self):

575 '''Get the I{polar} radius I{squared} (C{float}), M{b**2}.

576 '''

577 return Meter2(b2=self.b**2)

578

579 @Property_RO

580 def b2_a(self):

581 '''Get the I{equatorial} meridional radius of curvature (C{meter}), M{b**2 / a}, see C{rocMeridional}C{(0)}.

582

583 @see: U{Radii of Curvature<https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}.

584 '''

585 return Radius(b2_a=self.b2 / self.a if self.f else self.b)

586

587 @Property_RO

588 def b2_a2(self):

589 '''Get the ratio I{polar} over I{equatorial} radius I{squared} (C{float}), M{(b / a)**2}

590 == M{(1 - f)**2} == M{1 - e**2} == C{e21}.

591 '''

592 return Float(b2_a2=self.b_a**2 if self.f else _1_0)

593

594 def _beta(self, lat, inverse=False):

595 '''(INTERNAL) Get the I{parametric (or reduced) auxiliary latitude} or inverse thereof.

596 '''

597 s, c = sincos2d(lat) # like Karney's tand(lat)

598 s *= self.a_b if inverse else self.b_a

599 return atan1d(s, c)

600

601 @Property_RO

602 def BetaKs(self):

603 '''Get the I{Krüger} U{Beta series coefficients<https://GeographicLib.SourceForge.io/C++/doc/tmseries30.html>} (C{KsOrder}C{-tuple}).

604 '''

605 return self._Kseries( # XXX int/int quotients may require from __future__ import division as _; del _ # PYCHOK semicolon

606 # n n**2 n**3 n**4 n**5 n**6 n**7 n**8

607 _T(1/2, -2/3, 37/96, -1/360, -81/512, 96199/604800, -5406467/38707200, 7944359/67737600),

608 _T(1/48, 1/15, -437/1440, 46/105, -1118711/3870720, 51841/1209600, 24749483/348364800), # PYCHOK unaligned

609 _T(17/480, -37/840, -209/4480, 5569/90720, 9261899/58060800, -6457463/17740800), # PYCHOK unaligned

610 _T(4397/161280, -11/504, -830251/7257600, 466511/2494800, 324154477/7664025600), # PYCHOK unaligned

611 _T(4583/161280, -108847/3991680, -8005831/63866880, 22894433/124540416), # PYCHOK unaligned

612 _T(20648693/638668800, -16363163/518918400, -2204645983/12915302400), # PYCHOK unaligne

613 _T(219941297/5535129600, -497323811/12454041600), # PYCHOK unaligned

614 _T(191773887257/3719607091200)) # PYCHOK unaligned

615

616 @deprecated_Property_RO

617 def c(self): # PYCHOK no cover

618 '''DEPRECATED, use property C{R2} or C{Rauthalic}.'''

619 return self.R2

620

621 @Property_RO

622 def c2(self):

623 '''Get the I{authalic} earth radius I{squared} (C{meter} I{squared}).

624

625 @see: Properties L{c2x}, L{area}, L{R2}, L{Rauthalic}, I{Karney's} U{equation (60)

626 <https://Link.Springer.com/article/10.1007%2Fs00190-012-0578-z>} and C++ U{Ellipsoid.Area

627 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Ellipsoid.html>},

628 U{Authalic radius<https://WikiPedia.org/wiki/Earth_radius#Authalic_radius>}, U{Surface area

629 <https://WikiPedia.org/wiki/Ellipsoid>} and U{surface area

630 <https://www.Numericana.com/answer/geometry.htm#oblate>}.

631 '''

632 return self._c2f(False)

633

634 @Property_RO

635 def c2x(self):

636 '''Get the I{authalic} earth radius I{squared} (C{meter} I{squared}), more accurate for very I{oblate}

637 ellipsoids.

638

639 @see: Properties L{c2}, L{areax}, L{R2x}, L{Rauthalicx}, class L{GeodesicExact} and I{Karney}'s comments at C++

640 attribute U{GeodesicExact._c2<https://GeographicLib.SourceForge.io/C++/doc/GeodesicExact_8cpp_source.html>}.

641 '''

642 return self._c2f(True)

643

644 def _c2f(self, c2x):

645 '''(INTERNAL) Helper for C{.c2} and C{.c2x}.

646 '''

647 f, c2 = self.f, self.b2

648 if f:

649 e = self.e

650 if e > EPS0:

651 if f > 0: # .isOblate

652 c2 *= (asinh(sqrt(self.e22abs)) if c2x else atanh(e)) / e

653 elif f < 0: # .isProlate

654 c2 *= atan1(e) / e # XXX asin?

655 c2 = Meter2(c2=(self.a2 + c2) * _0_5)

656 return c2

657

658 def circle4(self, lat):

659 '''Get the equatorial or a parallel I{circle of latitude}.

660

661 @arg lat: Geodetic latitude (C{degrees90}, C{str}).

662

663 @return: A L{Circle4Tuple}C{(radius, height, lat, beta)}

664 instance.

665

666 @raise RangeError: Latitude B{C{lat}} outside valid range and

667 L{rangerrors<pygeodesy.rangerrors>} is C{True}.

668

669 @raise TypeError: Invalid B{C{lat}}.

670

671 @raise ValueError: Invalid B{C{lat}}.

672

673 @see: Definition of U{I{p} and I{z} under B{Parametric (or reduced) latitude}

674 <https://WikiPedia.org/wiki/Latitude>}, I{Karney's} C++ U{CircleRadius and CircleHeight

675 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Ellipsoid.html>}

676 and method C{Rlat}.

677 '''

678 lat = Lat(lat)

679 if lat:

680 b = lat

681 if fabs(lat) < _90_0:

682 if self.f:

683 b = self._beta(lat)

684 z, r = sincos2d(b)

685 r *= self.a

686 z *= self.b

687 else: # near-polar

688 r, z = _0_0, copysign0(self.b, lat)

689 else: # equator

690 r = self.a

691 z = lat = b = _0_0

692 return Circle4Tuple(r, z, lat, b)

693

694 def degrees2m(self, deg, lat=0):

695 '''Convert an angle to the distance along the equator or

696 along a parallel of (geodetic) latitude.

697

698 @arg deg: The angle (C{degrees}).

699 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

700

701 @return: Distance (C{meter}, same units as the equatorial

702 and polar radii) or C{0} for near-polar B{C{lat}}.

703

704 @raise RangeError: Latitude B{C{lat}} outside valid range and

705 L{rangerrors<pygeodesy.rangerrors>} is C{True}.

706

707 @raise ValueError: Invalid B{C{deg}} or B{C{lat}}.

708 '''

709 return self.radians2m(radians(deg), lat=lat)

710

711 def distance2(self, lat0, lon0, lat1, lon1):

712 '''I{Approximate} the distance and (initial) bearing between

713 two points based on the U{local, flat earth approximation

714 <https://www.EdWilliams.org/avform.htm#flat>} aka U{Hubeny

715 <https://www.OVG.AT/de/vgi/files/pdf/3781/>} formula.

716

717 I{Suitable only for distances of several hundred Km or Miles

718 and only between points not near-polar}.

719

720 @arg lat0: From latitude (C{degrees}).

721 @arg lon0: From longitude (C{degrees}).

722 @arg lat1: To latitude (C{degrees}).

723 @arg lon1: To longitude (C{degrees}).

724

725 @return: A L{Distance2Tuple}C{(distance, initial)} with C{distance}

726 in same units as this ellipsoid's axes.

727

728 @note: The meridional and prime_vertical radii of curvature are

729 taken and scaled I{at the initial latitude}, see C{roc2}.

730

731 @see: Function L{pygeodesy.flatLocal}/L{pygeodesy.hubeny}.

732 '''

733 phi0 = Phid(lat0=lat0)

734 m, n = self.roc2_(phi0, scaled=True)

735 m *= Phid(lat1=lat1) - phi0

736 n *= Lamd(lon1=lon1) - Lamd(lon0=lon0)

737 return Distance2Tuple(hypot(m, n), atan2b(n, m))

738

739 @Property_RO

740 def e(self):

741 '''Get the I{unsigned, (1st) eccentricity} (C{float}), M{sqrt(1 - (b / a)**2))}, see C{a_b2e}.

742

743 @see: Property L{es}.

744 '''

745 return Float(e=sqrt(self.e2abs) if self.e2 else _0_0)

746

747 @deprecated_Property_RO

748 def e12(self): # see property ._e12

749 '''DEPRECATED, use property C{e21}.'''

750 return self.e21

751

752# @Property_RO

753# def _e12(self): # see property ._elliptic_e12

754# # (INTERNAL) until e12 above can be replaced with e21.

755# return self.e2 / (_1_0 - self.e2) # see I{Karney}'s Ellipsoid._e12 = e2 / (1 - e2)

756

757 @Property_RO

758 def e2(self):

759 '''Get the I{signed, (1st) eccentricity squared} (C{float}), M{f * (2 - f)

760 == 1 - (b / a)**2}, see C{a_b2e2}.

761 '''

762 return self._assert(a_b2e2(self.a, self.b), e2=f2e2(self.f))

763

764 @Property_RO

765 def e2abs(self):

766 '''Get the I{unsigned, (1st) eccentricity squared} (C{float}).

767 '''

768 return fabs(self.e2)

769

770 @Property_RO

771 def e21(self):

772 '''Get 1 less I{1st eccentricity squared} (C{float}), M{1 - e**2}

773 == M{1 - e2} == M{(1 - f)**2} == M{b**2 / a**2}, see C{b2_a2}.

774 '''

775 return self._assert((_1_0 - self.f)**2, e21=_1_0 - self.e2, f0=_1_0)

776

777# _e2m = e21 # see I{Karney}'s Ellipsoid._e2m = 1 - _e2

778 _1_e21 = a2_b2 # == M{1 / e21} == M{1 / (1 - e**2)}

779

780 @Property_RO

781 def e22(self):

782 '''Get the I{signed, 2nd eccentricity squared} (C{float}), M{e2 / (1 - e2)

783 == e2 / (1 - f)**2 == (a / b)**2 - 1}, see C{a_b2e22}.

784 '''

785 return self._assert(a_b2e22(self.a, self.b), e22=f2e22(self.f))

786

787 @Property_RO

788 def e22abs(self):

789 '''Get the I{unsigned, 2nd eccentricity squared} (C{float}).

790 '''

791 return fabs(self.e22)

792

793 @Property_RO

794 def e32(self):

795 '''Get the I{signed, 3rd eccentricity squared} (C{float}), M{e2 / (2 - e2)

796 == (a**2 - b**2) / (a**2 + b**2)}, see C{a_b2e32}.

797 '''

798 return self._assert(a_b2e32(self.a, self.b), e32=f2e32(self.f))

799

800 @Property_RO

801 def e32abs(self):

802 '''Get the I{unsigned, 3rd eccentricity squared} (C{float}).

803 '''

804 return fabs(self.e32)

805

806 @Property_RO

807 def e4(self):

808 '''Get the I{unsignd, (1st) eccentricity} to 4th power (C{float}), M{e**4 == e2**2}.

809 '''

810 return Float(e4=self.e2**2 if self.e2 else _0_0)

811

812 eccentricity = e # eccentricity

813# eccentricity2 = e2 # eccentricity squared

814 eccentricity1st2 = e2 # first eccentricity squared

815 eccentricity2nd2 = e22 # second eccentricity squared

816 eccentricity3rd2 = e32 # third eccentricity squared

817

818 def ecef(self, Ecef=None):

819 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.

820

821 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.

822

823 @return: An ECEF converter for this C{ellipsoid}.

824

825 @raise TypeError: Invalid B{C{Ecef}}.

826

827 @see: Module L{pygeodesy.ecef}.

828 '''

829 return _MODS.ecef._4Ecef(self, Ecef)

830

831 @Property_RO

832 def _elliptic_e12(self): # see I{Karney}'s Ellipsoid._e12

833 '''(INTERNAL) Elliptic helper for C{Rhumb}.

834 '''

835 e12 = self.e2 / (self.e2 - _1_0) # NOT DEPRECATED .e12!

836 return _MODS.elliptic.Elliptic(e12)

837

838 @Property_RO

839 def _elliptic_e22(self): # aka ._elliptic_ep2

840 '''(INTERNAL) Elliptic helper for C{auxRectifying}, C{L}, C{Llat}.

841 '''

842 return _MODS.elliptic.Elliptic(-self.e22abs) # complex

843

844 equatoradius = a # Requatorial

845

846 def e2s(self, s):

847 '''Compute norm M{sqrt(1 - e2 * s**2)}.

848

849 @arg s: Sine value (C{scalar}).

850

851 @return: Norm (C{float}).

852

853 @raise ValueError: Invalid B{C{s}}.

854 '''

855 return sqrt(self.e2s2(s)) if self.e2 else _1_0

856

857 def e2s2(self, s):

858 '''Compute M{1 - e2 * s**2}.

859

860 @arg s: Sine value (C{scalar}).

861

862 @return: Result (C{float}).

863

864 @raise ValueError: Invalid B{C{s}}.

865 '''

866 r = _1_0

867 if self.e2:

868 try:

869 r -= self.e2 * Scalar(s=s)**2

870 if r < 0:

871 raise ValueError(_negative_)

872 except (TypeError, ValueError) as x:

873 t = self._DOT_(Ellipsoid.e2s2.__name__)

874 raise _ValueError(t, s, cause=x)

875 return r

876

877 @Property_RO

878 def es(self):

879 '''Get the I{signed (1st) eccentricity} (C{float}).

880

881 @see: Property L{e}.

882 '''

883 # note, self.e is always non-negative

884 return Float(es=copysign0(self.e, self.f)) # see .ups

885

886 def es_atanh(self, x):

887 '''Compute M{es * atanh(es * x)} or M{-es * atan(es * x)}

888 for I{oblate} respectively I{prolate} ellipsoids where

889 I{es} is the I{signed} (1st) eccentricity.

890

891 @raise ValueError: Invalid B{C{x}}.

892

893 @see: Function U{Math::eatanhe<https://GeographicLib.SourceForge.io/

894 C++/doc/classGeographicLib_1_1Math.html>}.

895 '''

896 return self._es_atanh(Scalar(x=x)) if self.f else _0_0

897

898 def _es_atanh(self, x): # see .albers._atanhee, .AuxLat._atanhee

899 '''(INTERNAL) Helper for .es_atanh, ._es_taupf2 and ._exp_es_atanh.

900 '''

901 es = self.es # signOf(es) == signOf(f)

902 return es * (atanh(es * x) if es > 0 else # .isOblate

903 (-atan(es * x) if es < 0 else # .isProlate

904 _0_0)) # .isSpherical

905

906 @Property_RO

907 def es_c(self):

908 '''Get M{(1 - f) * exp(es_atanh(1))} (C{float}), M{b_a * exp(es_atanh(1))}.

909 '''

910 return Float(es_c=(self._exp_es_atanh_1 * self.b_a) if self.f else _1_0)

911

912 def es_tauf(self, taup):

913 '''Compute I{Karney}'s U{equations (19), (20) and (21)

914 <https://ArXiv.org/abs/1002.1417>}.

915

916 @see: I{Karney}'s C++ method U{Math::tauf<https://GeographicLib.

917 SourceForge.io/C++/doc/classGeographicLib_1_1Math.html>} and

918 and I{Veness}' JavaScript method U{toLatLon<https://www.

919 Movable-Type.co.UK/scripts/latlong-utm-mgrs.html>}.

920 '''

921 t = Scalar(taup=taup)

922 if self.f: # .isEllipsoidal

923 a = fabs(t)

924 T = t * (self._exp_es_atanh_1 if a > 70 else self._1_e21)

925 if fabs(T * _EPSqrt) < _2_0: # handles +/- INF and NAN

926 s = (a * _TOL) if a > _1_0 else _TOL

927 for T, _, d in self._es_tauf3(t, T): # max 2

928 if fabs(d) < s:

929 break

930 t = Scalar(tauf=T)

931 return t

932

933 def _es_tauf3(self, taup, T, N=9): # in .utm.Utm._toLLEB

934 '''(INTERNAL) Yield a 3-tuple C{(τi, iteration, delta)} for at most

935 B{C{N}} Newton iterations, converging rapidly except when C{delta}

936 toggles on +/-1.12e-16 or +/-4.47e-16, see C{.utm.Utm._toLLEB}.

937 '''

938 e = self._1_e21

939 _F2_ = Fsum(T).fsum2f_ # τ0

940 _tf2 = self._es_taupf2

941 for i in range(1, N + 1):

942 a, h = _tf2(T)

943 d = (taup - a) * (e + T**2) / (hypot1(a) * h)

944 # = (taup - a) / hypot1(a) / ((e + T**2) / h)

945 T, d = _F2_(d) # τi, (τi - τi-1)

946 yield T, i, d

947

948 def es_taupf(self, tau):

949 '''Compute I{Karney}'s U{equations (7), (8) and (9)

950 <https://ArXiv.org/abs/1002.1417>}.

951

952 @see: I{Karney}'s C++ method U{Math::taupf<https://GeographicLib.

953 SourceForge.io/C++/doc/classGeographicLib_1_1Math.html>}.

954 '''

955 t = Scalar(tau=tau)

956 if self.f: # .isEllipsoidal

957 t, _ = self._es_taupf2(t)

958 t = Scalar(taupf=t)

959 return t

960

961 def _es_taupf2(self, tau):

962 '''(INTERNAL) Return 2-tuple C{(es_taupf(tau), hypot1(tau))}.

963 '''

964 if _isfinite(tau):

965 h = hypot1(tau)

966 s = sinh(self._es_atanh(tau / h))

967 a = hypot1(s) * tau - h * s

968 else:

969 a, h = tau, INF

970 return a, h

971

972 @Property_RO

973 def _exp_es_atanh_1(self):

974 '''(INTERNAL) Helper for .es_c and .es_tauf.

975 '''

976 return exp(self._es_atanh(_1_0)) if self.es else _1_0

977

978 @Property_RO

979 def f(self):

980 '''Get the I{flattening} (C{scalar}), M{(a - b) / a}, C{0} for spherical, negative for prolate.

981 '''

982 return self._f

983

984 @Property_RO

985 def f_(self):

986 '''Get the I{inverse flattening} (C{scalar}), M{1 / f} == M{a / (a - b)}, C{0} for spherical, see C{a_b2f_}.

987 '''

988 return self._f_

989

990 @Property_RO

991 def f1(self):

992 '''Get the I{1 - flattening} (C{float}), M{f1 == 1 - f == b / a}.

993

994 @see: Property L{b_a}.

995 '''

996 return Float(f1=_1_0 - self.f)

997

998 @Property_RO

999 def f2(self):

1000 '''Get the I{2nd flattening} (C{float}), M{(a - b) / b == f / (1 - f)}, C{0} for spherical, see C{a_b2f2}.

1001 '''

1002 return self._assert(self.a_b - _1_0, f2=f2f2(self.f))

1003

1004 @deprecated_Property_RO

1005 def geodesic(self):

1006 '''DEPRECATED, use property C{geodesicw}.'''

1007 return self.geodesicw

1008

1009 def geodesic_(self, exact=True):

1010 '''Get the an I{exact} C{Geodesic...} instance for this ellipsoid.

1011

1012 @kwarg exact: If C{bool} return L{GeodesicExact}C{(exact=B{exact}, ...)},

1013 otherwise a L{Geodesic}, L{GeodesicExact} or L{GeodesicSolve}

1014 instance for I{this} ellipsoid.

1015

1016 @return: The C{exact} geodesic (C{Geodesic...}).

1017

1018 @raise TypeError: Invalid B{C{exact}}.

1019

1020 @raise ValueError: Incompatible B{C{exact}} ellipsoid.

1021 '''

1022 if isbool(exact): # for consistenccy with C{.rhumb_}

1023 g = _MODS.geodesicx.GeodesicExact(self, C4order=30 if exact else 24,

1024 name=self.name)

1025 else:

1026 g = exact

1027 E = _xattr(g, ellipsoid=None)

1028 if not (E is self and isinstance(g, self._Geodesics)):

1029 raise _ValueError(exact=g, ellipsoid=E, txt_not_=self.name)

1030 return g

1031

1032 @property_RO

1033 def _Geodesics(self):

1034 '''(INTERNAL) Get all C{Geodesic...} classes, I{once}.

1035 '''

1036 try:

1037 t = _MODS.geodesicw._wrapped.Geodesic,

1038 except ImportError:

1039 t = ()

1040 Ellipsoid._Geodesics = t = (_MODS.geodesicx.GeodesicExact, # overwrite property_RO

1041 _MODS.geodsolve.GeodesicSolve) + t

1042 return t

1043

1044 @property_RO

1045 def geodesicw(self):

1046 '''Get this ellipsoid's I{wrapped} U{geodesicw.Geodesic

1047 <https://GeographicLib.SourceForge.io/Python/doc/code.html>}, provided

1048 I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}

1049 package is installed.

1050 '''

1051 # if not self.isEllipsoidal:

1052 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1053 return _MODS.geodesicw.Geodesic(self)

1054

1055 @property_RO

1056 def geodesicx(self):

1057 '''Get this ellipsoid's I{exact} L{GeodesicExact}.

1058 '''

1059 # if not self.isEllipsoidal:

1060 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1061 return _MODS.geodesicx.GeodesicExact(self, name=self.name)

1062

1063 @property

1064 def geodsolve(self):

1065 '''Get this ellipsoid's L{GeodesicSolve}, the I{wrapper} around utility

1066 U{GeodSolve<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>},

1067 provided the path to the C{GeodSolve} executable is specified with env

1068 variable C{PYGEODESY_GEODSOLVE} or re-/set with this property..

1069 '''

1070 # if not self.isEllipsoidal:

1071 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1072 return _MODS.geodsolve.GeodesicSolve(self, path=self._geodsolve, name=self.name)

1073

1074 @geodsolve.setter # PYCHOK setter!

1075 def geodsolve(self, path):

1076 '''Re-/set the (fully qualified) path to the U{GeodSolve

1077 <https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>} executable,

1078 overriding env variable C{PYGEODESY_GEODSOLVE} (C{str}).

1079 '''

1080 self._geodsolve = path

1081

1082 def hartzell4(self, pov, los=None):

1083 '''Compute the intersection of this ellipsoid's surface and a Line-Of-Sight

1084 from a Point-Of-View in space.

1085

1086 @arg pov: Point-Of-View outside this ellipsoid (C{Cartesian}, L{Ecef9Tuple}

1087 or L{Vector3d}).

1088 @kwarg los: Line-Of-Sight, I{direction} to this ellipsoid (L{Vector3d}) or

1089 C{None} to point to this ellipsoid's center.

1090

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

1092 C{y} and C{z} of the projection on or the intersection with this

1093 ellipsoid and the I{distance} C{h} from B{C{pov}} to C{(x, y, z)}

1094 along B{C{los}}, all in C{meter}, conventionally.

1095

1096 @raise IntersectionError: Null B{C{pov}} or B{C{los}} vector, or B{C{pov}}

1097 is inside this ellipsoid or B{C{los}} points

1098 outside this ellipsoid or points in an opposite

1099 direction.

1100

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

1102

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

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

1105 methods L{Ellipsoid.height4} and L{Triaxial.hartzell4}.

1106 '''

1107 try:

1108 v, d, _ = _MODS.triaxials._hartzell3(pov, los, self._triaxial)

1109 except Exception as x:

1110 raise IntersectionError(pov=pov, los=los, cause=x)

1111 return Vector4Tuple(v.x, v.y, v.z, d, name__=self.hartzell4)

1112

1113 @Property_RO

1114 def _hash(self):

1115 return hash((self.a, self.f))

1116

1117 def height4(self, xyz, normal=True):

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

1119 this ellipsoid's surface.

1120

1121 @arg xyz: The cartesian (C{Cartesian}, L{Ecef9Tuple}, L{Vector3d},

1122 L{Vector3Tuple} or L{Vector4Tuple}).

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

1124 point on) this ellipsoid's surface, otherwise the C{radial}

1125 line to this ellipsoid's center (C{bool}).

1126

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

1128 C{y} and C{z} of the projection on and the height C{h} above or

1129 below this ellipsoid's surface, all in C{meter}, conventionally.

1130

1131 @raise ValueError: Null B{C{xyz}}.

1132

1133 @raise TypeError: Non-cartesian B{C{xyz}}.

1134

1135 @see: U{Distance to<https://StackOverflow.com/questions/22959698/distance-from-given-point-to-given-ellipse>}

1136 and U{intersection with<https://MathWorld.wolfram.com/Ellipse-LineIntersection.html>} an ellipse and

1137 methods L{Ellipsoid.hartzell4} and L{Triaxial.height4}.

1138 '''

1139 v = _MODS.vector3d._otherV3d(xyz=xyz)

1140 r = v.length

1141

1142 a, b, i = self.a, self.b, None

1143 if r < EPS0: # EPS

1144 v = v.times(_0_0)

1145 h = -a

1146

1147 elif self.isSpherical:

1148 v = v.times(a / r)

1149 h = r - a

1150

1151 elif normal: # perpendicular to ellipsoid

1152 x, y = hypot(v.x, v.y), fabs(v.z)

1153 if x < EPS0: # PYCHOK no cover

1154 z = copysign0(b, v.z)

1155 v = Vector3Tuple(v.x, v.y, z)

1156 h = y - b # polar

1157 elif y < EPS0: # PYCHOK no cover

1158 t = a / r

1159 v = v.times_(t, t, 0) # force z=0.0

1160 h = x - a # equatorial

1161 else: # normal in 1st quadrant

1162 x, y, i = _plumbTo3(x, y, self)

1163 t, v = v, v.times_(x, x, y)

1164 h = t.minus(v).length

1165

1166 else: # radial to ellipsoid's center

1167 h = hypot_(a * v.z, b * v.x, b * v.y)

1168 t = (a * b / h) if h > EPS0 else _0_0 # EPS

1169 v = v.times(t)

1170 h = r * (_1_0 - t)

1171

1172 return Vector4Tuple(v.x, v.y, v.z, h, iteration=i, name__=self.height4)

1173

1174 def _hubeny_2(self, phi2, phi1, lam21, scaled=True, squared=True):

1175 '''(INTERNAL) like function C{pygeodesy.flatLocal_}/C{pygeodesy.hubeny_},

1176 returning the I{angular} distance in C{radians squared} or C{radians}

1177 '''

1178 m, n = self.roc2_((phi2 + phi1) * _0_5, scaled=scaled)

1179 h, r = (hypot2, self.a2_) if squared else (hypot, _1_0 / self.a)

1180 return h(m * (phi2 - phi1), n * lam21) * r

1181

1182 @Property_RO

1183 def isEllipsoidal(self):

1184 '''Is this model I{ellipsoidal} (C{bool})?

1185 '''

1186 return self.f != 0

1187

1188 @Property_RO

1189 def isOblate(self):

1190 '''Is this ellipsoid I{oblate} (C{bool})? I{Prolate} or

1191 spherical otherwise.

1192 '''

1193 return self.f > 0

1194

1195 @Property_RO

1196 def isProlate(self):

1197 '''Is this ellipsoid I{prolate} (C{bool})? I{Oblate} or

1198 spherical otherwise.

1199 '''

1200 return self.f < 0

1201

1202 @Property_RO

1203 def isSpherical(self):

1204 '''Is this ellipsoid I{spherical} (C{bool})?

1205 '''

1206 return self.f == 0

1207

1208 def _Kseries(self, *AB8Ks):

1209 '''(INTERNAL) Compute the 4-, 6- or 8-th order I{Krüger} Alpha

1210 or Beta series coefficients per I{Karney}'s U{equations (35)

1211 and (36)<https://ArXiv.org/pdf/1002.1417v3.pdf>}.

1212

1213 @arg AB8Ks: 8-Tuple of 8-th order I{Krüger} Alpha or Beta series

1214 coefficient tuples.

1215

1216 @return: I{Krüger} series coefficients (L{KsOrder}C{-tuple}).

1217

1218 @see: I{Karney}'s 30-th order U{TMseries30

1219 <https://GeographicLib.SourceForge.io/C++/doc/tmseries30.html>}.

1220 '''

1221 k = self.KsOrder

1222 if self.n:

1223 ns = fpowers(self.n, k)

1224 ks = tuple(fdot(AB8Ks[i][:k-i], *ns[i:]) for i in range(k))

1225 else:

1226 ks = _0_0s(k)

1227 return ks

1228

1229 @property_doc_(''' the I{Krüger} series' order (C{int}), see properties C{AlphaKs}, C{BetaKs}.''')

1230 def KsOrder(self):

1231 '''Get the I{Krüger} series' order (C{int} 4, 6 or 8).

1232 '''

1233 return self._KsOrder

1234

1235 @KsOrder.setter # PYCHOK setter!

1236 def KsOrder(self, order):

1237 '''Set the I{Krüger} series' order (C{int} 4, 6 or 8).

1238

1239 @raise ValueError: Invalid B{C{order}}.

1240 '''

1241 if not (isint(order) and order in (4, 6, 8)):

1242 raise _ValueError(order=order)

1243 if self._KsOrder != order:

1244 Ellipsoid.AlphaKs._update(self)

1245 Ellipsoid.BetaKs._update(self)

1246 self._KsOrder = order

1247

1248 @Property_RO

1249 def L(self):

1250 '''Get the I{quarter meridian} C{L}, aka the C{polar distance}

1251 along a meridian between the equator and a pole (C{meter}),

1252 M{b * Elliptic(-e2 / (1 - e2)).cE} or M{b * PI / 2}.

1253 '''

1254 r = self._elliptic_e22.cE if self.f else PI_2

1255 return Distance(L=self.b * r)

1256

1257 def Llat(self, lat):

1258 '''Return the I{meridional length}, the distance along a meridian

1259 between the equator and a (geodetic) latitude, see C{L}.

1260

1261 @arg lat: Geodetic latitude (C{degrees90}).

1262

1263 @return: The meridional length at B{C{lat}}, negative on southern

1264 hemisphere (C{meter}).

1265 '''

1266 r = self._elliptic_e22.fEd(self.auxParametric(lat)) if self.f else Phid(lat)

1267 return Distance(Llat=self.b * r)

1268

1269 Lmeridian = Llat # meridional distance

1270

1271 @property_RO

1272 def _Lpd(self):

1273 '''Get the I{quarter meridian} per degree (C{meter}), M{self.L / 90}.

1274 '''

1275 return Meter(_Lpd=self.L / _90_0)

1276

1277 @property_RO

1278 def _Lpr(self):

1279 '''Get the I{quarter meridian} per radian (C{meter}), M{self.L / PI_2}.

1280 '''

1281 return Meter(_Lpr=self.L / PI_2)

1282

1283 @deprecated_Property_RO

1284 def majoradius(self): # PYCHOK no cover

1285 '''DEPRECATED, use property C{a} or C{Requatorial}.'''

1286 return self.a

1287

1288 def m2degrees(self, distance, lat=0):

1289 '''Convert a distance to an angle along the equator or along

1290 a parallel of (geodetic) latitude.

1291

1292 @arg distance: Distance (C{meter}).

1293 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

1294

1295 @return: Angle (C{degrees}) or C{INF} for near-polar B{C{lat}}.

1296

1297 @raise RangeError: Latitude B{C{lat}} outside valid range and

1298 L{rangerrors<pygeodesy.rangerrors>} is C{True}.

1299

1300 @raise ValueError: Invalid B{C{distance}} or B{C{lat}}.

1301 '''

1302 return degrees(self.m2radians(distance, lat=lat))

1303

1304 def m2radians(self, distance, lat=0):

1305 '''Convert a distance to an angle along the equator or

1306 along a parallel of (geodetic) latitude.

1307

1308 @arg distance: Distance (C{meter}).

1309 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

1310

1311 @return: Angle (C{radians}) or C{INF} for near-polar B{C{lat}}.

1312

1313 @raise RangeError: Latitude B{C{lat}} outside valid range and

1314 L{rangerrors<pygeodesy.rangerrors>} is C{True}.

1315

1316 @raise ValueError: Invalid B{C{distance}} or B{C{lat}}.

1317 '''

1318 r = self.circle4(lat).radius if lat else self.a

1319 return m2radians(distance, radius=r, lat=0)

1320

1321 @deprecated_Property_RO

1322 def minoradius(self): # PYCHOK no cover

1323 '''DEPRECATED, use property C{b}, C{polaradius} or C{Rpolar}.'''

1324 return self.b

1325

1326 @Property_RO

1327 def n(self):

1328 '''Get the I{3rd flattening} (C{float}), M{f / (2 - f) == (a - b) / (a + b)}, see C{a_b2n}.

1329 '''

1330 return self._assert(a_b2n(self.a, self.b), n=f2n(self.f))

1331

1332 flattening = f

1333 flattening1st = f

1334 flattening2nd = f2

1335 flattening3rd = n

1336

1337 polaradius = b # Rpolar

1338

1339# @Property_RO

1340# def Q(self):

1341# '''Get the I{meridian arc unit} C{Q}, the mean, meridional length I{per radian} C({float}).

1342#

1343# @note: C{Q * PI / 2} ≈ C{L}, the I{quarter meridian}.

1344#

1345# @see: Property C{A} and U{Engsager, K., Poder, K.<https://StudyLib.net/doc/7443565/

1346# a-highly-accurate-world-wide-algorithm-for-the-transverse...>}.

1347# '''

1348# n = self.n

1349# d = (n + _1_0) / self.a

1350# return Float(Q=Fhorner(n**2, _1_0, _0_25, _1_16th, _0_25).fover(d) if d else self.b)

1351

1352# # Moritz, H. <https://Geodesy.Geology.Ohio-State.EDU/course/refpapers/00740128.pdf>

1353# # Q = (1 - 3/4 * e'2 + 45/64 * e'4 - 175/256 * e'6 + 11025/16384 * e'8) * rocPolar

1354# # = (4 + e'2 * (-3 + e'2 * (45/16 + e'2 * (-175/64 + e'2 * 11025/4096)))) * rocPolar / 4

1355# return Fhorner(self.e22, 4, -3, 45 / 16, -175 / 64, 11025 / 4096).fover(4 / self.rocPolar)

1356

1357 @deprecated_Property_RO

1358 def quarteradius(self): # PYCHOK no cover

1359 '''DEPRECATED, use property C{L} or method C{Llat}.'''

1360 return self.L

1361

1362 @Property_RO

1363 def R1(self):

1364 '''Get the I{mean} earth radius per I{IUGG} (C{meter}), M{(2 * a + b) / 3 == a * (1 - f / 3)}.

1365

1366 @see: U{Earth radius<https://WikiPedia.org/wiki/Earth_radius>}

1367 and method C{Rgeometric}.

1368 '''

1369 r = Fsum(self.a, self.a, self.b).fover(_3_0) if self.f else self.a

1370 return Radius(R1=r)

1371

1372 Rmean = R1

1373

1374 @Property_RO

1375 def R2(self):

1376 '''Get the I{authalic} earth radius (C{meter}), M{sqrt(c2)}.

1377

1378 @see: C{R2x}, C{c2}, C{area} and U{Earth radius

1379 <https://WikiPedia.org/wiki/Earth_radius>}.

1380 '''

1381 return Radius(R2=sqrt(self.c2) if self.f else self.a)

1382

1383 Rauthalic = R2

1384

1385# @Property_RO

1386# def R2(self):

1387# # Moritz, H. <https://Geodesy.Geology.Ohio-State.EDU/course/refpapers/00740128.pdf>

1388# # R2 = (1 - 2/3 * e'2 + 26/45 * e'4 - 100/189 * e'6 + 7034/14175 * e'8) * rocPolar

1389# # = (3 + e'2 * (-2 + e'2 * (26/15 + e'2 * (-100/63 + e'2 * 7034/4725)))) * rocPolar / 3

1390# return Fhorner(self.e22, 3, -2, 26 / 15, -100 / 63, 7034 / 4725).fover(3 / self.rocPolar)

1391

1392 @Property_RO

1393 def R2x(self):

1394 '''Get the I{authalic} earth radius (C{meter}), M{sqrt(c2x)}.

1395

1396 @see: C{R2}, C{c2x} and C{areax}.

1397 '''

1398 return Radius(R2x=sqrt(self.c2x) if self.f else self.a)

1399

1400 Rauthalicx = R2x

1401

1402 @Property_RO

1403 def R3(self):

1404 '''Get the I{volumetric} earth radius (C{meter}), M{(a * a * b)**(1/3)}.

1405

1406 @see: U{Earth radius<https://WikiPedia.org/wiki/Earth_radius>} and C{volume}.

1407 '''

1408 r = (cbrt(self.b_a) * self.a) if self.f else self.a

1409 return Radius(R3=r)

1410

1411 Rvolumetric = R3

1412

1413 def radians2m(self, rad, lat=0):

1414 '''Convert an angle to the distance along the equator or along

1415 a parallel of (geodetic) latitude.

1416

1417 @arg rad: The angle (C{radians}).

1418 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

1419

1420 @return: Distance (C{meter}, same units as the equatorial

1421 and polar radii) or C{0} for near-polar B{C{lat}}.

1422

1423 @raise RangeError: Latitude B{C{lat}} outside valid range and

1424 L{rangerrors<pygeodesy.rangerrors>} is C{True}.

1425

1426 @raise ValueError: Invalid B{C{rad}} or B{C{lat}}.

1427 '''

1428 r = self.circle4(lat).radius if lat else self.a

1429 return radians2m(rad, radius=r, lat=0)

1430

1431 @Property_RO

1432 def Rbiaxial(self):

1433 '''Get the I{biaxial, quadratic} mean earth radius (C{meter}), M{sqrt((a**2 + b**2) / 2)}.

1434

1435 @see: C{Rtriaxial}

1436 '''

1437 a, b = self.a, self.b

1438 if b < a:

1439 b = sqrt(_0_5 + self.b2_a2 * _0_5) * a

1440 elif b > a:

1441 b *= sqrt(_0_5 + self.a2_b2 * _0_5)

1442 return Radius(Rbiaxial=b)

1443

1444 Requatorial = a # for consistent naming

1445

1446 def Rgeocentric(self, lat):

1447 '''Compute the I{geocentric} earth radius of (geodetic) latitude.

1448

1449 @arg lat: Latitude (C{degrees90}).

1450

1451 @return: Geocentric earth radius (C{meter}).

1452

1453 @raise ValueError: Invalid B{C{lat}}.

1454

1455 @see: U{Geocentric Radius

1456 <https://WikiPedia.org/wiki/Earth_radius#Geocentric_radius>}

1457 '''

1458 r, a = self.a, Phid(lat)

1459 if a and self.f:

1460 if fabs(a) < PI_2:

1461 s2, c2 = _s2_c2(a)

1462 b2_a2_s2 = self.b2_a2 * s2

1463 # R == sqrt((a2**2 * c2 + b2**2 * s2) / (a2 * c2 + b2 * s2))

1464 # == sqrt(a2**2 * (c2 + (b2 / a2)**2 * s2) / (a2 * (c2 + b2 / a2 * s2)))

1465 # == sqrt(a2 * (c2 + (b2 / a2)**2 * s2) / (c2 + (b2 / a2) * s2))

1466 # == a * sqrt((c2 + b2_a2 * b2_a2 * s2) / (c2 + b2_a2 * s2))

1467 # == a * sqrt((c2 + b2_a2 * b2_a2_s2) / (c2 + b2_a2_s2))

1468 r *= sqrt((c2 + b2_a2_s2 * self.b2_a2) / (c2 + b2_a2_s2))

1469 else:

1470 r = self.b

1471 return Radius(Rgeocentric=r)

1472

1473 @Property_RO

1474 def Rgeometric(self):

1475 '''Get the I{geometric} mean earth radius (C{meter}), M{sqrt(a * b)}.

1476

1477 @see: C{R1}.

1478 '''

1479 g = sqrt(self.a * self.b) if self.f else self.a

1480 return Radius(Rgeometric=g)

1481

1482 def rhumb_(self, exact=True):

1483 '''Get the an I{exact} C{Rhumb...} instance for this ellipsoid.

1484

1485 @kwarg exact: If C{bool} or C{None} return L{Rhumb}C{(exact=B{exact}, ...)},

1486 otherwise a L{Rhumb}, L{RhumbAux} or L{RhumbSolve} instance

1487 for I{this} ellipsoid.

1488

1489 @return: The C{exact} rhumb (C{Rhumb...}).

1490

1491 @raise TypeError: Invalid B{C{exact}}.

1492

1493 @raise ValueError: Incompatible B{C{exact}} ellipsoid.

1494 '''

1495 if isbool(exact): # use Rhumb for backward compatibility

1496 r = _MODS.rhumb.ekx.Rhumb(self, exact=exact, name=self.name)

1497 else:

1498 r = exact

1499 E = _xattr(r, ellipsoid=None)

1500 if not (E is self and isinstance(r, self._Rhumbs)):

1501 raise _ValueError(exact=r, ellipsosid=E, txt_not_=self.name)

1502 return r

1503

1504 @property_RO

1505 def rhumbaux(self):

1506 '''Get this ellipsoid's I{Auxiliary} C{rhumb.RhumbAux}.

1507 '''

1508 # if not self.isEllipsoidal:

1509 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1510 return _MODS.rhumb.aux_.RhumbAux(self, name=self.name)

1511

1512 @property_RO

1513 def rhumbekx(self):

1514 '''Get this ellipsoid's I{Elliptic, Krüger} C{rhumb.Rhumb}.

1515 '''

1516 # if not self.isEllipsoidal:

1517 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1518 return _MODS.rhumb.ekx.Rhumb(self, name=self.name)

1519

1520 @property_RO

1521 def _Rhumbs(self):

1522 '''(INTERNAL) Get all C{Rhumb...} classes, I{once}.

1523 '''

1524 p = _MODS.rhumb

1525 Ellipsoid._Rhumbs = t = (p.aux_.RhumbAux, # overwrite property_RO

1526 p.ekx.Rhumb, p.solve.RhumbSolve)

1527 return t

1528

1529 @property

1530 def rhumbsolve(self):

1531 '''Get this ellipsoid's L{RhumbSolve}, the I{wrapper} around utility

1532 U{RhumbSolve<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>},

1533 provided the path to the C{RhumbSolve} executable is specified with env

1534 variable C{PYGEODESY_RHUMBSOLVE} or re-/set with this property.

1535 '''

1536 # if not self.isEllipsoidal:

1537 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1538 return _MODS.rhumb.solve.RhumbSolve(self, path=self._rhumbsolve, name=self.name)

1539

1540 @rhumbsolve.setter # PYCHOK setter!

1541 def rhumbsolve(self, path):

1542 '''Re-/set the (fully qualified) path to the U{RhumbSolve

1543 <https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>} executable,

1544 overriding env variable C{PYGEODESY_RHUMBSOLVE} (C{str}).

1545 '''

1546 self._rhumbsolve = path

1547

1548 @deprecated_property_RO

1549 def rhumbx(self):

1550 '''DEPRECATED on 2023.11.28, use property C{rhumbekx}. '''

1551 return self.rhumbekx

1552

1553 def Rlat(self, lat):

1554 '''I{Approximate} the earth radius of (geodetic) latitude.

1555

1556 @arg lat: Latitude (C{degrees90}).

1557

1558 @return: Approximate earth radius (C{meter}).

1559

1560 @raise RangeError: Latitude B{C{lat}} outside valid range and

1561 L{rangerrors<pygeodesy.rangerrors>} is C{True}.

1562

1563 @raise TypeError: Invalid B{C{lat}}.

1564

1565 @raise ValueError: Invalid B{C{lat}}.

1566

1567 @note: C{Rlat(B{90})} equals C{Rpolar}.

1568

1569 @see: Method C{circle4}.

1570 '''

1571 # r = a - (a - b) * |lat| / 90

1572 r = self.a

1573 if self.f and lat: # .isEllipsoidal

1574 r -= (r - self.b) * fabs(Lat(lat)) / _90_0

1575 r = Radius(Rlat=r)

1576 return r

1577

1578 Rpolar = b # for consistent naming

1579

1580 def roc1_(self, sa, ca=None):

1581 '''Compute the I{prime-vertical}, I{normal} radius of curvature

1582 of (geodetic) latitude, I{unscaled}.

1583

1584 @arg sa: Sine of the latitude (C{float}, [-1.0..+1.0]).

1585 @kwarg ca: Optional cosine of the latitude (C{float}, [-1.0..+1.0])

1586 to use an alternate formula.

1587

1588 @return: The prime-vertical radius of curvature (C{float}).

1589

1590 @note: The delta between both formulae with C{Ellipsoids.WGS84}

1591 is less than 2 nanometer over the entire latitude range.

1592

1593 @see: Method L{roc2_} and class L{EcefYou}.

1594 '''

1595 if not self.f: # .isSpherical

1596 n = self.a

1597 elif ca is None:

1598 r = self.e2s2(sa) # see .roc2_ and _EcefBase._forward

1599 n = sqrt(self.a2 / r) if r > EPS02 else _0_0

1600 elif ca: # derived from EcefYou.forward

1601 h = hypot(ca, self.b_a * sa) if sa else fabs(ca)

1602 n = self.a / h

1603 elif sa:

1604 n = self.a2_b / fabs(sa)

1605 else:

1606 n = self.a

1607 return n

1608

1609 def roc2(self, lat, scaled=False):

1610 '''Compute the I{meridional} and I{prime-vertical}, I{normal}

1611 radii of curvature of (geodetic) latitude.

1612

1613 @arg lat: Latitude (C{degrees90}).

1614 @kwarg scaled: Scale prime_vertical by C{cos(radians(B{lat}))} (C{bool}).

1615

1616 @return: An L{Curvature2Tuple}C{(meridional, prime_vertical)} with

1617 the radii of curvature.

1618

1619 @raise ValueError: Invalid B{C{lat}}.

1620

1621 @see: Methods L{roc2_} and L{roc1_}, U{Local, flat earth approximation

1622 <https://www.EdWilliams.org/avform.htm#flat>} and meridional and

1623 prime vertical U{Radii of Curvature<https://WikiPedia.org/wiki/

1624 Earth_radius#Radii_of_curvature>}.

1625 '''

1626 return self.roc2_(Phid(lat), scaled=scaled)

1627

1628 def roc2_(self, phi, scaled=False):

1629 '''Compute the I{meridional} and I{prime-vertical}, I{normal} radii of

1630 curvature of (geodetic) latitude.

1631

1632 @arg phi: Latitude (C{radians}).

1633 @kwarg scaled: Scale prime_vertical by C{cos(B{phi})} (C{bool}).

1634

1635 @return: An L{Curvature2Tuple}C{(meridional, prime_vertical)} with the

1636 radii of curvature.

1637

1638 @raise ValueError: Invalid B{C{phi}}.

1639

1640 @see: Methods L{roc2} and L{roc1_}, property L{rocEquatorial2}, U{Local,

1641 flat earth approximation<https://www.EdWilliams.org/avform.htm#flat>}

1642 and the meridional and prime vertical U{Radii of Curvature

1643 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}.

1644 '''

1645 a = fabs(Phi(phi))

1646 if self.f:

1647 r = self.e2s2(sin(a))

1648 if r > EPS02:

1649 n = self.a / sqrt(r)

1650 m = n * self.e21 / r

1651 else:

1652 m = n = _0_0

1653 else:

1654 m = n = self.a

1655 if scaled and a:

1656 n *= cos(a) if a < PI_2 else _0_0

1657 return Curvature2Tuple(Radius(rocMeridional=m),

1658 Radius(rocPrimeVertical=n))

1659

1660 def rocBearing(self, lat, bearing):

1661 '''Compute the I{directional} radius of curvature of (geodetic)

1662 latitude and compass direction.

1663

1664 @arg lat: Latitude (C{degrees90}).

1665 @arg bearing: Direction (compass C{degrees360}).

1666

1667 @return: Directional radius of curvature (C{meter}).

1668

1669 @raise RangeError: Latitude B{C{lat}} outside valid range and

1670 L{rangerrors<pygeodesy.rangerrors>} is C{True}.

1671

1672 @raise ValueError: Invalid B{C{lat}} or B{C{bearing}}.

1673

1674 @see: U{Radii of Curvature<https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}

1675 '''

1676 if self.f:

1677 s2, c2 = _s2_c2(Bearing_(bearing))

1678 m, n = self.roc2_(Phid(lat))

1679 if n < m: # == n / (c2 * n / m + s2)

1680 c2 *= n / m

1681 elif m < n: # == m / (c2 + s2 * m / n)

1682 s2 *= m / n

1683 n = m

1684 b = n / (c2 + s2) # == 1 / (c2 / m + s2 / n)

1685 else:

1686 b = self.b # == self.a

1687 return Radius(rocBearing=b)

1688

1689 @Property_RO

1690 def rocEquatorial2(self):

1691 '''Get the I{meridional} and I{prime-vertical}, I{normal} radii of curvature

1692 at the equator as L{Curvature2Tuple}C{(meridional, prime_vertical)}.

1693

1694 @see: Methods L{rocMeridional} and L{rocPrimeVertical}, properties L{b2_a},

1695 L{a2_b}, C{rocPolar} and polar and equatorial U{Radii of Curvature

1696 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}.

1697 '''

1698 return Curvature2Tuple(Radius(rocMeridional0=self.b2_a if self.f else self.a),

1699 Radius(rocPrimeVertical0=self.a))

1700

1701 def rocGauss(self, lat):

1702 '''Compute the I{Gaussian} radius of curvature of (geodetic) latitude.

1703

1704 @arg lat: Latitude (C{degrees90}).

1705

1706 @return: Gaussian radius of curvature (C{meter}).

1707

1708 @raise ValueError: Invalid B{C{lat}}.

1709

1710 @see: Non-directional U{Radii of Curvature<https://WikiPedia.org/wiki/

1711 Earth_radius#Radii_of_curvature>}

1712 '''

1713 # using ...

1714 # m, n = self.roc2_(Phid(lat))

1715 # return sqrt(m * n)

1716 # ... requires 1 or 2 sqrt

1717 g = self.b

1718 if self.f:

1719 s2, c2 = _s2_c2(Phid(lat))

1720 g = g / (c2 + self.b2_a2 * s2)

1721 return Radius(rocGauss=g)

1722

1723 def rocMean(self, lat):

1724 '''Compute the I{mean} radius of curvature of (geodetic) latitude.

1725

1726 @arg lat: Latitude (C{degrees90}).

1727

1728 @return: Mean radius of curvature (C{meter}).

1729

1730 @raise ValueError: Invalid B{C{lat}}.

1731

1732 @see: Non-directional U{Radii of Curvature<https://WikiPedia.org/wiki/

1733 Earth_radius#Radii_of_curvature>}

1734 '''

1735 if self.f:

1736 m, n = self.roc2_(Phid(lat))

1737 m *= n * _2_0 / (m + n) # == 2 / (1 / m + 1 / n)

1738 else:

1739 m = self.a

1740 return Radius(rocMean=m)

1741

1742 def rocMeridional(self, lat):

1743 '''Compute the I{meridional} radius of curvature of (geodetic) latitude.

1744

1745 @arg lat: Latitude (C{degrees90}).

1746

1747 @return: Meridional radius of curvature (C{meter}).

1748

1749 @raise ValueError: Invalid B{C{lat}}.

1750

1751 @see: Methods L{roc2} and L{roc2_}, U{Local, flat earth approximation

1752 <https://www.EdWilliams.org/avform.htm#flat>} and U{Radii of

1753 Curvature<https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}.

1754 '''

1755 return self.roc2_(Phid(lat)).meridional if lat else \

1756 self.rocEquatorial2.meridional

1757

1758 rocPolar = a2_b # synonymous

1759

1760 def rocPrimeVertical(self, lat):

1761 '''Compute the I{prime-vertical}, I{normal} radius of curvature of

1762 (geodetic) latitude, aka the I{transverse} radius of curvature.

1763

1764 @arg lat: Latitude (C{degrees90}).

1765

1766 @return: Prime-vertical radius of curvature (C{meter}).

1767

1768 @raise ValueError: Invalid B{C{lat}}.

1769

1770 @see: Methods L{roc2}, L{roc2_} and L{roc1_}, U{Local, flat earth

1771 approximation<https://www.EdWilliams.org/avform.htm#flat>} and

1772 U{Radii of Curvature<https://WikiPedia.org/wiki/

1773 Earth_radius#Radii_of_curvature>}.

1774 '''

1775 return self.roc2_(Phid(lat)).prime_vertical if lat else \

1776 self.rocEquatorial2.prime_vertical

1777

1778 rocTransverse = rocPrimeVertical # synonymous

1779

1780 @deprecated_Property_RO

1781 def Rquadratic(self): # PYCHOK no cover

1782 '''DEPRECATED, use property C{Rbiaxial} or C{Rtriaxial}.'''

1783 return self.Rbiaxial

1784

1785 @deprecated_Property_RO

1786 def Rr(self): # PYCHOK no cover

1787 '''DEPRECATED, use property C{Rrectifying}.'''

1788 return self.Rrectifying

1789

1790 @Property_RO

1791 def Rrectifying(self):

1792 '''Get the I{rectifying} earth radius (C{meter}), M{((a**(3/2) + b**(3/2)) / 2)**(2/3)}.

1793

1794 @see: U{Earth radius<https://WikiPedia.org/wiki/Earth_radius>}.

1795 '''

1796 r = (cbrt2((_1_0 + sqrt3(self.b_a)) * _0_5) * self.a) if self.f else self.a

1797 return Radius(Rrectifying=r)

1798

1799 @deprecated_Property_RO

1800 def Rs(self): # PYCHOK no cover

1801 '''DEPRECATED, use property C{Rgeometric}.'''

1802 return self.Rgeometric

1803

1804 @Property_RO

1805 def Rtriaxial(self):

1806 '''Get the I{triaxial, quadratic} mean earth radius (C{meter}), M{sqrt((3 * a**2 + b**2) / 4)}.

1807

1808 @see: C{Rbiaxial}

1809 '''

1810 a, b = self.a, self.b

1811 q = (sqrt((_3_0 + self.b2_a2) * _0_25) * a) if a > b else (

1812 (sqrt((_3_0 * self.a2_b2 + _1_0) * _0_25) * b) if a < b else a)

1813 return Radius(Rtriaxial=q)

1814

1815 def toEllipsoid2(self, **name):

1816 '''Get a copy of this ellipsoid as an L{Ellipsoid2}.

1817

1818 @kwarg name: Optional, unique C{B{name}=NN} (C{str}).

1819

1820 @see: Property C{a_f}.

1821 '''

1822 return Ellipsoid2(self, None, **name)

1823

1824 def toStr(self, prec=8, terse=4, **sep_name): # PYCHOK expected

1825 '''Return this ellipsoid as a text string.

1826

1827 @kwarg prec: Number of decimal digits, unstripped (C{int}).

1828 @kwarg terse: Limit the number of items (C{int}, 0...18),

1829 use C{B{terse}=0} or C{=None} for all.

1830 @kwarg sep_name: Optional C{B{name}=NN} (C{str}) or C{None}

1831 to exclude this ellipsoid's name and separator

1832 C{B{sep}=", "} to join the items (C{str}).

1833

1834 @return: This C{Ellipsoid}'s attributes (C{str}).

1835 '''

1836 E = Ellipsoid

1837 t = (E.a, E.f, E.f_, E.b, E.f2, E.n, E.e,

1838 E.e2, E.e21, E.e22, E.e32,

1839 E.A, E.L, E.R1, E.R2, E.R3,

1840 E.Rbiaxial, E.Rtriaxial)

1841 if terse:

1842 t = t[:terse]

1843 return self._instr(prec=prec, props=t, **sep_name)

1844

1845 def toTriaxial(self, **name):

1846 '''Convert this ellipsoid to a L{Triaxial_}.

1847

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

1849

1850 @return: A L{Triaxial_} or L{Triaxial} with the C{X} axis

1851 pointing east and C{Z} pointing north.

1852

1853 @see: Method L{Triaxial_.toEllipsoid}.

1854 '''

1855 T = self._triaxial

1856 return T.copy(**name) if name else T

1857

1858 @property_RO

1859 def _triaxial(self):

1860 '''(INTERNAL) Get this ellipsoid's un-/ordered C{Triaxial/_}.

1861 '''

1862 a, b, m = self.a, self.b, _MODS.triaxials

1863 T = m.Triaxial if a > b else m.Triaxial_

1864 return T(a, a, b, name=self.name)

1865

1866 @Property_RO

1867 def volume(self):

1868 '''Get the ellipsoid's I{volume} (C{meter**3}), M{4 / 3 * PI * R3**3}.

1869

1870 @see: C{R3}.

1871 '''

1872 return Meter3(volume=self.a2 * self.b * PI_3 * _4_0)

1873

1874

1875 class Ellipsoid2(Ellipsoid):

1876 '''An L{Ellipsoid} specified by I{equatorial} radius and I{flattening}.

1877 '''

1878 def __init__(self, a, f=None, **name):

1879 '''New L{Ellipsoid2}.

1880

1881 @arg a: Equatorial radius, semi-axis (C{meter}) or a previous

1882 L{Ellipsoid} instance.

1883 @arg f: Flattening: (C{float} < 1.0, negative for I{prolate}),

1884 if B{C{a}} is in C{meter}.

1885 @kwarg name: Optional, unique C{B{name}=NN} (C{str}).

1886

1887 @raise NameError: Ellipsoid with that B{C{name}} already exists.

1888

1889 @raise ValueError: Invalid B{C{a}} or B{C{f}}.

1890

1891 @note: C{abs(B{f}) < EPS} is forced to C{B{f}=0}, I{spherical}.

1892 Negative C{B{f}} produces a I{prolate} ellipsoid.

1893 '''

1894 if f is None and isinstance(a, Ellipsoid):

1895 Ellipsoid.__init__(self, a.a, f =a.f,

1896 b=a.b, f_=a.f_, **name)

1897 else:

1898 Ellipsoid.__init__(self, a, f=f, **name)

1899

1900

1901 def _spherical_a_b(a, b):

1902 '''(INTERNAL) C{True} for spherical or invalid C{a} or C{b}.

1903 '''

1904 return a < EPS0 or b < EPS0 or fabs(a - b) < EPS0

1905

1906

1907 def _spherical_f(f):

1908 '''(INTERNAL) C{True} for spherical or invalid C{f}.

1909 '''

1910 return fabs(f) < EPS or f > EPS1

1911

1912

1913 def _spherical_f_(f_):

1914 '''(INTERNAL) C{True} for spherical or invalid C{f_}.

1915 '''

1916 return fabs(f_) < EPS or fabs(f_) > _1_EPS

1917

1918

1919 def a_b2e(a, b):

1920 '''Return C{e}, the I{1st eccentricity} for a given I{equatorial} and I{polar} radius.

1921

1922 @arg a: Equatorial radius (C{scalar} > 0).

1923 @arg b: Polar radius (C{scalar} > 0).

1924

1925 @return: The I{unsigned}, (1st) eccentricity (C{float} or C{0}),

1926 M{sqrt(1 - (b / a)**2)}.

1927

1928 @note: The result is always I{non-negative} and C{0} for I{near-spherical} ellipsoids.

1929 '''

1930 return Float(e=sqrt(fabs(a_b2e2(a, b)))) # == sqrt(fabs(a - b) * (a + b)) / a)

1931

1932

1933 def a_b2e2(a, b):

1934 '''Return C{e2}, the I{1st eccentricity squared} for a given I{equatorial} and I{polar} radius.

1935

1936 @arg a: Equatorial radius (C{scalar} > 0).

1937 @arg b: Polar radius (C{scalar} > 0).

1938

1939 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or C{0}),

1940 M{1 - (b / a)**2}.

1941

1942 @note: The result is positive for I{oblate}, negative for I{prolate}

1943 or C{0} for I{near-spherical} ellipsoids.

1944 '''

1945 return Float(e2=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / a**2))

1946

1947

1948 def a_b2e22(a, b):

1949 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{equatorial} and I{polar} radius.

1950

1951 @arg a: Equatorial radius (C{scalar} > 0).

1952 @arg b: Polar radius (C{scalar} > 0).

1953

1954 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} or C{0}),

1955 M{(a / b)**2 - 1}.

1956

1957 @note: The result is positive for I{oblate}, negative for I{prolate}

1958 or C{0} for I{near-spherical} ellipsoids.

1959 '''

1960 return Float(e22=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / b**2))

1961

1962

1963 def a_b2e32(a, b):

1964 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{equatorial} and I{polar} radius.

1965

1966 @arg a: Equatorial radius (C{scalar} > 0).

1967 @arg b: Polar radius (C{scalar} > 0).

1968

1969 @return: The I{signed}, 3rd eccentricity I{squared} (C{float} or C{0}),

1970 M{(a**2 - b**2) / (a**2 + b**2)}.

1971

1972 @note: The result is positive for I{oblate}, negative for I{prolate}

1973 or C{0} for I{near-spherical} ellipsoids.

1974 '''

1975 a2, b2 = a**2, b**2

1976 return Float(e32=_0_0 if _spherical_a_b(a2, b2) else ((a2 - b2) / (a2 + b2)))

1977

1978

1979 def a_b2f(a, b):

1980 '''Return C{f}, the I{flattening} for a given I{equatorial} and I{polar} radius.

1981

1982 @arg a: Equatorial radius (C{scalar} > 0).

1983 @arg b: Polar radius (C{scalar} > 0).

1984

1985 @return: The flattening (C{scalar} or C{0}), M{(a - b) / a}.

1986

1987 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

1988 for I{near-spherical} ellipsoids.

1989 '''

1990 f = 0 if _spherical_a_b(a, b) else ((a - b) / a)

1991 return _f_0_0 if _spherical_f(f) else Float(f=f)

1992

1993

1994 def a_b2f_(a, b):

1995 '''Return C{f_}, the I{inverse flattening} for a given I{equatorial} and I{polar} radius.

1996

1997 @arg a: Equatorial radius (C{scalar} > 0).

1998 @arg b: Polar radius (C{scalar} > 0).

1999

2000 @return: The inverse flattening (C{scalar} or C{0}), M{a / (a - b)}.

2001

2002 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2003 for I{near-spherical} ellipsoids.

2004 '''

2005 f_ = 0 if _spherical_a_b(a, b) else (a / float(a - b))

2006 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_)

2007

2008

2009 def a_b2f2(a, b):

2010 '''Return C{f2}, the I{2nd flattening} for a given I{equatorial} and I{polar} radius.

2011

2012 @arg a: Equatorial radius (C{scalar} > 0).

2013 @arg b: Polar radius (C{scalar} > 0).

2014

2015 @return: The I{signed}, 2nd flattening (C{scalar} or C{0}), M{(a - b) / b}.

2016

2017 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2018 for I{near-spherical} ellipsoids.

2019 '''

2020 t = 0 if _spherical_a_b(a, b) else float(a - b)

2021 return Float(f2=_0_0 if fabs(t) < EPS0 else (t / b))

2022

2023

2024 def a_b2n(a, b):

2025 '''Return C{n}, the I{3rd flattening} for a given I{equatorial} and I{polar} radius.

2026

2027 @arg a: Equatorial radius (C{scalar} > 0).

2028 @arg b: Polar radius (C{scalar} > 0).

2029

2030 @return: The I{signed}, 3rd flattening (C{scalar} or C{0}), M{(a - b) / (a + b)}.

2031

2032 @note: The result is positive for I{oblate}, negative for I{prolate}

2033 or C{0} for I{near-spherical} ellipsoids.

2034 '''

2035 t = 0 if _spherical_a_b(a, b) else float(a - b)

2036 return Float(n=_0_0 if fabs(t) < EPS0 else (t / (a + b)))

2037

2038

2039 def a_f2b(a, f):

2040 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{flattening}.

2041

2042 @arg a: Equatorial radius (C{scalar} > 0).

2043 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2044

2045 @return: The polar radius (C{float}), M{a * (1 - f)}.

2046 '''

2047 b = a if _spherical_f(f) else (a * (_1_0 - f))

2048 return Radius_(b=a if _spherical_a_b(a, b) else b)

2049

2050

2051 def a_f_2b(a, f_):

2052 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{inverse flattening}.

2053

2054 @arg a: Equatorial radius (C{scalar} > 0).

2055 @arg f_: Inverse flattening (C{scalar} >>> 1).

2056

2057 @return: The polar radius (C{float}), M{a * (f_ - 1) / f_}.

2058 '''

2059 b = a if _spherical_f_(f_) else (a * (f_ - _1_0) / f_)

2060 return Radius_(b=a if _spherical_a_b(a, b) else b)

2061

2062

2063 def b_f2a(b, f):

2064 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{flattening}.

2065

2066 @arg b: Polar radius (C{scalar} > 0).

2067 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2068

2069 @return: The equatorial radius (C{float}), M{b / (1 - f)}.

2070 '''

2071 t = _1_0 - f

2072 a = b if fabs(t) < EPS0 else (b / t)

2073 return Radius_(a=b if _spherical_a_b(a, b) else a)

2074

2075

2076 def b_f_2a(b, f_):

2077 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{inverse flattening}.

2078

2079 @arg b: Polar radius (C{scalar} > 0).

2080 @arg f_: Inverse flattening (C{scalar} >>> 1).

2081

2082 @return: The equatorial radius (C{float}), M{b * f_ / (f_ - 1)}.

2083 '''

2084 t = f_ - _1_0

2085 a = b if _spherical_f_(f_) or fabs(t - f_) < EPS0 \

2086 or fabs(t) < EPS0 else (b * f_ / t)

2087 return Radius_(a=b if _spherical_a_b(a, b) else a)

2088

2089

2090 def e2f(e):

2091 '''Return C{f}, the I{flattening} for a given I{1st eccentricity}.

2092

2093 @arg e: The (1st) eccentricity (0 <= C{float} < 1)

2094

2095 @return: The flattening (C{scalar} or C{0}).

2096

2097 @see: Function L{e22f}.

2098 '''

2099 return e22f(e**2)

2100

2101

2102 def e22f(e2):

2103 '''Return C{f}, the I{flattening} for a given I{1st eccentricity squared}.

2104

2105 @arg e2: The (1st) eccentricity I{squared}, I{signed} (L{NINF} < C{float} < 1)

2106

2107 @return: The flattening (C{float} or C{0}), M{e2 / (sqrt(e2 - 1) + 1)}.

2108 '''

2109 return Float(f=e2 / (sqrt(_1_0 - e2) + _1_0)) if e2 else _f_0_0

2110

2111

2112 def f2e2(f):

2113 '''Return C{e2}, the I{1st eccentricity squared} for a given I{flattening}.

2114

2115 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2116

2117 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} < 1),

2118 M{f * (2 - f)}.

2119

2120 @note: The result is positive for I{oblate}, negative for I{prolate}

2121 or C{0} for I{near-spherical} ellipsoids.

2122

2123 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2124 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2125 <https://WikiPedia.org/wiki/Flattening>}.

2126 '''

2127 return Float(e2=_0_0 if _spherical_f(f) else (f * (_2_0 - f)))

2128

2129

2130 def f2e22(f):

2131 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{flattening}.

2132

2133 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2134

2135 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} > -1 or

2136 C{INF}), M{f * (2 - f) / (1 - f)**2}.

2137

2138 @note: The result is positive for I{oblate}, negative for I{prolate}

2139 or C{0} for near-spherical ellipsoids.

2140

2141 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2142 C++/doc/classGeographicLib_1_1Ellipsoid.html>}.

2143 '''

2144 # e2 / (1 - e2) == f * (2 - f) / (1 - f)**2

2145 t = (_1_0 - f)**2

2146 return Float(e22=INF if t < EPS0 else (f2e2(f) / t)) # PYCHOK type

2147

2148

2149 def f2e32(f):

2150 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{flattening}.

2151

2152 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2153

2154 @return: The I{signed}, 3rd eccentricity I{squared} (C{float}),

2155 M{f * (2 - f) / (1 + (1 - f)**2)}.

2156

2157 @note: The result is positive for I{oblate}, negative for I{prolate}

2158 or C{0} for I{near-spherical} ellipsoids.

2159

2160 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2161 C++/doc/classGeographicLib_1_1Ellipsoid.html>}.

2162 '''

2163 # e2 / (2 - e2) == f * (2 - f) / (1 + (1 - f)**2)

2164 e2 = f2e2(f)

2165 return Float(e32=e2 / (_2_0 - e2))

2166

2167

2168 def f_2f(f_):

2169 '''Return C{f}, the I{flattening} for a given I{inverse flattening}.

2170

2171 @arg f_: Inverse flattening (C{scalar} >>> 1).

2172

2173 @return: The flattening (C{scalar} or C{0}), M{1 / f_}.

2174

2175 @note: The result is positive for I{oblate}, negative for I{prolate}

2176 or C{0} for I{near-spherical} ellipsoids.

2177 '''

2178 f = 0 if _spherical_f_(f_) else _1_0 / f_

2179 return _f_0_0 if _spherical_f(f) else Float(f=f) # PYCHOK type

2180

2181

2182 def f2f_(f):

2183 '''Return C{f_}, the I{inverse flattening} for a given I{flattening}.

2184

2185 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2186

2187 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2188

2189 @note: The result is positive for I{oblate}, negative for I{prolate}

2190 or C{0} for I{near-spherical} ellipsoids.

2191 '''

2192 f_ = 0 if _spherical_f(f) else _1_0 / f

2193 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) # PYCHOK type

2194

2195

2196 def f2f2(f):

2197 '''Return C{f2}, the I{2nd flattening} for a given I{flattening}.

2198

2199 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2200

2201 @return: The I{signed}, 2nd flattening (C{scalar} or C{INF}), M{f / (1 - f)}.

2202

2203 @note: The result is positive for I{oblate}, negative for I{prolate}

2204 or C{0} for I{near-spherical} ellipsoids.

2205

2206 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2207 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2208 <https://WikiPedia.org/wiki/Flattening>}.

2209 '''

2210 t = _1_0 - f

2211 return Float(f2=_0_0 if _spherical_f(f) else (INF if fabs(t) < EPS

2212 else (f / t))) # PYCHOK type

2213

2214

2215 def f2n(f):

2216 '''Return C{n}, the I{3rd flattening} for a given I{flattening}.

2217

2218 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2219

2220 @return: The I{signed}, 3rd flattening (-1 <= C{float} < 1),

2221 M{f / (2 - f)}.

2222

2223 @note: The result is positive for I{oblate}, negative for I{prolate}

2224 or C{0} for I{near-spherical} ellipsoids.

2225

2226 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2227 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2228 <https://WikiPedia.org/wiki/Flattening>}.

2229 '''

2230 return Float(n=_0_0 if _spherical_f(f) else (f / float(_2_0 - f)))

2231

2232

2233 def n2e2(n):

2234 '''Return C{e2}, the I{1st eccentricity squared} for a given I{3rd flattening}.

2235

2236 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2237

2238 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or NINF),

2239 M{4 * n / (1 + n)**2}.

2240

2241 @note: The result is positive for I{oblate}, negative for I{prolate}

2242 or C{0} for I{near-spherical} ellipsoids.

2243

2244 @see: U{Flattening<https://WikiPedia.org/wiki/Flattening>}.

2245 '''

2246 t = (n + _1_0)**2

2247 return Float(e2=_0_0 if fabs(n) < EPS0 else

2248 (NINF if t < EPS0 else (_4_0 * n / t)))

2249

2250

2251 def n2f(n):

2252 '''Return C{f}, the I{flattening} for a given I{3rd flattening}.

2253

2254 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2255

2256 @return: The flattening (C{scalar} or NINF), M{2 * n / (1 + n)}.

2257

2258 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2259 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2260 <https://WikiPedia.org/wiki/Flattening>}.

2261 '''

2262 t = n + _1_0

2263 f = 0 if fabs(n) < EPS0 else (NINF if t < EPS0 else (_2_0 * n / t))

2264 return _f_0_0 if _spherical_f(f) else Float(f=f)

2265

2266

2267 def n2f_(n):

2268 '''Return C{f_}, the I{inverse flattening} for a given I{3rd flattening}.

2269

2270 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2271

2272 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2273

2274 @see: L{n2f} and L{f2f_}.

2275 '''

2276 return f2f_(n2f(n))

2277

2278

2279 def _plumbTo3(px, py, E, eps=EPS): # in .height4 above

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

2281

2282 @see: Functions C{pygeodesy.triaxial._plumbTo4} and C{-._plumbTo5}.

2283 '''

2284 a, b, e0 = E.a, E.b, EPS0

2285 if min(px, py, a, b) < e0:

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

2287

2288 a2 = a - b * E.b_a

2289 b2 = b - a * E.a_b

2290 tx = ty = _SQRT2_2

2291 _a, _h = fabs, hypot

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

2293 ex = a2 * tx**3

2294 ey = b2 * ty**3

2295

2296 qx = px - ex

2297 qy = py - ey

2298 q = _h(qx, qy)

2299 if q < e0: # PYCHOK no cover

2300 break

2301 r = _h(ex - tx * a,

2302 ey - ty * b) / q

2303

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

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

2306 t = _h(ty, tx)

2307 if t < e0: # PYCHOK no cover

2308 break

2309 tx = tx / t # /= chokes PyChecker

2310 ty = ty / t

2311 if _a(sx - tx) < eps and \

2312 _a(sy - ty) < eps:

2313 break

2314

2315 tx *= a / px

2316 ty *= b / py

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

2318

2319

2320 class Ellipsoids(_NamedEnum):

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

2322 to accommodate the L{_LazyNamedEnumItem} properties.

2323 '''

2324 def _Lazy(self, a, b, f_, **kwds):

2325 '''(INTERNAL) Instantiate the L{Ellipsoid}.

2326 '''

2327 return Ellipsoid(a, b=b, f_=f_, **kwds)

2328

2329 Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

2330'''Some pre-defined L{Ellipsoid}s, all I{lazily} instantiated.'''

2331# <https://www.GNU.org/software/gama/manual/html_node/Supported-ellipsoids.html>

2332# <https://GSSC.ESA.int/navipedia/index.php/Reference_Frames_in_GNSS>

2333# <https://kb.OSU.edu/dspace/handle/1811/77986>

2334# <https://www.IBM.com/docs/en/db2/11.5?topic=systems-supported-spheroids>

2335# <https://w3.Energistics.org/archive/Epicentre/Epicentre_v3.0/DataModel/LogicalDictionary/StandardValues/ellipsoid.html>

2336# <https://GitHub.com/locationtech/proj4j/blob/master/src/main/java/org/locationtech/proj4j/datum/Ellipsoid.java>

2337 Ellipsoids._assert( # <https://WikiPedia.org/wiki/Earth_ellipsoid>

2338 Airy1830 = _lazy(_Airy1830_, *_T(6377563.396, _0_0, 299.3249646)), # b=6356256.909

2339 AiryModified = _lazy(_AiryModified_, *_T(6377340.189, _0_0, 299.3249646)), # b=6356034.448

2340# APL4_9 = _lazy('APL4_9', *_T(6378137.0, _0_0, 298.24985392)), # Appl. Phys. Lab. 1965

2341# ANS = _lazy('ANS', *_T(6378160.0, _0_0, 298.25)), # Australian Nat. Spheroid

2342# AN_SA96 = _lazy('AN_SA96', *_T(6378160.0, _0_0, 298.24985392)), # Australian Nat. South America

2343 Australia1966 = _lazy('Australia1966', *_T(6378160.0, _0_0, 298.25)), # b=6356774.7192

2344 ATS1977 = _lazy('ATS1977', *_T(6378135.0, _0_0, 298.257)), # "Average Terrestrial System"

2345 Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.962818, 299.152812797)),

2346 BesselModified = _lazy('BesselModified', *_T(6377492.018, _0_0, 299.1528128)),

2347# BesselNamibia = _lazy('BesselNamibia', *_T(6377483.865, _0_0, 299.1528128)),

2348 CGCS2000 = _lazy('CGCS2000', *_T(R_MA, _0_0, 298.257222101)), # BeiDou Coord System (BDC)

2349# Clarke1858 = _lazy('Clarke1858', *_T(6378293.639, _0_0, 294.260676369)),

2350 Clarke1866 = _lazy(_Clarke1866_, *_T(6378206.4, 6356583.8, 294.978698214)),

2351 Clarke1880 = _lazy('Clarke1880', *_T(6378249.145, 6356514.86954978, 293.465)),

2352 Clarke1880IGN = _lazy(_Clarke1880IGN_, *_T(6378249.2, 6356515.0, 293.466021294)),

2353 Clarke1880Mod = _lazy('Clarke1880Mod', *_T(6378249.145, 6356514.96639549, 293.466307656)), # aka Clarke1880Arc

2354 CPM1799 = _lazy('CPM1799', *_T(6375738.7, 6356671.92557493, 334.39)), # Comm. des Poids et Mesures

2355 Delambre1810 = _lazy('Delambre1810', *_T(6376428.0, 6355957.92616372, 311.5)), # Belgium

2356 Engelis1985 = _lazy('Engelis1985', *_T(6378136.05, 6356751.32272154, 298.2566)),

2357# Everest1830 = _lazy('Everest1830', *_T(6377276.345, _0_0, 300.801699997)),

2358# Everest1948 = _lazy('Everest1948', *_T(6377304.063, _0_0, 300.801699997)),

2359# Everest1956 = _lazy('Everest1956', *_T(6377301.243, _0_0, 300.801699997)),

2360 Everest1969 = _lazy('Everest1969', *_T(6377295.664, 6356094.667915, 300.801699997)),

2361 Everest1975 = _lazy('Everest1975', *_T(6377299.151, 6356098.14512013, 300.8017255)),

2362 Fisher1968 = _lazy('Fisher1968', *_T(6378150.0, 6356768.33724438, 298.3)),

2363# Fisher1968Mod = _lazy('Fisher1968Mod', *_T(6378155.0, _0_0, 298.3)),

2364 GEM10C = _lazy('GEM10C', *_T(R_MA, 6356752.31424783, 298.2572236)),

2365 GPES = _lazy('GPES', *_T(6378135.0, 6356750.0, _0_0)), # "Gen. Purpose Earth Spheroid"

2366 GRS67 = _lazy('GRS67', *_T(6378160.0, _0_0, 298.247167427)), # Lucerne b=6356774.516

2367# GRS67Truncated = _lazy('GRS67Truncated', *_T(6378160.0, _0_0, 298.25)),

2368 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.25722210088)), # IUGG, ITRS, ETRS89

2369# Hayford1924 = _lazy('Hayford1924', *_T(6378388.0, 6356911.94612795, None)), # aka Intl1924 f_=297

2370 Helmert1906 = _lazy('Helmert1906', *_T(6378200.0, 6356818.16962789, 298.3)),

2371# Hough1960 = _lazy('Hough1960', *_T(6378270.0, _0_0, 297.0)),

2372 IAU76 = _lazy('IAU76', *_T(6378140.0, _0_0, 298.257)), # Int'l Astronomical Union

2373 IERS1989 = _lazy('IERS1989', *_T(6378136.0, _0_0, 298.257)), # b=6356751.302

2374 IERS1992TOPEX = _lazy('IERS1992TOPEX', *_T(6378136.3, 6356751.61659215, 298.257223563)), # IERS/TOPEX/Poseidon/McCarthy

2375 IERS2003 = _lazy('IERS2003', *_T(6378136.6, 6356751.85797165, 298.25642)),

2376 Intl1924 = _lazy(_Intl1924_, *_T(6378388.0, _0_0, 297.0)), # aka Hayford b=6356911.9462795

2377 Intl1967 = _lazy('Intl1967', *_T(6378157.5, 6356772.2, 298.24961539)), # New Int'l

2378 Krassovski1940 = _lazy(_Krassovski1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2379 Krassowsky1940 = _lazy(_Krassowsky1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2380# Kaula = _lazy('Kaula', *_T(6378163.0, _0_0, 298.24)), # Kaula 1961

2381# Lerch = _lazy('Lerch', *_T(6378139.0, _0_0, 298.257)), # Lerch 1979

2382 Maupertuis1738 = _lazy('Maupertuis1738', *_T(6397300.0, 6363806.28272251, 191.0)), # France

2383 Mercury1960 = _lazy('Mercury1960', *_T(6378166.0, 6356784.28360711, 298.3)),

2384 Mercury1968Mod = _lazy('Mercury1968Mod', *_T(6378150.0, 6356768.33724438, 298.3)),

2385# MERIT = _lazy('MERIT', *_T(6378137.0, _0_0, 298.257)), # MERIT 1983

2386# NWL10D = _lazy('NWL10D', *_T(6378135.0, _0_0, 298.26)), # Naval Weapons Lab.

2387 NWL1965 = _lazy('NWL1965', *_T(6378145.0, 6356759.76948868, 298.25)), # Naval Weapons Lab.

2388# NWL9D = _lazy('NWL9D', *_T(6378145.0, 6356759.76948868, 298.25)), # NWL1965

2389 OSU86F = _lazy('OSU86F', *_T(6378136.2, 6356751.51693008, 298.2572236)),

2390 OSU91A = _lazy('OSU91A', *_T(6378136.3, 6356751.6165948, 298.2572236)),

2391# Plessis1817 = _lazy('Plessis1817', *_T(6397523.0, 6355863.0, 153.56512242)), # XXX incorrect?

2392 Plessis1817 = _lazy('Plessis1817', *_T(6376523.0, 6355862.93325557, 308.64)), # XXX IGN France 1972

2393# Prolate = _lazy('Prolate', *_T(6356752.3, R_MA, _0_0)),

2394 PZ90 = _lazy('PZ90', *_T(6378136.0, _0_0, 298.257839303)), # GLOSNASS PZ-90 and PZ-90.11

2395# SEAsia = _lazy('SEAsia', *_T(6378155.0, _0_0, 298.3)), # SouthEast Asia

2396 SGS85 = _lazy('SGS85', *_T(6378136.0, 6356751.30156878, 298.257)), # Soviet Geodetic System

2397 SoAmerican1969 = _lazy('SoAmerican1969', *_T(6378160.0, 6356774.71919531, 298.25)), # South American

2398 Sphere = _lazy(_Sphere_, *_T(R_M, R_M, _0_0)), # pseudo

2399 SphereAuthalic = _lazy('SphereAuthalic', *_T(R_FM, R_FM, _0_0)), # pseudo

2400 SpherePopular = _lazy('SpherePopular', *_T(R_MA, R_MA, _0_0)), # EPSG:3857 Spheroid

2401 Struve1860 = _lazy('Struve1860', *_T(6378298.3, 6356657.14266956, 294.73)),

2402# Walbeck = _lazy('Walbeck', *_T(6376896.0, _0_0, 302.78)),

2403# WarOffice = _lazy('WarOffice', *_T(6378300.0, _0_0, 296.0)),

2404 WGS60 = _lazy('WGS60', *_T(6378165.0, 6356783.28695944, 298.3)),

2405 WGS66 = _lazy('WGS66', *_T(6378145.0, 6356759.76948868, 298.25)),

2406 WGS72 = _lazy(_WGS72_, *_T(6378135.0, _0_0, 298.26)), # b=6356750.52

2407 WGS84 = _lazy(_WGS84_, *_T(R_MA, _0_0, _f__WGS84)), # GPS b=6356752.3142451793

2408# U{NOAA/NOS/NGS/inverse<https://GitHub.com/noaa-ngs/inverse/blob/main/invers3d.f>}

2409 WGS84_NGS = _lazy('WGS84_NGS', *_T(R_MA, _0_0, 298.257222100882711243162836600094))

2410)

2411

2412_EWGS84 = Ellipsoids.WGS84 # (INTERNAL) shared

2413

2414 if __name__ == '__main__':

2415

2416 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

2417 from pygeodesy import nameof, printf

2418

2419 for E in (_EWGS84, Ellipsoids.GRS80, # NAD83,

2420 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

2421 Ellipsoid(_EWGS84.b, _EWGS84.a, name='_Prolate')):

2422 e = f2n(E.f) - E.n

2423 printf('# %s: %s', _DOT_('Ellipsoids', E.name), E.toStr(prec=10, terse=0), nl=1)

2424 printf('# e=%s, f_=%s, f=%s, n=%s (%s)', fstr(E.e, prec=13, fmt=Fmt.e),

2425 fstr(E.f_, prec=13, fmt=Fmt.e),

2426 fstr(E.f, prec=13, fmt=Fmt.e),

2427 fstr(E.n, prec=13, fmt=Fmt.e),

2428 fstr(e, prec=9, fmt=Fmt.e))

2429 printf('# %s %s', Ellipsoid.AlphaKs.name, fstr(E.AlphaKs, prec=20))

2430 printf('# %s %s', Ellipsoid.BetaKs.name, fstr(E.BetaKs, prec=20))

2431 printf('# %s %s', nameof(Ellipsoid.KsOrder), E.KsOrder) # property

2432

2433 # __doc__ of this file, force all into registry

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

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

2436

2437# % python3 -m pygeodesy.ellipsoids

2438

2439# Ellipsoids.WGS84: name='WGS84', a=6378137, f=0.0033528107, f_=298.257223563, b=6356752.3142451793, f2=0.0033640898, n=0.0016792204, e=0.0818191908, e2=0.00669438, e21=0.99330562, e22=0.0067394967, e32=0.0033584313, A=6367449.1458234144, L=10001965.7293127235, R1=6371008.7714150595, R2=6371007.1809184738, R3=6371000.7900091587, Rbiaxial=6367453.6345163295, Rtriaxial=6372797.5559594007

2440# e=8.1819190842622e-02, f_=2.98257223563e+02, f=3.3528106647475e-03, n=1.6792203863837e-03 (0.0e+00)

2441# AlphaKs 0.00083773182062446994, 0.00000076085277735725, 0.00000000119764550324, 0.00000000000242917068, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2442# BetaKs 0.00083773216405794875, 0.0000000590587015222, 0.00000000016734826653, 0.00000000000021647981, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2443# KsOrder 8

2444

2445# Ellipsoids.GRS80: name='GRS80', a=6378137, f=0.0033528107, f_=298.2572221009, b=6356752.3141403468, f2=0.0033640898, n=0.0016792204, e=0.081819191, e2=0.00669438, e21=0.99330562, e22=0.0067394968, e32=0.0033584313, A=6367449.1457710434, L=10001965.7292304561, R1=6371008.7713801153, R2=6371007.1808835147, R3=6371000.7899741363, Rbiaxial=6367453.6344640013, Rtriaxial=6372797.5559332585

2446# e=8.1819191042833e-02, f_=2.9825722210088e+02, f=3.3528106811837e-03, n=1.6792203946295e-03 (0.0e+00)

2447# AlphaKs 0.00083773182472890429, 0.00000076085278481561, 0.00000000119764552086, 0.00000000000242917073, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2448# BetaKs 0.0008377321681623882, 0.00000005905870210374, 0.000000000167348269, 0.00000000000021647982, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2449# KsOrder 8

2450

2451# Ellipsoids.Sphere: name='Sphere', a=6371008.7714149999, f=0, f_=0, b=6371008.7714149999, f2=0, n=0, e=0, e2=0, e21=1, e22=0, e32=0, A=6371008.7714149999, L=10007557.1761167478, R1=6371008.7714149999, R2=6371008.7714149999, R3=6371008.7714149999, Rbiaxial=6371008.7714149999, Rtriaxial=6371008.7714149999

2452# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2453# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2454# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2455# KsOrder 8

2456

2457# Ellipsoids.SpherePopular: name='SpherePopular', a=6378137, f=0, f_=0, b=6378137, f2=0, n=0, e=0, e2=0, e21=1, e22=0, e32=0, A=6378137, L=10018754.171394622, R1=6378137, R2=6378137, R3=6378137, Rbiaxial=6378137, Rtriaxial=6378137

2458# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2459# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2460# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2461# KsOrder 8

2462

2463# Ellipsoids._Prolate: name='_Prolate', a=6356752.3142451793, f=-0.0033640898, f_=-297.257223563, b=6378137, f2=-0.0033528107, n=-0.0016792204, e=0.0820944379, e2=-0.0067394967, e21=1.0067394967, e22=-0.00669438, e32=-0.0033584313, A=6367449.1458234144, L=10035500.5204500314, R1=6363880.5428301189, R2=6363878.9413582645, R3=6363872.5644020075, Rbiaxial=6367453.6345163295, Rtriaxial=6362105.2243882557

2464# e=8.2094437949696e-02, f_=-2.97257223563e+02, f=-3.3640898209765e-03, n=-1.6792203863837e-03 (0.0e+00)

2465# AlphaKs -0.00084149152514366627, 0.00000076653480614871, -0.00000000120934503389, 0.0000000000024576225, -0.00000000000000578863, 0.00000000000000001502, -0.00000000000000000004, 0.0

2466# BetaKs -0.00084149187224351817, 0.00000005842735196773, -0.0000000001680487236, 0.00000000000021706261, -0.00000000000000038002, 0.00000000000000000073, -0.0, 0.0

2467# KsOrder 8

2468

2469# **) MIT License

2470#

2471# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

2472#

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

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

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

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

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

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

2479#

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

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

2482#

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

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

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

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

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

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

2489# OTHER DEALINGS IN THE SOFTWARE.