1"""Line-like geometrical entities. 2 3Contains 4======== 5LinearEntity 6Line 7Ray 8Segment 9""" 10 11from ..core import Dummy, Eq, Integer, factor_terms, oo, pi 12from ..core.compatibility import is_sequence 13from ..core.sympify import sympify 14from ..functions import Piecewise, acos, sqrt, tan 15from ..functions.elementary.trigonometric import _pi_coeff as pi_coeff 16from ..logic import And, false, true 17from ..simplify import simplify 18from ..solvers import solve 19from .entity import GeometryEntity, GeometrySet 20from .exceptions import GeometryError 21from .point import Point 22from .util import _symbol 23 24 25# TODO: this should be placed elsewhere and reused in other modules 26 27 28class Undecidable(ValueError): 29 """Raised when can't decide on relation.""" 30 31 32class LinearEntity(GeometrySet): 33 """A base class for all linear entities (line, ray and segment) 34 in a 2-dimensional Euclidean space. 35 36 Notes 37 ===== 38 39 This is an abstract class and is not meant to be instantiated. 40 41 See Also 42 ======== 43 44 diofant.geometry.entity.GeometryEntity 45 46 """ 47 48 def __new__(cls, p1, p2, **kwargs): 49 p1 = Point(p1) 50 p2 = Point(p2) 51 if p1 == p2: 52 # sometimes we return a single point if we are not given two unique 53 # points. This is done in the specific subclass 54 raise ValueError( 55 f'{cls.__name__}.__new__ requires two unique Points.') 56 57 return GeometryEntity.__new__(cls, p1, p2, **kwargs) 58 59 @property 60 def p1(self): 61 """The first defining point of a linear entity. 62 63 See Also 64 ======== 65 66 diofant.geometry.point.Point 67 68 Examples 69 ======== 70 71 >>> p1, p2 = Point(0, 0), Point(5, 3) 72 >>> l = Line(p1, p2) 73 >>> l.p1 74 Point(0, 0) 75 76 """ 77 return self.args[0] 78 79 @property 80 def p2(self): 81 """The second defining point of a linear entity. 82 83 See Also 84 ======== 85 86 diofant.geometry.point.Point 87 88 Examples 89 ======== 90 91 >>> p1, p2 = Point(0, 0), Point(5, 3) 92 >>> l = Line(p1, p2) 93 >>> l.p2 94 Point(5, 3) 95 96 """ 97 return self.args[1] 98 99 @property 100 def coefficients(self): 101 """The coefficients (`a`, `b`, `c`) for `ax + by + c = 0`. 102 103 See Also 104 ======== 105 106 diofant.geometry.line.Line.equation 107 108 Examples 109 ======== 110 111 >>> p1, p2 = Point(0, 0), Point(5, 3) 112 >>> l = Line(p1, p2) 113 >>> l.coefficients 114 (-3, 5, 0) 115 116 >>> p3 = Point(x, y) 117 >>> l2 = Line(p1, p3) 118 >>> l2.coefficients 119 (-y, x, 0) 120 121 """ 122 p1, p2 = self.points 123 if p1.x == p2.x: 124 return Integer(1), Integer(0), -p1.x 125 elif p1.y == p2.y: 126 return Integer(0), Integer(1), -p1.y 127 return tuple(simplify(i) 128 for i in (self.p1.y - self.p2.y, 129 self.p2.x - self.p1.x, 130 self.p1.x*self.p2.y - self.p1.y*self.p2.x)) 131 132 @staticmethod 133 def are_concurrent(*lines): 134 """Is a sequence of linear entities concurrent? 135 136 Two or more linear entities are concurrent if they all 137 intersect at a single point. 138 139 Parameters 140 ========== 141 142 lines : a sequence of linear entities. 143 144 Returns 145 ======= 146 147 True : if the set of linear entities are concurrent, 148 False : otherwise. 149 150 Notes 151 ===== 152 153 Simply take the first two lines and find their intersection. 154 If there is no intersection, then the first two lines were 155 parallel and had no intersection so concurrency is impossible 156 amongst the whole set. Otherwise, check to see if the 157 intersection point of the first two lines is a member on 158 the rest of the lines. If so, the lines are concurrent. 159 160 See Also 161 ======== 162 163 diofant.geometry.util.intersection 164 165 Examples 166 ======== 167 168 >>> p1, p2 = Point(0, 0), Point(3, 5) 169 >>> p3, p4 = Point(-2, -2), Point(0, 2) 170 >>> l1, l2, l3 = Line(p1, p2), Line(p1, p3), Line(p1, p4) 171 >>> Line.are_concurrent(l1, l2, l3) 172 True 173 174 >>> l4 = Line(p2, p3) 175 >>> Line.are_concurrent(l2, l3, l4) 176 False 177 178 """ 179 # Concurrency requires intersection at a single point; One linear 180 # entity cannot be concurrent. 181 if len(lines) <= 1: 182 return False 183 184 # Get the intersection (if parallel) 185 p = lines[0].intersection(lines[1]) 186 187 # Make sure the intersection is on every linear entity 188 for line in lines[2:]: 189 if p[0] not in line: 190 return False 191 return True 192 193 def is_parallel(self, other): 194 """Are two linear entities parallel? 195 196 Parameters 197 ========== 198 199 self : LinearEntity 200 other : LinearEntity 201 202 Returns 203 ======= 204 205 True : if self and other are parallel, 206 False : otherwise. 207 208 See Also 209 ======== 210 211 coefficients 212 213 Examples 214 ======== 215 216 >>> p1, p2 = Point(0, 0), Point(1, 1) 217 >>> p3, p4 = Point(3, 4), Point(6, 7) 218 >>> l1, l2 = Line(p1, p2), Line(p3, p4) 219 >>> Line.is_parallel(l1, l2) 220 True 221 222 >>> p5 = Point(6, 6) 223 >>> l3 = Line(p3, p5) 224 >>> Line.is_parallel(l1, l3) 225 False 226 227 """ 228 a1, b1, _ = self.coefficients 229 a2, b2, _ = other.coefficients 230 return bool(simplify(a1*b2 - b1*a2) == 0) 231 232 def is_perpendicular(self, other): 233 """Are two linear entities perpendicular? 234 235 Parameters 236 ========== 237 238 self : LinearEntity 239 other : LinearEntity 240 241 Returns 242 ======= 243 244 True : if self and other are perpendicular, 245 False : otherwise. 246 247 See Also 248 ======== 249 250 coefficients 251 252 Examples 253 ======== 254 255 >>> p1, p2, p3 = Point(0, 0), Point(1, 1), Point(-1, 1) 256 >>> l1, l2 = Line(p1, p2), Line(p1, p3) 257 >>> l1.is_perpendicular(l2) 258 True 259 260 >>> p4 = Point(5, 3) 261 >>> l3 = Line(p1, p4) 262 >>> l1.is_perpendicular(l3) 263 False 264 265 """ 266 a1, b1, _ = self.coefficients 267 a2, b2, _ = other.coefficients 268 return bool(simplify(a1*a2 + b1*b2) == 0) 269 270 def angle_between(self, other): 271 """The angle formed between the two linear entities. 272 273 Parameters 274 ========== 275 276 self : LinearEntity 277 other : LinearEntity 278 279 Returns 280 ======= 281 282 angle : angle in radians 283 284 Notes 285 ===== 286 287 From the dot product of vectors v1 and v2 it is known that: 288 289 ``dot(v1, v2) = |v1|*|v2|*cos(A)`` 290 291 where A is the angle formed between the two vectors. We can 292 get the directional vectors of the two lines and readily 293 find the angle between the two using the above formula. 294 295 See Also 296 ======== 297 298 is_perpendicular 299 300 Examples 301 ======== 302 303 >>> p1, p2, p3 = Point(0, 0), Point(0, 4), Point(2, 0) 304 >>> l1, l2 = Line(p1, p2), Line(p1, p3) 305 >>> l1.angle_between(l2) 306 pi/2 307 308 """ 309 v1 = self.p2 - self.p1 310 v2 = other.p2 - other.p1 311 return acos(v1.dot(v2)/(abs(v1)*abs(v2))) 312 313 def parallel_line(self, p): 314 """Create a new Line parallel to this linear entity which passes 315 through the point `p`. 316 317 Parameters 318 ========== 319 320 p : diofant.geometry.point.Point 321 322 Returns 323 ======= 324 325 line : Line 326 327 See Also 328 ======== 329 330 is_parallel 331 332 Examples 333 ======== 334 335 >>> p1, p2, p3 = Point(0, 0), Point(2, 3), Point(-2, 2) 336 >>> l1 = Line(p1, p2) 337 >>> l2 = l1.parallel_line(p3) 338 >>> p3 in l2 339 True 340 >>> l1.is_parallel(l2) 341 True 342 343 """ 344 d = self.p1 - self.p2 345 p = Point(p) 346 return Line(p, p + d) 347 348 def perpendicular_line(self, p): 349 """Create a new Line perpendicular to this linear entity which passes 350 through the point `p`. 351 352 Parameters 353 ========== 354 355 p : diofant.geometry.point.Point 356 357 Returns 358 ======= 359 360 line : Line 361 362 See Also 363 ======== 364 365 is_perpendicular, perpendicular_segment 366 367 Examples 368 ======== 369 370 >>> p1, p2, p3 = Point(0, 0), Point(2, 3), Point(-2, 2) 371 >>> l1 = Line(p1, p2) 372 >>> l2 = l1.perpendicular_line(p3) 373 >>> p3 in l2 374 True 375 >>> l1.is_perpendicular(l2) 376 True 377 378 """ 379 p = Point(p) 380 d1, d2 = (self.p1 - self.p2).args 381 if d2 == 0: # If a horizontal line 382 if p.y == self.p1.y: # if p is on this linear entity 383 return Line(p, p + Point(0, 1)) 384 else: 385 p2 = Point(p.x, self.p1.y) 386 return Line(p, p2) 387 else: 388 p2 = Point(p.x - d2, p.y + d1) 389 return Line(p, p2) 390 391 def perpendicular_segment(self, p): 392 """Create a perpendicular line segment from `p` to this line. 393 394 The enpoints of the segment are ``p`` and the closest point in 395 the line containing self. (If self is not a line, the point might 396 not be in self.) 397 398 Parameters 399 ========== 400 401 p : diofant.geometry.point.Point 402 403 Returns 404 ======= 405 406 segment : Segment 407 408 Notes 409 ===== 410 411 Returns `p` itself if `p` is on this linear entity. 412 413 See Also 414 ======== 415 416 perpendicular_line 417 418 Examples 419 ======== 420 421 >>> p1, p2, p3 = Point(0, 0), Point(1, 1), Point(0, 2) 422 >>> l1 = Line(p1, p2) 423 >>> s1 = l1.perpendicular_segment(p3) 424 >>> l1.is_perpendicular(s1) 425 True 426 >>> p3 in s1 427 True 428 >>> l1.perpendicular_segment(Point(4, 0)) 429 Segment(Point(2, 2), Point(4, 0)) 430 431 """ 432 p = Point(p) 433 if p in self: 434 return p 435 a, b, c = self.coefficients 436 if a == 0: # horizontal 437 p2 = Point(p.x, self.p1.y) 438 elif b == 0: # vertical 439 p2 = Point(self.p1.x, p.y) 440 else: 441 # ax + by + c = 0 442 y = (-c - a*p.x)/b 443 m = self.slope 444 d2 = 1 + m**2 445 H = p.y - y 446 dx = m*H/d2 447 dy = m*dx 448 p2 = (p.x + dx, y + dy) 449 return Segment(p, p2) 450 451 @property 452 def length(self): 453 """ 454 The length of the line. 455 456 Examples 457 ======== 458 459 >>> p1, p2 = Point(0, 0), Point(3, 5) 460 >>> l1 = Line(p1, p2) 461 >>> l1.length 462 oo 463 464 """ 465 return oo 466 467 @property 468 def slope(self): 469 """The slope of this linear entity, or infinity if vertical. 470 471 Returns 472 ======= 473 474 slope : Expr 475 476 See Also 477 ======== 478 479 coefficients 480 481 Examples 482 ======== 483 484 >>> p1, p2 = Point(0, 0), Point(3, 5) 485 >>> l1 = Line(p1, p2) 486 >>> l1.slope 487 5/3 488 489 >>> p3 = Point(0, 4) 490 >>> l2 = Line(p1, p3) 491 >>> l2.slope 492 oo 493 494 """ 495 d1, d2 = (self.p1 - self.p2).args 496 if d1 == 0: 497 return oo 498 return simplify(d2/d1) 499 500 @property 501 def points(self): 502 """The two points used to define this linear entity. 503 504 Returns 505 ======= 506 507 points : tuple of Points 508 509 See Also 510 ======== 511 512 diofant.geometry.point.Point 513 514 Examples 515 ======== 516 517 >>> p1, p2 = Point(0, 0), Point(5, 11) 518 >>> l1 = Line(p1, p2) 519 >>> l1.points 520 (Point(0, 0), Point(5, 11)) 521 522 """ 523 return self.p1, self.p2 524 525 def projection(self, o): 526 """Project a point, line, ray, or segment onto this linear entity. 527 528 Parameters 529 ========== 530 531 other : diofant.geometry.point.Point or LinearEntity 532 533 Returns 534 ======= 535 536 projection : diofant.geometry.point.Point or LinearEntity 537 The return type matches the type of the parameter ``other``. 538 539 Raises 540 ====== 541 542 diofant.geometry.exceptions.GeometryError 543 When method is unable to perform projection. 544 545 Notes 546 ===== 547 548 A projection involves taking the two points that define 549 the linear entity and projecting those points onto a 550 Line and then reforming the linear entity using these 551 projections. 552 A point P is projected onto a line L by finding the point 553 on L that is closest to P. This point is the intersection 554 of L and the line perpendicular to L that passes through P. 555 556 See Also 557 ======== 558 559 diofant.geometry.point.Point, perpendicular_line 560 561 Examples 562 ======== 563 564 >>> p1, p2, p3 = Point(0, 0), Point(1, 1), Point(Rational(1, 2), 0) 565 >>> l1 = Line(p1, p2) 566 >>> l1.projection(p3) 567 Point(1/4, 1/4) 568 569 >>> p4, p5 = Point(10, 0), Point(12, 1) 570 >>> s1 = Segment(p4, p5) 571 >>> l1.projection(s1) 572 Segment(Point(5, 5), Point(13/2, 13/2)) 573 574 """ 575 tline = Line(self.p1, self.p2) 576 577 def _project(p): 578 """Project a point onto the line representing self.""" 579 if p in tline: 580 return p 581 l1 = tline.perpendicular_line(p) 582 return tline.intersection(l1)[0] 583 584 projected = None 585 if isinstance(o, Point): 586 return _project(o) 587 elif isinstance(o, LinearEntity): 588 n_p1 = _project(o.p1) 589 n_p2 = _project(o.p2) 590 if n_p1 == n_p2: 591 projected = n_p1 592 else: 593 projected = o.__class__(n_p1, n_p2) 594 595 # Didn't know how to project so raise an error 596 if projected is None: 597 n1 = self.__class__.__name__ 598 n2 = o.__class__.__name__ 599 raise GeometryError( 600 f'Do not know how to project {n2} onto {n1}') 601 602 return self.intersection(projected)[0] 603 604 def intersection(self, o): 605 """The intersection with another geometrical entity. 606 607 Parameters 608 ========== 609 610 o : diofant.geometry.point.Point or LinearEntity 611 612 Returns 613 ======= 614 615 intersection : list of geometrical entities 616 617 See Also 618 ======== 619 620 diofant.geometry.point.Point 621 622 Examples 623 ======== 624 625 >>> p1, p2, p3 = Point(0, 0), Point(1, 1), Point(7, 7) 626 >>> l1 = Line(p1, p2) 627 >>> l1.intersection(p3) 628 [Point(7, 7)] 629 630 >>> p4, p5 = Point(5, 0), Point(0, 3) 631 >>> l2 = Line(p4, p5) 632 >>> l1.intersection(l2) 633 [Point(15/8, 15/8)] 634 635 >>> p6, p7 = Point(0, 5), Point(2, 6) 636 >>> s1 = Segment(p6, p7) 637 >>> l1.intersection(s1) 638 [] 639 640 """ 641 if isinstance(o, Point): 642 if o in self: 643 return [o] 644 else: 645 return [] 646 647 elif isinstance(o, LinearEntity): 648 a1, b1, c1 = self.coefficients 649 a2, b2, c2 = o.coefficients 650 t = simplify(a1*b2 - a2*b1) 651 if t.equals(0) is not False: # assume they are parallel 652 if isinstance(self, Line): 653 if o.p1 in self: 654 return [o] 655 return [] 656 elif isinstance(o, Line): 657 if self.p1 in o: 658 return [self] 659 return [] 660 elif isinstance(self, Ray): 661 if isinstance(o, Ray): 662 # case 1, rays in the same direction 663 if self.xdirection == o.xdirection and \ 664 self.ydirection == o.ydirection: 665 return [self] if (self.source in o) else [o] 666 # case 2, rays in the opposite directions 667 else: 668 if o.source in self: 669 if self.source == o.source: 670 return [self.source] 671 return [Segment(o.source, self.source)] 672 return [] 673 elif isinstance(o, Segment): 674 if o.p1 in self: 675 if o.p2 in self: 676 return [o] 677 return [Segment(o.p1, self.source)] 678 elif o.p2 in self: 679 return [Segment(o.p2, self.source)] 680 return [] 681 else: 682 raise NotImplementedError 683 elif isinstance(self, Segment): 684 if isinstance(o, Ray): 685 return o.intersection(self) 686 elif isinstance(o, Segment): 687 # A reminder that the points of Segments are ordered 688 # in such a way that the following works. See 689 # Segment.__new__ for details on the ordering. 690 if self.p1 not in o: 691 if self.p2 not in o: 692 # Neither of the endpoints are in o so either 693 # o is contained in this segment or it isn't 694 if o in self: 695 return [self] 696 return [] 697 else: 698 # p1 not in o but p2 is. Either there is a 699 # segment as an intersection, or they only 700 # intersect at an endpoint 701 if self.p2 == o.p1: 702 return [o.p1] 703 return [Segment(o.p1, self.p2)] 704 elif self.p2 not in o: 705 # p2 not in o but p1 is. Either there is a 706 # segment as an intersection, or they only 707 # intersect at an endpoint 708 if self.p1 == o.p2: 709 return [o.p2] 710 return [Segment(o.p2, self.p1)] 711 712 # Both points of self in o so the whole segment 713 # is in o 714 return [self] 715 else: 716 raise NotImplementedError 717 else: 718 raise NotImplementedError 719 720 # Not parallel, so find the point of intersection 721 px = simplify((b1*c2 - c1*b2) / t) 722 py = simplify((a2*c1 - a1*c2) / t) 723 inter = Point(px, py) 724 # we do not use a simplistic 'inter in self and inter in o' 725 # because that requires an equality test that is fragile; 726 # instead we employ some diagnostics to see if the intersection 727 # is valid 728 729 def inseg(self): 730 def _between(a, b, c): 731 return c >= a and c <= b or c <= a and c >= b 732 if _between(self.p1.x, self.p2.x, inter.x) and \ 733 _between(self.p1.y, self.p2.y, inter.y): 734 return True 735 736 def inray(self): 737 if self.p1 == inter: 738 return True 739 sray = Ray(self.p1, inter) 740 if sray.xdirection == self.xdirection and \ 741 sray.ydirection == self.ydirection: 742 return True 743 744 prec = (Line, Ray, Segment) 745 expr = self 746 if prec.index(expr.func) > prec.index(o.func): 747 expr, o = o, expr 748 rv = [inter] 749 if isinstance(expr, Line): 750 if isinstance(o, Line): 751 return rv 752 elif isinstance(o, Ray) and inray(o): 753 return rv 754 elif isinstance(o, Segment) and inseg(o): 755 return rv 756 elif isinstance(expr, Ray) and inray(expr): 757 if isinstance(o, Ray) and inray(o): 758 return rv 759 elif isinstance(o, Segment) and inseg(o): 760 return rv 761 elif isinstance(expr, Segment) and inseg(expr): 762 if isinstance(o, Segment) and inseg(o): 763 return rv 764 return [] 765 766 return o.intersection(self) 767 768 def arbitrary_point(self, parameter='t'): 769 """A parameterized point on the Line. 770 771 Parameters 772 ========== 773 774 parameter : str, optional 775 The name of the parameter which will be used for the parametric 776 point. The default value is 't'. When this parameter is 0, the 777 first point used to define the line will be returned, and when 778 it is 1 the second point will be returned. 779 780 Returns 781 ======= 782 783 point : diofant.geometry.point.Point 784 785 Raises 786 ====== 787 788 ValueError 789 When ``parameter`` already appears in the Line's definition. 790 791 See Also 792 ======== 793 794 diofant.geometry.point.Point 795 796 Examples 797 ======== 798 799 >>> p1, p2 = Point(1, 0), Point(5, 3) 800 >>> l1 = Line(p1, p2) 801 >>> l1.arbitrary_point() 802 Point(4*t + 1, 3*t) 803 804 """ 805 t = _symbol(parameter) 806 if t.name in (f.name for f in self.free_symbols): 807 raise ValueError(f'Symbol {t.name} already appears in object ' 808 'and cannot be used as a parameter.') 809 # multiply on the right so the variable gets 810 # combined witht he coordinates of the point 811 return self.p1 + (self.p2 - self.p1)*t 812 813 def random_point(self): 814 """A random point on a LinearEntity. 815 816 Returns 817 ======= 818 819 point : diofant.geometry.point.Point 820 821 See Also 822 ======== 823 824 diofant.geometry.point.Point 825 826 Examples 827 ======== 828 829 >>> p1, p2 = Point(0, 0), Point(5, 3) 830 >>> l1 = Line(p1, p2) 831 >>> p3 = l1.random_point() 832 >>> # random point - don't know its coords in advance 833 >>> p3 834 Point(...) 835 >>> # point should belong to the line 836 >>> p3 in l1 837 True 838 839 """ 840 from random import randint 841 842 # The lower and upper 843 lower, upper = -2**32 - 1, 2**32 844 845 if self.slope is oo: 846 if isinstance(self, Ray): 847 if self.ydirection is oo: 848 lower = self.p1.y 849 else: 850 upper = self.p1.y 851 elif isinstance(self, Segment): 852 lower = self.p1.y 853 upper = self.p2.y 854 855 x = self.p1.x 856 y = randint(lower, upper) 857 else: 858 if isinstance(self, Ray): 859 if self.xdirection is oo: 860 lower = self.p1.x 861 else: 862 upper = self.p1.x 863 elif isinstance(self, Segment): 864 lower = self.p1.x 865 upper = self.p2.x 866 867 a, b, c = self.coefficients 868 x = randint(lower, upper) 869 y = (-c - a*x) / b 870 return Point(x, y) 871 872 def is_similar(self, other): 873 """ 874 Return True if self and other are contained in the same line. 875 876 Examples 877 ======== 878 879 >>> p1, p2, p3 = Point(0, 1), Point(3, 4), Point(2, 3) 880 >>> l1 = Line(p1, p2) 881 >>> l2 = Line(p1, p3) 882 >>> l1.is_similar(l2) 883 True 884 885 """ 886 def _norm(a, b, c): 887 if a != 0: 888 return 1, b/a, c/a 889 elif b != 0: 890 return a/b, 1, c/b 891 else: 892 raise NotImplementedError 893 return _norm(*self.coefficients) == _norm(*other.coefficients) 894 895 def __contains__(self, other): 896 """Return a definitive answer or else raise an error if it cannot 897 be determined that other is on the boundaries of self. 898 899 """ 900 result = self.contains(other) 901 902 if result is not None: 903 return result 904 else: 905 raise Undecidable( 906 f"can't decide whether '{self}' contains '{other}'") 907 908 def contains(self, other): 909 """Subclasses should implement this method and should return 910 True if other is on the boundaries of self; 911 False if not on the boundaries of self; 912 None if a determination cannot be made. 913 914 """ 915 raise NotImplementedError() 916 917 918class Line(LinearEntity): 919 """An infinite line in space. 920 921 A line is declared with two distinct points or a point and slope 922 as defined using keyword `slope`. 923 924 Notes 925 ===== 926 927 At the moment only lines in a 2D space can be declared, because 928 Points can be defined only for 2D spaces. 929 930 Parameters 931 ========== 932 933 p1 : diofant.geometry.point.Point 934 pt : diofant.geometry.point.Point 935 slope : Expr 936 937 See Also 938 ======== 939 940 diofant.geometry.point.Point 941 942 Examples 943 ======== 944 945 >>> L = Line(Point(2, 3), Point(3, 5)) 946 >>> L 947 Line(Point(2, 3), Point(3, 5)) 948 >>> L.points 949 (Point(2, 3), Point(3, 5)) 950 >>> L.equation() 951 -2*x + y + 1 952 >>> L.coefficients 953 (-2, 1, 1) 954 955 Instantiate with keyword ``slope``: 956 957 >>> Line(Point(0, 0), slope=0) 958 Line(Point(0, 0), Point(1, 0)) 959 960 Instantiate with another linear object 961 962 >>> s = Segment((0, 0), (0, 1)) 963 >>> Line(s).equation() 964 x 965 966 """ 967 968 def __new__(cls, p1, pt=None, slope=None, **kwargs): 969 if isinstance(p1, LinearEntity): 970 p1, pt = p1.args 971 else: 972 p1 = Point(p1) 973 if pt is not None and slope is None: 974 p2 = Point(pt) 975 elif slope is not None and pt is None: 976 slope = sympify(slope) 977 if slope.is_finite is False: 978 # when infinite slope, don't change x 979 dx = 0 980 dy = 1 981 else: 982 # go over 1 up slope 983 dx = 1 984 dy = slope 985 # XXX avoiding simplification by adding to coords directly 986 p2 = Point(p1.x + dx, p1.y + dy) 987 else: 988 raise ValueError('A 2nd Point or keyword "slope" must be used.') 989 990 return LinearEntity.__new__(cls, p1, p2, **kwargs) 991 992 def plot_interval(self, parameter='t'): 993 """The plot interval for the default geometric plot of line. Gives 994 values that will produce a line that is +/- 5 units long (where a 995 unit is the distance between the two points that define the line). 996 997 Parameters 998 ========== 999 1000 parameter : str, optional 1001 Default value is 't'. 1002 1003 Returns 1004 ======= 1005 1006 plot_interval : list (plot interval) 1007 [parameter, lower_bound, upper_bound] 1008 1009 Examples 1010 ======== 1011 1012 >>> p1, p2 = Point(0, 0), Point(5, 3) 1013 >>> l1 = Line(p1, p2) 1014 >>> l1.plot_interval() 1015 [t, -5, 5] 1016 1017 """ 1018 t = _symbol(parameter) 1019 return [t, -5, 5] 1020 1021 def equation(self, x='x', y='y'): 1022 """The equation of the line: ax + by + c. 1023 1024 Parameters 1025 ========== 1026 1027 x : str, optional 1028 The name to use for the x-axis, default value is 'x'. 1029 y : str, optional 1030 The name to use for the y-axis, default value is 'y'. 1031 1032 Returns 1033 ======= 1034 1035 equation : diofant expression 1036 1037 See Also 1038 ======== 1039 1040 LinearEntity.coefficients 1041 1042 Examples 1043 ======== 1044 1045 >>> p1, p2 = Point(1, 0), Point(5, 3) 1046 >>> l1 = Line(p1, p2) 1047 >>> l1.equation() 1048 -3*x + 4*y + 3 1049 1050 """ 1051 x, y = _symbol(x), _symbol(y) 1052 p1, p2 = self.points 1053 if p1.x == p2.x: 1054 return x - p1.x 1055 elif p1.y == p2.y: 1056 return y - p1.y 1057 1058 a, b, c = self.coefficients 1059 return a*x + b*y + c 1060 1061 def contains(self, other): 1062 """ 1063 Return True if o is on this Line, or False otherwise. 1064 1065 Examples 1066 ======== 1067 1068 >>> p1, p2 = Point(0, 1), Point(3, 4) 1069 >>> l = Line(p1, p2) 1070 >>> l.contains(p1) 1071 True 1072 >>> l.contains((0, 1)) 1073 True 1074 >>> l.contains((0, 0)) 1075 False 1076 1077 """ 1078 if is_sequence(other): 1079 other = Point(other) 1080 if isinstance(other, Point): 1081 other = other.func(*[simplify(i) for i in other.args]) 1082 x, y = Dummy(), Dummy() 1083 eq = self.equation(x, y) 1084 if not eq.has(y): 1085 return (solve(eq, x)[0][x] - other.x).equals(0) 1086 if not eq.has(x): 1087 return (solve(eq, y)[0][y] - other.y).equals(0) 1088 return (solve(eq.subs({x: other.x}), y)[0][y] - other.y).equals(0) 1089 elif not isinstance(other, LinearEntity): 1090 return False 1091 elif isinstance(other, Line): 1092 return self.equal(other) 1093 elif not self.is_similar(other): 1094 return False 1095 else: 1096 return other.p1 in self and other.p2 in self 1097 1098 def distance(self, o): 1099 """ 1100 Finds the shortest distance between a line and a point. 1101 1102 Raises 1103 ====== 1104 1105 NotImplementedError 1106 if o is not a Point 1107 1108 Examples 1109 ======== 1110 1111 >>> p1, p2 = Point(0, 0), Point(1, 1) 1112 >>> s = Line(p1, p2) 1113 >>> s.distance(Point(-1, 1)) 1114 sqrt(2) 1115 >>> s.distance((-1, 2)) 1116 3*sqrt(2)/2 1117 1118 """ 1119 if not isinstance(o, Point): 1120 o = Point(o) 1121 a, b, c = self.coefficients 1122 if 0 in (a, b): 1123 return self.perpendicular_segment(o).length 1124 m = self.slope 1125 x = o.x 1126 y = m*x - c/b 1127 return abs(factor_terms(o.y - y))/sqrt(1 + m**2) 1128 1129 def equal(self, other): 1130 """Returns True if self and other are the same mathematical entities.""" 1131 if not isinstance(other, Line): 1132 return False 1133 return Point.is_collinear(self.p1, other.p1, self.p2, other.p2) 1134 1135 1136class Ray(LinearEntity): 1137 """ 1138 A Ray is a semi-line in the space with a source point and a direction. 1139 1140 Parameters 1141 ========== 1142 1143 p1 : diofant.geometry.point.Point 1144 The source of the Ray 1145 p2 : diofant.geometry.point.Point or Expr 1146 This point determines the direction in which the Ray propagates. 1147 If given as an angle it is interpreted in radians with the positive 1148 direction being ccw. 1149 1150 See Also 1151 ======== 1152 1153 diofant.geometry.point.Point, Line 1154 1155 Notes 1156 ===== 1157 1158 At the moment only rays in a 2D space can be declared, because 1159 Points can be defined only for 2D spaces. 1160 1161 Examples 1162 ======== 1163 1164 >>> r = Ray(Point(2, 3), Point(3, 5)) 1165 >>> r = Ray(Point(2, 3), Point(3, 5)) 1166 >>> r 1167 Ray(Point(2, 3), Point(3, 5)) 1168 >>> r.points 1169 (Point(2, 3), Point(3, 5)) 1170 >>> r.source 1171 Point(2, 3) 1172 >>> r.xdirection 1173 oo 1174 >>> r.ydirection 1175 oo 1176 >>> r.slope 1177 2 1178 >>> Ray(Point(0, 0), angle=pi/4).slope 1179 1 1180 1181 """ 1182 1183 def __new__(cls, p1, pt=None, angle=None, **kwargs): 1184 p1 = Point(p1) 1185 if pt is not None and angle is None: 1186 try: 1187 p2 = Point(pt) 1188 except ValueError: 1189 from ..utilities import filldedent 1190 raise ValueError(filldedent(""" 1191 The 2nd argument was not a valid Point; if 1192 it was meant to be an angle it should be 1193 given with keyword "angle".""")) 1194 if p1 == p2: 1195 raise ValueError('A Ray requires two distinct points.') 1196 elif angle is not None and pt is None: 1197 # we need to know if the angle is an odd multiple of pi/2 1198 c = pi_coeff(sympify(angle)) 1199 p2 = None 1200 if c is not None: 1201 if c.is_Rational: 1202 if c.denominator == 2: 1203 if c.numerator == 1: 1204 p2 = p1 + Point(0, 1) 1205 else: 1206 assert c.numerator == 3 1207 p2 = p1 + Point(0, -1) 1208 elif c.denominator == 1: 1209 if c.numerator == 0: 1210 p2 = p1 + Point(1, 0) 1211 else: 1212 assert c.numerator == 1 1213 p2 = p1 + Point(-1, 0) 1214 if p2 is None: 1215 c *= pi 1216 else: 1217 c = angle % (2*pi) 1218 if not p2: 1219 m = 2*c/pi 1220 left = And(1 < m, m < 3) # is it in quadrant 2 or 3? 1221 x = Piecewise((-1, left), (Piecewise((0, Eq(m % 1, 0)), (1, True)), True)) 1222 y = Piecewise((-tan(c), left), (Piecewise((1, Eq(m, 1)), (-1, Eq(m, 3)), (tan(c), True)), True)) 1223 p2 = p1 + Point(x, y) 1224 else: 1225 raise ValueError('A 2nd point or keyword "angle" must be used.') 1226 1227 return LinearEntity.__new__(cls, p1, p2, **kwargs) 1228 1229 @property 1230 def source(self): 1231 """The point from which the ray emanates. 1232 1233 See Also 1234 ======== 1235 1236 diofant.geometry.point.Point 1237 1238 Examples 1239 ======== 1240 1241 >>> p1, p2 = Point(0, 0), Point(4, 1) 1242 >>> r1 = Ray(p1, p2) 1243 >>> r1.source 1244 Point(0, 0) 1245 1246 """ 1247 return self.p1 1248 1249 @property 1250 def direction(self): 1251 """The direction in which the ray emanates. 1252 1253 See Also 1254 ======== 1255 1256 diofant.geometry.point.Point 1257 1258 Examples 1259 ======== 1260 1261 >>> p1, p2 = Point(0, 0), Point(4, 1) 1262 >>> r1 = Ray(p1, p2) 1263 >>> r1.direction 1264 Point(4, 1) 1265 1266 """ 1267 return self.p2 - self.p1 1268 1269 @property 1270 def xdirection(self): 1271 """The x direction of the ray. 1272 1273 Positive infinity if the ray points in the positive x direction, 1274 negative infinity if the ray points in the negative x direction, 1275 or 0 if the ray is vertical. 1276 1277 See Also 1278 ======== 1279 1280 ydirection 1281 1282 Examples 1283 ======== 1284 1285 >>> p1, p2, p3 = Point(0, 0), Point(1, 1), Point(0, -1) 1286 >>> r1, r2 = Ray(p1, p2), Ray(p1, p3) 1287 >>> r1.xdirection 1288 oo 1289 >>> r2.xdirection 1290 0 1291 1292 """ 1293 if self.p1.x < self.p2.x: 1294 return oo 1295 elif self.p1.x == self.p2.x: 1296 return Integer(0) 1297 else: 1298 return -oo 1299 1300 @property 1301 def ydirection(self): 1302 """The y direction of the ray. 1303 1304 Positive infinity if the ray points in the positive y direction, 1305 negative infinity if the ray points in the negative y direction, 1306 or 0 if the ray is horizontal. 1307 1308 See Also 1309 ======== 1310 1311 xdirection 1312 1313 Examples 1314 ======== 1315 1316 >>> p1, p2, p3 = Point(0, 0), Point(-1, -1), Point(-1, 0) 1317 >>> r1, r2 = Ray(p1, p2), Ray(p1, p3) 1318 >>> r1.ydirection 1319 -oo 1320 >>> r2.ydirection 1321 0 1322 1323 """ 1324 if self.p1.y < self.p2.y: 1325 return oo 1326 elif self.p1.y == self.p2.y: 1327 return Integer(0) 1328 else: 1329 return -oo 1330 1331 def distance(self, o): 1332 """ 1333 Finds the shortest distance between the ray and a point. 1334 1335 Raises 1336 ====== 1337 1338 NotImplementedError 1339 if o is not a Point 1340 1341 Examples 1342 ======== 1343 1344 >>> p1, p2 = Point(0, 0), Point(1, 1) 1345 >>> s = Ray(p1, p2) 1346 >>> s.distance(Point(-1, -1)) 1347 sqrt(2) 1348 >>> s.distance((-1, 2)) 1349 3*sqrt(2)/2 1350 1351 """ 1352 if not isinstance(o, Point): 1353 o = Point(o) 1354 s = self.perpendicular_segment(o) 1355 if isinstance(s, Point): 1356 if self.contains(s): 1357 return Integer(0) 1358 else: 1359 # since arg-order is arbitrary, find the non-o point 1360 non_o = s.p1 if s.p1 != o else s.p2 1361 if self.contains(non_o): 1362 return Line(self).distance(o) # = s.length but simpler 1363 # the following applies when neither of the above apply 1364 return self.source.distance(o) 1365 1366 def plot_interval(self, parameter='t'): 1367 """The plot interval for the default geometric plot of the Ray. Gives 1368 values that will produce a ray that is 10 units long (where a unit is 1369 the distance between the two points that define the ray). 1370 1371 Parameters 1372 ========== 1373 1374 parameter : str, optional 1375 Default value is 't'. 1376 1377 Returns 1378 ======= 1379 1380 plot_interval : list 1381 [parameter, lower_bound, upper_bound] 1382 1383 Examples 1384 ======== 1385 1386 >>> r = Ray((0, 0), angle=pi/4) 1387 >>> r.plot_interval() 1388 [t, 0, 10] 1389 1390 """ 1391 t = _symbol(parameter) 1392 return [t, 0, 10] 1393 1394 def contains(self, other): 1395 """ 1396 Is other GeometryEntity contained in this Ray? 1397 1398 Examples 1399 ======== 1400 1401 >>> p1, p2 = Point(0, 0), Point(4, 4) 1402 >>> r = Ray(p1, p2) 1403 >>> r.contains(p1) 1404 True 1405 >>> r.contains((1, 1)) 1406 True 1407 >>> r.contains((1, 3)) 1408 False 1409 >>> s = Segment((1, 1), (2, 2)) 1410 >>> r.contains(s) 1411 True 1412 >>> s = Segment((1, 2), (2, 5)) 1413 >>> r.contains(s) 1414 False 1415 >>> r1 = Ray((2, 2), (3, 3)) 1416 >>> r.contains(r1) 1417 True 1418 >>> r1 = Ray((2, 2), (3, 5)) 1419 >>> r.contains(r1) 1420 False 1421 1422 """ 1423 if isinstance(other, Ray): 1424 return (Point.is_collinear(self.p1, self.p2, other.p1, other.p2) and 1425 self.xdirection == other.xdirection and 1426 self.ydirection == other.ydirection) 1427 elif isinstance(other, Segment): 1428 return other.p1 in self and other.p2 in self 1429 elif is_sequence(other): 1430 other = Point(other) 1431 if isinstance(other, Point): 1432 if Point.is_collinear(self.p1, self.p2, other): 1433 if self.xdirection is oo: 1434 rv = other.x >= self.source.x 1435 elif self.xdirection == -oo: 1436 rv = other.x <= self.source.x 1437 elif self.ydirection is oo: 1438 rv = other.y >= self.source.y 1439 else: 1440 rv = other.y <= self.source.y 1441 if rv in (true, false): 1442 return bool(rv) 1443 raise Undecidable( 1444 f'Cannot determine if {other} is in {self}') 1445 else: 1446 # Points are not collinear, so the rays are not parallel 1447 # and hence it is impossible for self to contain o 1448 return False 1449 1450 # No other known entity can be contained in a Ray 1451 return False 1452 1453 1454class Segment(LinearEntity): 1455 """A undirected line segment in space. 1456 1457 Parameters 1458 ========== 1459 1460 p1 : diofant.geometry.point.Point 1461 p2 : diofant.geometry.point.Point 1462 1463 See Also 1464 ======== 1465 1466 diofant.geometry.point.Point, Line 1467 1468 Notes 1469 ===== 1470 1471 At the moment only segments in a 2D space can be declared, because 1472 Points can be defined only for 2D spaces. 1473 1474 Examples 1475 ======== 1476 1477 >>> Segment((1, 0), (1, 1)) # tuples are interpreted as pts 1478 Segment(Point(1, 0), Point(1, 1)) 1479 >>> s = Segment(Point(4, 3), Point(1, 1)) 1480 >>> s 1481 Segment(Point(1, 1), Point(4, 3)) 1482 >>> s.points 1483 (Point(1, 1), Point(4, 3)) 1484 >>> s.slope 1485 2/3 1486 >>> s.length 1487 sqrt(13) 1488 >>> s.midpoint 1489 Point(5/2, 2) 1490 1491 """ 1492 1493 def __new__(cls, p1, p2, **kwargs): 1494 # Reorder the two points under the following ordering: 1495 # if p1.x != p2.x then p1.x < p2.x 1496 # if p1.x == p2.x then p1.y < p2.y 1497 p1 = Point(p1) 1498 p2 = Point(p2) 1499 if p1 == p2: 1500 return Point(p1) 1501 if (p1.x - p2.x).is_positive: 1502 p1, p2 = p2, p1 1503 elif (p1.x == p2.x) and (p1.y - p2.y).is_positive: 1504 p1, p2 = p2, p1 1505 return LinearEntity.__new__(cls, p1, p2, **kwargs) 1506 1507 def plot_interval(self, parameter='t'): 1508 """The plot interval for the default geometric plot of the Segment gives 1509 values that will produce the full segment in a plot. 1510 1511 Parameters 1512 ========== 1513 1514 parameter : str, optional 1515 Default value is 't'. 1516 1517 Returns 1518 ======= 1519 1520 plot_interval : list 1521 [parameter, lower_bound, upper_bound] 1522 1523 Examples 1524 ======== 1525 1526 >>> p1, p2 = Point(0, 0), Point(5, 3) 1527 >>> s1 = Segment(p1, p2) 1528 >>> s1.plot_interval() 1529 [t, 0, 1] 1530 1531 """ 1532 t = _symbol(parameter) 1533 return [t, 0, 1] 1534 1535 def perpendicular_bisector(self, p=None): 1536 """The perpendicular bisector of this segment. 1537 1538 If no point is specified or the point specified is not on the 1539 bisector then the bisector is returned as a Line. Otherwise a 1540 Segment is returned that joins the point specified and the 1541 intersection of the bisector and the segment. 1542 1543 Parameters 1544 ========== 1545 1546 p : diofant.geometry.point.Point 1547 1548 Returns 1549 ======= 1550 1551 bisector : Line or Segment 1552 1553 See Also 1554 ======== 1555 1556 LinearEntity.perpendicular_segment 1557 1558 Examples 1559 ======== 1560 1561 >>> p1, p2, p3 = Point(0, 0), Point(6, 6), Point(5, 1) 1562 >>> s1 = Segment(p1, p2) 1563 >>> s1.perpendicular_bisector() 1564 Line(Point(3, 3), Point(9, -3)) 1565 1566 >>> s1.perpendicular_bisector(p3) 1567 Segment(Point(3, 3), Point(5, 1)) 1568 1569 """ 1570 l = LinearEntity.perpendicular_line(self, self.midpoint) 1571 if p is None or Point(p) not in l: 1572 return l 1573 else: 1574 return Segment(self.midpoint, p) 1575 1576 @property 1577 def length(self): 1578 """The length of the line segment. 1579 1580 See Also 1581 ======== 1582 1583 diofant.geometry.point.Point.distance 1584 1585 Examples 1586 ======== 1587 1588 >>> p1, p2 = Point(0, 0), Point(4, 3) 1589 >>> s1 = Segment(p1, p2) 1590 >>> s1.length 1591 5 1592 1593 """ 1594 return Point.distance(self.p1, self.p2) 1595 1596 @property 1597 def midpoint(self): 1598 """The midpoint of the line segment. 1599 1600 See Also 1601 ======== 1602 1603 diofant.geometry.point.Point.midpoint 1604 1605 Examples 1606 ======== 1607 1608 >>> p1, p2 = Point(0, 0), Point(4, 3) 1609 >>> s1 = Segment(p1, p2) 1610 >>> s1.midpoint 1611 Point(2, 3/2) 1612 1613 """ 1614 return Point.midpoint(self.p1, self.p2) 1615 1616 def distance(self, o): 1617 """ 1618 Finds the shortest distance between a line segment and a point. 1619 1620 Raises 1621 ====== 1622 1623 NotImplementedError 1624 if o is not a Point 1625 1626 Examples 1627 ======== 1628 1629 >>> p1, p2 = Point(0, 1), Point(3, 4) 1630 >>> s = Segment(p1, p2) 1631 >>> s.distance(Point(10, 15)) 1632 sqrt(170) 1633 >>> s.distance((0, 12)) 1634 sqrt(73) 1635 1636 """ 1637 if is_sequence(o): 1638 o = Point(o) 1639 if isinstance(o, Point): 1640 seg_vector = self.p2 - self.p1 1641 pt_vector = o - self.p1 1642 t = seg_vector.dot(pt_vector)/self.length**2 1643 if t >= 1: 1644 distance = Point.distance(self.p2, o) 1645 elif t <= 0: 1646 distance = Point.distance(self.p1, o) 1647 else: 1648 distance = Point.distance( 1649 self.p1 + Point(t*seg_vector.x, t*seg_vector.y), o) 1650 return distance 1651 else: 1652 raise NotImplementedError 1653 1654 def contains(self, other): 1655 """ 1656 Is the other GeometryEntity contained within this Segment? 1657 1658 Examples 1659 ======== 1660 1661 >>> p1, p2 = Point(0, 1), Point(3, 4) 1662 >>> s = Segment(p1, p2) 1663 >>> s2 = Segment(p2, p1) 1664 >>> s.contains(s2) 1665 True 1666 1667 """ 1668 if isinstance(other, Segment): 1669 return other.p1 in self and other.p2 in self 1670 elif isinstance(other, Point): 1671 if Point.is_collinear(self.p1, self.p2, other): 1672 t = Dummy('t') 1673 x, y = self.arbitrary_point(t).args 1674 if self.p1.x != self.p2.x: 1675 ti = solve(x - other.x, t)[0][t] 1676 else: 1677 ti = solve(y - other.y, t)[0][t] 1678 if ti.is_number: 1679 return 0 <= ti <= 1 1680 return 1681 1682 return False 1683