plane.py 26 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878
  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. x, y, z, t = [Dummy('plane_dummy') for i in range(4)]
  21. class Plane(GeometryEntity):
  22. """
  23. A plane is a flat, two-dimensional surface. A plane is the two-dimensional
  24. analogue of a point (zero-dimensions), a line (one-dimension) and a solid
  25. (three-dimensions). A plane can generally be constructed by two types of
  26. inputs. They are:
  27. - three non-collinear points
  28. - a point and the plane's normal vector
  29. Attributes
  30. ==========
  31. p1
  32. normal_vector
  33. Examples
  34. ========
  35. >>> from sympy import Plane, Point3D
  36. >>> Plane(Point3D(1, 1, 1), Point3D(2, 3, 4), Point3D(2, 2, 2))
  37. Plane(Point3D(1, 1, 1), (-1, 2, -1))
  38. >>> Plane((1, 1, 1), (2, 3, 4), (2, 2, 2))
  39. Plane(Point3D(1, 1, 1), (-1, 2, -1))
  40. >>> Plane(Point3D(1, 1, 1), normal_vector=(1,4,7))
  41. Plane(Point3D(1, 1, 1), (1, 4, 7))
  42. """
  43. def __new__(cls, p1, a=None, b=None, **kwargs):
  44. p1 = Point3D(p1, dim=3)
  45. if a and b:
  46. p2 = Point(a, dim=3)
  47. p3 = Point(b, dim=3)
  48. if Point3D.are_collinear(p1, p2, p3):
  49. raise ValueError('Enter three non-collinear points')
  50. a = p1.direction_ratio(p2)
  51. b = p1.direction_ratio(p3)
  52. normal_vector = tuple(Matrix(a).cross(Matrix(b)))
  53. else:
  54. a = kwargs.pop('normal_vector', a)
  55. evaluate = kwargs.get('evaluate', True)
  56. if is_sequence(a) and len(a) == 3:
  57. normal_vector = Point3D(a).args if evaluate else a
  58. else:
  59. raise ValueError(filldedent('''
  60. Either provide 3 3D points or a point with a
  61. normal vector expressed as a sequence of length 3'''))
  62. if all(coord.is_zero for coord in normal_vector):
  63. raise ValueError('Normal vector cannot be zero vector')
  64. return GeometryEntity.__new__(cls, p1, normal_vector, **kwargs)
  65. def __contains__(self, o):
  66. k = self.equation(x, y, z)
  67. if isinstance(o, (LinearEntity, LinearEntity3D)):
  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. a = Point3D(o.arbitrary_point(t))
  334. p1, n = self.p1, Point3D(self.normal_vector)
  335. # TODO: Replace solve with solveset, when this line is tested
  336. c = solve((a - p1).dot(n), t)
  337. if not c:
  338. return []
  339. else:
  340. c = [i for i in c if i.is_real is not False]
  341. if len(c) > 1:
  342. c = [i for i in c if i.is_real]
  343. if len(c) != 1:
  344. raise Undecidable("not sure which point is real")
  345. p = a.subs(t, c[0])
  346. if p not in o:
  347. return [] # e.g. a segment might not intersect a plane
  348. return [p]
  349. if isinstance(o, Plane):
  350. if self.equals(o):
  351. return [self]
  352. if self.is_parallel(o):
  353. return []
  354. else:
  355. x, y, z = map(Dummy, 'xyz')
  356. a, b = Matrix([self.normal_vector]), Matrix([o.normal_vector])
  357. c = list(a.cross(b))
  358. d = self.equation(x, y, z)
  359. e = o.equation(x, y, z)
  360. result = list(linsolve([d, e], x, y, z))[0]
  361. for i in (x, y, z): result = result.subs(i, 0)
  362. return [Line3D(Point3D(result), direction_ratio=c)]
  363. def is_coplanar(self, o):
  364. """ Returns True if `o` is coplanar with self, else False.
  365. Examples
  366. ========
  367. >>> from sympy import Plane
  368. >>> o = (0, 0, 0)
  369. >>> p = Plane(o, (1, 1, 1))
  370. >>> p2 = Plane(o, (2, 2, 2))
  371. >>> p == p2
  372. False
  373. >>> p.is_coplanar(p2)
  374. True
  375. """
  376. if isinstance(o, Plane):
  377. return not cancel(self.equation(x, y, z)/o.equation(x, y, z)).has(x, y, z)
  378. if isinstance(o, Point3D):
  379. return o in self
  380. elif isinstance(o, LinearEntity3D):
  381. return all(i in self for i in self)
  382. elif isinstance(o, GeometryEntity): # XXX should only be handling 2D objects now
  383. return all(i == 0 for i in self.normal_vector[:2])
  384. def is_parallel(self, l):
  385. """Is the given geometric entity parallel to the plane?
  386. Parameters
  387. ==========
  388. LinearEntity3D or Plane
  389. Returns
  390. =======
  391. Boolean
  392. Examples
  393. ========
  394. >>> from sympy import Plane, Point3D
  395. >>> a = Plane(Point3D(1,4,6), normal_vector=(2, 4, 6))
  396. >>> b = Plane(Point3D(3,1,3), normal_vector=(4, 8, 12))
  397. >>> a.is_parallel(b)
  398. True
  399. """
  400. if isinstance(l, LinearEntity3D):
  401. a = l.direction_ratio
  402. b = self.normal_vector
  403. return sum(i*j for i, j in zip(a, b)) == 0
  404. if isinstance(l, Plane):
  405. a = Matrix(l.normal_vector)
  406. b = Matrix(self.normal_vector)
  407. return bool(a.cross(b).is_zero_matrix)
  408. def is_perpendicular(self, l):
  409. """Is the given geometric entity perpendicualar to the given plane?
  410. Parameters
  411. ==========
  412. LinearEntity3D or Plane
  413. Returns
  414. =======
  415. Boolean
  416. Examples
  417. ========
  418. >>> from sympy import Plane, Point3D
  419. >>> a = Plane(Point3D(1,4,6), normal_vector=(2, 4, 6))
  420. >>> b = Plane(Point3D(2, 2, 2), normal_vector=(-1, 2, -1))
  421. >>> a.is_perpendicular(b)
  422. True
  423. """
  424. if isinstance(l, LinearEntity3D):
  425. a = Matrix(l.direction_ratio)
  426. b = Matrix(self.normal_vector)
  427. if a.cross(b).is_zero_matrix:
  428. return True
  429. else:
  430. return False
  431. elif isinstance(l, Plane):
  432. a = Matrix(l.normal_vector)
  433. b = Matrix(self.normal_vector)
  434. if a.dot(b) == 0:
  435. return True
  436. else:
  437. return False
  438. else:
  439. return False
  440. @property
  441. def normal_vector(self):
  442. """Normal vector of the given plane.
  443. Examples
  444. ========
  445. >>> from sympy import Point3D, Plane
  446. >>> a = Plane(Point3D(1, 1, 1), Point3D(2, 3, 4), Point3D(2, 2, 2))
  447. >>> a.normal_vector
  448. (-1, 2, -1)
  449. >>> a = Plane(Point3D(1, 1, 1), normal_vector=(1, 4, 7))
  450. >>> a.normal_vector
  451. (1, 4, 7)
  452. """
  453. return self.args[1]
  454. @property
  455. def p1(self):
  456. """The only defining point of the plane. Others can be obtained from the
  457. arbitrary_point method.
  458. See Also
  459. ========
  460. sympy.geometry.point.Point3D
  461. Examples
  462. ========
  463. >>> from sympy import Point3D, Plane
  464. >>> a = Plane(Point3D(1, 1, 1), Point3D(2, 3, 4), Point3D(2, 2, 2))
  465. >>> a.p1
  466. Point3D(1, 1, 1)
  467. """
  468. return self.args[0]
  469. def parallel_plane(self, pt):
  470. """
  471. Plane parallel to the given plane and passing through the point pt.
  472. Parameters
  473. ==========
  474. pt: Point3D
  475. Returns
  476. =======
  477. Plane
  478. Examples
  479. ========
  480. >>> from sympy import Plane, Point3D
  481. >>> a = Plane(Point3D(1, 4, 6), normal_vector=(2, 4, 6))
  482. >>> a.parallel_plane(Point3D(2, 3, 5))
  483. Plane(Point3D(2, 3, 5), (2, 4, 6))
  484. """
  485. a = self.normal_vector
  486. return Plane(pt, normal_vector=a)
  487. def perpendicular_line(self, pt):
  488. """A line perpendicular to the given plane.
  489. Parameters
  490. ==========
  491. pt: Point3D
  492. Returns
  493. =======
  494. Line3D
  495. Examples
  496. ========
  497. >>> from sympy import Plane, Point3D
  498. >>> a = Plane(Point3D(1,4,6), normal_vector=(2, 4, 6))
  499. >>> a.perpendicular_line(Point3D(9, 8, 7))
  500. Line3D(Point3D(9, 8, 7), Point3D(11, 12, 13))
  501. """
  502. a = self.normal_vector
  503. return Line3D(pt, direction_ratio=a)
  504. def perpendicular_plane(self, *pts):
  505. """
  506. Return a perpendicular passing through the given points. If the
  507. direction ratio between the points is the same as the Plane's normal
  508. vector then, to select from the infinite number of possible planes,
  509. a third point will be chosen on the z-axis (or the y-axis
  510. if the normal vector is already parallel to the z-axis). If less than
  511. two points are given they will be supplied as follows: if no point is
  512. given then pt1 will be self.p1; if a second point is not given it will
  513. be a point through pt1 on a line parallel to the z-axis (if the normal
  514. is not already the z-axis, otherwise on the line parallel to the
  515. y-axis).
  516. Parameters
  517. ==========
  518. pts: 0, 1 or 2 Point3D
  519. Returns
  520. =======
  521. Plane
  522. Examples
  523. ========
  524. >>> from sympy import Plane, Point3D
  525. >>> a, b = Point3D(0, 0, 0), Point3D(0, 1, 0)
  526. >>> Z = (0, 0, 1)
  527. >>> p = Plane(a, normal_vector=Z)
  528. >>> p.perpendicular_plane(a, b)
  529. Plane(Point3D(0, 0, 0), (1, 0, 0))
  530. """
  531. if len(pts) > 2:
  532. raise ValueError('No more than 2 pts should be provided.')
  533. pts = list(pts)
  534. if len(pts) == 0:
  535. pts.append(self.p1)
  536. if len(pts) == 1:
  537. x, y, z = self.normal_vector
  538. if x == y == 0:
  539. dir = (0, 1, 0)
  540. else:
  541. dir = (0, 0, 1)
  542. pts.append(pts[0] + Point3D(*dir))
  543. p1, p2 = [Point(i, dim=3) for i in pts]
  544. l = Line3D(p1, p2)
  545. n = Line3D(p1, direction_ratio=self.normal_vector)
  546. if l in n: # XXX should an error be raised instead?
  547. # there are infinitely many perpendicular planes;
  548. x, y, z = self.normal_vector
  549. if x == y == 0:
  550. # the z axis is the normal so pick a pt on the y-axis
  551. p3 = Point3D(0, 1, 0) # case 1
  552. else:
  553. # else pick a pt on the z axis
  554. p3 = Point3D(0, 0, 1) # case 2
  555. # in case that point is already given, move it a bit
  556. if p3 in l:
  557. p3 *= 2 # case 3
  558. else:
  559. p3 = p1 + Point3D(*self.normal_vector) # case 4
  560. return Plane(p1, p2, p3)
  561. def projection_line(self, line):
  562. """Project the given line onto the plane through the normal plane
  563. containing the line.
  564. Parameters
  565. ==========
  566. LinearEntity or LinearEntity3D
  567. Returns
  568. =======
  569. Point3D, Line3D, Ray3D or Segment3D
  570. Notes
  571. =====
  572. For the interaction between 2D and 3D lines(segments, rays), you should
  573. convert the line to 3D by using this method. For example for finding the
  574. intersection between a 2D and a 3D line, convert the 2D line to a 3D line
  575. by projecting it on a required plane and then proceed to find the
  576. intersection between those lines.
  577. Examples
  578. ========
  579. >>> from sympy import Plane, Line, Line3D, Point3D
  580. >>> a = Plane(Point3D(1, 1, 1), normal_vector=(1, 1, 1))
  581. >>> b = Line(Point3D(1, 1), Point3D(2, 2))
  582. >>> a.projection_line(b)
  583. Line3D(Point3D(4/3, 4/3, 1/3), Point3D(5/3, 5/3, -1/3))
  584. >>> c = Line3D(Point3D(1, 1, 1), Point3D(2, 2, 2))
  585. >>> a.projection_line(c)
  586. Point3D(1, 1, 1)
  587. """
  588. if not isinstance(line, (LinearEntity, LinearEntity3D)):
  589. raise NotImplementedError('Enter a linear entity only')
  590. a, b = self.projection(line.p1), self.projection(line.p2)
  591. if a == b:
  592. # projection does not imply intersection so for
  593. # this case (line parallel to plane's normal) we
  594. # return the projection point
  595. return a
  596. if isinstance(line, (Line, Line3D)):
  597. return Line3D(a, b)
  598. if isinstance(line, (Ray, Ray3D)):
  599. return Ray3D(a, b)
  600. if isinstance(line, (Segment, Segment3D)):
  601. return Segment3D(a, b)
  602. def projection(self, pt):
  603. """Project the given point onto the plane along the plane normal.
  604. Parameters
  605. ==========
  606. Point or Point3D
  607. Returns
  608. =======
  609. Point3D
  610. Examples
  611. ========
  612. >>> from sympy import Plane, Point3D
  613. >>> A = Plane(Point3D(1, 1, 2), normal_vector=(1, 1, 1))
  614. The projection is along the normal vector direction, not the z
  615. axis, so (1, 1) does not project to (1, 1, 2) on the plane A:
  616. >>> b = Point3D(1, 1)
  617. >>> A.projection(b)
  618. Point3D(5/3, 5/3, 2/3)
  619. >>> _ in A
  620. True
  621. But the point (1, 1, 2) projects to (1, 1) on the XY-plane:
  622. >>> XY = Plane((0, 0, 0), (0, 0, 1))
  623. >>> XY.projection((1, 1, 2))
  624. Point3D(1, 1, 0)
  625. """
  626. rv = Point(pt, dim=3)
  627. if rv in self:
  628. return rv
  629. return self.intersection(Line3D(rv, rv + Point3D(self.normal_vector)))[0]
  630. def random_point(self, seed=None):
  631. """ Returns a random point on the Plane.
  632. Returns
  633. =======
  634. Point3D
  635. Examples
  636. ========
  637. >>> from sympy import Plane
  638. >>> p = Plane((1, 0, 0), normal_vector=(0, 1, 0))
  639. >>> r = p.random_point(seed=42) # seed value is optional
  640. >>> r.n(3)
  641. Point3D(2.29, 0, -1.35)
  642. The random point can be moved to lie on the circle of radius
  643. 1 centered on p1:
  644. >>> c = p.p1 + (r - p.p1).unit
  645. >>> c.distance(p.p1).equals(1)
  646. True
  647. """
  648. if seed is not None:
  649. rng = random.Random(seed)
  650. else:
  651. rng = random
  652. params = {
  653. x: 2*Rational(rng.gauss(0, 1)) - 1,
  654. y: 2*Rational(rng.gauss(0, 1)) - 1}
  655. return self.arbitrary_point(x, y).subs(params)
  656. def parameter_value(self, other, u, v=None):
  657. """Return the parameter(s) corresponding to the given point.
  658. Examples
  659. ========
  660. >>> from sympy import pi, Plane
  661. >>> from sympy.abc import t, u, v
  662. >>> p = Plane((2, 0, 0), (0, 0, 1), (0, 1, 0))
  663. By default, the parameter value returned defines a point
  664. that is a distance of 1 from the Plane's p1 value and
  665. in line with the given point:
  666. >>> on_circle = p.arbitrary_point(t).subs(t, pi/4)
  667. >>> on_circle.distance(p.p1)
  668. 1
  669. >>> p.parameter_value(on_circle, t)
  670. {t: pi/4}
  671. Moving the point twice as far from p1 does not change
  672. the parameter value:
  673. >>> off_circle = p.p1 + (on_circle - p.p1)*2
  674. >>> off_circle.distance(p.p1)
  675. 2
  676. >>> p.parameter_value(off_circle, t)
  677. {t: pi/4}
  678. If the 2-value parameter is desired, supply the two
  679. parameter symbols and a replacement dictionary will
  680. be returned:
  681. >>> p.parameter_value(on_circle, u, v)
  682. {u: sqrt(10)/10, v: sqrt(10)/30}
  683. >>> p.parameter_value(off_circle, u, v)
  684. {u: sqrt(10)/5, v: sqrt(10)/15}
  685. """
  686. if not isinstance(other, GeometryEntity):
  687. other = Point(other, dim=self.ambient_dimension)
  688. if not isinstance(other, Point):
  689. raise ValueError("other must be a point")
  690. if other == self.p1:
  691. return other
  692. if isinstance(u, Symbol) and v is None:
  693. delta = self.arbitrary_point(u) - self.p1
  694. eq = delta - (other - self.p1).unit
  695. sol = solve(eq, u, dict=True)
  696. elif isinstance(u, Symbol) and isinstance(v, Symbol):
  697. pt = self.arbitrary_point(u, v)
  698. sol = solve(pt - other, (u, v), dict=True)
  699. else:
  700. raise ValueError('expecting 1 or 2 symbols')
  701. if not sol:
  702. raise ValueError("Given point is not on %s" % func_name(self))
  703. return sol[0] # {t: tval} or {u: uval, v: vval}
  704. @property
  705. def ambient_dimension(self):
  706. return self.p1.ambient_dimension