Coverage for pygeodesy/ellipses.py: 92%

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

4u'''Class C{Ellipse} for 2-D ellipse attributes, like perimeter, area, etc.

5'''

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

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

9from pygeodesy.constants import EPS, EPS_2, INT0, PI, PI_2, PI2, \

10 _0_0, _1_0, _4_0, _over, _1_over

11from pygeodesy.constants import _0_5, _3_0, _10_0, MANT_DIG as _DIG53 # PYCHOK used!

12# from pygeodesy.ellipsoids import Ellipsoid # _MODS

13from pygeodesy.errors import _ConvergenceError, _ValueError

14from pygeodesy.fmath import fhorner, hypot

15from pygeodesy.fsums import _fsum # PYCHOK used!

16from pygeodesy.internals import typename, _DOT_

17# from pygeodesy.interns import _DOT_ # from .internals

18# from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS # from .utily

19from pygeodesy.named import _Named, unstr

20from pygeodesy.props import Property_RO, property_RO, property_ROnce

21# from pygeodesy.streprs import unstr # from .named

22# from pygeodesy.triaxials import Triaxial_ # _MODS

23from pygeodesy.utily import atan2, sincos2, _ALL_LAZY, _MODS

25from math import fabs, radians, sqrt

26# import operator as _operator # from .fmath

28__all__ = _ALL_LAZY.ellipses

29__version__ = '26.02.12'

31_TOL53 = sqrt(EPS_2) # sqrt(pow(_0_5, _DIG53))

32_TOL53_53 = _TOL53 / _DIG53 # "flat" b/a tolerance, 1.9e-10

33# assert _DIG53 == 53

36class Ellipse(_Named):

37 '''Class to compute attributes of a 2-D ellipse, like perimeter, area and arcs.

38 '''

39 _flat = False

40 _maxit = _DIG53

42 def __init__(self, a, b, **name):

43 '''New L{Ellipse} with semi-axes B{C{a}} and B{C{b}}.

45 The ellipse is oblate if C{a > b}, prolate if C{a < b}

46 circular if C{a == b} and "flat" if C{min(a, b) near 0}.

48 @arg a: X semi-axis length (C{meter}, conventionally).

49 @arg b: Y semi-axis length (C{meter}, conventionally).

51 @raise ValueError: Invalid B{C{a}} or B{C{b}}.

52 '''

53 if name:

54 self.name = name

55 self._a = a

56 self._b = b

57 r = a < b

58 if r: # prolate

59 a, b = b, a

60 if b < 0: # PYCHOK no cover

61 raise self._Error(None)

62 if a > b:

63 if _isFlat(a, b):

64 self._flat = True

65 P = a * _4_0

66 else: # pro-/oblate

67 P = None

68 else: # circle

69 P = a * PI2

70 self._Pab4 = r, P, a, b

72 @Property_RO

73 def a(self):

74 '''Return semi-axis C{B{a}} of this ellipse (C{meter}, conventionally).

75 '''

76 return self._a

78 def arc(self, deg2, deg1=0):

79 '''Compute the length of U{elliptic arc<https://www.JohnDCook.com/blog/2022/11/02/elliptic-arc-length/>}

80 C{(B{deg2} - B{deg1})}, both counter-clockwise from semi-axis B{C{a}} to B{C{b}} of the ellipse.

82 @arg deg2: End angle of the elliptic arc (C{degrees}).

83 @kwarg deg1: Start angle of the elliptic arc (C{degrees}).

85 @return: Arc length, signed (C{meter}, same units as semi-axes B{C{a}} and B{C{b}}).

86 '''

87 return self.arc_(radians(deg2), (radians(deg1) if deg1 else _0_0))

89 def arc_(self, rad2, rad1=0):

90 '''Compute the length of U{elliptic arc<https://www.JohnDCook.com/blog/2022/11/02/elliptic-arc-length/>}

91 C{(B{rad2} - B{rad1})}, both counter-clockwise from semi-axis B{C{a}} to B{C{b}} of the ellipse.

93 @arg rad2: End angle of the elliptic arc (C{radians}).

94 @kwarg rad1: Start angle of the elliptic arc (C{radians}).

96 @return: Arc length, signed (C{meter}, same units as semi-axes B{C{a}} and B{C{b}}).

97 '''

98 r, L, a, _ = self._Pab4

99 if L is None:

100 _e = self._ellipe or self._ellipE

101 k = self._k

102 r = PI_2 if r else _0_0

103 L = self._arc(_e, k, r + rad2)

104 r += rad1

105 if r:

106 L -= self._arc(_e, k, r)

107 L *= a

108 else:

109 L *= (rad2 - rad1) / PI2

110 return L

111

112 def _arc(self, _e, k, r):

113 '''(INTERNAL) Helper for method C{.arc_}.

114 '''

115 t, r = divmod(r, PI2)

116 L = _e(k, r) # phi=r

117 if t: # + t * perimeter

118 t *= _e(k) * _4_0

119 L += t

120 return L

121

122 @Property_RO

123 def area(self):

124 '''Return the area of this ellipse (C{meter**2}, conventionally).

125 '''

126 return self.a * self.b * PI

127

128 @Property_RO

129 def b(self):

130 '''Return semi-axis C{B{b}} of this ellipse (C{meter}, conventionally).

131 '''

132 return self._b

133

134 @Property_RO

135 def _Ek(self):

136 '''(INTERNAL) Get the C{Elliptic(k)} instance.

137 '''

138 return _MODS.elliptic._Ek(self._k)

139

140 def _ellipE(self, k, phi=None): # PYCHOK k

141 '''(INTERNAL) Get the in-/complete integral of the 2nd kind.

142 '''

143 # assert k == self._Ek.k2

144 return self._Ek.cE if phi is None else self._Ek.fE(phi)

145

146 @property_ROnce

147 def _ellipe(self):

148 '''(INTERNAL) Wrap functions C{scipy.special.ellipe} and C{-.ellipeinc}, I{once}.

149 '''

150 try:

151 from scipy.special import ellipe, ellipeinc

152

153 def _ellipe(k, phi=None):

154 r = ellipe(k) if phi is None else ellipeinc(phi, k)

155 return float(r)

156

157 except (AttributeError, ImportError):

158 _ellipe = None

159 return _ellipe # overwrite property_ROnce

160

161 def _Error(self, where, **cause):

162 '''(INTERNAL) Build an error.

163 '''

164 t = self.named3

165 u = unstr(t, a=self.a, b=self.b)

166 if where:

167 t = typename(where, where)

168 u = _DOT_(u, t)

169 return _ValueError(u, **cause)

170

171 @Property_RO

172 def foci(self):

173 '''Get the foci lengths (C{meter}, conventionally), C{positive} if

174 the ellipse is oblate with foci on semi-axis C{a}, C{negative}

175 if prolate with foci on semi-axis C{b} or C{0} if circular with

176 coincident foci in the ellipse' center.

177 '''

178 c = self._k

179 if c:

180 r, _, a, _ = self._Pab4

181 c = sqrt(c) * a

182 if r: # prolate

183 c = -c

184 return c

185

186 @property_ROnce

187 def _GKs(self):

188 '''(INTERNAL) Compute the coefficients for property C{.perimeterGK}, I{once}.

189 '''

190 # U{numerators<https://OEIS.org/A056981>}, U{denominators<https://OEIS.org/A056982>}

191 return (1, 1 / 4, 1 / 64, 1 / 256, 25 / 16384, 49 / 65536,

192 441 / 1048576, 1089 / 4194304) # overwrite property_ROnce

193

194 def _HGKs(self, h, maxit):

195 '''(INTERNAL) Yield the terms for property C{.perimeterHGK}.

196 '''

197 s = t = _1_0

198 yield s

199 for u in range(-1, maxit * 2, 2):

200 t *= u / (u + 3) * h

201 t2 = t**2

202 yield t2

203 p = s

204 s += t2

205 if s == p: # 44 trips

206 break

207 else: # PYCHOK no cover

208 raise _ConvergenceError(maxit, s, p)

209

210 @property_RO

211 def isCircular(self):

212 '''Is this ellipse circular? (C{bool})

213 '''

214 return self.a == self.b

215

216 @property_RO

217 def isFlat(self):

218 '''Is this ellipse "flat", too pro-/oblate? (C{bool})

219 '''

220 return self._flat

221

222 @property_RO

223 def isOblate(self):

224 '''Is this ellipse oblate? (C{bool})

225 '''

226 return self.a > self.b

227

228 @property_RO

229 def isProlate(self):

230 '''Is this ellipse prolate? (C{bool})

231 '''

232 return self.a < self.b

233

234 @Property_RO

235 def _k(self):

236 '''(INTERNAL) Get C{0 <= k <= 1}.

237 '''

238 # C{k} is aka Elliptic C{k2} and SciPy's C{m}

239 _, _, a, b = self._Pab4

240 return ((_1_0 - (b / a)**2) if 0 < b else _1_0) if b < a else _0_0

241

242 @Property_RO

243 def perimeterAGM(self):

244 '''Compute the perimeter of this ellipse using the U{Arithmetic-Geometric Mean

245 <https://PaulBourke.net/geometry/ellipsecirc>} formula (C{meter}, same units

246 as semi-axes B{C{a}} and B{C{b}}).

247 '''

248 _, P, a, b = self._Pab4

249 if P is None:

250 t = _TOL53

251 m = -1

252 c = a + b

253 ds = [c**2]

254 _d = ds.append

255 for _ in range(self._maxit): # 4..5 trips

256 b = sqrt(a * b)

257 a = c * _0_5

258 c = a + b

259 d = a - b

260 m *= 2

261 _d(d**2 * m)

262 if d <= (b * t):

263 break

264 else: # PYCHOK no cover

265 raise _ConvergenceError(self._maxit, _over(d, b), t)

266 P = _over(_fsum(ds) * PI, c) # nonfinites=True

267 return P

268

269 @Property_RO

270 def perimeter4Arc3(self):

271 '''Compute the perimeter (and arcs) of this ellipse using the U{4-Arc

272 <https://PaulBourke.net/geometry/ellipsecirc>} approximation as a

273 3-Tuple C{(P, Ra, Rb)} with perimeter C{P}, arc radii C{Ra} and

274 C{Rb} at the respective semi-axes (all in C{meter}, same units as

275 semi-axes B{C{a}} and B{C{b}}).

276 '''

277 r, P, a, b = self._Pab4

278 if P is None:

279 h = hypot(a, b)

280 t = atan2(b, a)

281 s, c = sincos2(t)

282 L = (h - (a - b)) * _0_5

283 Ra = _over(L, c)

284 Rb = _over(h - L, s)

285 P = (t * Rb + (PI_2 - t) * Ra) * _4_0

286 elif a > b: # flat

287 Ra, Rb = _0_0, _1_over(b) # INF

288 else: # circle

289 Ra, Rb = a, b

290 return (P, Rb, Ra) if r else (P, Ra, Rb)

291

292# @Property_RO

293# def perimeterCR(self):

294# '''Compute the perimeter of this ellipse using U{Rackauckas'

295# <https://www.ChrisRackauckas.com/assets/Papers/ChrisRackauckas-The_Circumference_of_an_Ellipse.pdf>}

296# approximation, also U{here<https://ExtremeLearning.com.AU/a-formula-for-the-perimeter-of-an-ellipse>}

297# (C{meter}, same units as semi-axes B{C{a}} and B{C{b}}).

298# '''

299# _, P, a, b = self._Pab4

300# if P is None:

301# P = a + b

302# h = ((a - b) / P)**2

303# P *= (fhorner(h, 135168, -85760, -5568, 3867) /

304# fhorner(h, 135168, -119552, 22208, 345)) * PI

305# return P

306

307 @Property_RO

308 def perimeterGK(self):

309 '''Compute the perimeter of this ellipse using the U{Gauss-Kummer

310 <https://www.JohnDCook.com/blog/2023/05/28/approximate-ellipse-perimeter>} series,

311 C{B{b / a} > 0.75} (C{meter}, same units as semi-axes B{C{a}} and B{C{b}}).

312 '''

313 P, h = self._Ph2

314 if h:

315 P *= fhorner(h**2, *self._GKs) * PI

316 return P

317

318 @Property_RO

319 def perimeter2k(self):

320 '''Compute the perimeter of this ellipse using the complete integral of the 2nd

321 kind, C{Elliptic.cE} (C{meter}, same units as semi-axes B{C{a}} and B{C{b}}).

322 '''

323 return self._perimeter2k(self._ellipE)

324

325 @Property_RO

326 def perimeter2k_(self):

327 '''Compute the perimeter of this ellipse using U{SciPy's ellipe

328 <https://www.JohnDCook.com/perimeter_ellipse.html>} function

329 if available, otherwise use property C{perimeter2k} (C{meter},

330 same units as semi-axes B{C{a}} and B{C{b}}).

331 '''

332 return self._perimeter2k(self._ellipe or self._ellipE)

333

334 def _perimeter2k(self, _ellip):

335 '''(INTERNAL) Helper for methods C{.PE2k} and C{.Pe2k}.

336 '''

337 _, P, a, _ = self._Pab4

338 if P is None: # see .ellipsoids.Ellipsoid.L

339 k = self._k

340 P = _ellip(k) * a * _4_0

341 return P

342

343 @Property_RO

344 def _Ph2(self):

345 _, P, a, b = self._Pab4

346 if P is None:

347 P = a + b

348 h = (a - b) / P

349 else:

350 h = None

351 return P, h

352

353 @Property_RO

354 def perimeterHGK(self):

355 '''Compute the perimeter of this ellipse using the U{Hypergeometric Gauss-Kummer

356 <https://web.Tecnico.ULisboa.PT/~mcasquilho/compute/com/,ellips/PerimeterOfEllipse.pdf>}

357 series (C{meter}, same units as semi-axes B{C{a}} and B{C{b}}).

358 '''

359 P, h = self._Ph2

360 if h:

361 hs = self._HGKs(h, self._maxit)

362 P *= _fsum(hs) * PI # nonfinites=True

363 return P

364

365# @Property_RO

366# def perimeterLS(self):

367# '''Compute the perimeter of this ellipse using the U{Linderholm-Segal

368# <https://www.JohnDCook.com/blog/2021/03/24/perimeter-of-an-ellipse>}

369# formula, aka C{3/2 norm} (C{meter}, same units as semi-axes B{C{a}} and B{C{b}}).

370# '''

371# _, P, a, b = self._Pab4

372# if P is None:

373# n = pow(a, _1_5) + pow(b, _1_5)

374# P = pow(n * _0_5, _2_3rd) * PI2

375# return P

376

377 @Property_RO

378 def perimeter2R(self):

379 '''Compute the perimeter of this ellipse using U{Ramanujan's 2nd

380 <https://PaulBourke.net/geometry/ellipsecirc>} approximation,

381 C{B{b / a} > 0.9} (C{meter}, same units as semi-axes B{C{a}} and B{C{b}}).

382 '''

383 P, h = self._Ph2

384 if h:

385 h *= _3_0 * h

386 h /= sqrt(_4_0 - h) + _10_0 # /= chokes PyChecker?

387 P *= (h + _1_0) * PI

388 return P

389

390 @Property_RO

391 def R2(self):

392 '''Compute the I{authalic} radius of this ellipse, C{sqrt(B{a} * B{b})}

393 (C{meter}, same units as semi-axes B{C{a}} and B{C{b}}).

394 '''

395 a, b = self.a, self.b

396 return sqrt(a * b) if a != b else float(a)

397

398 def Roc_(self, rad_x, *y):

399 '''Compute the U{radius of curvature<https://WikiPedia.org/wiki/Radius_of_curvature>}

400 at angle C{B{rad} radians} or at point C{(B{x}, B{y})} I{on} this ellipse.

401

402 @return: Curvature (C{meter}, same units as semi-axes B{C{a}} and B{C{b}}).

403

404 @raise ValueError: C{(B{x}, B{y})} not located on the ellipse.

405 '''

406 try:

407 a, b = self.a, self.b

408 if b != a:

409 r = float(rad_x)

410 if y:

411 x, y = r, float(y[0])

412 if self._sideOf(x, y, EPS):

413 raise _ValueError(x=x, y=y)

414 r = atan2(y, x)

415 s, c = sincos2(r)

416 r = _over(hypot(s * a, c * b)**3, a * b)

417 else: # circle

418 r = float(a)

419 except Exception as x:

420 raise self._Error(self.Roc_, cause=x)

421 return r

422

423 def _sideOf(self, x, y, eps):

424 '''(INTERNAL) Helper for methods C{.Roc} and C{.sideOf}.

425 '''

426 a, b = self.a, self.b

427 s = hypot(x * b, y * a) - (a * b)

428 return INT0 if fabs(s) < eps else s

429

430 def sideOf(self, x, y, eps=EPS):

431 '''Return a C{positive}, C{negative} or C{0} scalar if point C{(B{x}, B{y})}

432 is C{outside}, C{inside} respectively C{on} this ellipse.

433 '''

434 try:

435 return self._sideOf(x, y, eps)

436 except Exception as X:

437 raise self._Error(self.sideOf, x=x, y=y, cause=X)

438

439 def toEllipsoid(self):

440 '''Return an L{Ellipsoid<pygeodesy.Ellipsoid>} from this ellipse's

441 C{a} and C{b} semi-axes.

442 '''

443 return _MODS.ellipsoids.Ellipsoid(self.a, b=self.b, name=self.name)

444

445 def toTriaxial_(self, c=EPS):

446 '''Return a L{Triaxial_<pygeodesy.Triaxial_>} from this ellipse's

447 C{a} and C{b} semi-axes with B{C{c}} as minor semi-axis.

448 '''

449 return _MODS.triaxials.Triaxial_(self.a, b=self.b, c=c, name=self.name) # 'NN'

450

451

452def _isFlat(a, b): # in .triaxials.bases

453 '''(INTERNAL) Is C{b <<< a}?

454 '''

455 return b < (a * _TOL53_53)

456

457# **) MIT License

458#

460#

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

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

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

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

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

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

467#

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

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

470#

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

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

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

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

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

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

477# OTHER DEALINGS IN THE SOFTWARE.