Coverage for pygeodesy/ellipsoids.py: 96%

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

4u'''Ellipsoidal and spherical earth models.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

28@var Ellipsoids.Everest1967: Ellipsoid(name='Everest1967', a=6377298.556, f=0.00332445, e=0.08147298, b=6356097.55030068, f_=300.8017)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

64'''

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

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

68# from pygeodesy.albers import AlbersEqualAreaCylindrical # _MODS

69from pygeodesy.basics import copysign0, isbool, _isin, isint, typename

70from pygeodesy.constants import EPS, EPS_2, EPS0, EPS02, EPS1, INF, NINF, \

71 _EPSqrt, PI_2, PI_3, PI2, PI4, R_M, R_MA, R_FM, \

72 _EPStol as _TOL, _floatuple as _T, _isfinite, \

73 _1_over, _over, _SQRT2, _0_0s, _0_0, _0_5, _1_0, \

74 _1_EPS, _2_0, _4_0, _90_0

75from pygeodesy.constants import _0_25, _3_0 # PYCHOK used!

76# from pygeodesy.ellipses import Ellipse # _MODS

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

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

79 hypot1, hypot2, sqrt3, Fsum

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

81# from pygeodesy.internals import typename # from .basics

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

83 _Clarke1866_, _Clarke1880IGN_, _DMAIN_, _DOT_, _f_, _GRS80_, \

84 _Intl1924_, _incompatible_, _invalid_, _Krassovski1940_, \

85 _Krassowsky1940_, _lat_, _meridional_, _negative_, _not_finite_, \

86 _prime_vertical_, _Sphere_, _SPACE_, _vs_, _WGS72_, _WGS84_

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

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

89 _NamedEnumItem, _NamedTuple, _Pass, _ALL_LAZY, _MODS

90from pygeodesy.namedTuples import Circle4Tuple, Distance2Tuple, Vector3Tuple, Vector4Tuple

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

92 deprecated_property_RO, property_RO, property_ROver, \

93 deprecated_method

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

95# from pygeodesy.triaxials.triaxial5 import _hartzell3, _plumbTo3 # _MODS

96from pygeodesy.units import Azimuth, Bearing, Distance, Float, Float_, Lamd, Lat, _Lat0, \

97 Meter, Meter2, Meter3, Phi, Phid, Radius, Radius_, Scalar

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

100from math import asinh, atan, atanh, cos, degrees, exp, fabs, radians, sin, sinh, sqrt, tan # as _tan

101

102__all__ = _ALL_LAZY.ellipsoids

103__version__ = '26.10.30'

104

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

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

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

108# _f_WGS84 = Float(f =_1_0 / (298257223563 / 1000000000)) # 0.003_352_810_664_747_480_5

109_f__WGS84 = Float(f_=_1_0 / (1000000000 / 298257223563)) # 298.257223562_999_97 vs 298.257223563

110

111

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

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

114 '''

115 return Lat(lat, clip=clip, name=_lat_ if inverse else typename(auxLat))

116

117

118def _sin2cos2(rad):

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

120 '''

121 if rad:

122 s2 = sin(rad)**2

123 if s2 > EPS:

124 c2 = _1_0 - s2

125 if c2 > EPS:

126 if c2 < EPS1:

127 return s2, c2

128 else:

129 return _1_0, _0_0 # rad == PI_2

130 return _0_0, _1_0 # rad == 0

131

132

133class a_f2Tuple(_NamedTuple):

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

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

136

137 @see: Class L{Ellipsoid2}.

138 '''

139 _Names_ = (_a_, _f_) # name 'f' NOT 'f_'

140 _Units_ = (_Pass, _Pass)

141

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

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

144

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

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

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

148

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

150

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

152

153 @note: C{abs(B{f}) < }L{EPS<pygeodesy.constants.EPS>} is

154 forced to C{B{f}=0}, I{spherical}.

155

156 @note: Negative C{B{f}} produces a I{prolate} ellipsoid.

157 '''

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

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

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

161 f = _f_0_0

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

163

164 @Property_RO

165 def b(self):

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

167 '''

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

169

170 def ellipsoid(self, **name):

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

172

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

174

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

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

177 '''

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

179

180 @Property_RO

181 def f_(self):

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

183 '''

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

185

186

187class Curvature2Tuple(_NamedTuple):

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

189 C{meter}, conventionally.

190 '''

191 _Names_ = (_meridional_, _prime_vertical_)

192 _Units_ = ( Meter, Meter)

193

194 @property_RO

195 def transverse(self):

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

197 '''

198 return self.prime_vertical

199

200

201class Ellipsoid(_NamedEnumItem):

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

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

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

205 '''

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

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

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

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

210

211 _geodsolve = NN # means, use PYGEODESY_GEODSOLVE

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

213 _rhumbsolve = NN # means, use PYGEODESY_RHUMBSOLVE

214

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

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

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

218

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

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

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

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

223 spherical).

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

225

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

227

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

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

230

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

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

233 '''

234 ff_ = f, f_ # assertion below

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

236 try:

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

238 if not _isfinite(a):

239 raise ValueError(_SPACE_(_a_, _not_finite_))

240

241 if b: # not _isin(b, None, _0_0)

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

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

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

245 elif f is not None:

246 f = Float(f=f)

247 b = a_f2b(a, f)

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

249 elif f_:

250 f_ = Float(f_=f_)

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

252 f = f_2f(f_)

253 else: # only a, spherical

254 f_ = f = 0

255 b = a # superfluous

256

257 if not (f < _1_0): # sanity check, see .ecef.Ecef.__init__

258 raise ValueError(_SPACE_(_f_, _invalid_))

259 if not _isfinite(b):

260 raise ValueError(_SPACE_(_b_, _not_finite_))

261

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

263 b = a

264 f = _f_0_0

265 f_ = _f__0_0

266

267 except (TypeError, ValueError) as x:

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

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

270 raise _ValueError(t, cause=x)

271

272 self._a = a

273 self._b = b

274 self._f = f

275 self._f_ = f_

276

277 self._register(Ellipsoids, n)

278

279 if f and f_: # see test/testEllipsoidal

280 d = dict(eps=_TOL)

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

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

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

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

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

286

287 def __eq__(self, other):

288 '''Compare this and an other ellipsoid.

289

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

291

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

293 '''

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

295 self.a == other.a and

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

297

298 def __hash__(self):

299 return self._hash # memoized

300

301 @Property_RO

302 def a(self):

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

304 '''

305 return self._a

306

307 equatoradius = a # = Requatorial

308

309 @Property_RO

310 def a2(self):

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

312 '''

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

314

315 @Property_RO

316 def a2_(self):

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

318 '''

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

320

321 @Property_RO

322 def A(self):

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

324

325 @note: C{A * PI / 2} ≈= L{L<Ellipsoid.L>}, the I{quarter meridian}.

326

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

328 the mean, meridional length I{per radian}.

329 '''

330 A, n = self.a, self.n

331 if n:

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

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

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

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

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

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

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

339 return A

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

341# # q = 4 / self.rocPolar

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

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

344# return Radius(Q=Fhorner(self.e22, 4, -3, 45 / 16, -175 / 64, 11025 / 4096).fover(q))

345

346 @Property_RO

347 def a_b(self):

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

349 '''

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

351

352 @Property_RO

353 def a2_b(self):

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

355

356 @see: U{Radii of Curvature

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

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

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

360

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

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

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

364 '''

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

366

367 @Property_RO

368 def a2_b2(self):

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

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

371 '''

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

373

374 @Property_RO

375 def a_f(self):

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

377 '''

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

379

380 a_f2 = a_f # synonym

381

382 @Property_RO

383 def _albersCyl(self):

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

385 '''

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

387

388 @Property_RO

389 def AlphaKs(self):

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

391 '''

392 return self._Kseries( # XXX int/int quotients may require from __future__ import division as _; del _ # noqa: E702 ;

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

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

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

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

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

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

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

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

401 _T(1424729850961/743921418240)) # PYCHOK unaligned

402

403 @Property_RO

404 def area(self):

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

406

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

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

409 '''

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

411

412 @Property_RO

413 def areax(self):

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

415 accurate for very I{oblate} ellipsoids.

416

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

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

419 '''

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

421

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

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

424 '''

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

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

427 t = (v, _vs_, val)

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

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

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

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

432 raise Error(unstr(self._DOT_(typename(self._assert)), val,

433 eps=eps, f0=f0, **name_value))

434

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

436 '''Compute the I{authalic} auxiliary latitude (Xi) or the I{inverse} thereof.

437

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

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

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

441

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

443

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

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

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

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

448 '''

449 if self.f:

450 _f = self._albersCyl._tanf if inverse else \

451 self._albersCyl._txif # PYCHOK attr

452 lat = atan1d(_f(tan(Phid(lat)))) # PYCHOK attr

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

454

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

456 '''Compute the I{conformal} auxiliary latitude (Chi) or the I{inverse} thereof.

457

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

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

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

461

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

463

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

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

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

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

468 '''

469 if self.f:

470 _f = self.es_tauf if inverse else self.es_taupf # PYCHOK attr

471 lat = atan1d(_f(tan(Phid(lat)))) # PYCHOK attr

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

473

474 def auxGeocentric(self, lat, inverse=False, height=_0_0):

475 '''Compute the I{geocentric} auxiliary latitude (Theta) or the I{inverse} thereof.

476

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

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

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

480 @kwarg height: Optional, ellipsoidal height (C{meter}).

481

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

483

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

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

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

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

488 '''

489 if self.f: # and lat

490 t = tan(Phid(lat))

491 f = self.b2_a2

492 if height:

493 if inverse:

494 lat = atan1d(t * f) # geodetic n

495 d, f = f, _1_0

496 else:

497 d = _1_0

498 n = self.rocPrimeVertical(lat)

499 f = _over(n * f + height, n * d + height)

500 elif inverse:

501 f = self.a2_b2

502 lat = atan1d(t * f)

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

504

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

506 '''Compute the I{isometric} auxiliary latitude (Psi) or the I{inverse} thereof.

507

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

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

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

511

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

513

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

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

516 is the original geodetic latitude.

517

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

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

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

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

522 '''

523 if self.f:

524 r = Phid(lat, clip=0)

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

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

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

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

529

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

531 '''Compute the I{parametric} auxiliary latitude (Beta) or the I{inverse} thereof.

532

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

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

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

536

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

538

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

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

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

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

543 '''

544 if self.f:

545 lat = self._Betad(Lat(lat), inverse=inverse)

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

547

548 auxReduced = auxParametric # synonymous

549

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

551 '''Compute the I{rectifying} auxiliary latitude (Mu) or the I{inverse} thereof.

552

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

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

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

556

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

558

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

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

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

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

563 '''

564 if self.f:

565 lat = Lat(lat)

566 if 0 < fabs(lat) < _90_0:

567 if inverse:

568 e = self._elliptic_e22

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

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

571 else:

572 lat = _over(self.Llat(lat), self.L) * _90_0

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

574

575 @Property_RO

576 def b(self):

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

578 '''

579 return self._b

580

581 polaradius = b # = Rpolar

582

583 @Property_RO

584 def b_a(self):

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

586

587 @see: Property L{f1}.

588 '''

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

590

591 @Property_RO

592 def b2(self):

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

594 '''

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

596

597 @Property_RO

598 def b2_a(self):

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

600

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

602 '''

603 return Radius(b2_a=_over(self.b2, self.a) if self.f else self.b)

604

605 @Property_RO

606 def b2_a2(self):

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

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

609 '''

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

611

612 def _Betad(self, lat, inverse=False):

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

614 '''

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

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

617 return atan1d(s, c) # (s * a, c * b) if inverse else (s * b, c * a)

618

619 @Property_RO

620 def BetaKs(self):

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

622 '''

623 return self._Kseries( # XXX int/int quotients may require from __future__ import division as _; del _ # noqa: E702 ;

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

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

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

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

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

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

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

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

632 _T(191773887257/3719607091200)) # PYCHOK unaligned

633

634 @deprecated_Property_RO

635 def c(self): # PYCHOK no cover

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

637 return self.R2

638

639 @Property_RO

640 def c2(self):

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

642

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

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

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

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

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

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

649 '''

650 return self._c2f(False)

651

652 @Property_RO

653 def c2x(self):

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

655 ellipsoids.

656

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

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

659 '''

660 return self._c2f(True)

661

662 def _c2f(self, c2x):

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

664 '''

665 f, c2 = self.f, self.b2

666 if f:

667 e = self.e

668 if e > EPS0:

669 if f > 0: # .isOblate

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

671 elif f < 0: # .isProlate

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

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

674 return c2

675

676 def circle4(self, lat):

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

678

679 @arg lat: Geodetic latitude (C{degrees90} or C{str}).

680

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

682

683 @raise RangeError: Latitude B{C{lat}} outside valid range, only if

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

685

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

687

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

689

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

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

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

693 and method C{Triaxial.ellipse5}.

694 '''

695 lat = _Lat0(lat)

696 if lat:

697 B = lat # beta

698 if fabs(lat) < _90_0:

699 if self.f:

700 B = self._Betad(lat)

701 z, r = sincos2d(B)

702 r *= self.a

703 z *= self.b

704 else: # near-polar

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

706 else: # equatorial

707 r = self.a

708 z = lat = B = _0_0

709 return Circle4Tuple(r, z, lat, B)

710

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

712 '''Convert an arc in degrees to a distance along the equator or

713 along a parallel of (geodetic) latitude.

714

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

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

717

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

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

720

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

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

723

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

725 '''

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

727

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

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

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

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

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

733

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

735 and only between points not near-polar}.

736

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

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

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

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

741

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

743 in same units as this ellipsoid's axes.

744

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

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

747

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

749 '''

750 phi0 = Phid(lat0=lat0)

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

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

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

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

755

756 @Property_RO

757 def e(self):

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

759

760 @see: Property L{es}.

761 '''

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

763

764 @deprecated_Property_RO

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

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

767 return self.e21

768

769# @Property_RO

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

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

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

773

774 @Property_RO

775 def e2(self):

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

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

778 '''

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

780

781 @Property_RO

782 def e2abs(self):

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

784 '''

785 return fabs(self.e2)

786

787 @Property_RO

788 def e2not(self):

789 '''Get I{self.e**2 and NOT .e2 NOR .e2abs} (C{float}).

790 '''

791 return self.e**2

792

793 @Property_RO

794 def e21(self):

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

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

797 '''

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

799

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

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

802

803 @Property_RO

804 def e22(self):

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

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

807 '''

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

809

810 @Property_RO

811 def e22abs(self):

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

813 '''

814 return fabs(self.e22)

815

816 @Property_RO

817 def e32(self):

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

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

820 '''

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

822

823 @Property_RO

824 def e32abs(self):

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

826 '''

827 return fabs(self.e32)

828

829 @Property_RO

830 def e4(self):

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

832 '''

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

834

835 eccentricity = e # eccentricity

836# eccentricity2 = e2 # eccentricity squared

837 eccentricity1st2 = e2 # first eccentricity squared, signed

838 eccentricity2nd2 = e22 # second eccentricity squared, signed

839 eccentricity3rd2 = e32 # third eccentricity squared, signed

840

841 def ecef(self, Ecef=None):

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

843

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

845

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

847

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

849

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

851 '''

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

853

854 @Property_RO

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

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

857 '''

858 e12 = _over(self.e2, self.e2 - _1_0) # NOT DEPRECATED .e12!

859 return _MODS.elliptic.Elliptic(e12)

860

861 @Property_RO

862 def _elliptic_e22(self): # aka ._elliptic_ep2

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

864 '''

865 # courtesy gaoflow <https://GitHub.com/mrJean1/PyGeodesy/pull/86>

866 return _MODS.elliptic.Elliptic(-self.e22) # signed, like ._elliptic_e12

867

868 equatoradius = a # Requatorial

869

870 @Property_RO

871 def equatorimeter(self):

872 '''Get the ellipsoid's I{equatorial} perimeter (C{meter}).

873 '''

874 return Meter(equatorimeter=self.a * PI2)

875

876 def e2s(self, s):

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

878

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

880

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

882

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

884 '''

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

886

887 def e2s2(self, s):

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

889

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

891

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

893

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

895 '''

896 r, e2 = _1_0, self.e2

897 if e2: # and s

898 try:

899 r -= e2 * Scalar(s=s)**2

900 if r < 0:

901 raise ValueError(_negative_)

902 except (TypeError, ValueError) as x:

903 t = self._DOT_(typename(Ellipsoid.e2s2))

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

905 return r

906

907 @Property_RO

908 def es(self):

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

910

911 @see: Property L{e}.

912 '''

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

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

915

916 @deprecated_method

917 def es_atanh(self, x):

918 '''DEPRECATED on 2026.10.30, use method L{es_eatanhe}.'''

919 return self.es_eatanhe(x)

920

921 @Property_RO

922 def es_c(self):

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

924 '''

925 return Float(es_c=(self._es_Taus._exp_eatanhe_1 * self.b_a) if self.f else _1_0)

926

927 def es_eatanhe(self, x): # in .rhumb.ekx._DeatanhE, .webmercator

928 '''Compute M{es * atanh(es * x)} or M{-es * atan(es * x)} for I{oblate}

929 respectively I{prolate} ellipsoids where I{es} is the I{signed} (1st)

930 eccentricity.

931

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

933

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

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

936 '''

937 return _eatanhe(Scalar(x=x), self.es) if self.f else _0_0

938

939 def es_tauf(self, taup): # in .etm, .ktm, .ups

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

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

942

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

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

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

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

947 '''

948 return self._es_Taus.tauf(taup)

949

950 def es_taupf(self, tau): # in .etm, .ktm, .rhumb.ekx, .ups

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

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

953

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

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

956 '''

957 return self._es_Taus.taupf(tau)

958

959 @Property_RO

960 def _es_Taus(self): # in .utm.Utm._toLLEB

961 return _Taus(self.es, self._1_e21)

962

963 @Property_RO

964 def f(self):

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

966 '''

967 return self._f

968

969 @Property_RO

970 def f_(self):

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

972 '''

973 return self._f_

974

975 @Property_RO

976 def f1(self):

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

978

979 @see: Property L{b_a}.

980 '''

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

982

983 @Property_RO

984 def f2(self):

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

986 '''

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

988

989 @deprecated_Property_RO

990 def geodesic(self):

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

992 return self.geodesicw

993

994 def geodesic_(self, exact=True):

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

996

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

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

999 instance for I{this} ellipsoid.

1000

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

1002

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

1004

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

1006 '''

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

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

1009 name=self.name)

1010 else:

1011 g = exact

1012 E = _xattr(g, ellipsoid=None)

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

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

1015 return g

1016

1017 @property_ROver

1018 def _Geodesics(self):

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

1020 '''

1021 t = (_MODS.geodesicx.GeodesicExact,

1022 _MODS.geodsolve.GeodesicSolve)

1023 try:

1024 t += (_MODS.geodesicw.Geodesic,

1025 _MODS.geodesicw._wrapped.Geodesic)

1026 except ImportError:

1027 pass

1028 return t # overwrite property_ROver

1029

1030 @property_RO

1031 def geodesicw(self):

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

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

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

1035 package is installed.

1036 '''

1037 # if not self.isEllipsoidal:

1038 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1039 return _MODS.geodesicw.Geodesic(self)

1040

1041 @property_RO

1042 def geodesicx(self):

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

1044 '''

1045 # if not self.isEllipsoidal:

1046 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1048

1049 @property

1050 def geodsolve(self):

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

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

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

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

1055 '''

1056 # if not self.isEllipsoidal:

1057 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1059

1060 @geodsolve.setter # PYCHOK setter!

1061 def geodsolve(self, path):

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

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

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

1065 '''

1066 self._geodsolve = path

1067

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

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

1070 from a Point-Of-View in space.

1071

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

1073 or L{Vector3d}).

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

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

1076 of this ellipsoid or C{False} or C{None} to point to its center.

1077

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

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

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

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

1082

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

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

1085 outside this ellipsoid or in opposite direction.

1086

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

1088

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

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

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

1092 '''

1093 try:

1094 m = _MODS._triaxials_triaxial5

1095 v, d, i = m._hartzell3(pov, los, self._triaxial) # PYLINT E1101

1096 except Exception as x:

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

1098 return Vector4Tuple(v.x, v.y, v.z, d, iteration=i, name__=self.hartzell4)

1099

1100 @Property_RO

1101 def _hash(self):

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

1103

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

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

1106 this ellipsoid's surface.

1107

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

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

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

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

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

1113

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

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

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

1117

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

1119

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

1121

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

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

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

1125 '''

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

1127 r = v.length

1128

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

1130 if r < EPS0: # EPS

1131 v = v.times(_0_0)

1132 h = -a

1133

1134 elif self.isSpherical:

1135 v = v.times(a / r)

1136 h = r - a

1137

1138 elif normal: # perpendicular to ellipsoid

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

1140 if x < EPS0: # PYCHOK no cover

1141 z = copysign0(b, v.z)

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

1143 h = y - b # polar

1144 elif y < EPS0: # PYCHOK no cover

1145 t = a / r

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

1147 h = x - a # equatorial

1148 else: # normal in 1st quadrant

1149 m = _MODS._triaxials_triaxial5

1150 x, y, i = m._plumbTo3(x, y, self) # PYLINT E1101

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

1152 h = t.minus(v).length

1153

1154 else: # radial to ellipsoid's center

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

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

1157 v = v.times(t)

1158 h = r * (_1_0 - t)

1159

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

1161

1162 @Property_RO

1163 def heightMax(self):

1164 '''Get the height limit (C{meter}, conventionally).

1165 '''

1166 return self.a / EPS_2 # self.a * _2_EPS, about 12M lightyears

1167

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

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

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

1171 '''

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

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

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

1175

1176 @Property_RO

1177 def isEllipsoidal(self):

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

1179 '''

1180 return self.f != 0

1181

1182 @Property_RO

1183 def isOblate(self):

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

1185 spherical otherwise.

1186 '''

1187 return self.f > 0

1188

1189 @Property_RO

1190 def isProlate(self):

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

1192 spherical otherwise.

1193 '''

1194 return self.f < 0

1195

1196 @Property_RO

1197 def isSpherical(self):

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

1199 '''

1200 return self.f == 0

1201

1202 def _Kseries(self, *AB8Ks):

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

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

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

1206

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

1208 coefficient tuples.

1209

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

1211

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

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

1214 '''

1215 k = self.KsOrder

1216 if self.n:

1217 ns = fpowers(self.n, k)

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

1219 else:

1220 ks = _0_0s(k)

1221 return ks

1222

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

1224 def KsOrder(self):

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

1226 '''

1227 return self._KsOrder

1228

1229 @KsOrder.setter # PYCHOK setter!

1230 def KsOrder(self, order):

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

1232

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

1234 '''

1235 if not (isint(order) and _isin(order, 4, 6, 8)):

1236 raise _ValueError(order=order)

1237 if self._KsOrder != order:

1238 Ellipsoid.AlphaKs._update(self)

1239 Ellipsoid.BetaKs._update(self)

1240 self._KsOrder = order

1241

1242 @Property_RO

1243 def L(self):

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

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

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

1247 '''

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

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

1250

1251 def Llat(self, lat):

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

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

1254

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

1256

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

1258 hemisphere (C{meter}).

1259 '''

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

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

1262

1263 Lmeridian = Llat # meridional distance

1264

1265 @property_RO

1266 def _Lpd(self):

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

1268 '''

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

1270

1271 @property_RO

1272 def _Lpr(self):

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

1274 '''

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

1276

1277 @deprecated_Property_RO

1278 def majoradius(self): # PYCHOK no cover

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

1280 return self.a

1281

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

1283 '''Convert a distance to an arc in degrees along the equator or

1284 along a parallel of (geodetic) latitude.

1285

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

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

1288

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

1290

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

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

1293

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

1295 '''

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

1297

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

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

1300 a parallel of (geodetic) latitude.

1301

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

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

1304

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

1306

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

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

1309

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

1311 '''

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

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

1314

1315 @deprecated_Property_RO

1316 def minoradius(self): # PYCHOK no cover

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

1318 return self.b

1319

1320 @Property_RO

1321 def n(self):

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

1323 '''

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

1325

1326 flattening = f

1327 flattening1st = f

1328 flattening2nd = f2

1329 flattening3rd = n

1330

1331 polaradius = b # Rpolar

1332

1333 @property_RO

1334 def polarimeter(self):

1335 '''Get the ellipsoid's I{polar}, meridional perimeter (C{meter}).

1336 '''

1337 return Meter(polarimeter=self.L * _4_0)

1338

1339# Q = A # I{meridian arc unit} C{Q}, the mean, meridional length I{per radian}

1340

1341 @deprecated_Property_RO

1342 def quarteradius(self): # PYCHOK no cover

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

1344 return self.L

1345

1346 @Property_RO

1347 def R1(self):

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

1349

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

1351 and method C{Rgeometric}.

1352 '''

1353 return Radius(R1=fmean_(self.a, self.a, self.b) if self.f else self.a)

1354

1355 Rmean = R1

1356

1357 @Property_RO

1358 def R2(self):

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

1360

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

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

1363 '''

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

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

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

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

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

1369

1370 Rauthalic = R2

1371

1372 @Property_RO

1373 def R2x(self):

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

1375

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

1377 '''

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

1379

1380 Rauthalicx = R2x

1381

1382 @Property_RO

1383 def R3(self):

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

1385

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

1387 '''

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

1389 return Radius(R3=r)

1390

1391 Rvolumetric = R3

1392

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

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

1395 a parallel of (geodetic) latitude.

1396

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

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

1399

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

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

1402

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

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

1405

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

1407 '''

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

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

1410

1411 @Property_RO

1412 def Rbiaxial(self):

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

1414

1415 @see: C{Rtriaxial}

1416 '''

1417 a, b = self.a, self.b

1418 if b != a:

1419 b = hypot(a, b) / _SQRT2

1420 return Radius(Rbiaxial=b)

1421

1422 Requatorial = a # for consistent naming

1423

1424 def Rgeocentric(self, lat):

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

1426

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

1428

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

1430

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

1432

1433 @see: U{Geocentric Radius

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

1435 '''

1436 r, p = self.a, Phid(lat)

1437 if p and self.f:

1438 if fabs(p) < PI_2:

1439 s2, c2 = _sin2cos2(p)

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

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

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

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

1444 s2 *= self.b2_a2

1445 r *= sqrt((c2 + self.b2_a2 * s2) / (c2 + s2))

1446 else:

1447 r = self.b

1448 return Radius(Rgeocentric=r)

1449

1450 @Property_RO

1451 def Rgeometric(self):

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

1453

1454 @see: C{R1}.

1455 '''

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

1457 return Radius(Rgeometric=g)

1458

1459 def rhumb_(self, exact=True):

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

1461

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

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

1464 for I{this} ellipsoid.

1465

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

1467

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

1469

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

1471 '''

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

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

1474 else:

1475 r = exact

1476 E = _xattr(r, ellipsoid=None)

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

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

1479 return r

1480

1481 @property_RO

1482 def rhumbaux(self):

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

1484 '''

1485 # if not self.isEllipsoidal:

1486 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1488

1489 @property_RO

1490 def rhumbekx(self):

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

1492 '''

1493 # if not self.isEllipsoidal:

1494 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1496

1497 @property_ROver

1498 def _Rhumbs(self):

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

1500 '''

1501 r = _MODS.rhumb

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

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

1504

1505 @property

1506 def rhumbsolve(self):

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

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

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

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

1511 '''

1512 # if not self.isEllipsoidal:

1513 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1515

1516 @rhumbsolve.setter # PYCHOK setter!

1517 def rhumbsolve(self, path):

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

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

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

1521 '''

1522 self._rhumbsolve = path

1523

1524 @deprecated_property_RO

1525 def rhumbx(self):

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

1527 return self.rhumbekx

1528

1529 def Rlat(self, lat):

1530 '''I{Average} the earth radius between C{equatoradius} at C{0} and

1531 C{polaradius} at C{+/-90} degrees latitude.

1532

1533 @arg lat: Latitude (C{degrees90}, C{str} or C{Ang}).

1534

1535 @return: Averaged earth radius (C{meter}).

1536

1537 @raise RangeError: Latitude B{C{lat}} outside valid range, only if

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

1539

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

1541

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

1543 '''

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

1545 r = self.a

1546 if self.f: # .isEllipsoidal

1547 lat = _Lat0(lat)

1548 if lat:

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

1550 r = Radius(Rlat=r)

1551 return r

1552

1553 Rpolar = b # for consistent naming

1554

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

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

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

1558

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

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

1561 to use an alternate formula.

1562

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

1564

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

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

1567

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

1569 '''

1570 if sa and self.f: # .isEllipsoidal

1571 if ca is None:

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

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

1574 elif ca: # derived from EcefYou.forward

1575 h = hypot(ca, self.b_a * sa)

1576 n = self.a / h

1577 else:

1578 n = self.a2_b / fabs(sa)

1579 else:

1580 n = self.a

1581 return n

1582

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

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

1585 radii of curvature of (geodetic) latitude.

1586

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

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

1589

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

1591 the radii of curvature.

1592

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

1594

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

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

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

1598 Earth_radius#Radii_of_curvature>}.

1599 '''

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

1601

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

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

1604 curvature of (geodetic) latitude.

1605

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

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

1608

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

1610 radii of curvature.

1611

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

1613

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

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

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

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

1618 '''

1619 p = fabs(Phi(phi))

1620 if self.f:

1621 r = self.e2s2(sin(p))

1622 if r > EPS02:

1623 n = sqrt(self.a2 / r)

1624 m = n * self.e21 / r

1625 else:

1626 m = n = _0_0

1627 else:

1628 m = n = self.a

1629 if scaled and p:

1630 n *= cos(p) if p < PI_2 else _0_0

1631 return Curvature2Tuple(m, n)

1632

1633 def rocAzimuth(self, lat, azimuth):

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

1635 and C{azimuth} compass direction.

1636

1637 @see: Method L{rocBearing<Ellipsoid.rocBearing>} for details, using C{azimuth} for C{bearing}.

1638 '''

1639 return Radius(rocAzimuth=self._rocDirectional(lat, Azimuth(azimuth)))

1640

1641 def rocBearing(self, lat, bearing):

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

1643 and C{bearing} compass direction.

1644

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

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

1647

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

1649

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

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

1652

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

1654

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

1656 '''

1657 return Radius(rocBearing=self._rocDirectional(lat, Bearing(bearing)))

1658

1659 def _rocDirectional(self, lat, deg):

1660 '''(INTERNAL) Helper for C{rocAzimuth} and C{rocBearing}.

1661 '''

1662 if self.f:

1663 s2, c2 = _sin2cos2(radians(deg))

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

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

1666 c2 *= n / m

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

1668 s2 *= m / n

1669 n = m

1670 r = _over(n, c2 + s2) # == 1 / (c2 / m + s2 / n)

1671 else:

1672 r = self.b # == self.a

1673 return r

1674

1675 @Property_RO

1676 def rocEquatorial2(self):

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

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

1679

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

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

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

1683 '''

1684 m = self.b2_a if self.f else self.a

1685 return Curvature2Tuple(m, self.a)

1686

1687 def rocGauss(self, lat):

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

1689

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

1691

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

1693

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

1695

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

1697 Earth_radius#Radii_of_curvature>}

1698 '''

1699 # using ...

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

1701 # return sqrt(m * n)

1702 # ... requires 1 or 2 sqrt

1703 g = self.b

1704 if self.f:

1705 s2, c2 = _sin2cos2(Phid(lat))

1706 g = _over(g, c2 + self.b2_a2 * s2)

1707 return Radius(rocGauss=g)

1708

1709 def rocMean(self, lat):

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

1711

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

1713

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

1715

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

1717

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

1719 Earth_radius#Radii_of_curvature>}

1720 '''

1721 if self.f:

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

1723 m *= _over(n * _2_0, m + n) # == 2 / (1 / m + 1 / n)

1724 else:

1725 m = self.a

1726 return Radius(rocMean=m)

1727

1728 def rocMeridional(self, lat):

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

1730

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

1732

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

1734

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

1736

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

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

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

1740 '''

1741 r = self.roc2_(Phid(lat)) if lat else self.rocEquatorial2

1742 return Radius(rocMeridional=r.meridional)

1743

1744 rocPolar = a2_b # synonymous

1745

1746 def rocPrimeVertical(self, lat):

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

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

1749

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

1751

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

1753

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

1755

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

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

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

1759 Earth_radius#Radii_of_curvature>}.

1760 '''

1761 r = self.roc1_(sin(Phid(lat))) if lat else self.a

1762 return Radius(rocPrimeVertical=r)

1763

1764 rocTransverse = rocPrimeVertical # synonymous

1765

1766 @deprecated_Property_RO

1767 def Rquadratic(self): # PYCHOK no cover

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

1769 return self.Rbiaxial

1770

1771 @deprecated_Property_RO

1772 def Rr(self): # PYCHOK no cover

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

1774 return self.Rrectifying

1775

1776 @Property_RO

1777 def Rrectifying(self):

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

1779

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

1781 '''

1782 r = self.a

1783 if self.f:

1784 r *= cbrt2((sqrt3(self.b_a) + _1_0) * _0_5)

1785 return Radius(Rrectifying=r)

1786

1787 @deprecated_Property_RO

1788 def Rs(self): # PYCHOK no cover

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

1790 return self.Rgeometric

1791

1792 @Property_RO

1793 def Rtriaxial(self):

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

1795

1796 @see: C{Rbiaxial}

1797 '''

1798 q, b = self.a, self.b

1799 if b < q:

1800 q *= sqrt((self.b2_a2 + _3_0) * _0_25)

1801 elif b > q:

1802 q = sqrt((self.a2_b2 * _3_0 + _1_0) * _0_25) * b

1803 return Radius(Rtriaxial=q)

1804

1805 def toEllipse(self, **name):

1806 '''Get this ellipsoid's meridional ellipse L{Ellipse<pygeodesy.Ellipse>}.

1807

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

1809 '''

1810 return _MODS.ellipses.Ellipse(self.a, self.b, **name)

1811

1812 def toEllipsoid2(self, **name):

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

1814

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

1816

1817 @see: Property C{a_f}.

1818 '''

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

1820

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

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

1823

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

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

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

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

1828 to exclude this ellipsoid's name and separator

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

1830

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

1832 '''

1833 E = Ellipsoid

1834 t = (E.a, E.f, E.e, E.b, E.f_, E.e2,

1835 E.e21, E.e22, E.e32, E.n,

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

1837 E.Rbiaxial, E.Rtriaxial)

1838 if terse:

1839 t = t[:terse]

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

1841

1842 def toTriaxial(self, **name):

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

1844

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

1846

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

1848 pointing east and C{Z} pointing north.

1849

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

1851 '''

1852 T = self._triaxial

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

1854

1855 @property_RO

1856 def _triaxial(self):

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

1858 '''

1859 a, b, t = self.a, self.b, _MODS.triaxials

1860 T = t.Triaxial if a > b else t.Triaxial_

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

1862

1863 @Property_RO

1864 def volume(self):

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

1866

1867 @see: C{R3}.

1868 '''

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

1870

1871

1872class Ellipsoid2(Ellipsoid):

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

1874 '''

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

1876 '''New L{Ellipsoid2}.

1877

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

1879 L{Ellipsoid} instance.

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

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

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

1883

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

1885

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

1887

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

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

1890 '''

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

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

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

1894 else:

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

1896

1897

1898class _Taus(object): # in .karney

1899 '''(INTERNAL) Class providing C{tauf} and C{taupf} functions.

1900 '''

1901

1902 def __init__(self, es, _1_e21=None):

1903 self._es = es

1904 if _1_e21 is not None:

1905 self.__dict__.update(_1_e21=_1_e21)

1906

1907 @Property_RO

1908 def _1_e21(self): # == E.a2_b2 == 1 / e21 == 1 / e2m

1909 return _1_over(_1_0 - self.es**2) if self.es else _1_0

1910

1911 def eatanhe(self, x): # see .albers._atanheE, .AuxLat._atanhee

1912 '''See method L{Ellipsoid.es_eatanhe} for details.

1913 '''

1914 return _eatanhe(Scalar(x=x), self.es)

1915

1916 @Property_RO

1917 def es(self):

1918 return self._es

1919

1920 @Property_RO

1921 def _exp_eatanhe_1(self):

1922 '''(INTERNAL) Helper for Ellipsoid.es_c and .es_tauf.

1923 '''

1924 return exp(_eatanhe(_1_0, self.es)) if self.es else _1_0

1925

1926 def tauf(self, taup):

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

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

1929

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

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

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

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

1934 '''

1935 t = Scalar(taup=taup)

1936 if self.es: # .isEllipsoidal

1937 a = fabs(t)

1938 T = (self._exp_eatanhe_1 if a > 70 else self._1_e21) * t

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

1940 s = (a * _TOL) if a > 1 else _TOL

1941 for T, _, d in self._tauf3(t, T): # max 2

1942 if fabs(d) < s:

1943 break

1944 else: # PYCHOK no cover

1945 t = unstr(self.tauf, taup, es=self.es)

1946 raise _AssertionError(Fmt.no_convergence(d, s), txt=t)

1947 t = Scalar(tauf=T)

1948 return t

1949

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

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

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

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

1954 '''

1955 e = self._1_e21

1956 _f2 = Fsum(T).fsum2f_ # τ0

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

1958 a, h = self._taupf2(T)

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

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

1961 T, d = _f2(d) # τi, (τi - τi-1)

1962 yield T, i, d

1963

1964 def taupf(self, tau):

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

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

1967

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

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

1970 '''

1971 t = Scalar(tau=tau)

1972 if self.es: # .isEllipsoidal

1973 t, _ = self._taupf2(t)

1974 t = Scalar(taupf=t)

1975 return t

1976

1977 def _taupf2(self, tau):

1978 '''Return 2-tuple C{(taupf(tau), hypot1(tau))}.

1979 '''

1980 t = tau

1981 if _isfinite(t):

1982 h = hypot1(t)

1983 s = sinh(_eatanhe(t / h, self.es))

1984 t *= hypot1(s)

1985 t -= h * s

1986 else:

1987 h = INF # NAN?

1988 return t, h

1989

1990

1991def _eatanhe(x, es):

1992 '''(INTERNAL) Karney's GeographicLib 2.7+ C{Math.eatanhe}.

1993 '''

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

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

1996 _0_0)) # .isSpherical

1997

1998

1999def _ispherical_a_b(a, b):

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

2001 '''

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

2003

2004

2005def _ispherical_f(f):

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

2007 '''

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

2009

2010

2011def _ispherical_f_(f_):

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

2013 '''

2014 f_ = fabs(f_)

2015 return f_ < EPS or f_ > _1_EPS

2016

2017

2018def a_b2e(a, b):

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

2020

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

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

2023

2024 @return: The I{unsigned}, (1st) eccentricity (C{float} or C{0}), M{sqrt(1 - (b / a)**2)}.

2025

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

2027 '''

2028 e2 = _a2b2e2(a, b, b2=False)

2029 return Float(e=sqrt(fabs(e2)) if e2 else _0_0) # == sqrt(fabs((a - b) * (a + b))) / a

2030

2031

2032def a_b2e2(a, b):

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

2034

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

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

2037

2038 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or C{0}), M{1 - (b / a)**2}.

2039

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

2041 for I{near-spherical} ellipsoids.

2042 '''

2043 return Float(e2=_a2b2e2(a, b, b2=False))

2044

2045

2046def a_b2e22(a, b):

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

2048

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

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

2051

2052 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} or C{0}), M{(a / b)**2 - 1}.

2053

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

2055 for I{near-spherical} ellipsoids.

2056 '''

2057 return Float(e22=_a2b2e2(a, b, a2=False))

2058

2059

2060def a_b2e32(a, b):

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

2062

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

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

2065

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

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

2068

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

2070 for I{near-spherical} ellipsoids.

2071 '''

2072 return Float(e32=_a2b2e2(a, b))

2073

2074

2075def _a2b2e2(a, b, a2=True, b2=True):

2076 '''(INTERNAL) Helper for C{a_b2e}, C{a_b2e2}, C{a_b2e22} and C{a_b2e32}.

2077 '''

2078 if _ispherical_a_b(a, b):

2079 e2 = _0_0

2080 else: # a > 0, b > 0

2081 a, b = (_1_0, b / a) if a > b else (a / b, _1_0)

2082 a2b2 = float(a - b) * (a + b)

2083 e2 = _over(a2b2, (a**2 if a2 else _0_0) +

2084 (b**2 if b2 else _0_0)) if a2b2 else _0_0

2085 return e2

2086

2087

2088def a_b2f(a, b):

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

2090

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

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

2093

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

2095

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

2097 for I{near-spherical} ellipsoids.

2098 '''

2099 f = 0 if _ispherical_a_b(a, b) else _over(float(a - b), a)

2100 return _f_0_0 if _ispherical_f(f) else Float(f=f)

2101

2102

2103def a_b2f_(a, b):

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

2105

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

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

2108

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

2110

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

2112 for I{near-spherical} ellipsoids.

2113 '''

2114 f_ = 0 if _ispherical_a_b(a, b) else _over(a, float(a - b))

2115 return _f__0_0 if _ispherical_f_(f_) else Float(f_=f_)

2116

2117

2118def a_b2f2(a, b):

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

2120

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

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

2123

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

2125

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

2127 for I{near-spherical} ellipsoids.

2128 '''

2129 t = 0 if _ispherical_a_b(a, b) else float(a - b)

2130 return Float(f2=_0_0 if fabs(t) < EPS0 else _over(t, b))

2131

2132

2133def a_b2n(a, b):

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

2135

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

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

2138

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

2140

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

2142 for I{near-spherical} ellipsoids.

2143 '''

2144 t = 0 if _ispherical_a_b(a, b) else float(a - b)

2145 return Float(n=_0_0 if fabs(t) < EPS0 else _over(t, a + b))

2146

2147

2148def a_f2b(a, f):

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

2150

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

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

2153

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

2155 '''

2156 b = a if _ispherical_f(f) else (a * (_1_0 - f))

2157 return Radius_(b=a if _ispherical_a_b(a, b) else b)

2158

2159

2160def a_f_2b(a, f_):

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

2162

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

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

2165

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

2167 '''

2168 b = a if _ispherical_f_(f_) else _over(a * (f_ - _1_0), f_)

2169 return Radius_(b=a if _ispherical_a_b(a, b) else b)

2170

2171

2172def b_f2a(b, f):

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

2174

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

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

2177

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

2179 '''

2180 t = _1_0 - f

2181 a = b if fabs(t) < EPS0 else _over(b, t)

2182 return Radius_(a=b if _ispherical_a_b(a, b) else a)

2183

2184

2185def b_f_2a(b, f_):

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

2187

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

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

2190

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

2192 '''

2193 t = f_ - _1_0

2194 a = b if _ispherical_f_(f_) or fabs(t) < EPS0 \

2195 or fabs(t - f_) < EPS0 else _over(b * f_, t)

2196 return Radius_(a=b if _ispherical_a_b(a, b) else a)

2197

2198

2199def e2f(e):

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

2201

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

2203

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

2205

2206 @see: Function L{e22f}.

2207 '''

2208 return e22f(e**2)

2209

2210

2211def e22f(e2):

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

2213

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

2215

2216 @return: The flattening (C{float} or C{0}), M{e2 / (sqrt(1 - e2) + 1)}.

2217 '''

2218 return Float(f=_over(e2, sqrt(_1_0 - e2) + _1_0)) if e2 and e2 < _1_0 else _f_0_0

2219

2220

2221def f2e2(f):

2222 '''Return C{e2}, the I{1st eccentricity squared} for a given I{flattening}.

2223

2224 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2225

2226 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} < 1), M{f * (2 - f)}.

2227

2228 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2229 for I{near-spherical} ellipsoids.

2230

2231 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2232 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2233 <https://WikiPedia.org/wiki/Flattening>}.

2234 '''

2235 return Float(e2=_0_0 if _ispherical_f(f) else (f * (_2_0 - f)))

2236

2237

2238def f2e22(f):

2239 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{flattening}.

2240

2241 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2242

2243 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} > -1 or C{INF}),

2244 M{f * (2 - f) / (1 - f)**2}.

2245

2246 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2247 for near-spherical ellipsoids.

2248

2249 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2250 C++/doc/classGeographicLib_1_1Ellipsoid.html>}.

2251 '''

2252 # e2 / (1 - e2) == f * (2 - f) / (1 - f)**2

2253 t = (_1_0 - f)**2

2254 return Float(e22=INF if t < EPS0 else _over(f2e2(f), t)) # PYCHOK type

2255

2256

2257def f2e32(f):

2258 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{flattening}.

2259

2260 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2261

2262 @return: The I{signed}, 3rd eccentricity I{squared} (C{float}),

2263 M{f * (2 - f) / (1 + (1 - f)**2)}.

2264

2265 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2266 for I{near-spherical} ellipsoids.

2267

2268 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2269 C++/doc/classGeographicLib_1_1Ellipsoid.html>}.

2270 '''

2271 # e2 / (2 - e2) == f * (2 - f) / (1 + (1 - f)**2)

2272 e2 = f2e2(f)

2273 return Float(e32=_over(e2, _2_0 - e2))

2274

2275

2276def f_2f(f_):

2277 '''Return C{f}, the I{flattening} for a given I{inverse flattening}.

2278

2279 @arg f_: Inverse flattening (C{scalar} >>> 1).

2280

2281 @return: The flattening (C{scalar} or C{0}), M{1 / f_}.

2282

2283 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2284 for I{near-spherical} ellipsoids.

2285 '''

2286 f = 0 if _ispherical_f_(f_) else _over(_1_0, f_)

2287 return _f_0_0 if _ispherical_f(f) else Float(f=f) # PYCHOK type

2288

2289

2290def f2f_(f):

2291 '''Return C{f_}, the I{inverse flattening} for a given I{flattening}.

2292

2293 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2294

2295 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2296

2297 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2298 for I{near-spherical} ellipsoids.

2299 '''

2300 f_ = 0 if _ispherical_f(f) else _over(_1_0, f)

2301 return _f__0_0 if _ispherical_f_(f_) else Float(f_=f_) # PYCHOK type

2302

2303

2304def f2f2(f):

2305 '''Return C{f2}, the I{2nd flattening} for a given I{flattening}.

2306

2307 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2308

2309 @return: The I{signed}, 2nd flattening (C{scalar} or C{INF}), M{f / (1 - f)}.

2310

2311 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2312 for I{near-spherical} ellipsoids.

2313

2314 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2315 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2316 <https://WikiPedia.org/wiki/Flattening>}.

2317 '''

2318 t = _1_0 - f

2319 return Float(f2=_0_0 if _ispherical_f(f) else

2320 (INF if fabs(t) < EPS else _over(f, t))) # PYCHOK type

2321

2322

2323def f2n(f):

2324 '''Return C{n}, the I{3rd flattening} for a given I{flattening}.

2325

2326 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2327

2328 @return: The I{signed}, 3rd flattening (-1 <= C{float} < 1),

2329 M{f / (2 - f)}.

2330

2331 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2332 for I{near-spherical} ellipsoids.

2333

2334 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2335 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2336 <https://WikiPedia.org/wiki/Flattening>}.

2337 '''

2338 return Float(n=_0_0 if _ispherical_f(f) else _over(f, float(_2_0 - f)))

2339

2340

2341def n2e2(n):

2342 '''Return C{e2}, the I{1st eccentricity squared} for a given I{3rd flattening}.

2343

2344 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2345

2346 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or NINF),

2347 M{4 * n / (1 + n)**2}.

2348

2349 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2350 for I{near-spherical} ellipsoids.

2351

2352 @see: U{Flattening<https://WikiPedia.org/wiki/Flattening>}.

2353 '''

2354 t = (n + _1_0)**2

2355 return Float(e2=_0_0 if fabs(n) < EPS0 else

2356 (NINF if t < EPS0 else _over(_4_0 * n, t)))

2357

2358

2359def n2f(n):

2360 '''Return C{f}, the I{flattening} for a given I{3rd flattening}.

2361

2362 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2363

2364 @return: The flattening (C{scalar} or NINF), M{2 * n / (1 + n)}.

2365

2366 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2367 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2368 <https://WikiPedia.org/wiki/Flattening>}.

2369 '''

2370 t = n + _1_0

2371 f = 0 if fabs(n) < EPS0 else (NINF if t < EPS0 else _over(_2_0 * n, t))

2372 return _f_0_0 if _ispherical_f(f) else Float(f=f)

2373

2374

2375def n2f_(n):

2376 '''Return C{f_}, the I{inverse flattening} for a given I{3rd flattening}.

2377

2378 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2379

2380 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2381

2382 @see: L{n2f} and L{f2f_}.

2383 '''

2384 return f2f_(n2f(n))

2385

2386

2387class Ellipsoids(_NamedEnum):

2388 '''(INTERNAL) L{Ellipsoid} registry, I{must} be a sub-class

2389 to accommodate the L{_LazyNamedEnumItem} properties.

2390 '''

2391 def _Lazy(self, a, b, f_, **kwds):

2392 '''(INTERNAL) Instantiate the L{Ellipsoid}.

2393 '''

2394 return Ellipsoid(a, b=b, f_=f_, **kwds)

2395

2396Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

2397'''Some pre-defined L{Ellipsoid}s, all I{lazily} instantiated.'''

2398# <https://www.GNU.org/software/gama/manual/html_node/Supported-ellipsoids.html>

2399# <https://GSSC.ESA.int/navipedia/index.php/Reference_Frames_in_GNSS>

2400# <https://kb.OSU.edu/dspace/handle/1811/77986>

2401# <https://www.IBM.com/docs/en/db2/11.5?topic=systems-supported-spheroids>

2402# <https://w3.Energistics.org/archive/Epicentre/Epicentre_v3.0/DataModel/LogicalDictionary/StandardValues/ellipsoid.html>

2403# <https://GitHub.com/locationtech/proj4j/blob/master/src/main/java/org/locationtech/proj4j/datum/Ellipsoid.java>

2404Ellipsoids._assert( # <https://WikiPedia.org/wiki/Earth_ellipsoid>

2405 Airy1830 = _lazy(_Airy1830_, *_T(6377563.396, _0_0, 299.3249646)), # b=6356256.909

2406 AiryModified = _lazy(_AiryModified_, *_T(6377340.189, _0_0, 299.3249646)), # b=6356034.448

2407# APL4_9 = _lazy('APL4_9', *_T(6378137.0, _0_0, 298.24985392)), # Appl. Phys. Lab. 1965

2408# ANS = _lazy('ANS', *_T(6378160.0, _0_0, 298.25)), # Australian Nat. Spheroid

2409# AN_SA96 = _lazy('AN_SA96', *_T(6378160.0, _0_0, 298.24985392)), # Australian Nat. South America

2410 Australia1966 = _lazy('Australia1966', *_T(6378160.0, _0_0, 298.25)), # b=6356774.7192

2411 ATS1977 = _lazy('ATS1977', *_T(6378135.0, _0_0, 298.257)), # "Average Terrestrial System"

2412 Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.962818, 299.152812797)),

2413 BesselModified = _lazy('BesselModified', *_T(6377492.018, _0_0, 299.1528128)),

2414# BesselNamibia = _lazy('BesselNamibia', *_T(6377483.865, _0_0, 299.1528128)),

2415 CGCS2000 = _lazy('CGCS2000', *_T(R_MA, _0_0, 298.257222101)), # BeiDou Coord System (BDC)

2416# Clarke1858 = _lazy('Clarke1858', *_T(6378293.639, _0_0, 294.260676369)),

2417 Clarke1866 = _lazy(_Clarke1866_, *_T(6378206.4, 6356583.8, 294.978698214)),

2418 Clarke1880 = _lazy('Clarke1880', *_T(6378249.145, 6356514.86954978, 293.465)),

2419 Clarke1880IGN = _lazy(_Clarke1880IGN_, *_T(6378249.2, 6356515.0, 293.466021294)),

2420 Clarke1880Mod = _lazy('Clarke1880Mod', *_T(6378249.145, 6356514.96639549, 293.466307656)), # aka Clarke1880Arc

2421 CPM1799 = _lazy('CPM1799', *_T(6375738.7, 6356671.92557493, 334.39)), # Comm. des Poids et Mesures

2422 Delambre1810 = _lazy('Delambre1810', *_T(6376428.0, 6355957.92616372, 311.5)), # Belgium

2423 Engelis1985 = _lazy('Engelis1985', *_T(6378136.05, 6356751.32272154, 298.2566)),

2424# Everest1830 = _lazy('Everest1830', *_T(6377276.345, _0_0, 300.801699997)),

2425# Everest1948 = _lazy('Everest1948', *_T(6377304.063, _0_0, 300.801699997)),

2426# Everest1956 = _lazy('Everest1956', *_T(6377301.243, _0_0, 300.801699997)),

2427 Everest1967 = _lazy('Everest1967', *_T(6377298.556, _0_0, 300.801699997)),

2428 Everest1969 = _lazy('Everest1969', *_T(6377295.664, 6356094.667915, 300.801699997)),

2429 Everest1975 = _lazy('Everest1975', *_T(6377299.151, 6356098.14512013, 300.8017255)),

2430 Fisher1968 = _lazy('Fisher1968', *_T(6378150.0, 6356768.33724438, 298.3)),

2431# Fisher1968Mod = _lazy('Fisher1968Mod', *_T(6378155.0, _0_0, 298.3)),

2432 GEM10C = _lazy('GEM10C', *_T(R_MA, 6356752.31424783, 298.2572236)),

2433 GPES = _lazy('GPES', *_T(6378135.0, 6356750.0, _0_0)), # "Gen. Purpose Earth Spheroid"

2434 GRS67 = _lazy('GRS67', *_T(6378160.0, _0_0, 298.247167427)), # Lucerne b=6356774.516

2435# GRS67Truncated = _lazy('GRS67Truncated', *_T(6378160.0, _0_0, 298.25)),

2436 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.25722210088)), # IUGG, ITRS, ETRS89

2437# Hayford1924 = _lazy('Hayford1924', *_T(6378388.0, 6356911.94612795, None)), # aka Intl1924 f_=297

2438 Helmert1906 = _lazy('Helmert1906', *_T(6378200.0, 6356818.16962789, 298.3)),

2439# Hough1960 = _lazy('Hough1960', *_T(6378270.0, _0_0, 297.0)),

2440 IAU76 = _lazy('IAU76', *_T(6378140.0, _0_0, 298.257)), # Int'l Astronomical Union

2441 IERS1989 = _lazy('IERS1989', *_T(6378136.0, _0_0, 298.257)), # b=6356751.302

2442 IERS1992TOPEX = _lazy('IERS1992TOPEX', *_T(6378136.3, 6356751.61659215, 298.257223563)), # IERS/TOPEX/Poseidon/McCarthy

2443 IERS2003 = _lazy('IERS2003', *_T(6378136.6, 6356751.85797165, 298.25642)),

2444 Intl1924 = _lazy(_Intl1924_, *_T(6378388.0, _0_0, 297.0)), # aka Hayford b=6356911.9462795

2445 Intl1967 = _lazy('Intl1967', *_T(6378157.5, 6356772.2, 298.24961539)), # New Int'l

2446 Krassovski1940 = _lazy(_Krassovski1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2447 Krassowsky1940 = _lazy(_Krassowsky1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2448# Kaula = _lazy('Kaula', *_T(6378163.0, _0_0, 298.24)), # Kaula 1961

2449# Lerch = _lazy('Lerch', *_T(6378139.0, _0_0, 298.257)), # Lerch 1979

2450 Maupertuis1738 = _lazy('Maupertuis1738', *_T(6397300.0, 6363806.28272251, 191.0)), # France

2451 Mercury1960 = _lazy('Mercury1960', *_T(6378166.0, 6356784.28360711, 298.3)),

2452 Mercury1968Mod = _lazy('Mercury1968Mod', *_T(6378150.0, 6356768.33724438, 298.3)),

2453# MERIT = _lazy('MERIT', *_T(6378137.0, _0_0, 298.257)), # MERIT 1983

2454# NWL10D = _lazy('NWL10D', *_T(6378135.0, _0_0, 298.26)), # Naval Weapons Lab.

2455 NWL1965 = _lazy('NWL1965', *_T(6378145.0, 6356759.76948868, 298.25)), # Naval Weapons Lab.

2456# NWL9D = _lazy('NWL9D', *_T(6378145.0, 6356759.76948868, 298.25)), # NWL1965

2457 OSU86F = _lazy('OSU86F', *_T(6378136.2, 6356751.51693008, 298.2572236)),

2458 OSU91A = _lazy('OSU91A', *_T(6378136.3, 6356751.6165948, 298.2572236)),

2459# Plessis1817 = _lazy('Plessis1817', *_T(6397523.0, 6355863.0, 153.56512242)), # XXX incorrect?

2460 Plessis1817 = _lazy('Plessis1817', *_T(6376523.0, 6355862.93325557, 308.64)), # XXX IGN France 1972

2461# Prolate = _lazy('Prolate', *_T(6356752.3, R_MA, _0_0)),

2462 PZ90 = _lazy('PZ90', *_T(6378136.0, _0_0, 298.257839303)), # GLOSNASS PZ-90 and PZ-90.11

2463# SEAsia = _lazy('SEAsia', *_T(6378155.0, _0_0, 298.3)), # SouthEast Asia

2464 SGS85 = _lazy('SGS85', *_T(6378136.0, 6356751.30156878, 298.257)), # Soviet Geodetic System

2465 SoAmerican1969 = _lazy('SoAmerican1969', *_T(6378160.0, 6356774.71919531, 298.25)), # South American

2466 Sphere = _lazy(_Sphere_, *_T(R_M, R_M, _0_0)), # pseudo

2467 SphereAuthalic = _lazy('SphereAuthalic', *_T(R_FM, R_FM, _0_0)), # pseudo

2468 SpherePopular = _lazy('SpherePopular', *_T(R_MA, R_MA, _0_0)), # EPSG:3857 Spheroid

2469 Struve1860 = _lazy('Struve1860', *_T(6378298.3, 6356657.14266956, 294.73)),

2470# Walbeck = _lazy('Walbeck', *_T(6376896.0, _0_0, 302.78)),

2471# WarOffice = _lazy('WarOffice', *_T(6378300.0, _0_0, 296.0)),

2472 WGS60 = _lazy('WGS60', *_T(6378165.0, 6356783.28695944, 298.3)),

2473 WGS66 = _lazy('WGS66', *_T(6378145.0, 6356759.76948868, 298.25)),

2474 WGS72 = _lazy(_WGS72_, *_T(6378135.0, _0_0, 298.26)), # b=6356750.52

2475 WGS84 = _lazy(_WGS84_, *_T(R_MA, _0_0, _f__WGS84)), # b=6356752.3142451793

2476# U{NOAA/NOS/NGS/inverse<https://GitHub.com/noaa-ngs/inverse/blob/main/invers3d.f>}

2477 WGS84_NGS = _lazy('WGS84_NGS', *_T(R_MA, _0_0, 298.257222100882711243162836600094))

2478)

2479

2480_EWGS84 = Ellipsoids.WGS84 # (INTERNAL) shared

2481

2482if __name__ == _DMAIN_:

2483

2484 from pygeodesy import nameof, printf

2485

2486 def _e13(f):

2487 return fstr(f, prec=13, fmt=Fmt.e)

2488

2489 for E in (_EWGS84, Ellipsoids.GRS80, # NAD83,

2490 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

2491 Ellipsoid(_EWGS84.b, _EWGS84.a, name='_Prolate')):

2492 e = f2n(E.f) - E.n

2493 printf('# %s: %s', _DOT_('Ellipsoids', E.name), E.toStr(prec=10, terse=0), nl=1)

2494 printf('# (e-fmt) f_=%s, f=%s, e=%s, e2=%s, n=%s (%s)', _e13(E.f_), _e13(E.f),

2495 _e13(E.e), _e13(E.e2),

2496 _e13(E.n), fstr(e, prec=9, fmt=Fmt.e))

2497 printf('# %s %s', Ellipsoid.AlphaKs.name, fstr(E.AlphaKs, prec=20))

2498 printf('# %s %s', Ellipsoid.BetaKs.name, fstr(E.BetaKs, prec=20))

2499 printf('# %s %s', nameof(Ellipsoid.KsOrder), E.KsOrder) # property

2500

2501 from pygeodesy.internals import _pregistry

2502 # __doc__ of this file, force all into registry

2503 _pregistry(Ellipsoids)

2504

2505# % python3.13 -m pygeodesy.ellipsoids

2506

2507# Ellipsoids.WGS84: name='WGS84', a=6378137, f=0.0033528107, e=0.0818191908, b=6356752.3142451793, f_=298.257223563, e2=0.00669438, e21=0.99330562, e22=0.0067394967, e32=0.0033584313, n=0.0016792204, A=6367449.1458234144, L=10001965.7293127235, R1=6371008.7714150595, R2=6371007.1809184738, R3=6371000.7900091587, Rbiaxial=6367453.6345163295, Rtriaxial=6372797.5559594007

2508# (e-fmt) f_=2.98257223563e+02, f=3.3528106647475e-03, e=8.1819190842622e-02, e2=6.6943799901413e-03, n=1.6792203863837e-03 (0.0e+00)

2509# AlphaKs 0.00083773182062446994, 0.00000076085277735725, 0.00000000119764550324, 0.00000000000242917068, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2510# BetaKs 0.00083773216405794875, 0.0000000590587015222, 0.00000000016734826653, 0.00000000000021647981, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2511# KsOrder 8

2512

2513# Ellipsoids.GRS80: name='GRS80', a=6378137, f=0.0033528107, e=0.081819191, b=6356752.3141403468, f_=298.2572221009, e2=0.00669438, e21=0.99330562, e22=0.0067394968, e32=0.0033584313, n=0.0016792204, A=6367449.1457710434, L=10001965.7292304561, R1=6371008.7713801153, R2=6371007.1808835147, R3=6371000.7899741363, Rbiaxial=6367453.6344640013, Rtriaxial=6372797.5559332585

2514# (e-fmt) f_=2.9825722210088e+02, f=3.3528106811837e-03, e=8.1819191042833e-02, e2=6.6943800229036e-03, n=1.6792203946295e-03 (0.0e+00)

2515# AlphaKs 0.00083773182472890429, 0.00000076085278481561, 0.00000000119764552086, 0.00000000000242917073, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2516# BetaKs 0.0008377321681623882, 0.00000005905870210374, 0.000000000167348269, 0.00000000000021647982, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2517# KsOrder 8

2518

2519# Ellipsoids.Sphere: name='Sphere', a=6371008.7714149999, f=0, e=0, b=6371008.7714149999, f_=0, e2=0, e21=1, e22=0, e32=0, n=0, A=6371008.7714149999, L=10007557.1761167478, R1=6371008.7714149999, R2=6371008.7714149999, R3=6371008.7714149999, Rbiaxial=6371008.7714149999, Rtriaxial=6371008.7714149999

2520# (e-fmt) f_=0.0e+00, f=0.0e+00, e=0.0e+00, e2=0.0e+00, n=0.0e+00 (0.0e+00)

2521# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2522# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2523# KsOrder 8

2524

2525# Ellipsoids.SpherePopular: name='SpherePopular', a=6378137, f=0, e=0, b=6378137, f_=0, e2=0, e21=1, e22=0, e32=0, n=0, A=6378137, L=10018754.171394622, R1=6378137, R2=6378137, R3=6378137, Rbiaxial=6378137, Rtriaxial=6378137

2526# (e-fmt) f_=0.0e+00, f=0.0e+00, e=0.0e+00, e2=0.0e+00, n=0.0e+00 (0.0e+00)

2527# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2528# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2529# KsOrder 8

2530

2531# Ellipsoids._Prolate: name='_Prolate', a=6356752.3142451793, f=-0.0033640898, e=0.0820944379, b=6378137, f_=-297.257223563, e2=-0.0067394967, e21=1.0067394967, e22=-0.00669438, e32=-0.0033584313, n=-0.0016792204, A=6367449.1458234144, L=10001965.7293127235, R1=6363880.5428301189, R2=6363878.9413582645, R3=6363872.5644020075, Rbiaxial=6367453.6345163295, Rtriaxial=6362105.2243882557

2532# (e-fmt) f_=-2.97257223563e+02, f=-3.3640898209765e-03, e=8.2094437949696e-02, e2=-6.7394967422765e-03, n=-1.6792203863837e-03 (0.0e+00)

2533# AlphaKs -0.00084149152514366627, 0.00000076653480614871, -0.00000000120934503389, 0.0000000000024576225, -0.00000000000000578863, 0.00000000000000001502, -0.00000000000000000004, 0.0

2534# BetaKs -0.00084149187224351817, 0.00000005842735196773, -0.0000000001680487236, 0.00000000000021706261, -0.00000000000000038002, 0.00000000000000000073, -0.0, 0.0

2535# KsOrder 8

2536

2537# **) MIT License

2538#

2540#

2541# Permission is hereby granted, free of charge, to any person obtaining a

2542# copy of this software and associated documentation files (the "Software"),

2543# to deal in the Software without restriction, including without limitation

2544# the rights to use, copy, modify, merge, publish, distribute, sublicense,

2545# and/or sell copies of the Software, and to permit persons to whom the

2546# Software is furnished to do so, subject to the following conditions:

2547#

2548# The above copyright notice and this permission notice shall be included

2549# in all copies or substantial portions of the Software.

2550#

2551# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

2552# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

2553# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

2554# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

2555# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

2556# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

2557# OTHER DEALINGS IN THE SOFTWARE.