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

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, property_ROver

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 # as _tan

94

95__all__ = _ALL_LAZY.ellipsoids

96__version__ = '24.11.26'

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

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

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

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

667

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

669

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

671

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

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

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

675 and method C{Rlat}.

676 '''

677 lat = Lat(lat)

678 if lat:

679 b = lat

680 if fabs(lat) < _90_0:

681 if self.f:

682 b = self._beta(lat)

683 z, r = sincos2d(b)

684 r *= self.a

685 z *= self.b

686 else: # near-polar

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

688 else: # equator

689 r = self.a

690 z = lat = b = _0_0

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

692

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

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

695 along a parallel of (geodetic) latitude.

696

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

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

699

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

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

702

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

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

705

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

707 '''

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

709

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

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

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

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

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

715

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

717 and only between points not near-polar}.

718

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

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

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

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

723

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

725 in same units as this ellipsoid's axes.

726

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

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

729

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

731 '''

732 phi0 = Phid(lat0=lat0)

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

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

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

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

737

738 @Property_RO

739 def e(self):

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

741

742 @see: Property L{es}.

743 '''

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

745

746 @deprecated_Property_RO

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

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

749 return self.e21

750

751# @Property_RO

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

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

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

755

756 @Property_RO

757 def e2(self):

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

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

760 '''

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

762

763 @Property_RO

764 def e2abs(self):

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

766 '''

767 return fabs(self.e2)

768

769 @Property_RO

770 def e21(self):

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

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

773 '''

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

775

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

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

778

779 @Property_RO

780 def e22(self):

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

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

783 '''

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

785

786 @Property_RO

787 def e22abs(self):

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

789 '''

790 return fabs(self.e22)

791

792 @Property_RO

793 def e32(self):

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

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

796 '''

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

798

799 @Property_RO

800 def e32abs(self):

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

802 '''

803 return fabs(self.e32)

804

805 @Property_RO

806 def e4(self):

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

808 '''

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

810

811 eccentricity = e # eccentricity

812# eccentricity2 = e2 # eccentricity squared

813 eccentricity1st2 = e2 # first eccentricity squared, signed

814 eccentricity2nd2 = e22 # second eccentricity squared, signed

815 eccentricity3rd2 = e32 # third eccentricity squared, signed

816

817 def ecef(self, Ecef=None):

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

819

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

821

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

823

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

825

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

827 '''

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

829

830 @Property_RO

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

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

833 '''

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

835 return _MODS.elliptic.Elliptic(e12)

836

837 @Property_RO

838 def _elliptic_e22(self): # aka ._elliptic_ep2

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

840 '''

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

842

843 equatoradius = a # Requatorial

844

845 def e2s(self, s):

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

847

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

849

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

851

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

853 '''

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

855

856 def e2s2(self, s):

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

858

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

860

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

862

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

864 '''

865 r = _1_0

866 if self.e2:

867 try:

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

869 if r < 0:

870 raise ValueError(_negative_)

871 except (TypeError, ValueError) as x:

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

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

874 return r

875

876 @Property_RO

877 def es(self):

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

879

880 @see: Property L{e}.

881 '''

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

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

884

885 def es_atanh(self, x):

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

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

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

889

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

891

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

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

894 '''

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

896

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

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

899 '''

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

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

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

903 _0_0)) # .isSpherical

904

905 @Property_RO

906 def es_c(self):

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

908 '''

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

910

911 def es_tauf(self, taup):

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

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

914

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

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

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

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

919 '''

920 t = Scalar(taup=taup)

921 if self.f: # .isEllipsoidal

922 a = fabs(t)

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

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

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

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

927 if fabs(d) < s:

928 break

929 t = Scalar(tauf=T)

930 return t

931

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

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

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

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

936 '''

937 e = self._1_e21

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

939 _tf2 = self._es_taupf2

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

941 a, h = _tf2(T)

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

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

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

945 yield T, i, d

946

947 def es_taupf(self, tau):

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

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

950

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

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

953 '''

954 t = Scalar(tau=tau)

955 if self.f: # .isEllipsoidal

956 t, _ = self._es_taupf2(t)

957 t = Scalar(taupf=t)

958 return t

959

960 def _es_taupf2(self, tau):

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

962 '''

963 if _isfinite(tau):

964 h = hypot1(tau)

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

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

967 else:

968 a, h = tau, INF

969 return a, h

970

971 @Property_RO

972 def _exp_es_atanh_1(self):

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

974 '''

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

976

977 @Property_RO

978 def f(self):

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

980 '''

981 return self._f

982

983 @Property_RO

984 def f_(self):

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

986 '''

987 return self._f_

988

989 @Property_RO

990 def f1(self):

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

992

993 @see: Property L{b_a}.

994 '''

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

996

997 @Property_RO

998 def f2(self):

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

1000 '''

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

1002

1003 @deprecated_Property_RO

1004 def geodesic(self):

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

1006 return self.geodesicw

1007

1008 def geodesic_(self, exact=True):

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

1010

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

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

1013 instance for I{this} ellipsoid.

1014

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

1016

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

1018

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

1020 '''

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

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

1023 name=self.name)

1024 else:

1025 g = exact

1026 E = _xattr(g, ellipsoid=None)

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

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

1029 return g

1030

1031 @property_ROver

1032 def _Geodesics(self):

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

1034 '''

1035 t = (_MODS.geodesicx.GeodesicExact,

1036 _MODS.geodsolve.GeodesicSolve)

1037 try:

1038 t += (_MODS.geodesicw.Geodesic,

1039 _MODS.geodesicw._wrapped.Geodesic)

1040 except ImportError:

1041 pass

1042 return t # overwrite property_ROver

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 = _MODS.triaxials._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_ROver

1521 def _Rhumbs(self):

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

1523 '''

1524 r = _MODS.rhumb

1525 return (r.aux_.RhumbAux, # overwrite property_ROver

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

1527

1528 @property

1529 def rhumbsolve(self):

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

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

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

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

1534 '''

1535 # if not self.isEllipsoidal:

1536 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1538

1539 @rhumbsolve.setter # PYCHOK setter!

1540 def rhumbsolve(self, path):

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

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

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

1544 '''

1545 self._rhumbsolve = path

1546

1547 @deprecated_property_RO

1548 def rhumbx(self):

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

1550 return self.rhumbekx

1551

1552 def Rlat(self, lat):

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

1554

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

1556

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

1558

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

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

1561

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

1563

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

1565

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

1567

1568 @see: Method C{circle4}.

1569 '''

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

1571 r = self.a

1572 if self.f and lat: # .isEllipsoidal

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

1574 r = Radius(Rlat=r)

1575 return r

1576

1577 Rpolar = b # for consistent naming

1578

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

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

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

1582

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

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

1585 to use an alternate formula.

1586

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

1588

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

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

1591

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

1593 '''

1594 if not self.f: # .isSpherical

1595 n = self.a

1596 elif ca is None:

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

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

1599 elif ca: # derived from EcefYou.forward

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

1601 n = self.a / h

1602 elif sa:

1603 n = self.a2_b / fabs(sa)

1604 else:

1605 n = self.a

1606 return n

1607

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

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

1610 radii of curvature of (geodetic) latitude.

1611

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

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

1614

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

1616 the radii of curvature.

1617

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

1619

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

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

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

1623 Earth_radius#Radii_of_curvature>}.

1624 '''

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

1626

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

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

1629 curvature of (geodetic) latitude.

1630

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

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

1633

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

1635 radii of curvature.

1636

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

1638

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

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

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

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

1643 '''

1644 a = fabs(Phi(phi))

1645 if self.f:

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

1647 if r > EPS02:

1648 n = self.a / sqrt(r)

1649 m = n * self.e21 / r

1650 else:

1651 m = n = _0_0

1652 else:

1653 m = n = self.a

1654 if scaled and a:

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

1656 return Curvature2Tuple(Radius(rocMeridional=m),

1657 Radius(rocPrimeVertical=n))

1658

1659 def rocBearing(self, lat, bearing):

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

1661 latitude and compass direction.

1662

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

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

1665

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

1667

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

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

1670

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

1672

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

1674 '''

1675 if self.f:

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

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

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

1679 c2 *= n / m

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

1681 s2 *= m / n

1682 n = m

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

1684 else:

1685 b = self.b # == self.a

1686 return Radius(rocBearing=b)

1687

1688 @Property_RO

1689 def rocEquatorial2(self):

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

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

1692

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

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

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

1696 '''

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

1698 Radius(rocPrimeVertical0=self.a))

1699

1700 def rocGauss(self, lat):

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

1702

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

1704

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

1706

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

1708

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

1710 Earth_radius#Radii_of_curvature>}

1711 '''

1712 # using ...

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

1714 # return sqrt(m * n)

1715 # ... requires 1 or 2 sqrt

1716 g = self.b

1717 if self.f:

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

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

1720 return Radius(rocGauss=g)

1721

1722 def rocMean(self, lat):

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

1724

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

1726

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

1728

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

1730

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

1732 Earth_radius#Radii_of_curvature>}

1733 '''

1734 if self.f:

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

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

1737 else:

1738 m = self.a

1739 return Radius(rocMean=m)

1740

1741 def rocMeridional(self, lat):

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

1743

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

1745

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

1747

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

1749

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

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

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

1753 '''

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

1755 self.rocEquatorial2.meridional

1756

1757 rocPolar = a2_b # synonymous

1758

1759 def rocPrimeVertical(self, lat):

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

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

1762

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

1764

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

1766

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

1768

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

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

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

1772 Earth_radius#Radii_of_curvature>}.

1773 '''

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

1775 self.rocEquatorial2.prime_vertical

1776

1777 rocTransverse = rocPrimeVertical # synonymous

1778

1779 @deprecated_Property_RO

1780 def Rquadratic(self): # PYCHOK no cover

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

1782 return self.Rbiaxial

1783

1784 @deprecated_Property_RO

1785 def Rr(self): # PYCHOK no cover

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

1787 return self.Rrectifying

1788

1789 @Property_RO

1790 def Rrectifying(self):

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

1792

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

1794 '''

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

1796 return Radius(Rrectifying=r)

1797

1798 @deprecated_Property_RO

1799 def Rs(self): # PYCHOK no cover

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

1801 return self.Rgeometric

1802

1803 @Property_RO

1804 def Rtriaxial(self):

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

1806

1807 @see: C{Rbiaxial}

1808 '''

1809 a, b = self.a, self.b

1810 q = (sqrt((_3_0 + self.b2_a2) * _0_25) * a) if a > b else (

1811 (sqrt((_3_0 * self.a2_b2 + _1_0) * _0_25) * b) if a < b else a)

1812 return Radius(Rtriaxial=q)

1813

1814 def toEllipsoid2(self, **name):

1815 '''Get a copy of this ellipsoid as an L{Ellipsoid2}.

1816

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

1818

1819 @see: Property C{a_f}.

1820 '''

1821 return Ellipsoid2(self, None, **name)

1822

1823 def toStr(self, prec=8, terse=4, **sep_name): # PYCHOK expected

1824 '''Return this ellipsoid as a text string.

1825

1826 @kwarg prec: Number of decimal digits, unstripped (C{int}).

1827 @kwarg terse: Limit the number of items (C{int}, 0...18),

1828 use C{B{terse}=0} or C{=None} for all.

1829 @kwarg sep_name: Optional C{B{name}=NN} (C{str}) or C{None}

1830 to exclude this ellipsoid's name and separator

1831 C{B{sep}=", "} to join the items (C{str}).

1832

1833 @return: This C{Ellipsoid}'s attributes (C{str}).

1834 '''

1835 E = Ellipsoid

1836 t = (E.a, E.f, E.f_, E.b, E.f2, E.n, E.e,

1837 E.e2, E.e21, E.e22, E.e32,

1838 E.A, E.L, E.R1, E.R2, E.R3,

1839 E.Rbiaxial, E.Rtriaxial)

1840 if terse:

1841 t = t[:terse]

1842 return self._instr(prec=prec, props=t, **sep_name)

1843

1844 def toTriaxial(self, **name):

1845 '''Convert this ellipsoid to a L{Triaxial_}.

1846

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

1848

1849 @return: A L{Triaxial_} or L{Triaxial} with the C{X} axis

1850 pointing east and C{Z} pointing north.

1851

1852 @see: Method L{Triaxial_.toEllipsoid}.

1853 '''

1854 T = self._triaxial

1855 return T.copy(**name) if name else T

1856

1857 @property_RO

1858 def _triaxial(self):

1859 '''(INTERNAL) Get this ellipsoid's un-/ordered C{Triaxial/_}.

1860 '''

1861 a, b, m = self.a, self.b, _MODS.triaxials

1862 T = m.Triaxial if a > b else m.Triaxial_

1863 return T(a, a, b, name=self.name)

1864

1865 @Property_RO

1866 def volume(self):

1867 '''Get the ellipsoid's I{volume} (C{meter**3}), M{4 / 3 * PI * R3**3}.

1868

1869 @see: C{R3}.

1870 '''

1871 return Meter3(volume=self.a2 * self.b * PI_3 * _4_0)

1872

1873

1874 class Ellipsoid2(Ellipsoid):

1875 '''An L{Ellipsoid} specified by I{equatorial} radius and I{flattening}.

1876 '''

1877 def __init__(self, a, f=None, **name):

1878 '''New L{Ellipsoid2}.

1879

1880 @arg a: Equatorial radius, semi-axis (C{meter}) or a previous

1881 L{Ellipsoid} instance.

1882 @arg f: Flattening: (C{float} < 1.0, negative for I{prolate}),

1883 if B{C{a}} is in C{meter}.

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

1885

1886 @raise NameError: Ellipsoid with that B{C{name}} already exists.

1887

1888 @raise ValueError: Invalid B{C{a}} or B{C{f}}.

1889

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

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

1892 '''

1893 if f is None and isinstance(a, Ellipsoid):

1894 Ellipsoid.__init__(self, a.a, f =a.f,

1895 b=a.b, f_=a.f_, **name)

1896 else:

1897 Ellipsoid.__init__(self, a, f=f, **name)

1898

1899

1900 def _spherical_a_b(a, b):

1901 '''(INTERNAL) C{True} for spherical or invalid C{a} or C{b}.

1902 '''

1903 return a < EPS0 or b < EPS0 or fabs(a - b) < EPS0

1904

1905

1906 def _spherical_f(f):

1907 '''(INTERNAL) C{True} for spherical or invalid C{f}.

1908 '''

1909 return fabs(f) < EPS or f > EPS1

1910

1911

1912 def _spherical_f_(f_):

1913 '''(INTERNAL) C{True} for spherical or invalid C{f_}.

1914 '''

1915 return fabs(f_) < EPS or fabs(f_) > _1_EPS

1916

1917

1918 def a_b2e(a, b):

1919 '''Return C{e}, the I{1st eccentricity} for a given I{equatorial} and I{polar} radius.

1920

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

1922 @arg b: Polar radius (C{scalar} > 0).

1923

1924 @return: The I{unsigned}, (1st) eccentricity (C{float} or C{0}),

1925 M{sqrt(1 - (b / a)**2)}.

1926

1927 @note: The result is always I{non-negative} and C{0} for I{near-spherical} ellipsoids.

1928 '''

1929 return Float(e=sqrt(fabs(a_b2e2(a, b)))) # == sqrt(fabs(a - b) * (a + b)) / a)

1930

1931

1932 def a_b2e2(a, b):

1933 '''Return C{e2}, the I{1st eccentricity squared} for a given I{equatorial} and I{polar} radius.

1934

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

1936 @arg b: Polar radius (C{scalar} > 0).

1937

1938 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or C{0}),

1939 M{1 - (b / a)**2}.

1940

1941 @note: The result is positive for I{oblate}, negative for I{prolate}

1942 or C{0} for I{near-spherical} ellipsoids.

1943 '''

1944 return Float(e2=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / a**2))

1945

1946

1947 def a_b2e22(a, b):

1948 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{equatorial} and I{polar} radius.

1949

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

1951 @arg b: Polar radius (C{scalar} > 0).

1952

1953 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} or C{0}),

1954 M{(a / b)**2 - 1}.

1955

1956 @note: The result is positive for I{oblate}, negative for I{prolate}

1957 or C{0} for I{near-spherical} ellipsoids.

1958 '''

1959 return Float(e22=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / b**2))

1960

1961

1962 def a_b2e32(a, b):

1963 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{equatorial} and I{polar} radius.

1964

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

1966 @arg b: Polar radius (C{scalar} > 0).

1967

1968 @return: The I{signed}, 3rd eccentricity I{squared} (C{float} or C{0}),

1969 M{(a**2 - b**2) / (a**2 + b**2)}.

1970

1971 @note: The result is positive for I{oblate}, negative for I{prolate}

1972 or C{0} for I{near-spherical} ellipsoids.

1973 '''

1974 a2, b2 = a**2, b**2

1975 return Float(e32=_0_0 if _spherical_a_b(a2, b2) else ((a2 - b2) / (a2 + b2)))

1976

1977

1978 def a_b2f(a, b):

1979 '''Return C{f}, the I{flattening} for a given I{equatorial} and I{polar} radius.

1980

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

1982 @arg b: Polar radius (C{scalar} > 0).

1983

1984 @return: The flattening (C{scalar} or C{0}), M{(a - b) / a}.

1985

1986 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

1987 for I{near-spherical} ellipsoids.

1988 '''

1989 f = 0 if _spherical_a_b(a, b) else ((a - b) / a)

1990 return _f_0_0 if _spherical_f(f) else Float(f=f)

1991

1992

1993 def a_b2f_(a, b):

1994 '''Return C{f_}, the I{inverse flattening} for a given I{equatorial} and I{polar} radius.

1995

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

1997 @arg b: Polar radius (C{scalar} > 0).

1998

1999 @return: The inverse flattening (C{scalar} or C{0}), M{a / (a - b)}.

2000

2001 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2002 for I{near-spherical} ellipsoids.

2003 '''

2004 f_ = 0 if _spherical_a_b(a, b) else (a / float(a - b))

2005 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_)

2006

2007

2008 def a_b2f2(a, b):

2009 '''Return C{f2}, the I{2nd flattening} for a given I{equatorial} and I{polar} radius.

2010

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

2012 @arg b: Polar radius (C{scalar} > 0).

2013

2014 @return: The I{signed}, 2nd flattening (C{scalar} or C{0}), M{(a - b) / b}.

2015

2016 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2017 for I{near-spherical} ellipsoids.

2018 '''

2019 t = 0 if _spherical_a_b(a, b) else float(a - b)

2020 return Float(f2=_0_0 if fabs(t) < EPS0 else (t / b))

2021

2022

2023 def a_b2n(a, b):

2024 '''Return C{n}, the I{3rd flattening} for a given I{equatorial} and I{polar} radius.

2025

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

2027 @arg b: Polar radius (C{scalar} > 0).

2028

2029 @return: The I{signed}, 3rd flattening (C{scalar} or C{0}), M{(a - b) / (a + b)}.

2030

2031 @note: The result is positive for I{oblate}, negative for I{prolate}

2032 or C{0} for I{near-spherical} ellipsoids.

2033 '''

2034 t = 0 if _spherical_a_b(a, b) else float(a - b)

2035 return Float(n=_0_0 if fabs(t) < EPS0 else (t / (a + b)))

2036

2037

2038 def a_f2b(a, f):

2039 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{flattening}.

2040

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

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

2043

2044 @return: The polar radius (C{float}), M{a * (1 - f)}.

2045 '''

2046 b = a if _spherical_f(f) else (a * (_1_0 - f))

2047 return Radius_(b=a if _spherical_a_b(a, b) else b)

2048

2049

2050 def a_f_2b(a, f_):

2051 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{inverse flattening}.

2052

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

2054 @arg f_: Inverse flattening (C{scalar} >>> 1).

2055

2056 @return: The polar radius (C{float}), M{a * (f_ - 1) / f_}.

2057 '''

2058 b = a if _spherical_f_(f_) else (a * (f_ - _1_0) / f_)

2059 return Radius_(b=a if _spherical_a_b(a, b) else b)

2060

2061

2062 def b_f2a(b, f):

2063 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{flattening}.

2064

2065 @arg b: Polar radius (C{scalar} > 0).

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

2067

2068 @return: The equatorial radius (C{float}), M{b / (1 - f)}.

2069 '''

2070 t = _1_0 - f

2071 a = b if fabs(t) < EPS0 else (b / t)

2072 return Radius_(a=b if _spherical_a_b(a, b) else a)

2073

2074

2075 def b_f_2a(b, f_):

2076 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{inverse flattening}.

2077

2078 @arg b: Polar radius (C{scalar} > 0).

2079 @arg f_: Inverse flattening (C{scalar} >>> 1).

2080

2081 @return: The equatorial radius (C{float}), M{b * f_ / (f_ - 1)}.

2082 '''

2083 t = f_ - _1_0

2084 a = b if _spherical_f_(f_) or fabs(t - f_) < EPS0 \

2085 or fabs(t) < EPS0 else (b * f_ / t)

2086 return Radius_(a=b if _spherical_a_b(a, b) else a)

2087

2088

2089 def e2f(e):

2090 '''Return C{f}, the I{flattening} for a given I{1st eccentricity}.

2091

2092 @arg e: The (1st) eccentricity (0 <= C{float} < 1)

2093

2094 @return: The flattening (C{scalar} or C{0}).

2095

2096 @see: Function L{e22f}.

2097 '''

2098 return e22f(e**2)

2099

2100

2101 def e22f(e2):

2102 '''Return C{f}, the I{flattening} for a given I{1st eccentricity squared}.

2103

2104 @arg e2: The (1st) eccentricity I{squared}, I{signed} (L{NINF} < C{float} < 1)

2105

2106 @return: The flattening (C{float} or C{0}), M{e2 / (sqrt(e2 - 1) + 1)}.

2107 '''

2108 return Float(f=e2 / (sqrt(_1_0 - e2) + _1_0)) if e2 else _f_0_0

2109

2110

2111 def f2e2(f):

2112 '''Return C{e2}, the I{1st eccentricity squared} for a given I{flattening}.

2113

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

2115

2116 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} < 1),

2117 M{f * (2 - f)}.

2118

2119 @note: The result is positive for I{oblate}, negative for I{prolate}

2120 or C{0} for I{near-spherical} ellipsoids.

2121

2122 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2123 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2124 <https://WikiPedia.org/wiki/Flattening>}.

2125 '''

2126 return Float(e2=_0_0 if _spherical_f(f) else (f * (_2_0 - f)))

2127

2128

2129 def f2e22(f):

2130 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{flattening}.

2131

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

2133

2134 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} > -1 or

2135 C{INF}), M{f * (2 - f) / (1 - f)**2}.

2136

2137 @note: The result is positive for I{oblate}, negative for I{prolate}

2138 or C{0} for near-spherical ellipsoids.

2139

2140 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2141 C++/doc/classGeographicLib_1_1Ellipsoid.html>}.

2142 '''

2143 # e2 / (1 - e2) == f * (2 - f) / (1 - f)**2

2144 t = (_1_0 - f)**2

2145 return Float(e22=INF if t < EPS0 else (f2e2(f) / t)) # PYCHOK type

2146

2147

2148 def f2e32(f):

2149 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{flattening}.

2150

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

2152

2153 @return: The I{signed}, 3rd eccentricity I{squared} (C{float}),

2154 M{f * (2 - f) / (1 + (1 - f)**2)}.

2155

2156 @note: The result is positive for I{oblate}, negative for I{prolate}

2157 or C{0} for I{near-spherical} ellipsoids.

2158

2159 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2160 C++/doc/classGeographicLib_1_1Ellipsoid.html>}.

2161 '''

2162 # e2 / (2 - e2) == f * (2 - f) / (1 + (1 - f)**2)

2163 e2 = f2e2(f)

2164 return Float(e32=e2 / (_2_0 - e2))

2165

2166

2167 def f_2f(f_):

2168 '''Return C{f}, the I{flattening} for a given I{inverse flattening}.

2169

2170 @arg f_: Inverse flattening (C{scalar} >>> 1).

2171

2172 @return: The flattening (C{scalar} or C{0}), M{1 / f_}.

2173

2174 @note: The result is positive for I{oblate}, negative for I{prolate}

2175 or C{0} for I{near-spherical} ellipsoids.

2176 '''

2177 f = 0 if _spherical_f_(f_) else _1_0 / f_

2178 return _f_0_0 if _spherical_f(f) else Float(f=f) # PYCHOK type

2179

2180

2181 def f2f_(f):

2182 '''Return C{f_}, the I{inverse flattening} for a given I{flattening}.

2183

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

2185

2186 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2187

2188 @note: The result is positive for I{oblate}, negative for I{prolate}

2189 or C{0} for I{near-spherical} ellipsoids.

2190 '''

2191 f_ = 0 if _spherical_f(f) else _1_0 / f

2192 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) # PYCHOK type

2193

2194

2195 def f2f2(f):

2196 '''Return C{f2}, the I{2nd flattening} for a given I{flattening}.

2197

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

2199

2200 @return: The I{signed}, 2nd flattening (C{scalar} or C{INF}), M{f / (1 - f)}.

2201

2202 @note: The result is positive for I{oblate}, negative for I{prolate}

2203 or C{0} for I{near-spherical} ellipsoids.

2204

2205 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2206 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2207 <https://WikiPedia.org/wiki/Flattening>}.

2208 '''

2209 t = _1_0 - f

2210 return Float(f2=_0_0 if _spherical_f(f) else (INF if fabs(t) < EPS

2211 else (f / t))) # PYCHOK type

2212

2213

2214 def f2n(f):

2215 '''Return C{n}, the I{3rd flattening} for a given I{flattening}.

2216

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

2218

2219 @return: The I{signed}, 3rd flattening (-1 <= C{float} < 1),

2220 M{f / (2 - f)}.

2221

2222 @note: The result is positive for I{oblate}, negative for I{prolate}

2223 or C{0} for I{near-spherical} ellipsoids.

2224

2225 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2226 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2227 <https://WikiPedia.org/wiki/Flattening>}.

2228 '''

2229 return Float(n=_0_0 if _spherical_f(f) else (f / float(_2_0 - f)))

2230

2231

2232 def n2e2(n):

2233 '''Return C{e2}, the I{1st eccentricity squared} for a given I{3rd flattening}.

2234

2235 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2236

2237 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or NINF),

2238 M{4 * n / (1 + n)**2}.

2239

2240 @note: The result is positive for I{oblate}, negative for I{prolate}

2241 or C{0} for I{near-spherical} ellipsoids.

2242

2243 @see: U{Flattening<https://WikiPedia.org/wiki/Flattening>}.

2244 '''

2245 t = (n + _1_0)**2

2246 return Float(e2=_0_0 if fabs(n) < EPS0 else

2247 (NINF if t < EPS0 else (_4_0 * n / t)))

2248

2249

2250 def n2f(n):

2251 '''Return C{f}, the I{flattening} for a given I{3rd flattening}.

2252

2253 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2254

2255 @return: The flattening (C{scalar} or NINF), M{2 * n / (1 + n)}.

2256

2257 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2258 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2259 <https://WikiPedia.org/wiki/Flattening>}.

2260 '''

2261 t = n + _1_0

2262 f = 0 if fabs(n) < EPS0 else (NINF if t < EPS0 else (_2_0 * n / t))

2263 return _f_0_0 if _spherical_f(f) else Float(f=f)

2264

2265

2266 def n2f_(n):

2267 '''Return C{f_}, the I{inverse flattening} for a given I{3rd flattening}.

2268

2269 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2270

2271 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2272

2273 @see: L{n2f} and L{f2f_}.

2274 '''

2275 return f2f_(n2f(n))

2276

2277

2278 class Ellipsoids(_NamedEnum):

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

2280 to accommodate the L{_LazyNamedEnumItem} properties.

2281 '''

2282 def _Lazy(self, a, b, f_, **kwds):

2283 '''(INTERNAL) Instantiate the L{Ellipsoid}.

2284 '''

2285 return Ellipsoid(a, b=b, f_=f_, **kwds)

2286

2287 Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

2288'''Some pre-defined L{Ellipsoid}s, all I{lazily} instantiated.'''

2289# <https://www.GNU.org/software/gama/manual/html_node/Supported-ellipsoids.html>

2290# <https://GSSC.ESA.int/navipedia/index.php/Reference_Frames_in_GNSS>

2291# <https://kb.OSU.edu/dspace/handle/1811/77986>

2292# <https://www.IBM.com/docs/en/db2/11.5?topic=systems-supported-spheroids>

2293# <https://w3.Energistics.org/archive/Epicentre/Epicentre_v3.0/DataModel/LogicalDictionary/StandardValues/ellipsoid.html>

2294# <https://GitHub.com/locationtech/proj4j/blob/master/src/main/java/org/locationtech/proj4j/datum/Ellipsoid.java>

2295 Ellipsoids._assert( # <https://WikiPedia.org/wiki/Earth_ellipsoid>

2296 Airy1830 = _lazy(_Airy1830_, *_T(6377563.396, _0_0, 299.3249646)), # b=6356256.909

2297 AiryModified = _lazy(_AiryModified_, *_T(6377340.189, _0_0, 299.3249646)), # b=6356034.448

2298# APL4_9 = _lazy('APL4_9', *_T(6378137.0, _0_0, 298.24985392)), # Appl. Phys. Lab. 1965

2299# ANS = _lazy('ANS', *_T(6378160.0, _0_0, 298.25)), # Australian Nat. Spheroid

2300# AN_SA96 = _lazy('AN_SA96', *_T(6378160.0, _0_0, 298.24985392)), # Australian Nat. South America

2301 Australia1966 = _lazy('Australia1966', *_T(6378160.0, _0_0, 298.25)), # b=6356774.7192

2302 ATS1977 = _lazy('ATS1977', *_T(6378135.0, _0_0, 298.257)), # "Average Terrestrial System"

2303 Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.962818, 299.152812797)),

2304 BesselModified = _lazy('BesselModified', *_T(6377492.018, _0_0, 299.1528128)),

2305# BesselNamibia = _lazy('BesselNamibia', *_T(6377483.865, _0_0, 299.1528128)),

2306 CGCS2000 = _lazy('CGCS2000', *_T(R_MA, _0_0, 298.257222101)), # BeiDou Coord System (BDC)

2307# Clarke1858 = _lazy('Clarke1858', *_T(6378293.639, _0_0, 294.260676369)),

2308 Clarke1866 = _lazy(_Clarke1866_, *_T(6378206.4, 6356583.8, 294.978698214)),

2309 Clarke1880 = _lazy('Clarke1880', *_T(6378249.145, 6356514.86954978, 293.465)),

2310 Clarke1880IGN = _lazy(_Clarke1880IGN_, *_T(6378249.2, 6356515.0, 293.466021294)),

2311 Clarke1880Mod = _lazy('Clarke1880Mod', *_T(6378249.145, 6356514.96639549, 293.466307656)), # aka Clarke1880Arc

2312 CPM1799 = _lazy('CPM1799', *_T(6375738.7, 6356671.92557493, 334.39)), # Comm. des Poids et Mesures

2313 Delambre1810 = _lazy('Delambre1810', *_T(6376428.0, 6355957.92616372, 311.5)), # Belgium

2314 Engelis1985 = _lazy('Engelis1985', *_T(6378136.05, 6356751.32272154, 298.2566)),

2315# Everest1830 = _lazy('Everest1830', *_T(6377276.345, _0_0, 300.801699997)),

2316# Everest1948 = _lazy('Everest1948', *_T(6377304.063, _0_0, 300.801699997)),

2317# Everest1956 = _lazy('Everest1956', *_T(6377301.243, _0_0, 300.801699997)),

2318 Everest1969 = _lazy('Everest1969', *_T(6377295.664, 6356094.667915, 300.801699997)),

2319 Everest1975 = _lazy('Everest1975', *_T(6377299.151, 6356098.14512013, 300.8017255)),

2320 Fisher1968 = _lazy('Fisher1968', *_T(6378150.0, 6356768.33724438, 298.3)),

2321# Fisher1968Mod = _lazy('Fisher1968Mod', *_T(6378155.0, _0_0, 298.3)),

2322 GEM10C = _lazy('GEM10C', *_T(R_MA, 6356752.31424783, 298.2572236)),

2323 GPES = _lazy('GPES', *_T(6378135.0, 6356750.0, _0_0)), # "Gen. Purpose Earth Spheroid"

2324 GRS67 = _lazy('GRS67', *_T(6378160.0, _0_0, 298.247167427)), # Lucerne b=6356774.516

2325# GRS67Truncated = _lazy('GRS67Truncated', *_T(6378160.0, _0_0, 298.25)),

2326 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.25722210088)), # IUGG, ITRS, ETRS89

2327# Hayford1924 = _lazy('Hayford1924', *_T(6378388.0, 6356911.94612795, None)), # aka Intl1924 f_=297

2328 Helmert1906 = _lazy('Helmert1906', *_T(6378200.0, 6356818.16962789, 298.3)),

2329# Hough1960 = _lazy('Hough1960', *_T(6378270.0, _0_0, 297.0)),

2330 IAU76 = _lazy('IAU76', *_T(6378140.0, _0_0, 298.257)), # Int'l Astronomical Union

2331 IERS1989 = _lazy('IERS1989', *_T(6378136.0, _0_0, 298.257)), # b=6356751.302

2332 IERS1992TOPEX = _lazy('IERS1992TOPEX', *_T(6378136.3, 6356751.61659215, 298.257223563)), # IERS/TOPEX/Poseidon/McCarthy

2333 IERS2003 = _lazy('IERS2003', *_T(6378136.6, 6356751.85797165, 298.25642)),

2334 Intl1924 = _lazy(_Intl1924_, *_T(6378388.0, _0_0, 297.0)), # aka Hayford b=6356911.9462795

2335 Intl1967 = _lazy('Intl1967', *_T(6378157.5, 6356772.2, 298.24961539)), # New Int'l

2336 Krassovski1940 = _lazy(_Krassovski1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2337 Krassowsky1940 = _lazy(_Krassowsky1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2338# Kaula = _lazy('Kaula', *_T(6378163.0, _0_0, 298.24)), # Kaula 1961

2339# Lerch = _lazy('Lerch', *_T(6378139.0, _0_0, 298.257)), # Lerch 1979

2340 Maupertuis1738 = _lazy('Maupertuis1738', *_T(6397300.0, 6363806.28272251, 191.0)), # France

2341 Mercury1960 = _lazy('Mercury1960', *_T(6378166.0, 6356784.28360711, 298.3)),

2342 Mercury1968Mod = _lazy('Mercury1968Mod', *_T(6378150.0, 6356768.33724438, 298.3)),

2343# MERIT = _lazy('MERIT', *_T(6378137.0, _0_0, 298.257)), # MERIT 1983

2344# NWL10D = _lazy('NWL10D', *_T(6378135.0, _0_0, 298.26)), # Naval Weapons Lab.

2345 NWL1965 = _lazy('NWL1965', *_T(6378145.0, 6356759.76948868, 298.25)), # Naval Weapons Lab.

2346# NWL9D = _lazy('NWL9D', *_T(6378145.0, 6356759.76948868, 298.25)), # NWL1965

2347 OSU86F = _lazy('OSU86F', *_T(6378136.2, 6356751.51693008, 298.2572236)),

2348 OSU91A = _lazy('OSU91A', *_T(6378136.3, 6356751.6165948, 298.2572236)),

2349# Plessis1817 = _lazy('Plessis1817', *_T(6397523.0, 6355863.0, 153.56512242)), # XXX incorrect?

2350 Plessis1817 = _lazy('Plessis1817', *_T(6376523.0, 6355862.93325557, 308.64)), # XXX IGN France 1972

2351# Prolate = _lazy('Prolate', *_T(6356752.3, R_MA, _0_0)),

2352 PZ90 = _lazy('PZ90', *_T(6378136.0, _0_0, 298.257839303)), # GLOSNASS PZ-90 and PZ-90.11

2353# SEAsia = _lazy('SEAsia', *_T(6378155.0, _0_0, 298.3)), # SouthEast Asia

2354 SGS85 = _lazy('SGS85', *_T(6378136.0, 6356751.30156878, 298.257)), # Soviet Geodetic System

2355 SoAmerican1969 = _lazy('SoAmerican1969', *_T(6378160.0, 6356774.71919531, 298.25)), # South American

2356 Sphere = _lazy(_Sphere_, *_T(R_M, R_M, _0_0)), # pseudo

2357 SphereAuthalic = _lazy('SphereAuthalic', *_T(R_FM, R_FM, _0_0)), # pseudo

2358 SpherePopular = _lazy('SpherePopular', *_T(R_MA, R_MA, _0_0)), # EPSG:3857 Spheroid

2359 Struve1860 = _lazy('Struve1860', *_T(6378298.3, 6356657.14266956, 294.73)),

2360# Walbeck = _lazy('Walbeck', *_T(6376896.0, _0_0, 302.78)),

2361# WarOffice = _lazy('WarOffice', *_T(6378300.0, _0_0, 296.0)),

2362 WGS60 = _lazy('WGS60', *_T(6378165.0, 6356783.28695944, 298.3)),

2363 WGS66 = _lazy('WGS66', *_T(6378145.0, 6356759.76948868, 298.25)),

2364 WGS72 = _lazy(_WGS72_, *_T(6378135.0, _0_0, 298.26)), # b=6356750.52

2365 WGS84 = _lazy(_WGS84_, *_T(R_MA, _0_0, _f__WGS84)), # GPS b=6356752.3142451793

2366# U{NOAA/NOS/NGS/inverse<https://GitHub.com/noaa-ngs/inverse/blob/main/invers3d.f>}

2367 WGS84_NGS = _lazy('WGS84_NGS', *_T(R_MA, _0_0, 298.257222100882711243162836600094))

2368)

2369

2370_EWGS84 = Ellipsoids.WGS84 # (INTERNAL) shared

2371

2372 if __name__ == '__main__':

2373

2374 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

2375 from pygeodesy import nameof, printf

2376

2377 for E in (_EWGS84, Ellipsoids.GRS80, # NAD83,

2378 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

2379 Ellipsoid(_EWGS84.b, _EWGS84.a, name='_Prolate')):

2380 e = f2n(E.f) - E.n

2381 printf('# %s: %s', _DOT_('Ellipsoids', E.name), E.toStr(prec=10, terse=0), nl=1)

2382 printf('# e=%s, f_=%s, f=%s, n=%s (%s)', fstr(E.e, prec=13, fmt=Fmt.e),

2383 fstr(E.f_, prec=13, fmt=Fmt.e),

2384 fstr(E.f, prec=13, fmt=Fmt.e),

2385 fstr(E.n, prec=13, fmt=Fmt.e),

2386 fstr(e, prec=9, fmt=Fmt.e))

2387 printf('# %s %s', Ellipsoid.AlphaKs.name, fstr(E.AlphaKs, prec=20))

2388 printf('# %s %s', Ellipsoid.BetaKs.name, fstr(E.BetaKs, prec=20))

2389 printf('# %s %s', nameof(Ellipsoid.KsOrder), E.KsOrder) # property

2390

2391 # __doc__ of this file, force all into registry

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

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

2394

2395# % python3 -m pygeodesy.ellipsoids

2396

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

2398# e=8.1819190842622e-02, f_=2.98257223563e+02, f=3.3528106647475e-03, n=1.6792203863837e-03 (0.0e+00)

2399# AlphaKs 0.00083773182062446994, 0.00000076085277735725, 0.00000000119764550324, 0.00000000000242917068, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2400# BetaKs 0.00083773216405794875, 0.0000000590587015222, 0.00000000016734826653, 0.00000000000021647981, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2401# KsOrder 8

2402

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

2404# e=8.1819191042833e-02, f_=2.9825722210088e+02, f=3.3528106811837e-03, n=1.6792203946295e-03 (0.0e+00)

2405# AlphaKs 0.00083773182472890429, 0.00000076085278481561, 0.00000000119764552086, 0.00000000000242917073, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2406# BetaKs 0.0008377321681623882, 0.00000005905870210374, 0.000000000167348269, 0.00000000000021647982, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2407# KsOrder 8

2408

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

2410# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2411# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2412# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2413# KsOrder 8

2414

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

2416# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2417# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2418# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2419# KsOrder 8

2420

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

2422# e=8.2094437949696e-02, f_=-2.97257223563e+02, f=-3.3640898209765e-03, n=-1.6792203863837e-03 (0.0e+00)

2423# AlphaKs -0.00084149152514366627, 0.00000076653480614871, -0.00000000120934503389, 0.0000000000024576225, -0.00000000000000578863, 0.00000000000000001502, -0.00000000000000000004, 0.0

2424# BetaKs -0.00084149187224351817, 0.00000005842735196773, -0.0000000001680487236, 0.00000000000021706261, -0.00000000000000038002, 0.00000000000000000073, -0.0, 0.0

2425# KsOrder 8

2426

2427# **) MIT License

2428#

2429# Copyright (C) 2016-2025 -- mrJean1 at Gmail -- All Rights Reserved.

2430#

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

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

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

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

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

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

2437#

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

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

2440#

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

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

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

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

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

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

2447# OTHER DEALINGS IN THE SOFTWARE.