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

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, 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.27' 

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 = _4_0 * a 

76# b = _0_0 

77 else: # pro-/oblate 

78 P = None 

79 else: # circular 

80 P = a * PI2 

81# b = a 

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

83 

84 @Property_RO 

85 def a(self): 

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

87 ''' 

88 a, _, _ = self._ab3 

89 return Meter(a=a) 

90 

91 @Property_RO 

92 def apses2(self): 

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

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

95 ''' 

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

97 c = self.c 

98 if c: # a != p 

99 p = a - c 

100 a += c 

101 return a, p 

102 

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

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

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

106 

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

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

109 

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

111 ''' 

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

113 

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

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

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

117 

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

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

120 

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

122 ''' 

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

124 if R is None: 

125 _e = self._ellipe or self._ellipE 

126 k = self.e2 

127 r = PI_2 if r else _0_0 

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

129 r += rad1 

130 if r: 

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

132 else: 

133 a = (rad2 - rad1) / PI2 

134 return Meter(arc=R * a) 

135 

136 @Property_RO 

137 def area(self): 

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

139 ''' 

140 _, _, ab = self._ab3 

141 return Meter2(area=ab * PI) 

142 

143 @Property_RO 

144 def b(self): 

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

146 ''' 

147 _, b, _ = self._ab3 

148 return Meter(b=b) 

149 

150 @Property_RO 

151 def c(self): 

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

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

154 ''' 

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

156 

157 @Property_RO 

158 def e(self): 

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

160 ''' 

161 e2 = self.e2 

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

163 

164 @Property_RO 

165 def e2(self): 

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

167 ''' 

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

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

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

171 return Scalar(e2=e2) 

172 

173 @Property_RO 

174 def _Ek(self): 

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

176 ''' 

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

178 

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

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

181 ''' 

182 # assert k == self._Ek.k2 

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

184 

185 @property_ROnce 

186 def _ellipe(self): 

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

188 ''' 

189 try: 

190 from scipy.special import ellipe, ellipeinc 

191 

192 def _ellipe(k, phi=None): 

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

194 return float(r) 

195 

196 except (AttributeError, ImportError): 

197 _ellipe = None 

198 return _ellipe # overwrite property_ROnce 

199 

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

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

202 ''' 

203 t = self.named3 

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

205 if where: 

206 t = typename(where, where) 

207 u = _DOT_(u, t) 

208 return EllipseError(u, **cause) 

209 

210 @Property_RO 

211 def foci(self): 

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

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

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

215 ''' 

216 c = float(self.e) 

217 if c: 

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

219 c *= a 

220 if r: # prolate 

221 c = -c 

222 return Meter(foci=c) # signed 

223 

224 @property_ROnce 

225 def _GKs(self): 

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

227 ''' 

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

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

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

231 

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

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

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

235 

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

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

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

239 point to its center. 

240 

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

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

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

244 

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

246 outside or away from this ellipse. 

247 

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

249 ''' 

250 V3d = _MODS.vector3d.Vector3d 

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

252 try: 

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

254 except (AttributeError, TypeError): 

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

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

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

258 

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

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

261 in- or outside this ellipse. 

262 

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

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

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

266 use a negative value to skip validation. 

267 

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

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

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

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

272 negative C{h} means inside. 

273 

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

275 root finding or validation failed. 

276 

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

278 <triaxials.triaxial5.height4>}. 

279 ''' 

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

281 

282 def _HGKs(self, h, maxit): 

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

284 ''' 

285 s = t = _1_0 

286 yield s 

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

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

289 t2 = t**2 

290 yield t2 

291 p = s 

292 s += t2 

293 if s == p: # 44 trips 

294 break 

295 else: # PYCHOK no cover 

296 raise _ConvergenceError(maxit, s, p) 

297 

298 @property_RO 

299 def isCircular(self): 

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

301 ''' 

302 return self.a == self.b 

303 

304 @property_RO 

305 def isFlat(self): 

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

307 ''' 

308 return self._flat 

309 

310 @property_RO 

311 def isOblate(self): 

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

313 ''' 

314 return self.a > self.b 

315 

316 @property_RO 

317 def isProlate(self): 

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

319 ''' 

320 return self.a < self.b 

321 

322 @Property_RO 

323 def lati(self): 

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

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

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

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

328 ''' 

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

330 if 0 < p < a: 

331 p *= p / a 

332 if r: 

333 p = (-p) if p else NEG0 

334 return Meter(lati=p) # signed 

335 

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

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

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

339 

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

341 (C{scalar}). 

342 

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

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

345 

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

347 

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

349 C{Triaxial_.normal3d}. 

350 ''' 

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

352 

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

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

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

356 

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

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

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

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

361 

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

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

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

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

366 

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

368 

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

370 C{Triaxial_.normal4}. 

371 ''' 

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

373 

374 @Property_RO 

375 def p(self): 

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

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

378 ''' 

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

380 

381 @Property_RO 

382 def perimeterAGM(self): 

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

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

385 ''' 

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

387 if P is None: 

388 t = _TOL53 

389 m = -1 

390 c = a + b 

391 ds = [c**2] 

392 _d = ds.append 

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

394 b = sqrt(a * b) 

395 a = c * _0_5 

396 c = a + b 

397 d = a - b 

398 m *= 2 

399 _d(d**2 * m) 

400 if d <= (b * t): 

401 break 

402 else: # PYCHOK no cover 

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

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

405 return Meter(perimeterAGM=P) 

406 

407 @Property_RO 

408 def perimeter4Arc3(self): 

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

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

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

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

413 conventionally). 

414 ''' 

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

416 if P is None: 

417 h = hypot(a, b) 

418 t = atan2(b, a) 

419 s, c = sincos2(t) 

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

421 a = _over(L, c) 

422 b = _over(h - L, s) 

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

424 elif a > b: # flat 

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

426# else: # circular 

427# pass 

428 if r: 

429 a, b = b, a 

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

431 

432# @Property_RO 

433# def perimeterBPA(self): 

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

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

436# (C{meter}, conventionally). 

437# ''' 

438# P, h = self._Ph2 

439# if h: 

440# h *= h 

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

442# return Meter(perimeterBPA=P) 

443 

444 @Property_RO 

445 def perimeterCR(self): 

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

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

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

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

450 ''' 

451 P, h = self._Ph2 

452 if h: 

453 h *= h 

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

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

456 return Meter(perimeterCR=P) 

457 

458 @Property_RO 

459 def perimeterGK(self): 

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

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

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

463 ''' 

464 P, h = self._Ph2 

465 if h: 

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

467 return Meter(perimeterGK=P) 

468 

469 @Property_RO 

470 def perimeterHGK(self): 

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

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

473 series (C{meter}, conventionally). 

474 ''' 

475 P, h = self._Ph2 

476 if h: 

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

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

479 return Meter(perimeterHGK=P) 

480 

481# @Property_RO 

482# def perimeterJWPA(self): 

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

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

485# (C{meter}, conventionally). 

486# ''' 

487# P, h = self._Ph2 

488# if h: 

489# h *= h 

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

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

492# return Meter(perimeterJWPA=P) 

493 

494 @Property_RO 

495 def perimeter2k(self): 

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

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

498 ''' 

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

500 

501 @Property_RO 

502 def perimeter2k_(self): 

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

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

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

506 conventionally). 

507 ''' 

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

509 

510 def _perimeter2k(self, _ellip): 

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

512 ''' 

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

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

515 k = self.e2 

516 P = _ellip(k) * a * _4_0 

517 return P 

518 

519# @Property_RO 

520# def perimeterLS(self): 

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

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

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

524# ''' 

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

526# if P is None: 

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

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

529# return Meter(perimeterLS=P) 

530 

531 @Property_RO 

532 def perimeter2R(self): 

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

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

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

536 ''' 

537 P, h = self._Ph2 

538 if h: 

539 P *= _2RC(h, _1_0) 

540 return Meter(perimeter2R=P) 

541 

542 @Property_RO 

543 def perimeter2RC(self): 

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

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

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

547 ''' 

548 P, h = self._Ph2 

549 if h: 

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

551 return Meter(perimeter2RC=P) 

552 

553# @Property_RO 

554# def perimeterSPA(self): 

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

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

557# (C{meter}, conventionally). 

558# ''' 

559# P, h = self._Ph2 

560# if h: 

561# h *= h 

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

563# return Meter(perimeterSPA=P) 

564 

565 @Property_RO 

566 def _Ph2(self): 

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

568 if P is None: 

569 if b: 

570 P = a + b 

571 h = (a - b) / P 

572 else: 

573 P = a 

574 h = _1_0 

575 else: 

576 h = None 

577 return P, h 

578 

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

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

581 degrees} along this ellipse. 

582 

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

584 ''' 

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

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

587 

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

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

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

591 evenly along the minor semi-axis. 

592 

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

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

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

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

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

598 

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

600 ''' 

601 a, b, _ = self._ab3 

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

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

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

605 if not ccw: 

606 b = -b 

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

608 for p in ps: # 1st quadrant 

609 yield p 

610 p0 = ps.pop(0) # E 

611 pq = _0_0, b # N/S 

612 if q > 1: 

613 yield pq 

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

615 yield (-x), y 

616 pq = (-a), _0_0 # W 

617 if q > 2: 

618 yield pq 

619 for x, y in ps: # 3rd 

620 yield (-x), (-y) 

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

622 if q > 3: 

623 yield pq 

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

625 yield x, (-y) 

626 pq = p0 

627 if ended: 

628 yield pq 

629 else: # "flat" 

630 p0 = a, b 

631 yield p0 

632 if max(a, b) > eps: 

633 yield (-a), (-b) 

634 if ended: 

635 yield p0 

636 

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

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

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

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

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

642 ''' 

643 return polar2d(*self.point(deg_x, y)) 

644 

645 @Property_RO 

646 def R1(self): 

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

648 ''' 

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

650 if r: 

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

652 return Radius(R1=r or _0_0) 

653 

654 @Property_RO 

655 def R2(self): 

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

657 ''' 

658 a, b, ab = self._ab3 

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

660 

661 Rauthalic = Rgeometric = R2 

662 

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

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

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

666 

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

668 ''' 

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

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

671 

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

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

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

675 

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

677 

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

679 

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

681 ''' 

682 try: 

683 a, b, ab = self._ab3 

684 if b != a: 

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

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

687 else: # circular 

688 r = float(a) 

689 except Exception as X: 

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

691 return Radius(Roc=r) 

692 

693 @Property_RO 

694 def Rrectifying(self): 

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

696 ''' 

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

698 

699 def _sc2(self, x, y, eps): 

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

701 ''' 

702 if eps and eps > 0: 

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

704 if s: 

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

706 h = hypot(x, y) 

707 s = _over(y, h) 

708 c = _over(x, h) 

709 return s, c 

710 

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

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

713 ''' 

714 a, b, ab = self._ab3 

715 s = ab or max(a, b) 

716 if s: 

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

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

719 else: # dot 

720 s = _1_0 if x or y else _0_0 

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

722 

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

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

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

726 ''' 

727 try: 

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

729 except Exception as X: 

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

731 

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

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

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

735 

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

737 

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

739 ''' 

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

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

742 

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

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

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

746 

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

748 

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

750 

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

752 ''' 

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

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

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

756 if r >= PI3_2: 

757 r -= PI2 

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

759 

760 def toEllipsoid(self, **Ellipsoid_and_kwds): 

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

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

763 

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

765 and additional C{Ellipsoid} keyword arguments. 

766 ''' 

767 E, kwds = _xkwds_pop2(Ellipsoid_and_kwds, Ellipsoid= 

768 _MODS.ellipsoids.Ellipsoid) 

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

770 

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

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

773 

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

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

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

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

778 to exclude this ellipse' name and separator 

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

780 

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

782 ''' 

783 E = Ellipse 

784 t = E.a, E.b 

785 if (terse or 0) != 2: 

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

787 if terse: 

788 t = t[:terse] 

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

790 

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

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

793 

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

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

796 additional C{Triaxial} keyword arguments. 

797 ''' 

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

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

800 

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

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

803 ''' 

804 try: 

805 t = self._tEPS 

806 if t is None: 

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

808 _m = getattr(t, method.__name__) 

809 return _m(*args, **kwds) 

810 except Exception as x: 

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

812 

813 def _xy03(self, deg_x, y): 

814 if y is None: 

815 y, x = sincos2d(deg_x) 

816 y *= self.b 

817 x *= self.a 

818 else: 

819 x = float(deg_x) 

820 y = float(y) 

821 return x, y, 0 

822 

823 

824class EllipseError(_ValueError): 

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

826 ''' 

827 pass # ... 

828 

829 

830def _arc(_e, k, r): 

831 # in C{Ellipse.arc_} 

832 t, r = divmod(r, PI2) 

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

834 if t: # + t * perimeter 

835 t *= _e(k) * _4_0 

836 R += t 

837 return R 

838 

839 

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

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

842 return b < (a * _TOL53_53) 

843 

844 

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

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

847 if a > b: # oblate 

848 def _yx2(i): 

849 y = i / n 

850 return y, sqrt(_1_0 - y**2) 

851 

852 elif a < b: # prolate 

853 def _yx2(i): # PYCHOK redef