1from sympy.core import Expr, S, Symbol, oo, pi, sympify 2from sympy.core.compatibility import as_int, ordered 3from sympy.core.symbol import _symbol, Dummy, symbols 4from sympy.functions.elementary.complexes import sign 5from sympy.functions.elementary.piecewise import Piecewise 6from sympy.functions.elementary.trigonometric import cos, sin, tan 7from sympy.geometry.exceptions import GeometryError 8from sympy.logic import And 9from sympy.matrices import Matrix 10from sympy.simplify import simplify 11from sympy.utilities import default_sort_key 12from sympy.utilities.iterables import has_dups, has_variety, uniq, rotate_left, least_rotation 13from sympy.utilities.misc import func_name 14 15from .entity import GeometryEntity, GeometrySet 16from .point import Point 17from .ellipse import Circle 18from .line import Line, Segment, Ray 19 20import warnings 21 22 23class Polygon(GeometrySet): 24 """A two-dimensional polygon. 25 26 A simple polygon in space. Can be constructed from a sequence of points 27 or from a center, radius, number of sides and rotation angle. 28 29 Parameters 30 ========== 31 32 vertices : sequence of Points 33 34 Optional parameters 35 ========== 36 37 n : If > 0, an n-sided RegularPolygon is created. See below. 38 Default value is 0. 39 40 Attributes 41 ========== 42 43 area 44 angles 45 perimeter 46 vertices 47 centroid 48 sides 49 50 Raises 51 ====== 52 53 GeometryError 54 If all parameters are not Points. 55 56 See Also 57 ======== 58 59 sympy.geometry.point.Point, sympy.geometry.line.Segment, Triangle 60 61 Notes 62 ===== 63 64 Polygons are treated as closed paths rather than 2D areas so 65 some calculations can be be negative or positive (e.g., area) 66 based on the orientation of the points. 67 68 Any consecutive identical points are reduced to a single point 69 and any points collinear and between two points will be removed 70 unless they are needed to define an explicit intersection (see examples). 71 72 A Triangle, Segment or Point will be returned when there are 3 or 73 fewer points provided. 74 75 Examples 76 ======== 77 78 >>> from sympy import Polygon, pi 79 >>> p1, p2, p3, p4, p5 = [(0, 0), (1, 0), (5, 1), (0, 1), (3, 0)] 80 >>> Polygon(p1, p2, p3, p4) 81 Polygon(Point2D(0, 0), Point2D(1, 0), Point2D(5, 1), Point2D(0, 1)) 82 >>> Polygon(p1, p2) 83 Segment2D(Point2D(0, 0), Point2D(1, 0)) 84 >>> Polygon(p1, p2, p5) 85 Segment2D(Point2D(0, 0), Point2D(3, 0)) 86 87 The area of a polygon is calculated as positive when vertices are 88 traversed in a ccw direction. When the sides of a polygon cross the 89 area will have positive and negative contributions. The following 90 defines a Z shape where the bottom right connects back to the top 91 left. 92 93 >>> Polygon((0, 2), (2, 2), (0, 0), (2, 0)).area 94 0 95 96 When the the keyword `n` is used to define the number of sides of the 97 Polygon then a RegularPolygon is created and the other arguments are 98 interpreted as center, radius and rotation. The unrotated RegularPolygon 99 will always have a vertex at Point(r, 0) where `r` is the radius of the 100 circle that circumscribes the RegularPolygon. Its method `spin` can be 101 used to increment that angle. 102 103 >>> p = Polygon((0,0), 1, n=3) 104 >>> p 105 RegularPolygon(Point2D(0, 0), 1, 3, 0) 106 >>> p.vertices[0] 107 Point2D(1, 0) 108 >>> p.args[0] 109 Point2D(0, 0) 110 >>> p.spin(pi/2) 111 >>> p.vertices[0] 112 Point2D(0, 1) 113 114 """ 115 116 def __new__(cls, *args, n = 0, **kwargs): 117 if n: 118 args = list(args) 119 # return a virtual polygon with n sides 120 if len(args) == 2: # center, radius 121 args.append(n) 122 elif len(args) == 3: # center, radius, rotation 123 args.insert(2, n) 124 return RegularPolygon(*args, **kwargs) 125 126 vertices = [Point(a, dim=2, **kwargs) for a in args] 127 128 # remove consecutive duplicates 129 nodup = [] 130 for p in vertices: 131 if nodup and p == nodup[-1]: 132 continue 133 nodup.append(p) 134 if len(nodup) > 1 and nodup[-1] == nodup[0]: 135 nodup.pop() # last point was same as first 136 137 # remove collinear points 138 i = -3 139 while i < len(nodup) - 3 and len(nodup) > 2: 140 a, b, c = nodup[i], nodup[i + 1], nodup[i + 2] 141 if Point.is_collinear(a, b, c): 142 nodup.pop(i + 1) 143 if a == c: 144 nodup.pop(i) 145 else: 146 i += 1 147 148 vertices = list(nodup) 149 150 if len(vertices) > 3: 151 return GeometryEntity.__new__(cls, *vertices, **kwargs) 152 elif len(vertices) == 3: 153 return Triangle(*vertices, **kwargs) 154 elif len(vertices) == 2: 155 return Segment(*vertices, **kwargs) 156 else: 157 return Point(*vertices, **kwargs) 158 159 @property 160 def area(self): 161 """ 162 The area of the polygon. 163 164 Notes 165 ===== 166 167 The area calculation can be positive or negative based on the 168 orientation of the points. If any side of the polygon crosses 169 any other side, there will be areas having opposite signs. 170 171 See Also 172 ======== 173 174 sympy.geometry.ellipse.Ellipse.area 175 176 Examples 177 ======== 178 179 >>> from sympy import Point, Polygon 180 >>> p1, p2, p3, p4 = map(Point, [(0, 0), (1, 0), (5, 1), (0, 1)]) 181 >>> poly = Polygon(p1, p2, p3, p4) 182 >>> poly.area 183 3 184 185 In the Z shaped polygon (with the lower right connecting back 186 to the upper left) the areas cancel out: 187 188 >>> Z = Polygon((0, 1), (1, 1), (0, 0), (1, 0)) 189 >>> Z.area 190 0 191 192 In the M shaped polygon, areas do not cancel because no side 193 crosses any other (though there is a point of contact). 194 195 >>> M = Polygon((0, 0), (0, 1), (2, 0), (3, 1), (3, 0)) 196 >>> M.area 197 -3/2 198 199 """ 200 area = 0 201 args = self.args 202 for i in range(len(args)): 203 x1, y1 = args[i - 1].args 204 x2, y2 = args[i].args 205 area += x1*y2 - x2*y1 206 return simplify(area) / 2 207 208 @staticmethod 209 def _isright(a, b, c): 210 """Return True/False for cw/ccw orientation. 211 212 Examples 213 ======== 214 215 >>> from sympy import Point, Polygon 216 >>> a, b, c = [Point(i) for i in [(0, 0), (1, 1), (1, 0)]] 217 >>> Polygon._isright(a, b, c) 218 True 219 >>> Polygon._isright(a, c, b) 220 False 221 """ 222 ba = b - a 223 ca = c - a 224 t_area = simplify(ba.x*ca.y - ca.x*ba.y) 225 res = t_area.is_nonpositive 226 if res is None: 227 raise ValueError("Can't determine orientation") 228 return res 229 230 @property 231 def angles(self): 232 """The internal angle at each vertex. 233 234 Returns 235 ======= 236 237 angles : dict 238 A dictionary where each key is a vertex and each value is the 239 internal angle at that vertex. The vertices are represented as 240 Points. 241 242 See Also 243 ======== 244 245 sympy.geometry.point.Point, sympy.geometry.line.LinearEntity.angle_between 246 247 Examples 248 ======== 249 250 >>> from sympy import Point, Polygon 251 >>> p1, p2, p3, p4 = map(Point, [(0, 0), (1, 0), (5, 1), (0, 1)]) 252 >>> poly = Polygon(p1, p2, p3, p4) 253 >>> poly.angles[p1] 254 pi/2 255 >>> poly.angles[p2] 256 acos(-4*sqrt(17)/17) 257 258 """ 259 260 # Determine orientation of points 261 args = self.vertices 262 cw = self._isright(args[-1], args[0], args[1]) 263 264 ret = {} 265 for i in range(len(args)): 266 a, b, c = args[i - 2], args[i - 1], args[i] 267 ang = Ray(b, a).angle_between(Ray(b, c)) 268 if cw ^ self._isright(a, b, c): 269 ret[b] = 2*S.Pi - ang 270 else: 271 ret[b] = ang 272 return ret 273 274 @property 275 def ambient_dimension(self): 276 return self.vertices[0].ambient_dimension 277 278 @property 279 def perimeter(self): 280 """The perimeter of the polygon. 281 282 Returns 283 ======= 284 285 perimeter : number or Basic instance 286 287 See Also 288 ======== 289 290 sympy.geometry.line.Segment.length 291 292 Examples 293 ======== 294 295 >>> from sympy import Point, Polygon 296 >>> p1, p2, p3, p4 = map(Point, [(0, 0), (1, 0), (5, 1), (0, 1)]) 297 >>> poly = Polygon(p1, p2, p3, p4) 298 >>> poly.perimeter 299 sqrt(17) + 7 300 """ 301 p = 0 302 args = self.vertices 303 for i in range(len(args)): 304 p += args[i - 1].distance(args[i]) 305 return simplify(p) 306 307 @property 308 def vertices(self): 309 """The vertices of the polygon. 310 311 Returns 312 ======= 313 314 vertices : list of Points 315 316 Notes 317 ===== 318 319 When iterating over the vertices, it is more efficient to index self 320 rather than to request the vertices and index them. Only use the 321 vertices when you want to process all of them at once. This is even 322 more important with RegularPolygons that calculate each vertex. 323 324 See Also 325 ======== 326 327 sympy.geometry.point.Point 328 329 Examples 330 ======== 331 332 >>> from sympy import Point, Polygon 333 >>> p1, p2, p3, p4 = map(Point, [(0, 0), (1, 0), (5, 1), (0, 1)]) 334 >>> poly = Polygon(p1, p2, p3, p4) 335 >>> poly.vertices 336 [Point2D(0, 0), Point2D(1, 0), Point2D(5, 1), Point2D(0, 1)] 337 >>> poly.vertices[0] 338 Point2D(0, 0) 339 340 """ 341 return list(self.args) 342 343 @property 344 def centroid(self): 345 """The centroid of the polygon. 346 347 Returns 348 ======= 349 350 centroid : Point 351 352 See Also 353 ======== 354 355 sympy.geometry.point.Point, sympy.geometry.util.centroid 356 357 Examples 358 ======== 359 360 >>> from sympy import Point, Polygon 361 >>> p1, p2, p3, p4 = map(Point, [(0, 0), (1, 0), (5, 1), (0, 1)]) 362 >>> poly = Polygon(p1, p2, p3, p4) 363 >>> poly.centroid 364 Point2D(31/18, 11/18) 365 366 """ 367 A = 1/(6*self.area) 368 cx, cy = 0, 0 369 args = self.args 370 for i in range(len(args)): 371 x1, y1 = args[i - 1].args 372 x2, y2 = args[i].args 373 v = x1*y2 - x2*y1 374 cx += v*(x1 + x2) 375 cy += v*(y1 + y2) 376 return Point(simplify(A*cx), simplify(A*cy)) 377 378 379 def second_moment_of_area(self, point=None): 380 """Returns the second moment and product moment of area of a two dimensional polygon. 381 382 Parameters 383 ========== 384 385 point : Point, two-tuple of sympifyable objects, or None(default=None) 386 point is the point about which second moment of area is to be found. 387 If "point=None" it will be calculated about the axis passing through the 388 centroid of the polygon. 389 390 Returns 391 ======= 392 393 I_xx, I_yy, I_xy : number or sympy expression 394 I_xx, I_yy are second moment of area of a two dimensional polygon. 395 I_xy is product moment of area of a two dimensional polygon. 396 397 Examples 398 ======== 399 400 >>> from sympy import Polygon, symbols 401 >>> a, b = symbols('a, b') 402 >>> p1, p2, p3, p4, p5 = [(0, 0), (a, 0), (a, b), (0, b), (a/3, b/3)] 403 >>> rectangle = Polygon(p1, p2, p3, p4) 404 >>> rectangle.second_moment_of_area() 405 (a*b**3/12, a**3*b/12, 0) 406 >>> rectangle.second_moment_of_area(p5) 407 (a*b**3/9, a**3*b/9, a**2*b**2/36) 408 409 References 410 ========== 411 412 https://en.wikipedia.org/wiki/Second_moment_of_area 413 414 """ 415 416 I_xx, I_yy, I_xy = 0, 0, 0 417 args = self.vertices 418 for i in range(len(args)): 419 x1, y1 = args[i-1].args 420 x2, y2 = args[i].args 421 v = x1*y2 - x2*y1 422 I_xx += (y1**2 + y1*y2 + y2**2)*v 423 I_yy += (x1**2 + x1*x2 + x2**2)*v 424 I_xy += (x1*y2 + 2*x1*y1 + 2*x2*y2 + x2*y1)*v 425 A = self.area 426 c_x = self.centroid[0] 427 c_y = self.centroid[1] 428 # parallel axis theorem 429 I_xx_c = (I_xx/12) - (A*(c_y**2)) 430 I_yy_c = (I_yy/12) - (A*(c_x**2)) 431 I_xy_c = (I_xy/24) - (A*(c_x*c_y)) 432 if point is None: 433 return I_xx_c, I_yy_c, I_xy_c 434 435 I_xx = (I_xx_c + A*((point[1]-c_y)**2)) 436 I_yy = (I_yy_c + A*((point[0]-c_x)**2)) 437 I_xy = (I_xy_c + A*((point[0]-c_x)*(point[1]-c_y))) 438 439 return I_xx, I_yy, I_xy 440 441 442 def first_moment_of_area(self, point=None): 443 """ 444 Returns the first moment of area of a two-dimensional polygon with 445 respect to a certain point of interest. 446 447 First moment of area is a measure of the distribution of the area 448 of a polygon in relation to an axis. The first moment of area of 449 the entire polygon about its own centroid is always zero. Therefore, 450 here it is calculated for an area, above or below a certain point 451 of interest, that makes up a smaller portion of the polygon. This 452 area is bounded by the point of interest and the extreme end 453 (top or bottom) of the polygon. The first moment for this area is 454 is then determined about the centroidal axis of the initial polygon. 455 456 References 457 ========== 458 459 https://skyciv.com/docs/tutorials/section-tutorials/calculating-the-statical-or-first-moment-of-area-of-beam-sections/?cc=BMD 460 https://mechanicalc.com/reference/cross-sections 461 462 Parameters 463 ========== 464 465 point: Point, two-tuple of sympifyable objects, or None (default=None) 466 point is the point above or below which the area of interest lies 467 If ``point=None`` then the centroid acts as the point of interest. 468 469 Returns 470 ======= 471 472 Q_x, Q_y: number or sympy expressions 473 Q_x is the first moment of area about the x-axis 474 Q_y is the first moment of area about the y-axis 475 A negative sign indicates that the section modulus is 476 determined for a section below (or left of) the centroidal axis 477 478 Examples 479 ======== 480 481 >>> from sympy import Point, Polygon 482 >>> a, b = 50, 10 483 >>> p1, p2, p3, p4 = [(0, b), (0, 0), (a, 0), (a, b)] 484 >>> p = Polygon(p1, p2, p3, p4) 485 >>> p.first_moment_of_area() 486 (625, 3125) 487 >>> p.first_moment_of_area(point=Point(30, 7)) 488 (525, 3000) 489 """ 490 if point: 491 xc, yc = self.centroid 492 else: 493 point = self.centroid 494 xc, yc = point 495 496 h_line = Line(point, slope=0) 497 v_line = Line(point, slope=S.Infinity) 498 499 h_poly = self.cut_section(h_line) 500 v_poly = self.cut_section(v_line) 501 502 poly_1 = h_poly[0] if h_poly[0].area <= h_poly[1].area else h_poly[1] 503 poly_2 = v_poly[0] if v_poly[0].area <= v_poly[1].area else v_poly[1] 504 505 Q_x = (poly_1.centroid.y - yc)*poly_1.area 506 Q_y = (poly_2.centroid.x - xc)*poly_2.area 507 508 return Q_x, Q_y 509 510 511 def polar_second_moment_of_area(self): 512 """Returns the polar modulus of a two-dimensional polygon 513 514 It is a constituent of the second moment of area, linked through 515 the perpendicular axis theorem. While the planar second moment of 516 area describes an object's resistance to deflection (bending) when 517 subjected to a force applied to a plane parallel to the central 518 axis, the polar second moment of area describes an object's 519 resistance to deflection when subjected to a moment applied in a 520 plane perpendicular to the object's central axis (i.e. parallel to 521 the cross-section) 522 523 References 524 ========== 525 526 https://en.wikipedia.org/wiki/Polar_moment_of_inertia 527 528 Examples 529 ======== 530 531 >>> from sympy import Polygon, symbols 532 >>> a, b = symbols('a, b') 533 >>> rectangle = Polygon((0, 0), (a, 0), (a, b), (0, b)) 534 >>> rectangle.polar_second_moment_of_area() 535 a**3*b/12 + a*b**3/12 536 """ 537 second_moment = self.second_moment_of_area() 538 return second_moment[0] + second_moment[1] 539 540 541 def section_modulus(self, point=None): 542 """Returns a tuple with the section modulus of a two-dimensional 543 polygon. 544 545 Section modulus is a geometric property of a polygon defined as the 546 ratio of second moment of area to the distance of the extreme end of 547 the polygon from the centroidal axis. 548 549 References 550 ========== 551 552 https://en.wikipedia.org/wiki/Section_modulus 553 554 Parameters 555 ========== 556 557 point : Point, two-tuple of sympifyable objects, or None(default=None) 558 point is the point at which section modulus is to be found. 559 If "point=None" it will be calculated for the point farthest from the 560 centroidal axis of the polygon. 561 562 Returns 563 ======= 564 565 S_x, S_y: numbers or SymPy expressions 566 S_x is the section modulus with respect to the x-axis 567 S_y is the section modulus with respect to the y-axis 568 A negative sign indicates that the section modulus is 569 determined for a point below the centroidal axis 570 571 Examples 572 ======== 573 574 >>> from sympy import symbols, Polygon, Point 575 >>> a, b = symbols('a, b', positive=True) 576 >>> rectangle = Polygon((0, 0), (a, 0), (a, b), (0, b)) 577 >>> rectangle.section_modulus() 578 (a*b**2/6, a**2*b/6) 579 >>> rectangle.section_modulus(Point(a/4, b/4)) 580 (-a*b**2/3, -a**2*b/3) 581 """ 582 x_c, y_c = self.centroid 583 if point is None: 584 # taking x and y as maximum distances from centroid 585 x_min, y_min, x_max, y_max = self.bounds 586 y = max(y_c - y_min, y_max - y_c) 587 x = max(x_c - x_min, x_max - x_c) 588 else: 589 # taking x and y as distances of the given point from the centroid 590 y = point.y - y_c 591 x = point.x - x_c 592 593 second_moment= self.second_moment_of_area() 594 S_x = second_moment[0]/y 595 S_y = second_moment[1]/x 596 597 return S_x, S_y 598 599 600 @property 601 def sides(self): 602 """The directed line segments that form the sides of the polygon. 603 604 Returns 605 ======= 606 607 sides : list of sides 608 Each side is a directed Segment. 609 610 See Also 611 ======== 612 613 sympy.geometry.point.Point, sympy.geometry.line.Segment 614 615 Examples 616 ======== 617 618 >>> from sympy import Point, Polygon 619 >>> p1, p2, p3, p4 = map(Point, [(0, 0), (1, 0), (5, 1), (0, 1)]) 620 >>> poly = Polygon(p1, p2, p3, p4) 621 >>> poly.sides 622 [Segment2D(Point2D(0, 0), Point2D(1, 0)), 623 Segment2D(Point2D(1, 0), Point2D(5, 1)), 624 Segment2D(Point2D(5, 1), Point2D(0, 1)), Segment2D(Point2D(0, 1), Point2D(0, 0))] 625 626 """ 627 res = [] 628 args = self.vertices 629 for i in range(-len(args), 0): 630 res.append(Segment(args[i], args[i + 1])) 631 return res 632 633 @property 634 def bounds(self): 635 """Return a tuple (xmin, ymin, xmax, ymax) representing the bounding 636 rectangle for the geometric figure. 637 638 """ 639 640 verts = self.vertices 641 xs = [p.x for p in verts] 642 ys = [p.y for p in verts] 643 return (min(xs), min(ys), max(xs), max(ys)) 644 645 def is_convex(self): 646 """Is the polygon convex? 647 648 A polygon is convex if all its interior angles are less than 180 649 degrees and there are no intersections between sides. 650 651 Returns 652 ======= 653 654 is_convex : boolean 655 True if this polygon is convex, False otherwise. 656 657 See Also 658 ======== 659 660 sympy.geometry.util.convex_hull 661 662 Examples 663 ======== 664 665 >>> from sympy import Point, Polygon 666 >>> p1, p2, p3, p4 = map(Point, [(0, 0), (1, 0), (5, 1), (0, 1)]) 667 >>> poly = Polygon(p1, p2, p3, p4) 668 >>> poly.is_convex() 669 True 670 671 """ 672 # Determine orientation of points 673 args = self.vertices 674 cw = self._isright(args[-2], args[-1], args[0]) 675 for i in range(1, len(args)): 676 if cw ^ self._isright(args[i - 2], args[i - 1], args[i]): 677 return False 678 # check for intersecting sides 679 sides = self.sides 680 for i, si in enumerate(sides): 681 pts = si.args 682 # exclude the sides connected to si 683 for j in range(1 if i == len(sides) - 1 else 0, i - 1): 684 sj = sides[j] 685 if sj.p1 not in pts and sj.p2 not in pts: 686 hit = si.intersection(sj) 687 if hit: 688 return False 689 return True 690 691 def encloses_point(self, p): 692 """ 693 Return True if p is enclosed by (is inside of) self. 694 695 Notes 696 ===== 697 698 Being on the border of self is considered False. 699 700 Parameters 701 ========== 702 703 p : Point 704 705 Returns 706 ======= 707 708 encloses_point : True, False or None 709 710 See Also 711 ======== 712 713 sympy.geometry.point.Point, sympy.geometry.ellipse.Ellipse.encloses_point 714 715 Examples 716 ======== 717 718 >>> from sympy import Polygon, Point 719 >>> p = Polygon((0, 0), (4, 0), (4, 4)) 720 >>> p.encloses_point(Point(2, 1)) 721 True 722 >>> p.encloses_point(Point(2, 2)) 723 False 724 >>> p.encloses_point(Point(5, 5)) 725 False 726 727 References 728 ========== 729 730 [1] http://paulbourke.net/geometry/polygonmesh/#insidepoly 731 732 """ 733 p = Point(p, dim=2) 734 if p in self.vertices or any(p in s for s in self.sides): 735 return False 736 737 # move to p, checking that the result is numeric 738 lit = [] 739 for v in self.vertices: 740 lit.append(v - p) # the difference is simplified 741 if lit[-1].free_symbols: 742 return None 743 744 poly = Polygon(*lit) 745 746 # polygon closure is assumed in the following test but Polygon removes duplicate pts so 747 # the last point has to be added so all sides are computed. Using Polygon.sides is 748 # not good since Segments are unordered. 749 args = poly.args 750 indices = list(range(-len(args), 1)) 751 752 if poly.is_convex(): 753 orientation = None 754 for i in indices: 755 a = args[i] 756 b = args[i + 1] 757 test = ((-a.y)*(b.x - a.x) - (-a.x)*(b.y - a.y)).is_negative 758 if orientation is None: 759 orientation = test 760 elif test is not orientation: 761 return False 762 return True 763 764 hit_odd = False 765 p1x, p1y = args[0].args 766 for i in indices[1:]: 767 p2x, p2y = args[i].args 768 if 0 > min(p1y, p2y): 769 if 0 <= max(p1y, p2y): 770 if 0 <= max(p1x, p2x): 771 if p1y != p2y: 772 xinters = (-p1y)*(p2x - p1x)/(p2y - p1y) + p1x 773 if p1x == p2x or 0 <= xinters: 774 hit_odd = not hit_odd 775 p1x, p1y = p2x, p2y 776 return hit_odd 777 778 def arbitrary_point(self, parameter='t'): 779 """A parameterized point on the polygon. 780 781 The parameter, varying from 0 to 1, assigns points to the position on 782 the perimeter that is that fraction of the total perimeter. So the 783 point evaluated at t=1/2 would return the point from the first vertex 784 that is 1/2 way around the polygon. 785 786 Parameters 787 ========== 788 789 parameter : str, optional 790 Default value is 't'. 791 792 Returns 793 ======= 794 795 arbitrary_point : Point 796 797 Raises 798 ====== 799 800 ValueError 801 When `parameter` already appears in the Polygon's definition. 802 803 See Also 804 ======== 805 806 sympy.geometry.point.Point 807 808 Examples 809 ======== 810 811 >>> from sympy import Polygon, Symbol 812 >>> t = Symbol('t', real=True) 813 >>> tri = Polygon((0, 0), (1, 0), (1, 1)) 814 >>> p = tri.arbitrary_point('t') 815 >>> perimeter = tri.perimeter 816 >>> s1, s2 = [s.length for s in tri.sides[:2]] 817 >>> p.subs(t, (s1 + s2/2)/perimeter) 818 Point2D(1, 1/2) 819 820 """ 821 t = _symbol(parameter, real=True) 822 if t.name in (f.name for f in self.free_symbols): 823 raise ValueError('Symbol %s already appears in object and cannot be used as a parameter.' % t.name) 824 sides = [] 825 perimeter = self.perimeter 826 perim_fraction_start = 0 827 for s in self.sides: 828 side_perim_fraction = s.length/perimeter 829 perim_fraction_end = perim_fraction_start + side_perim_fraction 830 pt = s.arbitrary_point(parameter).subs( 831 t, (t - perim_fraction_start)/side_perim_fraction) 832 sides.append( 833 (pt, (And(perim_fraction_start <= t, t < perim_fraction_end)))) 834 perim_fraction_start = perim_fraction_end 835 return Piecewise(*sides) 836 837 def parameter_value(self, other, t): 838 from sympy.solvers.solvers import solve 839 if not isinstance(other,GeometryEntity): 840 other = Point(other, dim=self.ambient_dimension) 841 if not isinstance(other,Point): 842 raise ValueError("other must be a point") 843 if other.free_symbols: 844 raise NotImplementedError('non-numeric coordinates') 845 unknown = False 846 T = Dummy('t', real=True) 847 p = self.arbitrary_point(T) 848 for pt, cond in p.args: 849 sol = solve(pt - other, T, dict=True) 850 if not sol: 851 continue 852 value = sol[0][T] 853 if simplify(cond.subs(T, value)) == True: 854 return {t: value} 855 unknown = True 856 if unknown: 857 raise ValueError("Given point may not be on %s" % func_name(self)) 858 raise ValueError("Given point is not on %s" % func_name(self)) 859 860 def plot_interval(self, parameter='t'): 861 """The plot interval for the default geometric plot of the polygon. 862 863 Parameters 864 ========== 865 866 parameter : str, optional 867 Default value is 't'. 868 869 Returns 870 ======= 871 872 plot_interval : list (plot interval) 873 [parameter, lower_bound, upper_bound] 874 875 Examples 876 ======== 877 878 >>> from sympy import Polygon 879 >>> p = Polygon((0, 0), (1, 0), (1, 1)) 880 >>> p.plot_interval() 881 [t, 0, 1] 882 883 """ 884 t = Symbol(parameter, real=True) 885 return [t, 0, 1] 886 887 def intersection(self, o): 888 """The intersection of polygon and geometry entity. 889 890 The intersection may be empty and can contain individual Points and 891 complete Line Segments. 892 893 Parameters 894 ========== 895 896 other: GeometryEntity 897 898 Returns 899 ======= 900 901 intersection : list 902 The list of Segments and Points 903 904 See Also 905 ======== 906 907 sympy.geometry.point.Point, sympy.geometry.line.Segment 908 909 Examples 910 ======== 911 912 >>> from sympy import Point, Polygon, Line 913 >>> p1, p2, p3, p4 = map(Point, [(0, 0), (1, 0), (5, 1), (0, 1)]) 914 >>> poly1 = Polygon(p1, p2, p3, p4) 915 >>> p5, p6, p7 = map(Point, [(3, 2), (1, -1), (0, 2)]) 916 >>> poly2 = Polygon(p5, p6, p7) 917 >>> poly1.intersection(poly2) 918 [Point2D(1/3, 1), Point2D(2/3, 0), Point2D(9/5, 1/5), Point2D(7/3, 1)] 919 >>> poly1.intersection(Line(p1, p2)) 920 [Segment2D(Point2D(0, 0), Point2D(1, 0))] 921 >>> poly1.intersection(p1) 922 [Point2D(0, 0)] 923 """ 924 intersection_result = [] 925 k = o.sides if isinstance(o, Polygon) else [o] 926 for side in self.sides: 927 for side1 in k: 928 intersection_result.extend(side.intersection(side1)) 929 930 intersection_result = list(uniq(intersection_result)) 931 points = [entity for entity in intersection_result if isinstance(entity, Point)] 932 segments = [entity for entity in intersection_result if isinstance(entity, Segment)] 933 934 if points and segments: 935 points_in_segments = list(uniq([point for point in points for segment in segments if point in segment])) 936 if points_in_segments: 937 for i in points_in_segments: 938 points.remove(i) 939 return list(ordered(segments + points)) 940 else: 941 return list(ordered(intersection_result)) 942 943 944 def cut_section(self, line): 945 """ 946 Returns a tuple of two polygon segments that lie above and below 947 the intersecting line respectively. 948 949 Parameters 950 ========== 951 952 line: Line object of geometry module 953 line which cuts the Polygon. The part of the Polygon that lies 954 above and below this line is returned. 955 956 Returns 957 ======= 958 959 upper_polygon, lower_polygon: Polygon objects or None 960 upper_polygon is the polygon that lies above the given line. 961 lower_polygon is the polygon that lies below the given line. 962 upper_polygon and lower polygon are ``None`` when no polygon 963 exists above the line or below the line. 964 965 Raises 966 ====== 967 968 ValueError: When the line does not intersect the polygon 969 970 References 971 ========== 972 973 https://github.com/sympy/sympy/wiki/A-method-to-return-a-cut-section-of-any-polygon-geometry 974 975 Examples 976 ======== 977 978 >>> from sympy import Polygon, Line 979 >>> a, b = 20, 10 980 >>> p1, p2, p3, p4 = [(0, b), (0, 0), (a, 0), (a, b)] 981 >>> rectangle = Polygon(p1, p2, p3, p4) 982 >>> t = rectangle.cut_section(Line((0, 5), slope=0)) 983 >>> t 984 (Polygon(Point2D(0, 10), Point2D(0, 5), Point2D(20, 5), Point2D(20, 10)), 985 Polygon(Point2D(0, 5), Point2D(0, 0), Point2D(20, 0), Point2D(20, 5))) 986 >>> upper_segment, lower_segment = t 987 >>> upper_segment.area 988 100 989 >>> upper_segment.centroid 990 Point2D(10, 15/2) 991 >>> lower_segment.centroid 992 Point2D(10, 5/2) 993 """ 994 intersection_points = self.intersection(line) 995 if not intersection_points: 996 raise ValueError("This line does not intersect the polygon") 997 998 points = list(self.vertices) 999 points.append(points[0]) 1000 1001 x, y = symbols('x, y', real=True, cls=Dummy) 1002 eq = line.equation(x, y) 1003 1004 # considering equation of line to be `ax +by + c` 1005 a = eq.coeff(x) 1006 b = eq.coeff(y) 1007 1008 upper_vertices = [] 1009 lower_vertices = [] 1010 # prev is true when previous point is above the line 1011 prev = True 1012 prev_point = None 1013 for point in points: 1014 # when coefficient of y is 0, right side of the line is 1015 # considered 1016 compare = eq.subs({x: point.x, y: point.y})/b if b \ 1017 else eq.subs(x, point.x)/a 1018 1019 # if point lies above line 1020 if compare > 0: 1021 if not prev: 1022 # if previous point lies below the line, the intersection 1023 # point of the polygon egde and the line has to be included 1024 edge = Line(point, prev_point) 1025 new_point = edge.intersection(line) 1026 upper_vertices.append(new_point[0]) 1027 lower_vertices.append(new_point[0]) 1028 1029 upper_vertices.append(point) 1030 prev = True 1031 else: 1032 if prev and prev_point: 1033 edge = Line(point, prev_point) 1034 new_point = edge.intersection(line) 1035 upper_vertices.append(new_point[0]) 1036 lower_vertices.append(new_point[0]) 1037 lower_vertices.append(point) 1038 prev = False 1039 prev_point = point 1040 1041 upper_polygon, lower_polygon = None, None 1042 if upper_vertices and isinstance(Polygon(*upper_vertices), Polygon): 1043 upper_polygon = Polygon(*upper_vertices) 1044 if lower_vertices and isinstance(Polygon(*lower_vertices), Polygon): 1045 lower_polygon = Polygon(*lower_vertices) 1046 1047 return upper_polygon, lower_polygon 1048 1049 1050 def distance(self, o): 1051 """ 1052 Returns the shortest distance between self and o. 1053 1054 If o is a point, then self does not need to be convex. 1055 If o is another polygon self and o must be convex. 1056 1057 Examples 1058 ======== 1059 1060 >>> from sympy import Point, Polygon, RegularPolygon 1061 >>> p1, p2 = map(Point, [(0, 0), (7, 5)]) 1062 >>> poly = Polygon(*RegularPolygon(p1, 1, 3).vertices) 1063 >>> poly.distance(p2) 1064 sqrt(61) 1065 """ 1066 if isinstance(o, Point): 1067 dist = oo 1068 for side in self.sides: 1069 current = side.distance(o) 1070 if current == 0: 1071 return S.Zero 1072 elif current < dist: 1073 dist = current 1074 return dist 1075 elif isinstance(o, Polygon) and self.is_convex() and o.is_convex(): 1076 return self._do_poly_distance(o) 1077 raise NotImplementedError() 1078 1079 def _do_poly_distance(self, e2): 1080 """ 1081 Calculates the least distance between the exteriors of two 1082 convex polygons e1 and e2. Does not check for the convexity 1083 of the polygons as this is checked by Polygon.distance. 1084 1085 Notes 1086 ===== 1087 1088 - Prints a warning if the two polygons possibly intersect as the return 1089 value will not be valid in such a case. For a more through test of 1090 intersection use intersection(). 1091 1092 See Also 1093 ======== 1094 1095 sympy.geometry.point.Point.distance 1096 1097 Examples 1098 ======== 1099 1100 >>> from sympy.geometry import Point, Polygon 1101 >>> square = Polygon(Point(0, 0), Point(0, 1), Point(1, 1), Point(1, 0)) 1102 >>> triangle = Polygon(Point(1, 2), Point(2, 2), Point(2, 1)) 1103 >>> square._do_poly_distance(triangle) 1104 sqrt(2)/2 1105 1106 Description of method used 1107 ========================== 1108 1109 Method: 1110 [1] http://cgm.cs.mcgill.ca/~orm/mind2p.html 1111 Uses rotating calipers: 1112 [2] https://en.wikipedia.org/wiki/Rotating_calipers 1113 and antipodal points: 1114 [3] https://en.wikipedia.org/wiki/Antipodal_point 1115 """ 1116 e1 = self 1117 1118 '''Tests for a possible intersection between the polygons and outputs a warning''' 1119 e1_center = e1.centroid 1120 e2_center = e2.centroid 1121 e1_max_radius = S.Zero 1122 e2_max_radius = S.Zero 1123 for vertex in e1.vertices: 1124 r = Point.distance(e1_center, vertex) 1125 if e1_max_radius < r: 1126 e1_max_radius = r 1127 for vertex in e2.vertices: 1128 r = Point.distance(e2_center, vertex) 1129 if e2_max_radius < r: 1130 e2_max_radius = r 1131 center_dist = Point.distance(e1_center, e2_center) 1132 if center_dist <= e1_max_radius + e2_max_radius: 1133 warnings.warn("Polygons may intersect producing erroneous output") 1134 1135 ''' 1136 Find the upper rightmost vertex of e1 and the lowest leftmost vertex of e2 1137 ''' 1138 e1_ymax = Point(0, -oo) 1139 e2_ymin = Point(0, oo) 1140 1141 for vertex in e1.vertices: 1142 if vertex.y > e1_ymax.y or (vertex.y == e1_ymax.y and vertex.x > e1_ymax.x): 1143 e1_ymax = vertex 1144 for vertex in e2.vertices: 1145 if vertex.y < e2_ymin.y or (vertex.y == e2_ymin.y and vertex.x < e2_ymin.x): 1146 e2_ymin = vertex 1147 min_dist = Point.distance(e1_ymax, e2_ymin) 1148 1149 ''' 1150 Produce a dictionary with vertices of e1 as the keys and, for each vertex, the points 1151 to which the vertex is connected as its value. The same is then done for e2. 1152 ''' 1153 e1_connections = {} 1154 e2_connections = {} 1155 1156 for side in e1.sides: 1157 if side.p1 in e1_connections: 1158 e1_connections[side.p1].append(side.p2) 1159 else: 1160 e1_connections[side.p1] = [side.p2] 1161 1162 if side.p2 in e1_connections: 1163 e1_connections[side.p2].append(side.p1) 1164 else: 1165 e1_connections[side.p2] = [side.p1] 1166 1167 for side in e2.sides: 1168 if side.p1 in e2_connections: 1169 e2_connections[side.p1].append(side.p2) 1170 else: 1171 e2_connections[side.p1] = [side.p2] 1172 1173 if side.p2 in e2_connections: 1174 e2_connections[side.p2].append(side.p1) 1175 else: 1176 e2_connections[side.p2] = [side.p1] 1177 1178 e1_current = e1_ymax 1179 e2_current = e2_ymin 1180 support_line = Line(Point(S.Zero, S.Zero), Point(S.One, S.Zero)) 1181 1182 ''' 1183 Determine which point in e1 and e2 will be selected after e2_ymin and e1_ymax, 1184 this information combined with the above produced dictionaries determines the 1185 path that will be taken around the polygons 1186 ''' 1187 point1 = e1_connections[e1_ymax][0] 1188 point2 = e1_connections[e1_ymax][1] 1189 angle1 = support_line.angle_between(Line(e1_ymax, point1)) 1190 angle2 = support_line.angle_between(Line(e1_ymax, point2)) 1191 if angle1 < angle2: 1192 e1_next = point1 1193 elif angle2 < angle1: 1194 e1_next = point2 1195 elif Point.distance(e1_ymax, point1) > Point.distance(e1_ymax, point2): 1196 e1_next = point2 1197 else: 1198 e1_next = point1 1199 1200 point1 = e2_connections[e2_ymin][0] 1201 point2 = e2_connections[e2_ymin][1] 1202 angle1 = support_line.angle_between(Line(e2_ymin, point1)) 1203 angle2 = support_line.angle_between(Line(e2_ymin, point2)) 1204 if angle1 > angle2: 1205 e2_next = point1 1206 elif angle2 > angle1: 1207 e2_next = point2 1208 elif Point.distance(e2_ymin, point1) > Point.distance(e2_ymin, point2): 1209 e2_next = point2 1210 else: 1211 e2_next = point1 1212 1213 ''' 1214 Loop which determines the distance between anti-podal pairs and updates the 1215 minimum distance accordingly. It repeats until it reaches the starting position. 1216 ''' 1217 while True: 1218 e1_angle = support_line.angle_between(Line(e1_current, e1_next)) 1219 e2_angle = pi - support_line.angle_between(Line( 1220 e2_current, e2_next)) 1221 1222 if (e1_angle < e2_angle) is True: 1223 support_line = Line(e1_current, e1_next) 1224 e1_segment = Segment(e1_current, e1_next) 1225 min_dist_current = e1_segment.distance(e2_current) 1226 1227 if min_dist_current.evalf() < min_dist.evalf(): 1228 min_dist = min_dist_current 1229 1230 if e1_connections[e1_next][0] != e1_current: 1231 e1_current = e1_next 1232 e1_next = e1_connections[e1_next][0] 1233 else: 1234 e1_current = e1_next 1235 e1_next = e1_connections[e1_next][1] 1236 elif (e1_angle > e2_angle) is True: 1237 support_line = Line(e2_next, e2_current) 1238 e2_segment = Segment(e2_current, e2_next) 1239 min_dist_current = e2_segment.distance(e1_current) 1240 1241 if min_dist_current.evalf() < min_dist.evalf(): 1242 min_dist = min_dist_current 1243 1244 if e2_connections[e2_next][0] != e2_current: 1245 e2_current = e2_next 1246 e2_next = e2_connections[e2_next][0] 1247 else: 1248 e2_current = e2_next 1249 e2_next = e2_connections[e2_next][1] 1250 else: 1251 support_line = Line(e1_current, e1_next) 1252 e1_segment = Segment(e1_current, e1_next) 1253 e2_segment = Segment(e2_current, e2_next) 1254 min1 = e1_segment.distance(e2_next) 1255 min2 = e2_segment.distance(e1_next) 1256 1257 min_dist_current = min(min1, min2) 1258 if min_dist_current.evalf() < min_dist.evalf(): 1259 min_dist = min_dist_current 1260 1261 if e1_connections[e1_next][0] != e1_current: 1262 e1_current = e1_next 1263 e1_next = e1_connections[e1_next][0] 1264 else: 1265 e1_current = e1_next 1266 e1_next = e1_connections[e1_next][1] 1267 1268 if e2_connections[e2_next][0] != e2_current: 1269 e2_current = e2_next 1270 e2_next = e2_connections[e2_next][0] 1271 else: 1272 e2_current = e2_next 1273 e2_next = e2_connections[e2_next][1] 1274 if e1_current == e1_ymax and e2_current == e2_ymin: 1275 break 1276 return min_dist 1277 1278 def _svg(self, scale_factor=1., fill_color="#66cc99"): 1279 """Returns SVG path element for the Polygon. 1280 1281 Parameters 1282 ========== 1283 1284 scale_factor : float 1285 Multiplication factor for the SVG stroke-width. Default is 1. 1286 fill_color : str, optional 1287 Hex string for fill color. Default is "#66cc99". 1288 """ 1289 1290 from sympy.core.evalf import N 1291 1292 verts = map(N, self.vertices) 1293 coords = ["{},{}".format(p.x, p.y) for p in verts] 1294 path = "M {} L {} z".format(coords[0], " L ".join(coords[1:])) 1295 return ( 1296 '<path fill-rule="evenodd" fill="{2}" stroke="#555555" ' 1297 'stroke-width="{0}" opacity="0.6" d="{1}" />' 1298 ).format(2. * scale_factor, path, fill_color) 1299 1300 def _hashable_content(self): 1301 1302 D = {} 1303 def ref_list(point_list): 1304 kee = {} 1305 for i, p in enumerate(ordered(set(point_list))): 1306 kee[p] = i 1307 D[i] = p 1308 return [kee[p] for p in point_list] 1309 1310 S1 = ref_list(self.args) 1311 r_nor = rotate_left(S1, least_rotation(S1)) 1312 S2 = ref_list(list(reversed(self.args))) 1313 r_rev = rotate_left(S2, least_rotation(S2)) 1314 if r_nor < r_rev: 1315 r = r_nor 1316 else: 1317 r = r_rev 1318 canonical_args = [ D[order] for order in r ] 1319 return tuple(canonical_args) 1320 1321 def __contains__(self, o): 1322 """ 1323 Return True if o is contained within the boundary lines of self.altitudes 1324 1325 Parameters 1326 ========== 1327 1328 other : GeometryEntity 1329 1330 Returns 1331 ======= 1332 1333 contained in : bool 1334 The points (and sides, if applicable) are contained in self. 1335 1336 See Also 1337 ======== 1338 1339 sympy.geometry.entity.GeometryEntity.encloses 1340 1341 Examples 1342 ======== 1343 1344 >>> from sympy import Line, Segment, Point 1345 >>> p = Point(0, 0) 1346 >>> q = Point(1, 1) 1347 >>> s = Segment(p, q*2) 1348 >>> l = Line(p, q) 1349 >>> p in q 1350 False 1351 >>> p in s 1352 True 1353 >>> q*3 in s 1354 False 1355 >>> s in l 1356 True 1357 1358 """ 1359 1360 if isinstance(o, Polygon): 1361 return self == o 1362 elif isinstance(o, Segment): 1363 return any(o in s for s in self.sides) 1364 elif isinstance(o, Point): 1365 if o in self.vertices: 1366 return True 1367 for side in self.sides: 1368 if o in side: 1369 return True 1370 1371 return False 1372 1373 def bisectors(p, prec=None): 1374 """Returns angle bisectors of a polygon. If prec is given 1375 then approximate the point defining the ray to that precision. 1376 1377 The distance between the points defining the bisector ray is 1. 1378 1379 Examples 1380 ======== 1381 1382 >>> from sympy import Polygon, Point 1383 >>> p = Polygon(Point(0, 0), Point(2, 0), Point(1, 1), Point(0, 3)) 1384 >>> p.bisectors(2) 1385 {Point2D(0, 0): Ray2D(Point2D(0, 0), Point2D(0.71, 0.71)), 1386 Point2D(0, 3): Ray2D(Point2D(0, 3), Point2D(0.23, 2.0)), 1387 Point2D(1, 1): Ray2D(Point2D(1, 1), Point2D(0.19, 0.42)), 1388 Point2D(2, 0): Ray2D(Point2D(2, 0), Point2D(1.1, 0.38))} 1389 """ 1390 b = {} 1391 pts = list(p.args) 1392 pts.append(pts[0]) # close it 1393 cw = Polygon._isright(*pts[:3]) 1394 if cw: 1395 pts = list(reversed(pts)) 1396 for v, a in p.angles.items(): 1397 i = pts.index(v) 1398 p1, p2 = Point._normalize_dimension(pts[i], pts[i + 1]) 1399 ray = Ray(p1, p2).rotate(a/2, v) 1400 dir = ray.direction 1401 ray = Ray(ray.p1, ray.p1 + dir/dir.distance((0, 0))) 1402 if prec is not None: 1403 ray = Ray(ray.p1, ray.p2.n(prec)) 1404 b[v] = ray 1405 return b 1406 1407 1408class RegularPolygon(Polygon): 1409 """ 1410 A regular polygon. 1411 1412 Such a polygon has all internal angles equal and all sides the same length. 1413 1414 Parameters 1415 ========== 1416 1417 center : Point 1418 radius : number or Basic instance 1419 The distance from the center to a vertex 1420 n : int 1421 The number of sides 1422 1423 Attributes 1424 ========== 1425 1426 vertices 1427 center 1428 radius 1429 rotation 1430 apothem 1431 interior_angle 1432 exterior_angle 1433 circumcircle 1434 incircle 1435 angles 1436 1437 Raises 1438 ====== 1439 1440 GeometryError 1441 If the `center` is not a Point, or the `radius` is not a number or Basic 1442 instance, or the number of sides, `n`, is less than three. 1443 1444 Notes 1445 ===== 1446 1447 A RegularPolygon can be instantiated with Polygon with the kwarg n. 1448 1449 Regular polygons are instantiated with a center, radius, number of sides 1450 and a rotation angle. Whereas the arguments of a Polygon are vertices, the 1451 vertices of the RegularPolygon must be obtained with the vertices method. 1452 1453 See Also 1454 ======== 1455 1456 sympy.geometry.point.Point, Polygon 1457 1458 Examples 1459 ======== 1460 1461 >>> from sympy.geometry import RegularPolygon, Point 1462 >>> r = RegularPolygon(Point(0, 0), 5, 3) 1463 >>> r 1464 RegularPolygon(Point2D(0, 0), 5, 3, 0) 1465 >>> r.vertices[0] 1466 Point2D(5, 0) 1467 1468 """ 1469 1470 __slots__ = ('_n', '_center', '_radius', '_rot') 1471 1472 def __new__(self, c, r, n, rot=0, **kwargs): 1473 r, n, rot = map(sympify, (r, n, rot)) 1474 c = Point(c, dim=2, **kwargs) 1475 if not isinstance(r, Expr): 1476 raise GeometryError("r must be an Expr object, not %s" % r) 1477 if n.is_Number: 1478 as_int(n) # let an error raise if necessary 1479 if n < 3: 1480 raise GeometryError("n must be a >= 3, not %s" % n) 1481 1482 obj = GeometryEntity.__new__(self, c, r, n, **kwargs) 1483 obj._n = n 1484 obj._center = c 1485 obj._radius = r 1486 obj._rot = rot % (2*S.Pi/n) if rot.is_number else rot 1487 return obj 1488 1489 @property 1490 def args(self): 1491 """ 1492 Returns the center point, the radius, 1493 the number of sides, and the orientation angle. 1494 1495 Examples 1496 ======== 1497 1498 >>> from sympy import RegularPolygon, Point 1499 >>> r = RegularPolygon(Point(0, 0), 5, 3) 1500 >>> r.args 1501 (Point2D(0, 0), 5, 3, 0) 1502 """ 1503 return self._center, self._radius, self._n, self._rot 1504 1505 def __str__(self): 1506 return 'RegularPolygon(%s, %s, %s, %s)' % tuple(self.args) 1507 1508 def __repr__(self): 1509 return 'RegularPolygon(%s, %s, %s, %s)' % tuple(self.args) 1510 1511 @property 1512 def area(self): 1513 """Returns the area. 1514 1515 Examples 1516 ======== 1517 1518 >>> from sympy.geometry import RegularPolygon 1519 >>> square = RegularPolygon((0, 0), 1, 4) 1520 >>> square.area 1521 2 1522 >>> _ == square.length**2 1523 True 1524 """ 1525 c, r, n, rot = self.args 1526 return sign(r)*n*self.length**2/(4*tan(pi/n)) 1527 1528 @property 1529 def length(self): 1530 """Returns the length of the sides. 1531 1532 The half-length of the side and the apothem form two legs 1533 of a right triangle whose hypotenuse is the radius of the 1534 regular polygon. 1535 1536 Examples 1537 ======== 1538 1539 >>> from sympy.geometry import RegularPolygon 1540 >>> from sympy import sqrt 1541 >>> s = square_in_unit_circle = RegularPolygon((0, 0), 1, 4) 1542 >>> s.length 1543 sqrt(2) 1544 >>> sqrt((_/2)**2 + s.apothem**2) == s.radius 1545 True 1546 1547 """ 1548 return self.radius*2*sin(pi/self._n) 1549 1550 @property 1551 def center(self): 1552 """The center of the RegularPolygon 1553 1554 This is also the center of the circumscribing circle. 1555 1556 Returns 1557 ======= 1558 1559 center : Point 1560 1561 See Also 1562 ======== 1563 1564 sympy.geometry.point.Point, sympy.geometry.ellipse.Ellipse.center 1565 1566 Examples 1567 ======== 1568 1569 >>> from sympy.geometry import RegularPolygon, Point 1570 >>> rp = RegularPolygon(Point(0, 0), 5, 4) 1571 >>> rp.center 1572 Point2D(0, 0) 1573 """ 1574 return self._center 1575 1576 centroid = center 1577 1578 @property 1579 def circumcenter(self): 1580 """ 1581 Alias for center. 1582 1583 Examples 1584 ======== 1585 1586 >>> from sympy.geometry import RegularPolygon, Point 1587 >>> rp = RegularPolygon(Point(0, 0), 5, 4) 1588 >>> rp.circumcenter 1589 Point2D(0, 0) 1590 """ 1591 return self.center 1592 1593 @property 1594 def radius(self): 1595 """Radius of the RegularPolygon 1596 1597 This is also the radius of the circumscribing circle. 1598 1599 Returns 1600 ======= 1601 1602 radius : number or instance of Basic 1603 1604 See Also 1605 ======== 1606 1607 sympy.geometry.line.Segment.length, sympy.geometry.ellipse.Circle.radius 1608 1609 Examples 1610 ======== 1611 1612 >>> from sympy import Symbol 1613 >>> from sympy.geometry import RegularPolygon, Point 1614 >>> radius = Symbol('r') 1615 >>> rp = RegularPolygon(Point(0, 0), radius, 4) 1616 >>> rp.radius 1617 r 1618 1619 """ 1620 return self._radius 1621 1622 @property 1623 def circumradius(self): 1624 """ 1625 Alias for radius. 1626 1627 Examples 1628 ======== 1629 1630 >>> from sympy import Symbol 1631 >>> from sympy.geometry import RegularPolygon, Point 1632 >>> radius = Symbol('r') 1633 >>> rp = RegularPolygon(Point(0, 0), radius, 4) 1634 >>> rp.circumradius 1635 r 1636 """ 1637 return self.radius 1638 1639 @property 1640 def rotation(self): 1641 """CCW angle by which the RegularPolygon is rotated 1642 1643 Returns 1644 ======= 1645 1646 rotation : number or instance of Basic 1647 1648 Examples 1649 ======== 1650 1651 >>> from sympy import pi 1652 >>> from sympy.abc import a 1653 >>> from sympy.geometry import RegularPolygon, Point 1654 >>> RegularPolygon(Point(0, 0), 3, 4, pi/4).rotation 1655 pi/4 1656 1657 Numerical rotation angles are made canonical: 1658 1659 >>> RegularPolygon(Point(0, 0), 3, 4, a).rotation 1660 a 1661 >>> RegularPolygon(Point(0, 0), 3, 4, pi).rotation 1662 0 1663 1664 """ 1665 return self._rot 1666 1667 @property 1668 def apothem(self): 1669 """The inradius of the RegularPolygon. 1670 1671 The apothem/inradius is the radius of the inscribed circle. 1672 1673 Returns 1674 ======= 1675 1676 apothem : number or instance of Basic 1677 1678 See Also 1679 ======== 1680 1681 sympy.geometry.line.Segment.length, sympy.geometry.ellipse.Circle.radius 1682 1683 Examples 1684 ======== 1685 1686 >>> from sympy import Symbol 1687 >>> from sympy.geometry import RegularPolygon, Point 1688 >>> radius = Symbol('r') 1689 >>> rp = RegularPolygon(Point(0, 0), radius, 4) 1690 >>> rp.apothem 1691 sqrt(2)*r/2 1692 1693 """ 1694 return self.radius * cos(S.Pi/self._n) 1695 1696 @property 1697 def inradius(self): 1698 """ 1699 Alias for apothem. 1700 1701 Examples 1702 ======== 1703 1704 >>> from sympy import Symbol 1705 >>> from sympy.geometry import RegularPolygon, Point 1706 >>> radius = Symbol('r') 1707 >>> rp = RegularPolygon(Point(0, 0), radius, 4) 1708 >>> rp.inradius 1709 sqrt(2)*r/2 1710 """ 1711 return self.apothem 1712 1713 @property 1714 def interior_angle(self): 1715 """Measure of the interior angles. 1716 1717 Returns 1718 ======= 1719 1720 interior_angle : number 1721 1722 See Also 1723 ======== 1724 1725 sympy.geometry.line.LinearEntity.angle_between 1726 1727 Examples 1728 ======== 1729 1730 >>> from sympy.geometry import RegularPolygon, Point 1731 >>> rp = RegularPolygon(Point(0, 0), 4, 8) 1732 >>> rp.interior_angle 1733 3*pi/4 1734 1735 """ 1736 return (self._n - 2)*S.Pi/self._n 1737 1738 @property 1739 def exterior_angle(self): 1740 """Measure of the exterior angles. 1741 1742 Returns 1743 ======= 1744 1745 exterior_angle : number 1746 1747 See Also 1748 ======== 1749 1750 sympy.geometry.line.LinearEntity.angle_between 1751 1752 Examples 1753 ======== 1754 1755 >>> from sympy.geometry import RegularPolygon, Point 1756 >>> rp = RegularPolygon(Point(0, 0), 4, 8) 1757 >>> rp.exterior_angle 1758 pi/4 1759 1760 """ 1761 return 2*S.Pi/self._n 1762 1763 @property 1764 def circumcircle(self): 1765 """The circumcircle of the RegularPolygon. 1766 1767 Returns 1768 ======= 1769 1770 circumcircle : Circle 1771 1772 See Also 1773 ======== 1774 1775 circumcenter, sympy.geometry.ellipse.Circle 1776 1777 Examples 1778 ======== 1779 1780 >>> from sympy.geometry import RegularPolygon, Point 1781 >>> rp = RegularPolygon(Point(0, 0), 4, 8) 1782 >>> rp.circumcircle 1783 Circle(Point2D(0, 0), 4) 1784 1785 """ 1786 return Circle(self.center, self.radius) 1787 1788 @property 1789 def incircle(self): 1790 """The incircle of the RegularPolygon. 1791 1792 Returns 1793 ======= 1794 1795 incircle : Circle 1796 1797 See Also 1798 ======== 1799 1800 inradius, sympy.geometry.ellipse.Circle 1801 1802 Examples 1803 ======== 1804 1805 >>> from sympy.geometry import RegularPolygon, Point 1806 >>> rp = RegularPolygon(Point(0, 0), 4, 7) 1807 >>> rp.incircle 1808 Circle(Point2D(0, 0), 4*cos(pi/7)) 1809 1810 """ 1811 return Circle(self.center, self.apothem) 1812 1813 @property 1814 def angles(self): 1815 """ 1816 Returns a dictionary with keys, the vertices of the Polygon, 1817 and values, the interior angle at each vertex. 1818 1819 Examples 1820 ======== 1821 1822 >>> from sympy import RegularPolygon, Point 1823 >>> r = RegularPolygon(Point(0, 0), 5, 3) 1824 >>> r.angles 1825 {Point2D(-5/2, -5*sqrt(3)/2): pi/3, 1826 Point2D(-5/2, 5*sqrt(3)/2): pi/3, 1827 Point2D(5, 0): pi/3} 1828 """ 1829 ret = {} 1830 ang = self.interior_angle 1831 for v in self.vertices: 1832 ret[v] = ang 1833 return ret 1834 1835 def encloses_point(self, p): 1836 """ 1837 Return True if p is enclosed by (is inside of) self. 1838 1839 Notes 1840 ===== 1841 1842 Being on the border of self is considered False. 1843 1844 The general Polygon.encloses_point method is called only if 1845 a point is not within or beyond the incircle or circumcircle, 1846 respectively. 1847 1848 Parameters 1849 ========== 1850 1851 p : Point 1852 1853 Returns 1854 ======= 1855 1856 encloses_point : True, False or None 1857 1858 See Also 1859 ======== 1860 1861 sympy.geometry.ellipse.Ellipse.encloses_point 1862 1863 Examples 1864 ======== 1865 1866 >>> from sympy import RegularPolygon, S, Point, Symbol 1867 >>> p = RegularPolygon((0, 0), 3, 4) 1868 >>> p.encloses_point(Point(0, 0)) 1869 True 1870 >>> r, R = p.inradius, p.circumradius 1871 >>> p.encloses_point(Point((r + R)/2, 0)) 1872 True 1873 >>> p.encloses_point(Point(R/2, R/2 + (R - r)/10)) 1874 False 1875 >>> t = Symbol('t', real=True) 1876 >>> p.encloses_point(p.arbitrary_point().subs(t, S.Half)) 1877 False 1878 >>> p.encloses_point(Point(5, 5)) 1879 False 1880 1881 """ 1882 1883 c = self.center 1884 d = Segment(c, p).length 1885 if d >= self.radius: 1886 return False 1887 elif d < self.inradius: 1888 return True 1889 else: 1890 # now enumerate the RegularPolygon like a general polygon. 1891 return Polygon.encloses_point(self, p) 1892 1893 def spin(self, angle): 1894 """Increment *in place* the virtual Polygon's rotation by ccw angle. 1895 1896 See also: rotate method which moves the center. 1897 1898 >>> from sympy import Polygon, Point, pi 1899 >>> r = Polygon(Point(0,0), 1, n=3) 1900 >>> r.vertices[0] 1901 Point2D(1, 0) 1902 >>> r.spin(pi/6) 1903 >>> r.vertices[0] 1904 Point2D(sqrt(3)/2, 1/2) 1905 1906 See Also 1907 ======== 1908 1909 rotation 1910 rotate : Creates a copy of the RegularPolygon rotated about a Point 1911 1912 """ 1913 self._rot += angle 1914 1915 def rotate(self, angle, pt=None): 1916 """Override GeometryEntity.rotate to first rotate the RegularPolygon 1917 about its center. 1918 1919 >>> from sympy import Point, RegularPolygon, pi 1920 >>> t = RegularPolygon(Point(1, 0), 1, 3) 1921 >>> t.vertices[0] # vertex on x-axis 1922 Point2D(2, 0) 1923 >>> t.rotate(pi/2).vertices[0] # vertex on y axis now 1924 Point2D(0, 2) 1925 1926 See Also 1927 ======== 1928 1929 rotation 1930 spin : Rotates a RegularPolygon in place 1931 1932 """ 1933 1934 r = type(self)(*self.args) # need a copy or else changes are in-place 1935 r._rot += angle 1936 return GeometryEntity.rotate(r, angle, pt) 1937 1938 def scale(self, x=1, y=1, pt=None): 1939 """Override GeometryEntity.scale since it is the radius that must be 1940 scaled (if x == y) or else a new Polygon must be returned. 1941 1942 >>> from sympy import RegularPolygon 1943 1944 Symmetric scaling returns a RegularPolygon: 1945 1946 >>> RegularPolygon((0, 0), 1, 4).scale(2, 2) 1947 RegularPolygon(Point2D(0, 0), 2, 4, 0) 1948 1949 Asymmetric scaling returns a kite as a Polygon: 1950 1951 >>> RegularPolygon((0, 0), 1, 4).scale(2, 1) 1952 Polygon(Point2D(2, 0), Point2D(0, 1), Point2D(-2, 0), Point2D(0, -1)) 1953 1954 """ 1955 if pt: 1956 pt = Point(pt, dim=2) 1957 return self.translate(*(-pt).args).scale(x, y).translate(*pt.args) 1958 if x != y: 1959 return Polygon(*self.vertices).scale(x, y) 1960 c, r, n, rot = self.args 1961 r *= x 1962 return self.func(c, r, n, rot) 1963 1964 def reflect(self, line): 1965 """Override GeometryEntity.reflect since this is not made of only 1966 points. 1967 1968 Examples 1969 ======== 1970 1971 >>> from sympy import RegularPolygon, Line 1972 1973 >>> RegularPolygon((0, 0), 1, 4).reflect(Line((0, 1), slope=-2)) 1974 RegularPolygon(Point2D(4/5, 2/5), -1, 4, atan(4/3)) 1975 1976 """ 1977 c, r, n, rot = self.args 1978 v = self.vertices[0] 1979 d = v - c 1980 cc = c.reflect(line) 1981 vv = v.reflect(line) 1982 dd = vv - cc 1983 # calculate rotation about the new center 1984 # which will align the vertices 1985 l1 = Ray((0, 0), dd) 1986 l2 = Ray((0, 0), d) 1987 ang = l1.closing_angle(l2) 1988 rot += ang 1989 # change sign of radius as point traversal is reversed 1990 return self.func(cc, -r, n, rot) 1991 1992 @property 1993 def vertices(self): 1994 """The vertices of the RegularPolygon. 1995 1996 Returns 1997 ======= 1998 1999 vertices : list 2000 Each vertex is a Point. 2001 2002 See Also 2003 ======== 2004 2005 sympy.geometry.point.Point 2006 2007 Examples 2008 ======== 2009 2010 >>> from sympy.geometry import RegularPolygon, Point 2011 >>> rp = RegularPolygon(Point(0, 0), 5, 4) 2012 >>> rp.vertices 2013 [Point2D(5, 0), Point2D(0, 5), Point2D(-5, 0), Point2D(0, -5)] 2014 2015 """ 2016 c = self._center 2017 r = abs(self._radius) 2018 rot = self._rot 2019 v = 2*S.Pi/self._n 2020 2021 return [Point(c.x + r*cos(k*v + rot), c.y + r*sin(k*v + rot)) 2022 for k in range(self._n)] 2023 2024 def __eq__(self, o): 2025 if not isinstance(o, Polygon): 2026 return False 2027 elif not isinstance(o, RegularPolygon): 2028 return Polygon.__eq__(o, self) 2029 return self.args == o.args 2030 2031 def __hash__(self): 2032 return super().__hash__() 2033 2034 2035class Triangle(Polygon): 2036 """ 2037 A polygon with three vertices and three sides. 2038 2039 Parameters 2040 ========== 2041 2042 points : sequence of Points 2043 keyword: asa, sas, or sss to specify sides/angles of the triangle 2044 2045 Attributes 2046 ========== 2047 2048 vertices 2049 altitudes 2050 orthocenter 2051 circumcenter 2052 circumradius 2053 circumcircle 2054 inradius 2055 incircle 2056 exradii 2057 medians 2058 medial 2059 nine_point_circle 2060 2061 Raises 2062 ====== 2063 2064 GeometryError 2065 If the number of vertices is not equal to three, or one of the vertices 2066 is not a Point, or a valid keyword is not given. 2067 2068 See Also 2069 ======== 2070 2071 sympy.geometry.point.Point, Polygon 2072 2073 Examples 2074 ======== 2075 2076 >>> from sympy.geometry import Triangle, Point 2077 >>> Triangle(Point(0, 0), Point(4, 0), Point(4, 3)) 2078 Triangle(Point2D(0, 0), Point2D(4, 0), Point2D(4, 3)) 2079 2080 Keywords sss, sas, or asa can be used to give the desired 2081 side lengths (in order) and interior angles (in degrees) that 2082 define the triangle: 2083 2084 >>> Triangle(sss=(3, 4, 5)) 2085 Triangle(Point2D(0, 0), Point2D(3, 0), Point2D(3, 4)) 2086 >>> Triangle(asa=(30, 1, 30)) 2087 Triangle(Point2D(0, 0), Point2D(1, 0), Point2D(1/2, sqrt(3)/6)) 2088 >>> Triangle(sas=(1, 45, 2)) 2089 Triangle(Point2D(0, 0), Point2D(2, 0), Point2D(sqrt(2)/2, sqrt(2)/2)) 2090 2091 """ 2092 2093 def __new__(cls, *args, **kwargs): 2094 if len(args) != 3: 2095 if 'sss' in kwargs: 2096 return _sss(*[simplify(a) for a in kwargs['sss']]) 2097 if 'asa' in kwargs: 2098 return _asa(*[simplify(a) for a in kwargs['asa']]) 2099 if 'sas' in kwargs: 2100 return _sas(*[simplify(a) for a in kwargs['sas']]) 2101 msg = "Triangle instantiates with three points or a valid keyword." 2102 raise GeometryError(msg) 2103 2104 vertices = [Point(a, dim=2, **kwargs) for a in args] 2105 2106 # remove consecutive duplicates 2107 nodup = [] 2108 for p in vertices: 2109 if nodup and p == nodup[-1]: 2110 continue 2111 nodup.append(p) 2112 if len(nodup) > 1 and nodup[-1] == nodup[0]: 2113 nodup.pop() # last point was same as first 2114 2115 # remove collinear points 2116 i = -3 2117 while i < len(nodup) - 3 and len(nodup) > 2: 2118 a, b, c = sorted( 2119 [nodup[i], nodup[i + 1], nodup[i + 2]], key=default_sort_key) 2120 if Point.is_collinear(a, b, c): 2121 nodup[i] = a 2122 nodup[i + 1] = None 2123 nodup.pop(i + 1) 2124 i += 1 2125 2126 vertices = list(filter(lambda x: x is not None, nodup)) 2127 2128 if len(vertices) == 3: 2129 return GeometryEntity.__new__(cls, *vertices, **kwargs) 2130 elif len(vertices) == 2: 2131 return Segment(*vertices, **kwargs) 2132 else: 2133 return Point(*vertices, **kwargs) 2134 2135 @property 2136 def vertices(self): 2137 """The triangle's vertices 2138 2139 Returns 2140 ======= 2141 2142 vertices : tuple 2143 Each element in the tuple is a Point 2144 2145 See Also 2146 ======== 2147 2148 sympy.geometry.point.Point 2149 2150 Examples 2151 ======== 2152 2153 >>> from sympy.geometry import Triangle, Point 2154 >>> t = Triangle(Point(0, 0), Point(4, 0), Point(4, 3)) 2155 >>> t.vertices 2156 (Point2D(0, 0), Point2D(4, 0), Point2D(4, 3)) 2157 2158 """ 2159 return self.args 2160 2161 def is_similar(t1, t2): 2162 """Is another triangle similar to this one. 2163 2164 Two triangles are similar if one can be uniformly scaled to the other. 2165 2166 Parameters 2167 ========== 2168 2169 other: Triangle 2170 2171 Returns 2172 ======= 2173 2174 is_similar : boolean 2175 2176 See Also 2177 ======== 2178 2179 sympy.geometry.entity.GeometryEntity.is_similar 2180 2181 Examples 2182 ======== 2183 2184 >>> from sympy.geometry import Triangle, Point 2185 >>> t1 = Triangle(Point(0, 0), Point(4, 0), Point(4, 3)) 2186 >>> t2 = Triangle(Point(0, 0), Point(-4, 0), Point(-4, -3)) 2187 >>> t1.is_similar(t2) 2188 True 2189 2190 >>> t2 = Triangle(Point(0, 0), Point(-4, 0), Point(-4, -4)) 2191 >>> t1.is_similar(t2) 2192 False 2193 2194 """ 2195 if not isinstance(t2, Polygon): 2196 return False 2197 2198 s1_1, s1_2, s1_3 = [side.length for side in t1.sides] 2199 s2 = [side.length for side in t2.sides] 2200 2201 def _are_similar(u1, u2, u3, v1, v2, v3): 2202 e1 = simplify(u1/v1) 2203 e2 = simplify(u2/v2) 2204 e3 = simplify(u3/v3) 2205 return bool(e1 == e2) and bool(e2 == e3) 2206 2207 # There's only 6 permutations, so write them out 2208 return _are_similar(s1_1, s1_2, s1_3, *s2) or \ 2209 _are_similar(s1_1, s1_3, s1_2, *s2) or \ 2210 _are_similar(s1_2, s1_1, s1_3, *s2) or \ 2211 _are_similar(s1_2, s1_3, s1_1, *s2) or \ 2212 _are_similar(s1_3, s1_1, s1_2, *s2) or \ 2213 _are_similar(s1_3, s1_2, s1_1, *s2) 2214 2215 def is_equilateral(self): 2216 """Are all the sides the same length? 2217 2218 Returns 2219 ======= 2220 2221 is_equilateral : boolean 2222 2223 See Also 2224 ======== 2225 2226 sympy.geometry.entity.GeometryEntity.is_similar, RegularPolygon 2227 is_isosceles, is_right, is_scalene 2228 2229 Examples 2230 ======== 2231 2232 >>> from sympy.geometry import Triangle, Point 2233 >>> t1 = Triangle(Point(0, 0), Point(4, 0), Point(4, 3)) 2234 >>> t1.is_equilateral() 2235 False 2236 2237 >>> from sympy import sqrt 2238 >>> t2 = Triangle(Point(0, 0), Point(10, 0), Point(5, 5*sqrt(3))) 2239 >>> t2.is_equilateral() 2240 True 2241 2242 """ 2243 return not has_variety(s.length for s in self.sides) 2244 2245 def is_isosceles(self): 2246 """Are two or more of the sides the same length? 2247 2248 Returns 2249 ======= 2250 2251 is_isosceles : boolean 2252 2253 See Also 2254 ======== 2255 2256 is_equilateral, is_right, is_scalene 2257 2258 Examples 2259 ======== 2260 2261 >>> from sympy.geometry import Triangle, Point 2262 >>> t1 = Triangle(Point(0, 0), Point(4, 0), Point(2, 4)) 2263 >>> t1.is_isosceles() 2264 True 2265 2266 """ 2267 return has_dups(s.length for s in self.sides) 2268 2269 def is_scalene(self): 2270 """Are all the sides of the triangle of different lengths? 2271 2272 Returns 2273 ======= 2274 2275 is_scalene : boolean 2276 2277 See Also 2278 ======== 2279 2280 is_equilateral, is_isosceles, is_right 2281 2282 Examples 2283 ======== 2284 2285 >>> from sympy.geometry import Triangle, Point 2286 >>> t1 = Triangle(Point(0, 0), Point(4, 0), Point(1, 4)) 2287 >>> t1.is_scalene() 2288 True 2289 2290 """ 2291 return not has_dups(s.length for s in self.sides) 2292 2293 def is_right(self): 2294 """Is the triangle right-angled. 2295 2296 Returns 2297 ======= 2298 2299 is_right : boolean 2300 2301 See Also 2302 ======== 2303 2304 sympy.geometry.line.LinearEntity.is_perpendicular 2305 is_equilateral, is_isosceles, is_scalene 2306 2307 Examples 2308 ======== 2309 2310 >>> from sympy.geometry import Triangle, Point 2311 >>> t1 = Triangle(Point(0, 0), Point(4, 0), Point(4, 3)) 2312 >>> t1.is_right() 2313 True 2314 2315 """ 2316 s = self.sides 2317 return Segment.is_perpendicular(s[0], s[1]) or \ 2318 Segment.is_perpendicular(s[1], s[2]) or \ 2319 Segment.is_perpendicular(s[0], s[2]) 2320 2321 @property 2322 def altitudes(self): 2323 """The altitudes of the triangle. 2324 2325 An altitude of a triangle is a segment through a vertex, 2326 perpendicular to the opposite side, with length being the 2327 height of the vertex measured from the line containing the side. 2328 2329 Returns 2330 ======= 2331 2332 altitudes : dict 2333 The dictionary consists of keys which are vertices and values 2334 which are Segments. 2335 2336 See Also 2337 ======== 2338 2339 sympy.geometry.point.Point, sympy.geometry.line.Segment.length 2340 2341 Examples 2342 ======== 2343 2344 >>> from sympy.geometry import Point, Triangle 2345 >>> p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, 1) 2346 >>> t = Triangle(p1, p2, p3) 2347 >>> t.altitudes[p1] 2348 Segment2D(Point2D(0, 0), Point2D(1/2, 1/2)) 2349 2350 """ 2351 s = self.sides 2352 v = self.vertices 2353 return {v[0]: s[1].perpendicular_segment(v[0]), 2354 v[1]: s[2].perpendicular_segment(v[1]), 2355 v[2]: s[0].perpendicular_segment(v[2])} 2356 2357 @property 2358 def orthocenter(self): 2359 """The orthocenter of the triangle. 2360 2361 The orthocenter is the intersection of the altitudes of a triangle. 2362 It may lie inside, outside or on the triangle. 2363 2364 Returns 2365 ======= 2366 2367 orthocenter : Point 2368 2369 See Also 2370 ======== 2371 2372 sympy.geometry.point.Point 2373 2374 Examples 2375 ======== 2376 2377 >>> from sympy.geometry import Point, Triangle 2378 >>> p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, 1) 2379 >>> t = Triangle(p1, p2, p3) 2380 >>> t.orthocenter 2381 Point2D(0, 0) 2382 2383 """ 2384 a = self.altitudes 2385 v = self.vertices 2386 return Line(a[v[0]]).intersection(Line(a[v[1]]))[0] 2387 2388 @property 2389 def circumcenter(self): 2390 """The circumcenter of the triangle 2391 2392 The circumcenter is the center of the circumcircle. 2393 2394 Returns 2395 ======= 2396 2397 circumcenter : Point 2398 2399 See Also 2400 ======== 2401 2402 sympy.geometry.point.Point 2403 2404 Examples 2405 ======== 2406 2407 >>> from sympy.geometry import Point, Triangle 2408 >>> p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, 1) 2409 >>> t = Triangle(p1, p2, p3) 2410 >>> t.circumcenter 2411 Point2D(1/2, 1/2) 2412 """ 2413 a, b, c = [x.perpendicular_bisector() for x in self.sides] 2414 if not a.intersection(b): 2415 print(a,b,a.intersection(b)) 2416 return a.intersection(b)[0] 2417 2418 @property 2419 def circumradius(self): 2420 """The radius of the circumcircle of the triangle. 2421 2422 Returns 2423 ======= 2424 2425 circumradius : number of Basic instance 2426 2427 See Also 2428 ======== 2429 2430 sympy.geometry.ellipse.Circle.radius 2431 2432 Examples 2433 ======== 2434 2435 >>> from sympy import Symbol 2436 >>> from sympy.geometry import Point, Triangle 2437 >>> a = Symbol('a') 2438 >>> p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, a) 2439 >>> t = Triangle(p1, p2, p3) 2440 >>> t.circumradius 2441 sqrt(a**2/4 + 1/4) 2442 """ 2443 return Point.distance(self.circumcenter, self.vertices[0]) 2444 2445 @property 2446 def circumcircle(self): 2447 """The circle which passes through the three vertices of the triangle. 2448 2449 Returns 2450 ======= 2451 2452 circumcircle : Circle 2453 2454 See Also 2455 ======== 2456 2457 sympy.geometry.ellipse.Circle 2458 2459 Examples 2460 ======== 2461 2462 >>> from sympy.geometry import Point, Triangle 2463 >>> p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, 1) 2464 >>> t = Triangle(p1, p2, p3) 2465 >>> t.circumcircle 2466 Circle(Point2D(1/2, 1/2), sqrt(2)/2) 2467 2468 """ 2469 return Circle(self.circumcenter, self.circumradius) 2470 2471 def bisectors(self): 2472 """The angle bisectors of the triangle. 2473 2474 An angle bisector of a triangle is a straight line through a vertex 2475 which cuts the corresponding angle in half. 2476 2477 Returns 2478 ======= 2479 2480 bisectors : dict 2481 Each key is a vertex (Point) and each value is the corresponding 2482 bisector (Segment). 2483 2484 See Also 2485 ======== 2486 2487 sympy.geometry.point.Point, sympy.geometry.line.Segment 2488 2489 Examples 2490 ======== 2491 2492 >>> from sympy.geometry import Point, Triangle, Segment 2493 >>> p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, 1) 2494 >>> t = Triangle(p1, p2, p3) 2495 >>> from sympy import sqrt 2496 >>> t.bisectors()[p2] == Segment(Point(1, 0), Point(0, sqrt(2) - 1)) 2497 True 2498 2499 """ 2500 # use lines containing sides so containment check during 2501 # intersection calculation can be avoided, thus reducing 2502 # the processing time for calculating the bisectors 2503 s = [Line(l) for l in self.sides] 2504 v = self.vertices 2505 c = self.incenter 2506 l1 = Segment(v[0], Line(v[0], c).intersection(s[1])[0]) 2507 l2 = Segment(v[1], Line(v[1], c).intersection(s[2])[0]) 2508 l3 = Segment(v[2], Line(v[2], c).intersection(s[0])[0]) 2509 return {v[0]: l1, v[1]: l2, v[2]: l3} 2510 2511 @property 2512 def incenter(self): 2513 """The center of the incircle. 2514 2515 The incircle is the circle which lies inside the triangle and touches 2516 all three sides. 2517 2518 Returns 2519 ======= 2520 2521 incenter : Point 2522 2523 See Also 2524 ======== 2525 2526 incircle, sympy.geometry.point.Point 2527 2528 Examples 2529 ======== 2530 2531 >>> from sympy.geometry import Point, Triangle 2532 >>> p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, 1) 2533 >>> t = Triangle(p1, p2, p3) 2534 >>> t.incenter 2535 Point2D(1 - sqrt(2)/2, 1 - sqrt(2)/2) 2536 2537 """ 2538 s = self.sides 2539 l = Matrix([s[i].length for i in [1, 2, 0]]) 2540 p = sum(l) 2541 v = self.vertices 2542 x = simplify(l.dot(Matrix([vi.x for vi in v]))/p) 2543 y = simplify(l.dot(Matrix([vi.y for vi in v]))/p) 2544 return Point(x, y) 2545 2546 @property 2547 def inradius(self): 2548 """The radius of the incircle. 2549 2550 Returns 2551 ======= 2552 2553 inradius : number of Basic instance 2554 2555 See Also 2556 ======== 2557 2558 incircle, sympy.geometry.ellipse.Circle.radius 2559 2560 Examples 2561 ======== 2562 2563 >>> from sympy.geometry import Point, Triangle 2564 >>> p1, p2, p3 = Point(0, 0), Point(4, 0), Point(0, 3) 2565 >>> t = Triangle(p1, p2, p3) 2566 >>> t.inradius 2567 1 2568 2569 """ 2570 return simplify(2 * self.area / self.perimeter) 2571 2572 @property 2573 def incircle(self): 2574 """The incircle of the triangle. 2575 2576 The incircle is the circle which lies inside the triangle and touches 2577 all three sides. 2578 2579 Returns 2580 ======= 2581 2582 incircle : Circle 2583 2584 See Also 2585 ======== 2586 2587 sympy.geometry.ellipse.Circle 2588 2589 Examples 2590 ======== 2591 2592 >>> from sympy.geometry import Point, Triangle 2593 >>> p1, p2, p3 = Point(0, 0), Point(2, 0), Point(0, 2) 2594 >>> t = Triangle(p1, p2, p3) 2595 >>> t.incircle 2596 Circle(Point2D(2 - sqrt(2), 2 - sqrt(2)), 2 - sqrt(2)) 2597 2598 """ 2599 return Circle(self.incenter, self.inradius) 2600 2601 @property 2602 def exradii(self): 2603 """The radius of excircles of a triangle. 2604 2605 An excircle of the triangle is a circle lying outside the triangle, 2606 tangent to one of its sides and tangent to the extensions of the 2607 other two. 2608 2609 Returns 2610 ======= 2611 2612 exradii : dict 2613 2614 See Also 2615 ======== 2616 2617 sympy.geometry.polygon.Triangle.inradius 2618 2619 Examples 2620 ======== 2621 2622 The exradius touches the side of the triangle to which it is keyed, e.g. 2623 the exradius touching side 2 is: 2624 2625 >>> from sympy.geometry import Point, Triangle 2626 >>> p1, p2, p3 = Point(0, 0), Point(6, 0), Point(0, 2) 2627 >>> t = Triangle(p1, p2, p3) 2628 >>> t.exradii[t.sides[2]] 2629 -2 + sqrt(10) 2630 2631 References 2632 ========== 2633 2634 [1] http://mathworld.wolfram.com/Exradius.html 2635 [2] http://mathworld.wolfram.com/Excircles.html 2636 2637 """ 2638 2639 side = self.sides 2640 a = side[0].length 2641 b = side[1].length 2642 c = side[2].length 2643 s = (a+b+c)/2 2644 area = self.area 2645 exradii = {self.sides[0]: simplify(area/(s-a)), 2646 self.sides[1]: simplify(area/(s-b)), 2647 self.sides[2]: simplify(area/(s-c))} 2648 2649 return exradii 2650 2651 @property 2652 def excenters(self): 2653 """Excenters of the triangle. 2654 2655 An excenter is the center of a circle that is tangent to a side of the 2656 triangle and the extensions of the other two sides. 2657 2658 Returns 2659 ======= 2660 2661 excenters : dict 2662 2663 2664 Examples 2665 ======== 2666 2667 The excenters are keyed to the side of the triangle to which their corresponding 2668 excircle is tangent: The center is keyed, e.g. the excenter of a circle touching 2669 side 0 is: 2670 2671 >>> from sympy.geometry import Point, Triangle 2672 >>> p1, p2, p3 = Point(0, 0), Point(6, 0), Point(0, 2) 2673 >>> t = Triangle(p1, p2, p3) 2674 >>> t.excenters[t.sides[0]] 2675 Point2D(12*sqrt(10), 2/3 + sqrt(10)/3) 2676 2677 See Also 2678 ======== 2679 2680 sympy.geometry.polygon.Triangle.exradii 2681 2682 References 2683 ========== 2684 2685 .. [1] http://mathworld.wolfram.com/Excircles.html 2686 2687 """ 2688 2689 s = self.sides 2690 v = self.vertices 2691 a = s[0].length 2692 b = s[1].length 2693 c = s[2].length 2694 x = [v[0].x, v[1].x, v[2].x] 2695 y = [v[0].y, v[1].y, v[2].y] 2696 2697 exc_coords = { 2698 "x1": simplify(-a*x[0]+b*x[1]+c*x[2]/(-a+b+c)), 2699 "x2": simplify(a*x[0]-b*x[1]+c*x[2]/(a-b+c)), 2700 "x3": simplify(a*x[0]+b*x[1]-c*x[2]/(a+b-c)), 2701 "y1": simplify(-a*y[0]+b*y[1]+c*y[2]/(-a+b+c)), 2702 "y2": simplify(a*y[0]-b*y[1]+c*y[2]/(a-b+c)), 2703 "y3": simplify(a*y[0]+b*y[1]-c*y[2]/(a+b-c)) 2704 } 2705 2706 excenters = { 2707 s[0]: Point(exc_coords["x1"], exc_coords["y1"]), 2708 s[1]: Point(exc_coords["x2"], exc_coords["y2"]), 2709 s[2]: Point(exc_coords["x3"], exc_coords["y3"]) 2710 } 2711 2712 return excenters 2713 2714 @property 2715 def medians(self): 2716 """The medians of the triangle. 2717 2718 A median of a triangle is a straight line through a vertex and the 2719 midpoint of the opposite side, and divides the triangle into two 2720 equal areas. 2721 2722 Returns 2723 ======= 2724 2725 medians : dict 2726 Each key is a vertex (Point) and each value is the median (Segment) 2727 at that point. 2728 2729 See Also 2730 ======== 2731 2732 sympy.geometry.point.Point.midpoint, sympy.geometry.line.Segment.midpoint 2733 2734 Examples 2735 ======== 2736 2737 >>> from sympy.geometry import Point, Triangle 2738 >>> p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, 1) 2739 >>> t = Triangle(p1, p2, p3) 2740 >>> t.medians[p1] 2741 Segment2D(Point2D(0, 0), Point2D(1/2, 1/2)) 2742 2743 """ 2744 s = self.sides 2745 v = self.vertices 2746 return {v[0]: Segment(v[0], s[1].midpoint), 2747 v[1]: Segment(v[1], s[2].midpoint), 2748 v[2]: Segment(v[2], s[0].midpoint)} 2749 2750 @property 2751 def medial(self): 2752 """The medial triangle of the triangle. 2753 2754 The triangle which is formed from the midpoints of the three sides. 2755 2756 Returns 2757 ======= 2758 2759 medial : Triangle 2760 2761 See Also 2762 ======== 2763 2764 sympy.geometry.line.Segment.midpoint 2765 2766 Examples 2767 ======== 2768 2769 >>> from sympy.geometry import Point, Triangle 2770 >>> p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, 1) 2771 >>> t = Triangle(p1, p2, p3) 2772 >>> t.medial 2773 Triangle(Point2D(1/2, 0), Point2D(1/2, 1/2), Point2D(0, 1/2)) 2774 2775 """ 2776 s = self.sides 2777 return Triangle(s[0].midpoint, s[1].midpoint, s[2].midpoint) 2778 2779 @property 2780 def nine_point_circle(self): 2781 """The nine-point circle of the triangle. 2782 2783 Nine-point circle is the circumcircle of the medial triangle, which 2784 passes through the feet of altitudes and the middle points of segments 2785 connecting the vertices and the orthocenter. 2786 2787 Returns 2788 ======= 2789 2790 nine_point_circle : Circle 2791 2792 See also 2793 ======== 2794 2795 sympy.geometry.line.Segment.midpoint 2796 sympy.geometry.polygon.Triangle.medial 2797 sympy.geometry.polygon.Triangle.orthocenter 2798 2799 Examples 2800 ======== 2801 2802 >>> from sympy.geometry import Point, Triangle 2803 >>> p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, 1) 2804 >>> t = Triangle(p1, p2, p3) 2805 >>> t.nine_point_circle 2806 Circle(Point2D(1/4, 1/4), sqrt(2)/4) 2807 2808 """ 2809 return Circle(*self.medial.vertices) 2810 2811 @property 2812 def eulerline(self): 2813 """The Euler line of the triangle. 2814 2815 The line which passes through circumcenter, centroid and orthocenter. 2816 2817 Returns 2818 ======= 2819 2820 eulerline : Line (or Point for equilateral triangles in which case all 2821 centers coincide) 2822 2823 Examples 2824 ======== 2825 2826 >>> from sympy.geometry import Point, Triangle 2827 >>> p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, 1) 2828 >>> t = Triangle(p1, p2, p3) 2829 >>> t.eulerline 2830 Line2D(Point2D(0, 0), Point2D(1/2, 1/2)) 2831 2832 """ 2833 if self.is_equilateral(): 2834 return self.orthocenter 2835 return Line(self.orthocenter, self.circumcenter) 2836 2837def rad(d): 2838 """Return the radian value for the given degrees (pi = 180 degrees).""" 2839 return d*pi/180 2840 2841 2842def deg(r): 2843 """Return the degree value for the given radians (pi = 180 degrees).""" 2844 return r/pi*180 2845 2846 2847def _slope(d): 2848 rv = tan(rad(d)) 2849 return rv 2850 2851 2852def _asa(d1, l, d2): 2853 """Return triangle having side with length l on the x-axis.""" 2854 xy = Line((0, 0), slope=_slope(d1)).intersection( 2855 Line((l, 0), slope=_slope(180 - d2)))[0] 2856 return Triangle((0, 0), (l, 0), xy) 2857 2858 2859def _sss(l1, l2, l3): 2860 """Return triangle having side of length l1 on the x-axis.""" 2861 c1 = Circle((0, 0), l3) 2862 c2 = Circle((l1, 0), l2) 2863 inter = [a for a in c1.intersection(c2) if a.y.is_nonnegative] 2864 if not inter: 2865 return None 2866 pt = inter[0] 2867 return Triangle((0, 0), (l1, 0), pt) 2868 2869 2870def _sas(l1, d, l2): 2871 """Return triangle having side with length l2 on the x-axis.""" 2872 p1 = Point(0, 0) 2873 p2 = Point(l2, 0) 2874 p3 = Point(cos(rad(d))*l1, sin(rad(d))*l1) 2875 return Triangle(p1, p2, p3) 2876