orienters.py 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398
  1. from sympy.core.basic import Basic
  2. from sympy.core.sympify import sympify
  3. from sympy.functions.elementary.trigonometric import (cos, sin)
  4. from sympy.matrices.dense import (eye, rot_axis1, rot_axis2, rot_axis3)
  5. from sympy.matrices.immutable import ImmutableDenseMatrix as Matrix
  6. from sympy.core.cache import cacheit
  7. from sympy.core.symbol import Str
  8. import sympy.vector
  9. class Orienter(Basic):
  10. """
  11. Super-class for all orienter classes.
  12. """
  13. def rotation_matrix(self):
  14. """
  15. The rotation matrix corresponding to this orienter
  16. instance.
  17. """
  18. return self._parent_orient
  19. class AxisOrienter(Orienter):
  20. """
  21. Class to denote an axis orienter.
  22. """
  23. def __new__(cls, angle, axis):
  24. if not isinstance(axis, sympy.vector.Vector):
  25. raise TypeError("axis should be a Vector")
  26. angle = sympify(angle)
  27. obj = super().__new__(cls, angle, axis)
  28. obj._angle = angle
  29. obj._axis = axis
  30. return obj
  31. def __init__(self, angle, axis):
  32. """
  33. Axis rotation is a rotation about an arbitrary axis by
  34. some angle. The angle is supplied as a SymPy expr scalar, and
  35. the axis is supplied as a Vector.
  36. Parameters
  37. ==========
  38. angle : Expr
  39. The angle by which the new system is to be rotated
  40. axis : Vector
  41. The axis around which the rotation has to be performed
  42. Examples
  43. ========
  44. >>> from sympy.vector import CoordSys3D
  45. >>> from sympy import symbols
  46. >>> q1 = symbols('q1')
  47. >>> N = CoordSys3D('N')
  48. >>> from sympy.vector import AxisOrienter
  49. >>> orienter = AxisOrienter(q1, N.i + 2 * N.j)
  50. >>> B = N.orient_new('B', (orienter, ))
  51. """
  52. # Dummy initializer for docstrings
  53. pass
  54. @cacheit
  55. def rotation_matrix(self, system):
  56. """
  57. The rotation matrix corresponding to this orienter
  58. instance.
  59. Parameters
  60. ==========
  61. system : CoordSys3D
  62. The coordinate system wrt which the rotation matrix
  63. is to be computed
  64. """
  65. axis = sympy.vector.express(self.axis, system).normalize()
  66. axis = axis.to_matrix(system)
  67. theta = self.angle
  68. parent_orient = ((eye(3) - axis * axis.T) * cos(theta) +
  69. Matrix([[0, -axis[2], axis[1]],
  70. [axis[2], 0, -axis[0]],
  71. [-axis[1], axis[0], 0]]) * sin(theta) +
  72. axis * axis.T)
  73. parent_orient = parent_orient.T
  74. return parent_orient
  75. @property
  76. def angle(self):
  77. return self._angle
  78. @property
  79. def axis(self):
  80. return self._axis
  81. class ThreeAngleOrienter(Orienter):
  82. """
  83. Super-class for Body and Space orienters.
  84. """
  85. def __new__(cls, angle1, angle2, angle3, rot_order):
  86. if isinstance(rot_order, Str):
  87. rot_order = rot_order.name
  88. approved_orders = ('123', '231', '312', '132', '213',
  89. '321', '121', '131', '212', '232',
  90. '313', '323', '')
  91. original_rot_order = rot_order
  92. rot_order = str(rot_order).upper()
  93. if not (len(rot_order) == 3):
  94. raise TypeError('rot_order should be a str of length 3')
  95. rot_order = [i.replace('X', '1') for i in rot_order]
  96. rot_order = [i.replace('Y', '2') for i in rot_order]
  97. rot_order = [i.replace('Z', '3') for i in rot_order]
  98. rot_order = ''.join(rot_order)
  99. if rot_order not in approved_orders:
  100. raise TypeError('Invalid rot_type parameter')
  101. a1 = int(rot_order[0])
  102. a2 = int(rot_order[1])
  103. a3 = int(rot_order[2])
  104. angle1 = sympify(angle1)
  105. angle2 = sympify(angle2)
  106. angle3 = sympify(angle3)
  107. if cls._in_order:
  108. parent_orient = (_rot(a1, angle1) *
  109. _rot(a2, angle2) *
  110. _rot(a3, angle3))
  111. else:
  112. parent_orient = (_rot(a3, angle3) *
  113. _rot(a2, angle2) *
  114. _rot(a1, angle1))
  115. parent_orient = parent_orient.T
  116. obj = super().__new__(
  117. cls, angle1, angle2, angle3, Str(rot_order))
  118. obj._angle1 = angle1
  119. obj._angle2 = angle2
  120. obj._angle3 = angle3
  121. obj._rot_order = original_rot_order
  122. obj._parent_orient = parent_orient
  123. return obj
  124. @property
  125. def angle1(self):
  126. return self._angle1
  127. @property
  128. def angle2(self):
  129. return self._angle2
  130. @property
  131. def angle3(self):
  132. return self._angle3
  133. @property
  134. def rot_order(self):
  135. return self._rot_order
  136. class BodyOrienter(ThreeAngleOrienter):
  137. """
  138. Class to denote a body-orienter.
  139. """
  140. _in_order = True
  141. def __new__(cls, angle1, angle2, angle3, rot_order):
  142. obj = ThreeAngleOrienter.__new__(cls, angle1, angle2, angle3,
  143. rot_order)
  144. return obj
  145. def __init__(self, angle1, angle2, angle3, rot_order):
  146. """
  147. Body orientation takes this coordinate system through three
  148. successive simple rotations.
  149. Body fixed rotations include both Euler Angles and
  150. Tait-Bryan Angles, see https://en.wikipedia.org/wiki/Euler_angles.
  151. Parameters
  152. ==========
  153. angle1, angle2, angle3 : Expr
  154. Three successive angles to rotate the coordinate system by
  155. rotation_order : string
  156. String defining the order of axes for rotation
  157. Examples
  158. ========
  159. >>> from sympy.vector import CoordSys3D, BodyOrienter
  160. >>> from sympy import symbols
  161. >>> q1, q2, q3 = symbols('q1 q2 q3')
  162. >>> N = CoordSys3D('N')
  163. A 'Body' fixed rotation is described by three angles and
  164. three body-fixed rotation axes. To orient a coordinate system D
  165. with respect to N, each sequential rotation is always about
  166. the orthogonal unit vectors fixed to D. For example, a '123'
  167. rotation will specify rotations about N.i, then D.j, then
  168. D.k. (Initially, D.i is same as N.i)
  169. Therefore,
  170. >>> body_orienter = BodyOrienter(q1, q2, q3, '123')
  171. >>> D = N.orient_new('D', (body_orienter, ))
  172. is same as
  173. >>> from sympy.vector import AxisOrienter
  174. >>> axis_orienter1 = AxisOrienter(q1, N.i)
  175. >>> D = N.orient_new('D', (axis_orienter1, ))
  176. >>> axis_orienter2 = AxisOrienter(q2, D.j)
  177. >>> D = D.orient_new('D', (axis_orienter2, ))
  178. >>> axis_orienter3 = AxisOrienter(q3, D.k)
  179. >>> D = D.orient_new('D', (axis_orienter3, ))
  180. Acceptable rotation orders are of length 3, expressed in XYZ or
  181. 123, and cannot have a rotation about about an axis twice in a row.
  182. >>> body_orienter1 = BodyOrienter(q1, q2, q3, '123')
  183. >>> body_orienter2 = BodyOrienter(q1, q2, 0, 'ZXZ')
  184. >>> body_orienter3 = BodyOrienter(0, 0, 0, 'XYX')
  185. """
  186. # Dummy initializer for docstrings
  187. pass
  188. class SpaceOrienter(ThreeAngleOrienter):
  189. """
  190. Class to denote a space-orienter.
  191. """
  192. _in_order = False
  193. def __new__(cls, angle1, angle2, angle3, rot_order):
  194. obj = ThreeAngleOrienter.__new__(cls, angle1, angle2, angle3,
  195. rot_order)
  196. return obj
  197. def __init__(self, angle1, angle2, angle3, rot_order):
  198. """
  199. Space rotation is similar to Body rotation, but the rotations
  200. are applied in the opposite order.
  201. Parameters
  202. ==========
  203. angle1, angle2, angle3 : Expr
  204. Three successive angles to rotate the coordinate system by
  205. rotation_order : string
  206. String defining the order of axes for rotation
  207. See Also
  208. ========
  209. BodyOrienter : Orienter to orient systems wrt Euler angles.
  210. Examples
  211. ========
  212. >>> from sympy.vector import CoordSys3D, SpaceOrienter
  213. >>> from sympy import symbols
  214. >>> q1, q2, q3 = symbols('q1 q2 q3')
  215. >>> N = CoordSys3D('N')
  216. To orient a coordinate system D with respect to N, each
  217. sequential rotation is always about N's orthogonal unit vectors.
  218. For example, a '123' rotation will specify rotations about
  219. N.i, then N.j, then N.k.
  220. Therefore,
  221. >>> space_orienter = SpaceOrienter(q1, q2, q3, '312')
  222. >>> D = N.orient_new('D', (space_orienter, ))
  223. is same as
  224. >>> from sympy.vector import AxisOrienter
  225. >>> axis_orienter1 = AxisOrienter(q1, N.i)
  226. >>> B = N.orient_new('B', (axis_orienter1, ))
  227. >>> axis_orienter2 = AxisOrienter(q2, N.j)
  228. >>> C = B.orient_new('C', (axis_orienter2, ))
  229. >>> axis_orienter3 = AxisOrienter(q3, N.k)
  230. >>> D = C.orient_new('C', (axis_orienter3, ))
  231. """
  232. # Dummy initializer for docstrings
  233. pass
  234. class QuaternionOrienter(Orienter):
  235. """
  236. Class to denote a quaternion-orienter.
  237. """
  238. def __new__(cls, q0, q1, q2, q3):
  239. q0 = sympify(q0)
  240. q1 = sympify(q1)
  241. q2 = sympify(q2)
  242. q3 = sympify(q3)
  243. parent_orient = (Matrix([[q0 ** 2 + q1 ** 2 - q2 ** 2 -
  244. q3 ** 2,
  245. 2 * (q1 * q2 - q0 * q3),
  246. 2 * (q0 * q2 + q1 * q3)],
  247. [2 * (q1 * q2 + q0 * q3),
  248. q0 ** 2 - q1 ** 2 +
  249. q2 ** 2 - q3 ** 2,
  250. 2 * (q2 * q3 - q0 * q1)],
  251. [2 * (q1 * q3 - q0 * q2),
  252. 2 * (q0 * q1 + q2 * q3),
  253. q0 ** 2 - q1 ** 2 -
  254. q2 ** 2 + q3 ** 2]]))
  255. parent_orient = parent_orient.T
  256. obj = super().__new__(cls, q0, q1, q2, q3)
  257. obj._q0 = q0
  258. obj._q1 = q1
  259. obj._q2 = q2
  260. obj._q3 = q3
  261. obj._parent_orient = parent_orient
  262. return obj
  263. def __init__(self, angle1, angle2, angle3, rot_order):
  264. """
  265. Quaternion orientation orients the new CoordSys3D with
  266. Quaternions, defined as a finite rotation about lambda, a unit
  267. vector, by some amount theta.
  268. This orientation is described by four parameters:
  269. q0 = cos(theta/2)
  270. q1 = lambda_x sin(theta/2)
  271. q2 = lambda_y sin(theta/2)
  272. q3 = lambda_z sin(theta/2)
  273. Quaternion does not take in a rotation order.
  274. Parameters
  275. ==========
  276. q0, q1, q2, q3 : Expr
  277. The quaternions to rotate the coordinate system by
  278. Examples
  279. ========
  280. >>> from sympy.vector import CoordSys3D
  281. >>> from sympy import symbols
  282. >>> q0, q1, q2, q3 = symbols('q0 q1 q2 q3')
  283. >>> N = CoordSys3D('N')
  284. >>> from sympy.vector import QuaternionOrienter
  285. >>> q_orienter = QuaternionOrienter(q0, q1, q2, q3)
  286. >>> B = N.orient_new('B', (q_orienter, ))
  287. """
  288. # Dummy initializer for docstrings
  289. pass
  290. @property
  291. def q0(self):
  292. return self._q0
  293. @property
  294. def q1(self):
  295. return self._q1
  296. @property
  297. def q2(self):
  298. return self._q2
  299. @property
  300. def q3(self):
  301. return self._q3
  302. def _rot(axis, angle):
  303. """DCM for simple axis 1, 2 or 3 rotations. """
  304. if axis == 1:
  305. return Matrix(rot_axis1(angle).T)
  306. elif axis == 2:
  307. return Matrix(rot_axis2(angle).T)
  308. elif axis == 3:
  309. return Matrix(rot_axis3(angle).T)