Coverage for pygeodesy/ellipses.py: 93%

1 

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

3 

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 ; 

8 

9# from pygeodesy.basics import islistuple # _MODS 

10from pygeodesy.constants import EPS, EPS_2, INT0, NEG0, PI, PI_2, PI3_2, PI2, \ 

11 _0_0, _1_0, _4_0, _isfinite, _over, _1_over # _N_1_0 

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

13# from pygeodesy.ellipsoids import Ellipsoid # _MODS 

14from pygeodesy.errors import _ConvergenceError, _ValueError, _xkwds, _xkwds_pop2 

15from pygeodesy.fmath import euclid, fhorner, fmean_, hypot 

16from pygeodesy.fsums import _fsum # PYCHOK used! 

17from pygeodesy.internals import typename, _DOT_, _UNDER_ 

18# from pygeodesy.interns import _DOT_, _UNDER_ # from .internals 

19from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

20from pygeodesy.named import _NamedBase, unstr 

21from pygeodesy.props import Property_RO, property_RO, property_ROnce 

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

23# from pygeodesy.triaxials import Triaxial_ , TriaxialError # _MODS 

24from pygeodesy.units import Degrees, Meter, Meter2, Radians, Radius, Scalar 

25from pygeodesy.utily import atan2, atan2b, atan2p, sincos2, sincos2d 

26# from pygeodesy.vector3d import Vector3d # _MODS 

27 

28from math import degrees, fabs, radians, sqrt 

29# import operator as _operator # from .fmath 

30 

31__all__ = _ALL_LAZY.ellipses 

32__version__ = '26.03.25' 

33 

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

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

36# assert _DIG53 == 53 

37 

38 

39class Ellipse(_NamedBase): 

40 '''Class to compute various attributes of a 2-D ellipse. 

41 ''' 

42# _ab3 = (a, b, a * b) # unordered 

43 _flat = False 

44 _maxit = _DIG53 

45# _Pab4 = (r, P, a, b) # a >= b, ordered 

46 _4_PI_ = _4_0 / PI - (14 / 11) 

47 _tEPS = None 

48 

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

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

51 

52 The ellipse is C{oblate} if C{B{a} > B{b}}, C{prolate} if 

53 C{B{a} < B{b}}, C{circular} if C{B{a} == B{b}} and C{"flat"} 

54 if C{min(B{a}, B{b}) <<< max(B{a}, B{b})}. 

55 

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

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

58 

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

60 

61 @see: U{Ellipse<https://MathWorld.Wolfram.com/Ellipse.html>}. 

62 ''' 

63 if name: 

64 self.name = name 

65 self._ab3 = a, b, (a * b) # unordered 

66 

67 r = a < b 

68 if r: # prolate 

69 a, b = b, a 

70 if b < 0 or not _isfinite(a): # PYCHOK no cover 

71 raise self._Error(None) 

72 if a > b: 

73 if _isFlat(a, b): 

74 self._flat = True 

75 P = a * _4_0 

76 else: # pro-/oblate 

77 P = None 

78 else: # circular 

79 P = a * PI2 

80# b = a 

81 self._Pab4 = r, P, a, b # ordered 

82 

83 @Property_RO 

84 def a(self): 

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

86 ''' 

87 a, _, _ = self._ab3 

88 return Meter(a=a) 

89 

90 @Property_RO 

91 def apses2(self): 

92 '''Get 2-tuple C{(apoapsis, periapsis)} with the U{apo-<https://MathWorld.Wolfram.com/Apoapsis.html>} 

93 and U{periapsis<https://MathWorld.Wolfram.com/Periapsis.html>} of this ellipse, both C{meter}. 

94 ''' 

95 _, _, a, p = self._Pab4 

96 c = self.c 

97 if c: # a != p 

98 p = a - c 

99 a += c 

100 return a, p 

101 

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

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

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

105 

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

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

108 

109 @return: Arc length, signed (C{meter}, conventionally). 

110 ''' 

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

112 

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

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

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

116 

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

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

119 

120 @return: Arc length, signed (C{meter}, conventionally). 

121 ''' 

122 r, R, a, _ = self._Pab4 

123 if R is None: 

124 _e = self._ellipe or self._ellipE 

125 k = self.e2 

126 r = PI_2 if r else _0_0 

127 R = _arc(_e, k, r + rad2) 

128 r += rad1 

129 if r: 

130 R -= _arc(_e, k, r) 

131 else: 

132 a = (rad2 - rad1) / PI2 

133 return Meter(arc=R * a) 

134 

135 @Property_RO 

136 def area(self): 

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

138 ''' 

139 _, _, ab = self._ab3 

140 return Meter2(area=ab * PI) 

141 

142 @Property_RO 

143 def b(self): 

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

145 ''' 

146 _, b, _ = self._ab3 

147 return Meter(b=b) 

148 

149 @Property_RO 

150 def c(self): 

151 '''Get the U{linear eccentricity<https://WikiPedia.org/wiki/Ellipse#Linear_eccentricity>} 

152 C{c}, I{unsigned} (C{meter}, conventionally). 

153 ''' 

154 return Meter(c=fabs(self.foci)) 

155 

156 @Property_RO 

157 def e(self): 

158 '''Get the eccentricity (C{scalar, 0 <= B{e} <= 1}). 

159 ''' 

160 e2 = self.e2 

161 return Scalar(e=sqrt(e2) if 0 < e2 < 1 else e2) 

162 

163 @Property_RO 

164 def e2(self): 

165 '''Get the eccentricity I{squared} (C{scalar, 0 <= B{e2} <= 1}). 

166 ''' 

167 # C{e2} is aka C{k}, Elliptic C{k2} and SciPy's C{m} 

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

169 e2 = ((_1_0 - (b / a)**2) if 0 < b else _1_0) if b < a else _0_0 

170 return Scalar(e2=e2) 

171 

172 @Property_RO 

173 def _Ek(self): 

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

175 ''' 

176 return _MODS.elliptic._Ek(self.e2) 

177 

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

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

180 ''' 

181 # assert k == self._Ek.k2 

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

183 

184 @property_ROnce 

185 def _ellipe(self): 

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

187 ''' 

188 try: 

189 from scipy.special import ellipe, ellipeinc 

190 

191 def _ellipe(k, phi=None): 

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

193 return float(r) 

194 

195 except (AttributeError, ImportError): 

196 _ellipe = None 

197 return _ellipe # overwrite property_ROnce 

198 

199 def _Error(self, where, **cause): # PYCHOK no cover 

200 '''(INTERNAL) Build an L{EllipseError}. 

201 ''' 

202 t = self.named3 

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

204 if where: 

205 t = typename(where, where) 

206 u = _DOT_(u, t) 

207 return EllipseError(u, **cause) 

208 

209 @Property_RO 

210 def foci(self): 

211 '''Get the U{linear eccentricity<https://WikiPedia.org/wiki/Ellipse#Linear_eccentricity>}, 

212 I{signed} (C{meter}, conventionally), C{positive} if this ellipse is oblate, C{negative} 

213 if prolate or C{0} if circular. See also property L{Ellipse.c}. 

214 ''' 

215 c = float(self.e) 

216 if c: 

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

218 c *= a 

219 if r: # prolate 

220 c = -c 

221 return Meter(foci=c) # signed 

222 

223 @property_ROnce 

224 def _GKs(self): 

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

226 ''' 

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

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

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

230 

231 def hartzell4(self, x, y, los=False): 

232 '''Compute the intersection of this ellipse with a Line-Of-Sight from Point-Of-View 

233 C{(B{x}, B{y})} I{outside} this ellipse. 

234 

235 @kwarg los: Line-Of-Sight, I{direction} to the ellipse (L{Los}, L{Vector3d}, 

236 L{Vector2Tuple} or 2-tuple C{(dx, dy)}) or C{True} for the I{normal, 

237 perpendicular, plumb} to this ellipse or C{False} or C{None} to 

238 point to its center. 

239 

240 @return: L{Vector4Tuple}C{(x, y, z, h)} with coordinates C{x}, C{y} and C{z=0} 

241 of the intersection and C{h} the distance to "Point-Of-View" C{(B{x}, 

242 B{y})} I{along the} B{C{los}}, all in C{meter}, conventionally. 

243 

244 @raise EllipseError: Invalid B{C{x}}, B{C{y}} or B{C{los}} or B{C{los}} points 

245 outside or away from this ellipse. 

246 

247 @see: Function L{hartzell4<triaxials.triaxial5.hartzell4>} for further details. 

248 ''' 

249 V3d = _MODS.vector3d.Vector3d 

250 if los not in (True, False, None): 

251 try: 

252 los = V3d(los.x, los.y, 0) 

253 except (AttributeError, TypeError): 

254 if _MODS.basics.islistuple(los, minum=2): 

255 los = V3d(*map(float, los[:2])) 

256 return self._triaxialX(self.hartzell4, V3d(x, y, 0), los=los) 

257 

258 def height4(self, x, y, **normal_eps): 

259 '''Compute the projection on and distance to this ellipse from a point C{(B{x}, B{y})} 

260 in- or outside this ellipse. 

261 

262 @kwarg normal_eps: With default C{B{normal}=True} the projection is I{perpendicular, 

263 plumb} to this ellipse, otherwise C{radially} to its center (C{bool}). 

264 Tolerance C{B{eps}=EPS} for root finding and validation (C{scalar}), 

265 use a negative value to skip validation. 

266 

267 @return: L{Vector4Tuple}C{(x, y, z, h)} with coordinates C{x}, C{y} and C{z=0} 

268 of the projection on or the intersection with the ellipse and C{h} the 

269 I{signed, normal distance} to the ellipse in C{meter}, conventionally. 

270 Positive C{h} indicates, C{x} and/or C{y} are outside the ellipse, 

271 negative C{h} means inside. 

272 

273 @raise EllipseError: Invalid B{C{x}}, B{C{y}} or B{C{eps}}, no convergence in 

274 root finding or validation failed. 

275 

276 @see: Methods L{Ellipse.normal3d}, L{Ellipse.normal4} and function L{height4 

277 <triaxials.triaxial5.height4>}. 

278 ''' 

279 return self._triaxialX(self.height4, x, y, 0, **normal_eps) 

280 

281 def _HGKs(self, h, maxit): 

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

283 ''' 

284 s = t = _1_0 

285 yield s 

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

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

288 t2 = t**2 

289 yield t2 

290 p = s 

291 s += t2 

292 if s == p: # 44 trips 

293 break 

294 else: # PYCHOK no cover 

295 raise _ConvergenceError(maxit, s, p) 

296 

297 @property_RO 

298 def isCircular(self): 

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

300 ''' 

301 return self.a == self.b 

302 

303 @property_RO 

304 def isFlat(self): 

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

306 ''' 

307 return self._flat 

308 

309 @property_RO 

310 def isOblate(self): 

311 '''Is this ellipse oblate (foci on semi-axis C{a})? (C{bool}) 

312 ''' 

313 return self.a > self.b 

314 

315 @property_RO 

316 def isProlate(self): 

317 '''Is this ellipse prolate (foci on semi-axis C{b})? (C{bool}) 

318 ''' 

319 return self.a < self.b 

320 

321 @Property_RO 

322 def lati(self): 

323 '''Get the U{semi-latus rectum<https://WikiPedia.org/wiki/Ellipse#Semi-latus_rectum>}, 

324 I{signed} (C{meter}, conventionally), C{positive} if this ellipse is oblate or 

325 circular, C{0} if "flat" and oblate, C{negative} if prolate or C{NEG0} if "flat" 

326 and prolate. See also property L{Ellipse.p}. 

327 ''' 

328 r, _, a, p = self._Pab4 

329 if 0 < p < a: 

330 p *= p / a 

331 if r: 

332 p = -p if p else NEG0 

333 return Meter(lati=p) # signed 

334 

335 def normal3d(self, deg_x, y=None, **length): 

336 '''Get a 3-D vector I{perpendicular to} this ellipse from point C{(B{x}, B{y})} 

337 C{on} this ellipse or at C{B{deg} degrees} along this ellipse. 

338 

339 @kwarg length: Optional, signed C{B{length}=1} in out-/inward direction 

340 (C{scalar}). 

341 

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

343 out- or inward for postive respectively negative B{C{length}}. 

344 

345 @raise EllipseError: Invalid B{C{x}} and/or B{C{y}}. 

346 

347 @see: Methods L{Ellipse.height4}, L{Ellipse.normal4}, L{Ellipse.sideOf} and 

348 C{Triaxial_.normal3d}. 

349 ''' 

350 return self._triaxialX(self.normal3d, *self._xy03(deg_x, y), **length) 

351 

352 def normal4(self, deg_x, y=None, **height_normal): 

353 '''Compute a point at B{C{height}} above or below this ellipse point C{(B{x}, 

354 B{y})} C{on} this ellipse or at C{B{deg} degrees} along this ellipse. 

355 

356 @kwarg height_normal: The desired distance C{B{height}=0} in- or outside this 

357 ellipse (C{meter}, conventionally) and C{B{normal}=True}, If 

358 C{B{normal}=True}, the B{C{height}} is I{perpendicular, plumb} 

359 to this ellipse, otherwise C{radially} to its center (C{bool}). 

360 

361 @return: L{Vector4Tuple}C{(x, y, z, h)} with coordinates C{x}, C{y} and C{z=0} 

362 and C{h} the I{signed, normal distance} to the ellipse in C{meter}, 

363 conventionally. Positive C{h} indicates, C{x} and/or C{y} are outside 

364 the ellipse, negative C{h} means inside. 

365 

366 @raise EllipseError: Invalid B{C{x}} and/or B{C{y}}. 

367 

368 @see: Methods L{Ellipse.height4}, L{Ellipse.normal3d}, L{Ellipse.sideOf} and 

369 C{Triaxial_.normal4}. 

370 ''' 

371 return self._triaxialX(self.normal4, *self._xy03(deg_x, y), **height_normal) 

372 

373 @Property_RO 

374 def p(self): 

375 '''Get the U{semi-latus rectum<https://WikiPedia.org/wiki/Ellipse#Semi-latus_rectum>} 

376 C{p (aka B{𝓁}, script-small-l)}, I{unsigned} (C{meter}, conventionally). 

377 ''' 

378 return Meter(p=fabs(self.lati)) 

379 

380 @Property_RO 

381 def perimeterAGM(self): 

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

383 <https://PaulBourke.net/geometry/ellipsecirc>} formula (C{meter}, conventionally). 

384 ''' 

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

386 if P is None: 

387 t = _TOL53 

388 m = -1 

389 c = a + b 

390 ds = [c**2] 

391 _d = ds.append 

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

393 b = sqrt(a * b) 

394 a = c * _0_5 

395 c = a + b 

396 d = a - b 

397 m *= 2 

398 _d(d**2 * m) 

399 if d <= (b * t): 

400 break 

401 else: # PYCHOK no cover 

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

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

404 return Meter(perimeterAGM=P) 

405 

406 @Property_RO 

407 def perimeter4Arc3(self): 

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

409 <https://PaulBourke.net/geometry/ellipsecirc>} (aka 4-Center) 

410 approximation as a 3-Tuple C{(P, Ra, Rb)} with perimeter C{P}, arc 

411 radii C{Ra} and C{Rb} at the respective semi-axes (all in C{meter}, 

412 conventionally). 

413 ''' 

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

415 if P is None: 

416 h = hypot(a, b) 

417 t = atan2(b, a) 

418 s, c = sincos2(t) 

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

420 a = _over(L, c) 

421 b = _over(h - L, s) 

422 P = (t * b + (PI_2 - t) * a) * _4_0 

423 elif a > b: # flat 

424 a, b = _0_0, _1_over(b) # INF 

425# else: # circular 

426# pass 

427 if r: 

428 a, b = b, a 

429 return Meter(perimeter4Arc=P), Radius(Ra=a), Radius(Rb=b) 

430 

431# @Property_RO 

432# def perimeterBPA(self): 

433# '''Compute the perimeter of this ellipse using the U{Bronshtein Padé 

434# Approximant <https://www.math.TTU.edu/~pearce/papers/schov.pdf>} 

435# (C{meter}, conventionally). 

436# ''' 

437# P, h = self._Ph2 

438# if h: 

439# h *= h 

440# P *= _over(h**2 * _3_0 - _64_0, h * _16_0 - _64_0) * PI 

441# return Meter(perimeterBPA=P) 

442 

443 @Property_RO 

444 def perimeterCR(self): 

445 '''Compute the perimeter of this ellipse using U{Rackauckas' 

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

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

448 U{here<http://www.EByte.IT/library/docs/math05a/EllipsePerimeterApprox05.html>} (C{meter}, conventionally). 

449 ''' 

450 P, h = self._Ph2 

451 if h: 

452 h *= h 

453 P *= _over(fhorner(h, 135168, -85760, -5568, 3867), 

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

455 return Meter(perimeterCR=P) 

456 

457 @Property_RO 

458 def perimeterGK(self): 

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

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

461 series, C{B{b / a} > 0.75} (C{meter}, conventionally). 

462 ''' 

463 P, h = self._Ph2 

464 if h: 

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

466 return Meter(perimeterGK=P) 

467 

468 @Property_RO 

469 def perimeterHGK(self): 

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

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

472 series (C{meter}, conventionally). 

473 ''' 

474 P, h = self._Ph2 

475 if h: 

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

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

478 return Meter(perimeterHGK=P) 

479 

480# @Property_RO 

481# def perimeterJWPA(self): 

482# '''Compute the perimeter of this ellipse using the U{Jacobson-Waadeland 

483# Padé Approximant <https://www.math.TTU.edu/~pearce/papers/schov.pdf>} 

484# (C{meter}, conventionally). 

485# ''' 

486# P, h = self._Ph2 

487# if h: 

488# h *= h 

489# P *= _over(fhorner(h, 256, -48, -21), 

490# fhorner(h, 256, -112, 3)) * PI 

491# return Meter(perimeterJWPA=P) 

492 

493 @Property_RO 

494 def perimeter2k(self): 

495 '''Compute the perimeter of this ellipse using the complete integral 

496 of the 2nd kind, C{Elliptic.cE} (C{meter}, conventionally). 

497 ''' 

498 return Meter(perimeter2k=self._perimeter2k(self._ellipE)) 

499 

500 @Property_RO 

501 def perimeter2k_(self): 

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

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

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

505 conventionally). 

506 ''' 

507 return Meter(perimeter2k_=self._perimeter2k(self._ellipe or self._ellipE)) 

508 

509 def _perimeter2k(self, _ellip): 

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

511 ''' 

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

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

514 k = self.e2 

515 P = _ellip(k) * a * _4_0 

516 return P 

517 

518# @Property_RO 

519# def perimeterLS(self): 

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

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

522# formula, aka C{3/2 norm} (C{meter}, conventionally). 

523# ''' 

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

525# if P is None: 

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

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

528# return Meter(perimeterLS=P) 

529 

530 @Property_RO 

531 def perimeter2R(self): 

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

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

534 C{B{b / a} > 0.9} (C{meter}, conventionally). 

535 ''' 

536 P, h = self._Ph2 

537 if h: 

538 P *= _2RC(h, _1_0) 

539 return Meter(perimeter2R=P) 

540 

541 @Property_RO 

542 def perimeter2RC(self): 

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

544 <http://www.EByte.IT/library/docs/math05a/EllipsePerimeterApprox05.html>} 

545 approximation, C{B{b / a} > 0.9} (C{meter}, conventionally). 

546 ''' 

547 P, h = self._Ph2 

548 if h: 

549 P *= _2RC(h, pow(h, 24) * self._4_PI_ + _1_0) 

550 return Meter(perimeter2RC=P) 

551 

552# @Property_RO 

553# def perimeterSPA(self): 

554# '''Compute the perimeter of this ellipse using the U{Selmer Padé Approximant 

555# <http://www.EByte.IT/library/docs/math05a/EllipsePerimeterApprox05.html>} 

556# (C{meter}, conventionally). 

557# ''' 

558# P, h = self._Ph2 

559# if h: 

560# h *= h 

561# P *= _over(h * _3_0 + _16_0, _16_0 - h) * PI 

562# return Meter(perimeterSPA=P) 

563 

564 @Property_RO 

565 def _Ph2(self): 

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

567 if P is None: 

568 if b: 

569 P = a + b 

570 h = (a - b) / P 

571 else: 

572 P = a 

573 h = _1_0 

574 else: 

575 h = None 

576 return P, h 

577 

578 def point(self, deg_x, y=None): 

579 '''Return the point I{on} this ellipse at C{B{deg}} or C{atan2d(B{y}, B{x}) 

580 degrees} along this ellipse. 

581 

582 @return: A 2-tuple C{(x, y)}. 

583 ''' 

584 s, c = sincos2d(deg_x) if y is None else self._sc2(deg_x, y, eps=None) 

585 return (self.a * c), (self.b * s) 

586 

587 def points(self, np, nq=4, ccw=False, ended=False, eps=EPS): # MCCABE 13 

588 '''Yield up to C{np} points along this ellipse, each a 2-tuple C{(x, y)}, 

589 starting at semi-axis C{+a}, in (counter-)clockwise order and distributed 

590 evenly along the minor semi-axis. 

591 

592 @arg np: Number of points to generate (C{int}). 

593 @kwarg nq: Number of quarters to cover (C{int}, 1..4). 

594 @kwarg ccw: Use C{B{ccw}=True} for counter-clockwise order (C{bool}). 

595 @kwarg ended: If C{True}, include the last quadrant's end point (C{bool}). 

596 @kwarg eps: Tolerance for duplicate points (C{meter}, conventionally). 

597 

598 @see: U{Directrix<https://MathWorld.Wolfram.com/ConicSectionDirectrix.html>}. 

599 ''' 

600 a, b, _ = self._ab3 

601 if min(a, b) > eps and not self.isFlat: 

602 q = max(min(int(nq), 4), 1) 

603 n = max( int(np) // q, 1) 

604 if not ccw: 

605 b = -b 

606 ps = list(_q1ps(a, b, n, eps)) 

607 for p in ps: # 1st quadrant 

608 yield p 

609 p0 = ps.pop(0) # E 

610 pq = _0_0, b # N/S 

611 if q > 1: 

612 yield pq 

613 for x, y in reversed(ps): # 2nd 

614 yield (-x), y 

615 pq = (-a), _0_0 # W 

616 if q > 2: 

617 yield pq 

618 for x, y in ps: # 3rd 

619 yield (-x), (-y) 

620 pq = _0_0, (-b) # S/N 

621 if q > 3: 

622 yield pq 

623 for x, y in reversed(ps): # 4th 

624 yield x, (-y) 

625 pq = p0 

626 if ended: 

627 yield pq 

628 else: # "flat" 

629 p0 = a, b 

630 yield p0 

631 if max(a, b) > eps: 

632 yield (-a), (-b) 

633 if ended: 

634 yield p0 

635 

636 def polar2d(self, deg_x, y=None): 

637 '''For a point at C{B{deg}} or C{atan2d(B{y}, B{x}) degrees} along this 

638 ellipse, return 2-tuple C{(radius, angle)} with the polar U{radius 

639 <https://WikiPedia.org/wiki/Ellipse#Polar_form_relative_to_center>} 

640 from the center (C{meter}, conventionally) and C{angle} in C{degrees}. 

641 ''' 

642 s, c = sincos2d(deg_x) if y is None else self._sc2(deg_x, y, eps=None) 

643 a, b, ab = self._ab3 

644 r = (_over(ab, hypot(a * s, b * c)) if c else b) if s else a 

645 return r, atan2b(s, c) 

646 

647 @Property_RO 

648 def R1(self): 

649 '''Get this ellipse' I{arithmetic mean} radius, C{(2 * a + b) / 3} (C{meter}, conventionally). 

650 ''' 

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

652 if r: 

653 r = fmean_(a, a, r) if a > r else a 

654 return Radius(R1=r or _0_0) 

655 

656 @Property_RO 

657 def R2(self): 

658 '''Get this ellipse' I{authalic} radius, C{sqrt(B{a} * B{b})} (C{meter}, conventionally). 

659 ''' 

660 a, b, ab = self._ab3 

661 return Radius(R2=(sqrt(ab) if a != b else float(a)) if ab else _0_0) 

662 

663 Rauthalic = Rgeometric = R2 

664 

665 def Roc(self, deg_x, y=None, eps=None): 

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

667 at a point C{B{deg}} or C{atan2d(B{y}, B{x}) degrees} along this ellipse. 

668 

669 @see: Method L{Roc_<Ellipse.Roc_>} for ruther details. 

670 ''' 

671 x = radians(deg_x) if y is None else deg_x 

672 return self.Roc_(x, y, eps=eps) 

673 

674 def Roc_(self, rad_x, y=None, eps=None): 

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

676 at a point C{B{rad}} or C{atan2(B{y}, B{x}) radians} along this ellipse. 

677 

678 @kwarg eps: See method C{sideOf}, use C{B{eps}=0} to permit any points. 

679 

680 @return: Curvature (C{meter}, conventionally). 

681 

682 @raise ValueError: Point C{(B{x}, B{y})} off this ellipse, unless C{B{eps}=0}. 

683 ''' 

684 try: 

685 a, b, ab = self._ab3 

686 if b != a: 

687 s, c = sincos2(rad_x) if y is None else self._sc2(rad_x, y, eps=eps) 

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

689 else: # circular 

690 r = float(a) 

691 except Exception as X: 

692 raise self._Error(self.Roc_, cause=X) 

693 return Radius(Roc=r) 

694 

695 @Property_RO 

696 def Rrectifying(self): 

697 '''Get this ellipse' I{rectifying} radius, C{perimeter2k_ / PI2} (C{meter}, conventionally). 

698 ''' 

699 return Radius(Rrectifying=self.perimeter2k_ / PI2) 

700 

701 def _sc2(self, x, y, eps=EPS): 

702 '''(INTERNAL) Helper for methods C{.point}, C{.polar}, C{.Roc_} and C{.slope_}. 

703 ''' 

704 if eps and eps > 0: 

705 s = self._sideOf(x, y, eps) 

706 if s: 

707 raise _ValueError(x=x, y=y, eps=eps, sideOf=s) 

708 h = hypot(x, y) 

709 s = _over(y, h) 

710 c = _over(x, h) 

711 return s, c 

712 

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

714 '''(INTERNAL) Helper for methods C{._sc2} and C{.sideOf}. 

715 ''' 

716 a, b, ab = self._ab3 

717 s = ab or max(a, b) 

718 if s: 

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

720# s = max(_N_1_0, min(_1_0, s)) 

721 else: # dot 

722 s = _1_0 if x or y else _0_0 

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

724 

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

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

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

728 ''' 

729 try: 

730 return Scalar(sideOf=self._sideOf(x, y, eps)) 

731 except Exception as X: 

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

733 

734 def slope(self, deg_x, y=None, eps=None): 

735 '''Compute the U{tangent slope<https://WikiPedia.org/wiki/Ellipse#Tangent_slope_as_parameter>} 

736 at a point C{B{deg}} or C{atan2d(B{y}, B{x}) degrees} along this ellipse. 

737 

738 @return: Slope (C{degrees}), negative for C{0 <= B{deg} < 90}. 

739 

740 @see: Method L{slope_<Ellipse.slope_>} for further details. 

741 ''' 

742 x = radians(deg_x) if y is None else deg_x 

743 return Degrees(slope=degrees(self.slope_(x, y, eps=eps))) 

744 

745 def slope_(self, rad_x, y=None, eps=None): 

746 '''Compute the U{tangent slope<https://WikiPedia.org/wiki/Ellipse#Tangent_slope_as_parameter>} 

747 at a point C{B{rad}} or C{atan2(B{y}, B{x}) radians} along this ellipse. 

748 

749 @kwarg eps: See method C{sideOf}, use C{B{eps}=0} to permit any points. 

750 

751 @return: Slope (C{radians}), negative for C{0 <= B{rad} < PI/2}. 

752 

753 @raise ValueError: C{(B{x}, B{y})} off this ellipse, unless C{B{eps}=0}. 

754 ''' 

755 # <https://UNacademy.com/content/jee/study-material/mathematics/equation-of-a-tangent-to-the-ellipse/> 

756 s, c = sincos2(rad_x) if y is None else self._sc2(rad_x, y, eps=eps) 

757 r = atan2p(-self.b * c, self.a * s) 

758 if r >= PI3_2: 

759 r -= PI2 

760 return Radians(slope=r or _0_0) # no -0.0 

761 

762 def toEllipsoid(self, **Ellipsoid_and_kwds): 

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

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

765 

766 @kwarg Ellipsoid_and_kwds: Optional C{B{Ellipsoid}=Ellipsoid} class 

767 and additional C{Ellipsoid} keyword arguments. 

768 ''' 

769 E, kwds = _xkwds_pop2(Ellipsoid_and_kwds, Ellipsoid= 

770 _MODS.ellipsoids.Ellipsoid) 

771 return E(self.a, b=self.b, **_xkwds(kwds, name=self.name)) 

772 

773 def toStr(self, prec=8, terse=2, **sep_name): # PYCHOK signature 

774 '''Return this ellipse as a text string. 

775 

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

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

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

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

780 to exclude this ellipse' name and separator 

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

782 

783 @return: This C{Ellipse}' attributes (C{str}). 

784 ''' 

785 E = Ellipse 

786 t = E.a, E.b 

787 if (terse or 0) != 2: 

788 t += E.c, E.e, E.e2, E.p, E.area, E.perimeter2k, E.R2 

789 if terse: 

790 t = t[:terse] 

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

792 

793 def toTriaxial_(self, c=EPS, **Triaxial_and_kwds): # like .Ellipse5Tuple.toTriaxial_ 

794 '''Return a L{Triaxial_<pygeodesy.Triaxial_>} from this ellipse' semi-axes. 

795 

796 @kwarg c: Near-zero, minor semi-axis (C{meter}, conventionally). 

797 @kwarg Triaxial_and_kwds: Optional C{B{Triaxial}=Triaxial_} class and 

798 additional C{Triaxial} keyword arguments. 

799 ''' 

800 T, kwds = _xkwds_pop2(Triaxial_and_kwds, Triaxial=_MODS.triaxials.Triaxial_) 

801 return T(self.a, b=self.b, c=c, **_xkwds(kwds, name=self.name or _UNDER_)) # 'NN' 

802 

803 def _triaxialX(self, method, *args, **kwds): 

804 '''(INTERNAL) Invoke a triaxial method and map exceptions to L{EllipseError}s. 

805 ''' 

806 try: 

807 t = self._tEPS 

808 if t is None: 

809 self._tEPS = t = self.toTriaxial_(EPS) 

810 _m = getattr(t, method.__name__) 

811 return _m(*args, **kwds) 

812 except Exception as x: 

813 raise self._Error(method, Triaxial_=t, cause=x) 

814 

815 def _xy03(self, deg_x, y): 

816 if y is None: 

817 y, x = sincos2d(deg_x) 

818 y *= self.b 

819 x *= self.a 

820 else: 

821 x = float(deg_x) 

822 y = float(y) 

823 return x, y, 0 

824 

825 

826class EllipseError(_ValueError): 

827 '''Raised for any L{Ellipse} or C{ellipses} issue. 

828 ''' 

829 pass # ... 

830 

831 

832def _arc(_e, k, r): 

833 # in C{Ellipse.arc_} 

834 t, r = divmod(r, PI2) 

835 R = _e(k, r) # phi=r 

836 if t: # + t * perimeter 

837 t *= _e(k) * _4_0 

838 R += t 

839 return R 

840 

841 

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

843 # is C{b <<< a}? 

844 return b < (a * _TOL53_53) 

845 

846 

847def _q1ps(a, b, n, eps): 

848 # yield the 1st quadrant C{Ellipse.points} 

849 if a > b: # oblate