sympy/geometry/point.py

   1 from sympy.core.basic import Basic, S, sympify
   2 from sympy.simplify import simplify
   3 from sympy.geometry.exceptions import GeometryError
   4 from sympy.functions.elementary.miscellaneous import sqrt
   5 from entity import GeometryEntity
   6
   7
   8 class Point(GeometryEntity):
   9     """
  10     A point in Euclidean N-space defined by a sequence of values. Can be
  11     constructed from a sequence of points or a list of points.
  12
  13     Examples:
  14     ======
  15         >>> Point(1, 2)
  16         Point(1, 2)
  17         >>> Point([1, 2])
  18         Point(1, 2)
  19
  20     Notes:
  21     ======
  22         - Currently only 2-dimensional points are supported.
  23     """
  24     def __new__(cls, *args, **kwargs):
  25         if isinstance(args[0], (tuple, list, set)):
  26             coords = tuple([sympify(x) for x in args[0]])
  27         else:
  28             coords = tuple([sympify(x) for x in args])
  29
  30         if len(coords) != 2:
  31             raise NotImplementedError("Only two dimensional points currently supported")
  32
  33         return GeometryEntity.__new__(cls, *coords)
  34
  35     def is_collinear(*points):
  36         """
  37         Test whether or not a set of points are collinear. Returns True if
  38         the set of points are collinear, or False otherwise.
  39
  40         Examples:
  41         =========
  42             >>> from sympy import *
  43             >>> x = Symbol('x')
  44             >>> p1,p2 = Point(0, 0), Point(1, 1)
  45             >>> p3,p4,p5 = Point(2, 2), Point(x, x), Point(1, 2)
  46             >>> Point.is_collinear(p1, p2, p3, p4)
  47             True
  48             >>> Point.is_collinear(p1, p2, p3, p5)
  49             False
  50
  51         Description of method used:
  52         ===========================
  53             Slope is preserved everywhere on a line, so the slope between
  54             any two points on the line should be the same. Take the first
  55             two points, p1 and p2, and create a translated point v1
  56             with p1 as the origin. Now for every other point we create
  57             a translated point, vi with p1 also as the origin. Note that
  58             these translations preserve slope since everything is
  59             consistently translated to a new origin of p1. Since slope
  60             is preserved then we have the following equality:
  61                     v1_slope = vi_slope
  62               =>    v1.y/v1.x = vi.y/vi.x (due to translation)
  63               =>    v1.y*vi.x = vi.y*v1.x
  64               =>    v1.y*vi.x - vi.y*v1.x = 0           (*)
  65             Hence, if we have a vi such that the equality in (*) is False
  66             then the points are not collinear. We do this test for every
  67             point in the list, and if all pass then they are collinear.
  68         """
  69         points = GeometryEntity.extract_entities(points)
  70         if len(points) == 0: return False
  71         if len(points) <= 2: return True # two points always form a line
  72
  73         # XXX Cross product is used now, but that only extends to three
  74         #     dimensions. If the concept needs to extend to greater
  75         #     dimensions then another method would have to be used
  76         p1 = points[0]
  77         p2 = points[1]
  78         v1 = p2 - p1
  79         for p3 in points[2:]:
  80             v2 = p3 - p1
  81             test = simplify(v1[0]*v2[1] - v1[1]*v2[0])
  82             if simplify(test) != 0:
  83                 return False
  84         return True
  85
  86     def is_concyclic(*points):
  87         """
  88         Test whether or not a set of points are concyclic (i.e., on the same
  89         circle). Returns True if they are concyclic, or False otherwise.
  90
  91         Example:
  92         ========
  93             >>> p1,p2 = Point(-1, 0), Point(1, 0)
  94             >>> p3,p4 = Point(0, 1), Point(-1, 2)
  95             >>> Point.is_concyclic(p1, p2, p3)
  96             True
  97             >>> Point.is_concyclic(p1, p2, p3, p4)
  98             False
  99
 100         Description of method used:
 101         ===========================
 102             No points are not considered to be concyclic. One or two points
 103             are definitely concyclic and three points are conyclic iff they
 104             are not collinear.
 105
 106             For more than three points, we pick the first three points and
 107             attempt to create a circle. If the circle cannot be created
 108             (i.e., they are collinear) then all of the points cannot be
 109             concyclic. If the circle is created successfully then simply
 110             check all of the other points for containment in the circle.
 111         """
 112         points = GeometryEntity.extract_entities(points)
 113         if len(points) == 0: return False
 114         if len(points) <= 2: return True
 115         if len(points) == 3: return (not Point.is_collinear(*points))
 116
 117         try:
 118             from ellipse import Circle
 119             c = Circle(points[0], points[1], points[2])
 120             for point in points[3:]:
 121                 if point not in c:
 122                     return False
 123             return True
 124         except GeometryError,e:
 125             # Circle could not be created, because of collinearity of the
 126             # three points passed in, hence they are not concyclic.
 127             return False
 128         """
 129         # This code is from Maple
 130         def f(u):
 131             dd = u[0]**2 + u[1]**2 + 1
 132             u1 = 2*u[0] / dd
 133             u2 = 2*u[1] / dd
 134             u3 = (dd - 2) / dd
 135             return u1,u2,u3
 136
 137         u1,u2,u3 = f(points[0])
 138         v1,v2,v3 = f(points[1])
 139         w1,w2,w3 = f(points[2])
 140         p = [v1 - u1, v2 - u2, v3 - u3]
 141         q = [w1 - u1, w2 - u2, w3 - u3]
 142         r = [p[1]*q[2] - p[2]*q[1], p[2]*q[0] - p[0]*q[2], p[0]*q[1] - p[1]*q[0]]
 143         for ind in xrange(3, len(points)):
 144             s1,s2,s3 = f(points[ind])
 145             test = simplify(r[0]*(s1-u1) + r[1]*(s2-u2) + r[2]*(s3-u3))
 146             if test != 0:
 147                 return False
 148         return True
 149         """
 150
 151     @staticmethod
 152     def distance(p1, p2):
 153         """
 154         Get the Euclidean distance between two points.
 155
 156         Example:
 157         ========
 158             >>> p1,p2 = Point(1, 1), Point(4, 5)
 159             >>> Point.distance(p1, p2)
 160             5
 161         """
 162         return sqrt( sum([(a-b)**2 for a,b in zip(p1,p2)]) )
 163
 164     @staticmethod
 165     def midpoint(p1, p2):
 166         """
 167         Get the midpoint of two points.
 168
 169         Example:
 170         ========
 171             >>> p1,p2 = Point(1, 1), Point(13, 5)
 172             >>> Point.midpoint(p1, p2)
 173             Point(7, 3)
 174         """
 175         return Point( [simplify((a + b)*S.Half) for a,b in zip(p1,p2)] )
 176
 177     def evalf(self):
 178         """
 179         Evaluate and return a Point where every coordinate is evaluated to
 180         a floating point number.
 181
 182         Example:
 183         ========
 184             >>> from sympy import *
 185             >>> Point(Rational(1,2), Rational(3,2)).evalf()
 186             Point(0.5, 1.5)
 187         """
 188         return Point([x.evalf() for x in self])
 189
 190     def intersection(self, o):
 191         GeometryEntity.intersection.__doc__
 192         if isinstance(o, Point):
 193             if self == o:
 194                 return [self]
 195             return []
 196         raise NotImplementedError()
 197
 198     def __add__(self, other):
 199         """
 200         Create a new point where each coordinate in this point is
 201         increased by the corresponding coordinate in other.
 202         """
 203         if isinstance(other, Point):
 204             if len(other) == len(self):
 205                 return Point( [simplify(a+b) for a,b in zip(self, other)] )
 206             else:
 207                 raise Exception("Points must have the same number of dimensions")
 208         else:
 209             other = sympify(other)
 210             return Point( [simplify(a+other) for a in self] )
 211
 212     def __sub__(self, other):
 213         """
 214         Create a new point where each coordinate in this point is
 215         decreased by the corresponding coordinate in other.
 216         """
 217         return self + (-other)
 218
 219     def __mul__(self, factor):
 220         """
 221         Create a new point where each coordinate in this point is
 222         multiplied by factor.
 223         """
 224         factor = sympify(factor)
 225         return Point( [x*factor for x in self] )
 226
 227     def __div__(self, divisor):
 228         """
 229         Create a new point where each coordinate in this point is
 230         divided by factor.
 231         """
 232         divisor = sympify(divisor)
 233         return Point( [x/divisor for x in self] )
 234
 235     def __neg__(self):
 236         """
 237         Create a new point where each oordinate in this point is negated.
 238         """
 239         return Point( [-x for x in self] )
 240
 241     def __abs__(self):
 242         """Returns the distance between this point and the origin."""
 243         origin = Point([0] * len(self))
 244         return Point.distance(origin, self)