plane.py 26 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887
  1. """Geometrical Planes.
  2. Contains
  3. ========
  4. Plane
  5. """
  6. from sympy.core import Dummy, Rational, S, Symbol
  7. from sympy.core.symbol import _symbol
  8. from sympy.functions.elementary.trigonometric import cos, sin, acos, asin, sqrt
  9. from .entity import GeometryEntity
  10. from .line import (Line, Ray, Segment, Line3D, LinearEntity, LinearEntity3D,
  11. Ray3D, Segment3D)
  12. from .point import Point, Point3D
  13. from sympy.matrices import Matrix
  14. from sympy.polys.polytools import cancel
  15. from sympy.solvers import solve, linsolve
  16. from sympy.utilities.iterables import uniq, is_sequence
  17. from sympy.utilities.misc import filldedent, func_name, Undecidable
  18. from mpmath.libmp.libmpf import prec_to_dps
  19. import random
  20. class Plane(GeometryEntity):
  21. """
  22. A plane is a flat, two-dimensional surface. A plane is the two-dimensional
  23. analogue of a point (zero-dimensions), a line (one-dimension) and a solid
  24. (three-dimensions). A plane can generally be constructed by two types of
  25. inputs. They are three non-collinear points and a point and the plane's
  26. normal vector.
  27. Attributes
  28. ==========
  29. p1
  30. normal_vector
  31. Examples
  32. ========
  33. >>> from sympy import Plane, Point3D
  34. >>> Plane(Point3D(1, 1, 1), Point3D(2, 3, 4), Point3D(2, 2, 2))
  35. Plane(Point3D(1, 1, 1), (-1, 2, -1))
  36. >>> Plane((1, 1, 1), (2, 3, 4), (2, 2, 2))
  37. Plane(Point3D(1, 1, 1), (-1, 2, -1))
  38. >>> Plane(Point3D(1, 1, 1), normal_vector=(1,4,7))
  39. Plane(Point3D(1, 1, 1), (1, 4, 7))
  40. """
  41. def __new__(cls, p1, a=None, b=None, **kwargs):
  42. p1 = Point3D(p1, dim=3)
  43. if a and b:
  44. p2 = Point(a, dim=3)
  45. p3 = Point(b, dim=3)
  46. if Point3D.are_collinear(p1, p2, p3):
  47. raise ValueError('Enter three non-collinear points')
  48. a = p1.direction_ratio(p2)
  49. b = p1.direction_ratio(p3)
  50. normal_vector = tuple(Matrix(a).cross(Matrix(b)))
  51. else:
  52. a = kwargs.pop('normal_vector', a)
  53. evaluate = kwargs.get('evaluate', True)
  54. if is_sequence(a) and len(a) == 3:
  55. normal_vector = Point3D(a).args if evaluate else a
  56. else:
  57. raise ValueError(filldedent('''
  58. Either provide 3 3D points or a point with a
  59. normal vector expressed as a sequence of length 3'''))
  60. if all(coord.is_zero for coord in normal_vector):
  61. raise ValueError('Normal vector cannot be zero vector')
  62. return GeometryEntity.__new__(cls, p1, normal_vector, **kwargs)
  63. def __contains__(self, o):
  64. x, y, z = map(Dummy, 'xyz')
  65. k = self.equation(x, y, z)
  66. if isinstance(o, (LinearEntity, LinearEntity3D)):
  67. t = Dummy()
  68. d = Point3D(o.arbitrary_point(t))
  69. e = k.subs([(x, d.x), (y, d.y), (z, d.z)])
  70. return e.equals(0)
  71. try:
  72. o = Point(o, dim=3, strict=True)
  73. d = k.xreplace(dict(zip((x, y, z), o.args)))
  74. return d.equals(0)
  75. except TypeError:
  76. return False
  77. def _eval_evalf(self, prec=15, **options):
  78. pt, tup = self.args
  79. dps = prec_to_dps(prec)
  80. pt = pt.evalf(n=dps, **options)
  81. tup = tuple([i.evalf(n=dps, **options) for i in tup])
  82. return self.func(pt, normal_vector=tup, evaluate=False)
  83. def angle_between(self, o):
  84. """Angle between the plane and other geometric entity.
  85. Parameters
  86. ==========
  87. LinearEntity3D, Plane.
  88. Returns
  89. =======
  90. angle : angle in radians
  91. Notes
  92. =====
  93. This method accepts only 3D entities as it's parameter, but if you want
  94. to calculate the angle between a 2D entity and a plane you should
  95. first convert to a 3D entity by projecting onto a desired plane and
  96. then proceed to calculate the angle.
  97. Examples
  98. ========
  99. >>> from sympy import Point3D, Line3D, Plane
  100. >>> a = Plane(Point3D(1, 2, 2), normal_vector=(1, 2, 3))
  101. >>> b = Line3D(Point3D(1, 3, 4), Point3D(2, 2, 2))
  102. >>> a.angle_between(b)
  103. -asin(sqrt(21)/6)
  104. """
  105. if isinstance(o, LinearEntity3D):
  106. a = Matrix(self.normal_vector)
  107. b = Matrix(o.direction_ratio)
  108. c = a.dot(b)
  109. d = sqrt(sum([i**2 for i in self.normal_vector]))
  110. e = sqrt(sum([i**2 for i in o.direction_ratio]))
  111. return asin(c/(d*e))
  112. if isinstance(o, Plane):
  113. a = Matrix(self.normal_vector)
  114. b = Matrix(o.normal_vector)
  115. c = a.dot(b)
  116. d = sqrt(sum([i**2 for i in self.normal_vector]))
  117. e = sqrt(sum([i**2 for i in o.normal_vector]))
  118. return acos(c/(d*e))
  119. def arbitrary_point(self, u=None, v=None):
  120. """ Returns an arbitrary point on the Plane. If given two
  121. parameters, the point ranges over the entire plane. If given 1
  122. or no parameters, returns a point with one parameter which,
  123. when varying from 0 to 2*pi, moves the point in a circle of
  124. radius 1 about p1 of the Plane.
  125. Examples
  126. ========
  127. >>> from sympy import Plane, Ray
  128. >>> from sympy.abc import u, v, t, r
  129. >>> p = Plane((1, 1, 1), normal_vector=(1, 0, 0))
  130. >>> p.arbitrary_point(u, v)
  131. Point3D(1, u + 1, v + 1)
  132. >>> p.arbitrary_point(t)
  133. Point3D(1, cos(t) + 1, sin(t) + 1)
  134. While arbitrary values of u and v can move the point anywhere in
  135. the plane, the single-parameter point can be used to construct a
  136. ray whose arbitrary point can be located at angle t and radius
  137. r from p.p1:
  138. >>> Ray(p.p1, _).arbitrary_point(r)
  139. Point3D(1, r*cos(t) + 1, r*sin(t) + 1)
  140. Returns
  141. =======
  142. Point3D
  143. """
  144. circle = v is None
  145. if circle:
  146. u = _symbol(u or 't', real=True)
  147. else:
  148. u = _symbol(u or 'u', real=True)
  149. v = _symbol(v or 'v', real=True)
  150. x, y, z = self.normal_vector
  151. a, b, c = self.p1.args
  152. # x1, y1, z1 is a nonzero vector parallel to the plane
  153. if x.is_zero and y.is_zero:
  154. x1, y1, z1 = S.One, S.Zero, S.Zero
  155. else:
  156. x1, y1, z1 = -y, x, S.Zero
  157. # x2, y2, z2 is also parallel to the plane, and orthogonal to x1, y1, z1
  158. x2, y2, z2 = tuple(Matrix((x, y, z)).cross(Matrix((x1, y1, z1))))
  159. if circle:
  160. x1, y1, z1 = (w/sqrt(x1**2 + y1**2 + z1**2) for w in (x1, y1, z1))
  161. x2, y2, z2 = (w/sqrt(x2**2 + y2**2 + z2**2) for w in (x2, y2, z2))
  162. p = Point3D(a + x1*cos(u) + x2*sin(u), \
  163. b + y1*cos(u) + y2*sin(u), \
  164. c + z1*cos(u) + z2*sin(u))
  165. else:
  166. p = Point3D(a + x1*u + x2*v, b + y1*u + y2*v, c + z1*u + z2*v)
  167. return p
  168. @staticmethod
  169. def are_concurrent(*planes):
  170. """Is a sequence of Planes concurrent?
  171. Two or more Planes are concurrent if their intersections
  172. are a common line.
  173. Parameters
  174. ==========
  175. planes: list
  176. Returns
  177. =======
  178. Boolean
  179. Examples
  180. ========
  181. >>> from sympy import Plane, Point3D
  182. >>> a = Plane(Point3D(5, 0, 0), normal_vector=(1, -1, 1))
  183. >>> b = Plane(Point3D(0, -2, 0), normal_vector=(3, 1, 1))
  184. >>> c = Plane(Point3D(0, -1, 0), normal_vector=(5, -1, 9))
  185. >>> Plane.are_concurrent(a, b)
  186. True
  187. >>> Plane.are_concurrent(a, b, c)
  188. False
  189. """
  190. planes = list(uniq(planes))
  191. for i in planes:
  192. if not isinstance(i, Plane):
  193. raise ValueError('All objects should be Planes but got %s' % i.func)
  194. if len(planes) < 2:
  195. return False
  196. planes = list(planes)
  197. first = planes.pop(0)
  198. sol = first.intersection(planes[0])
  199. if sol == []:
  200. return False
  201. else:
  202. line = sol[0]
  203. for i in planes[1:]:
  204. l = first.intersection(i)
  205. if not l or l[0] not in line:
  206. return False
  207. return True
  208. def distance(self, o):
  209. """Distance between the plane and another geometric entity.
  210. Parameters
  211. ==========
  212. Point3D, LinearEntity3D, Plane.
  213. Returns
  214. =======
  215. distance
  216. Notes
  217. =====
  218. This method accepts only 3D entities as it's parameter, but if you want
  219. to calculate the distance between a 2D entity and a plane you should
  220. first convert to a 3D entity by projecting onto a desired plane and
  221. then proceed to calculate the distance.
  222. Examples
  223. ========
  224. >>> from sympy import Point3D, Line3D, Plane
  225. >>> a = Plane(Point3D(1, 1, 1), normal_vector=(1, 1, 1))
  226. >>> b = Point3D(1, 2, 3)
  227. >>> a.distance(b)
  228. sqrt(3)
  229. >>> c = Line3D(Point3D(2, 3, 1), Point3D(1, 2, 2))
  230. >>> a.distance(c)
  231. 0
  232. """
  233. if self.intersection(o) != []:
  234. return S.Zero
  235. if isinstance(o, (Segment3D, Ray3D)):
  236. a, b = o.p1, o.p2
  237. pi, = self.intersection(Line3D(a, b))
  238. if pi in o:
  239. return self.distance(pi)
  240. elif a in Segment3D(pi, b):
  241. return self.distance(a)
  242. else:
  243. assert isinstance(o, Segment3D) is True
  244. return self.distance(b)
  245. # following code handles `Point3D`, `LinearEntity3D`, `Plane`
  246. a = o if isinstance(o, Point3D) else o.p1
  247. n = Point3D(self.normal_vector).unit
  248. d = (a - self.p1).dot(n)
  249. return abs(d)
  250. def equals(self, o):
  251. """
  252. Returns True if self and o are the same mathematical entities.
  253. Examples
  254. ========
  255. >>> from sympy import Plane, Point3D
  256. >>> a = Plane(Point3D(1, 2, 3), normal_vector=(1, 1, 1))
  257. >>> b = Plane(Point3D(1, 2, 3), normal_vector=(2, 2, 2))
  258. >>> c = Plane(Point3D(1, 2, 3), normal_vector=(-1, 4, 6))
  259. >>> a.equals(a)
  260. True
  261. >>> a.equals(b)
  262. True
  263. >>> a.equals(c)
  264. False
  265. """
  266. if isinstance(o, Plane):
  267. a = self.equation()
  268. b = o.equation()
  269. return cancel(a/b).is_constant()
  270. else:
  271. return False
  272. def equation(self, x=None, y=None, z=None):
  273. """The equation of the Plane.
  274. Examples
  275. ========
  276. >>> from sympy import Point3D, Plane
  277. >>> a = Plane(Point3D(1, 1, 2), Point3D(2, 4, 7), Point3D(3, 5, 1))
  278. >>> a.equation()
  279. -23*x + 11*y - 2*z + 16
  280. >>> a = Plane(Point3D(1, 4, 2), normal_vector=(6, 6, 6))
  281. >>> a.equation()
  282. 6*x + 6*y + 6*z - 42
  283. """
  284. x, y, z = [i if i else Symbol(j, real=True) for i, j in zip((x, y, z), 'xyz')]
  285. a = Point3D(x, y, z)
  286. b = self.p1.direction_ratio(a)
  287. c = self.normal_vector
  288. return (sum(i*j for i, j in zip(b, c)))
  289. def intersection(self, o):
  290. """ The intersection with other geometrical entity.
  291. Parameters
  292. ==========
  293. Point, Point3D, LinearEntity, LinearEntity3D, Plane
  294. Returns
  295. =======
  296. List
  297. Examples
  298. ========
  299. >>> from sympy import Point3D, Line3D, Plane
  300. >>> a = Plane(Point3D(1, 2, 3), normal_vector=(1, 1, 1))
  301. >>> b = Point3D(1, 2, 3)
  302. >>> a.intersection(b)
  303. [Point3D(1, 2, 3)]
  304. >>> c = Line3D(Point3D(1, 4, 7), Point3D(2, 2, 2))
  305. >>> a.intersection(c)
  306. [Point3D(2, 2, 2)]
  307. >>> d = Plane(Point3D(6, 0, 0), normal_vector=(2, -5, 3))
  308. >>> e = Plane(Point3D(2, 0, 0), normal_vector=(3, 4, -3))
  309. >>> d.intersection(e)
  310. [Line3D(Point3D(78/23, -24/23, 0), Point3D(147/23, 321/23, 23))]
  311. """
  312. if not isinstance(o, GeometryEntity):
  313. o = Point(o, dim=3)
  314. if isinstance(o, Point):
  315. if o in self:
  316. return [o]
  317. else:
  318. return []
  319. if isinstance(o, (LinearEntity, LinearEntity3D)):
  320. # recast to 3D
  321. p1, p2 = o.p1, o.p2
  322. if isinstance(o, Segment):
  323. o = Segment3D(p1, p2)
  324. elif isinstance(o, Ray):
  325. o = Ray3D(p1, p2)
  326. elif isinstance(o, Line):
  327. o = Line3D(p1, p2)
  328. else:
  329. raise ValueError('unhandled linear entity: %s' % o.func)
  330. if o in self:
  331. return [o]
  332. else:
  333. t = Dummy() # unnamed else it may clash with a symbol in o
  334. a = Point3D(o.arbitrary_point(t))
  335. p1, n = self.p1, Point3D(self.normal_vector)
  336. # TODO: Replace solve with solveset, when this line is tested
  337. c = solve((a - p1).dot(n), t)
  338. if not c:
  339. return []
  340. else:
  341. c = [i for i in c if i.is_real is not False]
  342. if len(c) > 1:
  343. c = [i for i in c if i.is_real]
  344. if len(c) != 1:
  345. raise Undecidable("not sure which point is real")
  346. p = a.subs(t, c[0])
  347. if p not in o:
  348. return [] # e.g. a segment might not intersect a plane
  349. return [p]
  350. if isinstance(o, Plane):
  351. if self.equals(o):
  352. return [self]
  353. if self.is_parallel(o):
  354. return []
  355. else:
  356. x, y, z = map(Dummy, 'xyz')
  357. a, b = Matrix([self.normal_vector]), Matrix([o.normal_vector])
  358. c = list(a.cross(b))
  359. d = self.equation(x, y, z)
  360. e = o.equation(x, y, z)
  361. result = list(linsolve([d, e], x, y, z))[0]
  362. for i in (x, y, z): result = result.subs(i, 0)
  363. return [Line3D(Point3D(result), direction_ratio=c)]
  364. def is_coplanar(self, o):
  365. """ Returns True if `o` is coplanar with self, else False.
  366. Examples
  367. ========
  368. >>> from sympy import Plane
  369. >>> o = (0, 0, 0)
  370. >>> p = Plane(o, (1, 1, 1))
  371. >>> p2 = Plane(o, (2, 2, 2))
  372. >>> p == p2
  373. False
  374. >>> p.is_coplanar(p2)
  375. True
  376. """
  377. if isinstance(o, Plane):
  378. x, y, z = map(Dummy, 'xyz')
  379. return not cancel(self.equation(x, y, z)/o.equation(x, y, z)).has(x, y, z)
  380. if isinstance(o, Point3D):
  381. return o in self
  382. elif isinstance(o, LinearEntity3D):
  383. return all(i in self for i in self)
  384. elif isinstance(o, GeometryEntity): # XXX should only be handling 2D objects now
  385. return all(i == 0 for i in self.normal_vector[:2])
  386. def is_parallel(self, l):
  387. """Is the given geometric entity parallel to the plane?
  388. Parameters
  389. ==========
  390. LinearEntity3D or Plane
  391. Returns
  392. =======
  393. Boolean
  394. Examples
  395. ========
  396. >>> from sympy import Plane, Point3D
  397. >>> a = Plane(Point3D(1,4,6), normal_vector=(2, 4, 6))
  398. >>> b = Plane(Point3D(3,1,3), normal_vector=(4, 8, 12))
  399. >>> a.is_parallel(b)
  400. True
  401. """
  402. if isinstance(l, LinearEntity3D):
  403. a = l.direction_ratio
  404. b = self.normal_vector
  405. c = sum([i*j for i, j in zip(a, b)])
  406. if c == 0:
  407. return True
  408. else:
  409. return False
  410. elif isinstance(l, Plane):
  411. a = Matrix(l.normal_vector)
  412. b = Matrix(self.normal_vector)
  413. if a.cross(b).is_zero_matrix:
  414. return True
  415. else:
  416. return False
  417. def is_perpendicular(self, l):
  418. """is the given geometric entity perpendicualar to the given plane?
  419. Parameters
  420. ==========
  421. LinearEntity3D or Plane
  422. Returns
  423. =======
  424. Boolean
  425. Examples
  426. ========
  427. >>> from sympy import Plane, Point3D
  428. >>> a = Plane(Point3D(1,4,6), normal_vector=(2, 4, 6))
  429. >>> b = Plane(Point3D(2, 2, 2), normal_vector=(-1, 2, -1))
  430. >>> a.is_perpendicular(b)
  431. True
  432. """
  433. if isinstance(l, LinearEntity3D):
  434. a = Matrix(l.direction_ratio)
  435. b = Matrix(self.normal_vector)
  436. if a.cross(b).is_zero_matrix:
  437. return True
  438. else:
  439. return False
  440. elif isinstance(l, Plane):
  441. a = Matrix(l.normal_vector)
  442. b = Matrix(self.normal_vector)
  443. if a.dot(b) == 0:
  444. return True
  445. else:
  446. return False
  447. else:
  448. return False
  449. @property
  450. def normal_vector(self):
  451. """Normal vector of the given plane.
  452. Examples
  453. ========
  454. >>> from sympy import Point3D, Plane
  455. >>> a = Plane(Point3D(1, 1, 1), Point3D(2, 3, 4), Point3D(2, 2, 2))
  456. >>> a.normal_vector
  457. (-1, 2, -1)
  458. >>> a = Plane(Point3D(1, 1, 1), normal_vector=(1, 4, 7))
  459. >>> a.normal_vector
  460. (1, 4, 7)
  461. """
  462. return self.args[1]
  463. @property
  464. def p1(self):
  465. """The only defining point of the plane. Others can be obtained from the
  466. arbitrary_point method.
  467. See Also
  468. ========
  469. sympy.geometry.point.Point3D
  470. Examples
  471. ========
  472. >>> from sympy import Point3D, Plane
  473. >>> a = Plane(Point3D(1, 1, 1), Point3D(2, 3, 4), Point3D(2, 2, 2))
  474. >>> a.p1
  475. Point3D(1, 1, 1)
  476. """
  477. return self.args[0]
  478. def parallel_plane(self, pt):
  479. """
  480. Plane parallel to the given plane and passing through the point pt.
  481. Parameters
  482. ==========
  483. pt: Point3D
  484. Returns
  485. =======
  486. Plane
  487. Examples
  488. ========
  489. >>> from sympy import Plane, Point3D
  490. >>> a = Plane(Point3D(1, 4, 6), normal_vector=(2, 4, 6))
  491. >>> a.parallel_plane(Point3D(2, 3, 5))
  492. Plane(Point3D(2, 3, 5), (2, 4, 6))
  493. """
  494. a = self.normal_vector
  495. return Plane(pt, normal_vector=a)
  496. def perpendicular_line(self, pt):
  497. """A line perpendicular to the given plane.
  498. Parameters
  499. ==========
  500. pt: Point3D
  501. Returns
  502. =======
  503. Line3D
  504. Examples
  505. ========
  506. >>> from sympy import Plane, Point3D
  507. >>> a = Plane(Point3D(1,4,6), normal_vector=(2, 4, 6))
  508. >>> a.perpendicular_line(Point3D(9, 8, 7))
  509. Line3D(Point3D(9, 8, 7), Point3D(11, 12, 13))
  510. """
  511. a = self.normal_vector
  512. return Line3D(pt, direction_ratio=a)
  513. def perpendicular_plane(self, *pts):
  514. """
  515. Return a perpendicular passing through the given points. If the
  516. direction ratio between the points is the same as the Plane's normal
  517. vector then, to select from the infinite number of possible planes,
  518. a third point will be chosen on the z-axis (or the y-axis
  519. if the normal vector is already parallel to the z-axis). If less than
  520. two points are given they will be supplied as follows: if no point is
  521. given then pt1 will be self.p1; if a second point is not given it will
  522. be a point through pt1 on a line parallel to the z-axis (if the normal
  523. is not already the z-axis, otherwise on the line parallel to the
  524. y-axis).
  525. Parameters
  526. ==========
  527. pts: 0, 1 or 2 Point3D
  528. Returns
  529. =======
  530. Plane
  531. Examples
  532. ========
  533. >>> from sympy import Plane, Point3D
  534. >>> a, b = Point3D(0, 0, 0), Point3D(0, 1, 0)
  535. >>> Z = (0, 0, 1)
  536. >>> p = Plane(a, normal_vector=Z)
  537. >>> p.perpendicular_plane(a, b)
  538. Plane(Point3D(0, 0, 0), (1, 0, 0))
  539. """
  540. if len(pts) > 2:
  541. raise ValueError('No more than 2 pts should be provided.')
  542. pts = list(pts)
  543. if len(pts) == 0:
  544. pts.append(self.p1)
  545. if len(pts) == 1:
  546. x, y, z = self.normal_vector
  547. if x == y == 0:
  548. dir = (0, 1, 0)
  549. else:
  550. dir = (0, 0, 1)
  551. pts.append(pts[0] + Point3D(*dir))
  552. p1, p2 = [Point(i, dim=3) for i in pts]
  553. l = Line3D(p1, p2)
  554. n = Line3D(p1, direction_ratio=self.normal_vector)
  555. if l in n: # XXX should an error be raised instead?
  556. # there are infinitely many perpendicular planes;
  557. x, y, z = self.normal_vector
  558. if x == y == 0:
  559. # the z axis is the normal so pick a pt on the y-axis
  560. p3 = Point3D(0, 1, 0) # case 1
  561. else:
  562. # else pick a pt on the z axis
  563. p3 = Point3D(0, 0, 1) # case 2
  564. # in case that point is already given, move it a bit
  565. if p3 in l:
  566. p3 *= 2 # case 3
  567. else:
  568. p3 = p1 + Point3D(*self.normal_vector) # case 4
  569. return Plane(p1, p2, p3)
  570. def projection_line(self, line):
  571. """Project the given line onto the plane through the normal plane
  572. containing the line.
  573. Parameters
  574. ==========
  575. LinearEntity or LinearEntity3D
  576. Returns
  577. =======
  578. Point3D, Line3D, Ray3D or Segment3D
  579. Notes
  580. =====
  581. For the interaction between 2D and 3D lines(segments, rays), you should
  582. convert the line to 3D by using this method. For example for finding the
  583. intersection between a 2D and a 3D line, convert the 2D line to a 3D line
  584. by projecting it on a required plane and then proceed to find the
  585. intersection between those lines.
  586. Examples
  587. ========
  588. >>> from sympy import Plane, Line, Line3D, Point3D
  589. >>> a = Plane(Point3D(1, 1, 1), normal_vector=(1, 1, 1))
  590. >>> b = Line(Point3D(1, 1), Point3D(2, 2))
  591. >>> a.projection_line(b)
  592. Line3D(Point3D(4/3, 4/3, 1/3), Point3D(5/3, 5/3, -1/3))
  593. >>> c = Line3D(Point3D(1, 1, 1), Point3D(2, 2, 2))
  594. >>> a.projection_line(c)
  595. Point3D(1, 1, 1)
  596. """
  597. if not isinstance(line, (LinearEntity, LinearEntity3D)):
  598. raise NotImplementedError('Enter a linear entity only')
  599. a, b = self.projection(line.p1), self.projection(line.p2)
  600. if a == b:
  601. # projection does not imply intersection so for
  602. # this case (line parallel to plane's normal) we
  603. # return the projection point
  604. return a
  605. if isinstance(line, (Line, Line3D)):
  606. return Line3D(a, b)
  607. if isinstance(line, (Ray, Ray3D)):
  608. return Ray3D(a, b)
  609. if isinstance(line, (Segment, Segment3D)):
  610. return Segment3D(a, b)
  611. def projection(self, pt):
  612. """Project the given point onto the plane along the plane normal.
  613. Parameters
  614. ==========
  615. Point or Point3D
  616. Returns
  617. =======
  618. Point3D
  619. Examples
  620. ========
  621. >>> from sympy import Plane, Point3D
  622. >>> A = Plane(Point3D(1, 1, 2), normal_vector=(1, 1, 1))
  623. The projection is along the normal vector direction, not the z
  624. axis, so (1, 1) does not project to (1, 1, 2) on the plane A:
  625. >>> b = Point3D(1, 1)
  626. >>> A.projection(b)
  627. Point3D(5/3, 5/3, 2/3)
  628. >>> _ in A
  629. True
  630. But the point (1, 1, 2) projects to (1, 1) on the XY-plane:
  631. >>> XY = Plane((0, 0, 0), (0, 0, 1))
  632. >>> XY.projection((1, 1, 2))
  633. Point3D(1, 1, 0)
  634. """
  635. rv = Point(pt, dim=3)
  636. if rv in self:
  637. return rv
  638. return self.intersection(Line3D(rv, rv + Point3D(self.normal_vector)))[0]
  639. def random_point(self, seed=None):
  640. """ Returns a random point on the Plane.
  641. Returns
  642. =======
  643. Point3D
  644. Examples
  645. ========
  646. >>> from sympy import Plane
  647. >>> p = Plane((1, 0, 0), normal_vector=(0, 1, 0))
  648. >>> r = p.random_point(seed=42) # seed value is optional
  649. >>> r.n(3)
  650. Point3D(2.29, 0, -1.35)
  651. The random point can be moved to lie on the circle of radius
  652. 1 centered on p1:
  653. >>> c = p.p1 + (r - p.p1).unit
  654. >>> c.distance(p.p1).equals(1)
  655. True
  656. """
  657. if seed is not None:
  658. rng = random.Random(seed)
  659. else:
  660. rng = random
  661. u, v = Dummy('u'), Dummy('v')
  662. params = {
  663. u: 2*Rational(rng.gauss(0, 1)) - 1,
  664. v: 2*Rational(rng.gauss(0, 1)) - 1}
  665. return self.arbitrary_point(u, v).subs(params)
  666. def parameter_value(self, other, u, v=None):
  667. """Return the parameter(s) corresponding to the given point.
  668. Examples
  669. ========
  670. >>> from sympy import pi, Plane
  671. >>> from sympy.abc import t, u, v
  672. >>> p = Plane((2, 0, 0), (0, 0, 1), (0, 1, 0))
  673. By default, the parameter value returned defines a point
  674. that is a distance of 1 from the Plane's p1 value and
  675. in line with the given point:
  676. >>> on_circle = p.arbitrary_point(t).subs(t, pi/4)
  677. >>> on_circle.distance(p.p1)
  678. 1
  679. >>> p.parameter_value(on_circle, t)
  680. {t: pi/4}
  681. Moving the point twice as far from p1 does not change
  682. the parameter value:
  683. >>> off_circle = p.p1 + (on_circle - p.p1)*2
  684. >>> off_circle.distance(p.p1)
  685. 2
  686. >>> p.parameter_value(off_circle, t)
  687. {t: pi/4}
  688. If the 2-value parameter is desired, supply the two
  689. parameter symbols and a replacement dictionary will
  690. be returned:
  691. >>> p.parameter_value(on_circle, u, v)
  692. {u: sqrt(10)/10, v: sqrt(10)/30}
  693. >>> p.parameter_value(off_circle, u, v)
  694. {u: sqrt(10)/5, v: sqrt(10)/15}
  695. """
  696. from sympy.geometry.point import Point
  697. if not isinstance(other, GeometryEntity):
  698. other = Point(other, dim=self.ambient_dimension)
  699. if not isinstance(other, Point):
  700. raise ValueError("other must be a point")
  701. if other == self.p1:
  702. return other
  703. if isinstance(u, Symbol) and v is None:
  704. delta = self.arbitrary_point(u) - self.p1
  705. eq = delta - (other - self.p1).unit
  706. sol = solve(eq, u, dict=True)
  707. elif isinstance(u, Symbol) and isinstance(v, Symbol):
  708. pt = self.arbitrary_point(u, v)
  709. sol = solve(pt - other, (u, v), dict=True)
  710. else:
  711. raise ValueError('expecting 1 or 2 symbols')
  712. if not sol:
  713. raise ValueError("Given point is not on %s" % func_name(self))
  714. return sol[0] # {t: tval} or {u: uval, v: vval}
  715. @property
  716. def ambient_dimension(self):
  717. return self.p1.ambient_dimension