polygon.py 80 KB

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