coordsysrect.py 36 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031
  1. from collections.abc import Callable
  2. from sympy.core.basic import Basic
  3. from sympy.core.cache import cacheit
  4. from sympy.core import S, Dummy, Lambda
  5. from sympy.core.symbol import Str
  6. from sympy.core.symbol import symbols
  7. from sympy.matrices.immutable import ImmutableDenseMatrix as Matrix
  8. from sympy.matrices.matrixbase import MatrixBase
  9. from sympy.solvers import solve
  10. from sympy.vector.scalar import BaseScalar
  11. from sympy.core.containers import Tuple
  12. from sympy.core.function import diff
  13. from sympy.functions.elementary.miscellaneous import sqrt
  14. from sympy.functions.elementary.trigonometric import (acos, atan2, cos, sin)
  15. from sympy.matrices.dense import eye
  16. from sympy.matrices.immutable import ImmutableDenseMatrix
  17. from sympy.simplify.simplify import simplify
  18. from sympy.simplify.trigsimp import trigsimp
  19. import sympy.vector
  20. from sympy.vector.orienters import (Orienter, AxisOrienter, BodyOrienter,
  21. SpaceOrienter, QuaternionOrienter)
  22. class CoordSys3D(Basic):
  23. """
  24. Represents a coordinate system in 3-D space.
  25. """
  26. def __new__(cls, name, transformation=None, parent=None, location=None,
  27. rotation_matrix=None, vector_names=None, variable_names=None):
  28. """
  29. The orientation/location parameters are necessary if this system
  30. is being defined at a certain orientation or location wrt another.
  31. Parameters
  32. ==========
  33. name : str
  34. The name of the new CoordSys3D instance.
  35. transformation : Lambda, Tuple, str
  36. Transformation defined by transformation equations or chosen
  37. from predefined ones.
  38. location : Vector
  39. The position vector of the new system's origin wrt the parent
  40. instance.
  41. rotation_matrix : SymPy ImmutableMatrix
  42. The rotation matrix of the new coordinate system with respect
  43. to the parent. In other words, the output of
  44. new_system.rotation_matrix(parent).
  45. parent : CoordSys3D
  46. The coordinate system wrt which the orientation/location
  47. (or both) is being defined.
  48. vector_names, variable_names : iterable(optional)
  49. Iterables of 3 strings each, with custom names for base
  50. vectors and base scalars of the new system respectively.
  51. Used for simple str printing.
  52. """
  53. name = str(name)
  54. Vector = sympy.vector.Vector
  55. Point = sympy.vector.Point
  56. if not isinstance(name, str):
  57. raise TypeError("name should be a string")
  58. if transformation is not None:
  59. if (location is not None) or (rotation_matrix is not None):
  60. raise ValueError("specify either `transformation` or "
  61. "`location`/`rotation_matrix`")
  62. if isinstance(transformation, (Tuple, tuple, list)):
  63. if isinstance(transformation[0], MatrixBase):
  64. rotation_matrix = transformation[0]
  65. location = transformation[1]
  66. else:
  67. transformation = Lambda(transformation[0],
  68. transformation[1])
  69. elif isinstance(transformation, Callable):
  70. x1, x2, x3 = symbols('x1 x2 x3', cls=Dummy)
  71. transformation = Lambda((x1, x2, x3),
  72. transformation(x1, x2, x3))
  73. elif isinstance(transformation, str):
  74. transformation = Str(transformation)
  75. elif isinstance(transformation, (Str, Lambda)):
  76. pass
  77. else:
  78. raise TypeError("transformation: "
  79. "wrong type {}".format(type(transformation)))
  80. # If orientation information has been provided, store
  81. # the rotation matrix accordingly
  82. if rotation_matrix is None:
  83. rotation_matrix = ImmutableDenseMatrix(eye(3))
  84. else:
  85. if not isinstance(rotation_matrix, MatrixBase):
  86. raise TypeError("rotation_matrix should be an Immutable" +
  87. "Matrix instance")
  88. rotation_matrix = rotation_matrix.as_immutable()
  89. # If location information is not given, adjust the default
  90. # location as Vector.zero
  91. if parent is not None:
  92. if not isinstance(parent, CoordSys3D):
  93. raise TypeError("parent should be a " +
  94. "CoordSys3D/None")
  95. if location is None:
  96. location = Vector.zero
  97. else:
  98. if not isinstance(location, Vector):
  99. raise TypeError("location should be a Vector")
  100. # Check that location does not contain base
  101. # scalars
  102. for x in location.free_symbols:
  103. if isinstance(x, BaseScalar):
  104. raise ValueError("location should not contain" +
  105. " BaseScalars")
  106. origin = parent.origin.locate_new(name + '.origin',
  107. location)
  108. else:
  109. location = Vector.zero
  110. origin = Point(name + '.origin')
  111. if transformation is None:
  112. transformation = Tuple(rotation_matrix, location)
  113. if isinstance(transformation, Tuple):
  114. lambda_transformation = CoordSys3D._compose_rotation_and_translation(
  115. transformation[0],
  116. transformation[1],
  117. parent
  118. )
  119. r, l = transformation
  120. l = l._projections
  121. lambda_lame = CoordSys3D._get_lame_coeff('cartesian')
  122. lambda_inverse = lambda x, y, z: r.inv()*Matrix(
  123. [x-l[0], y-l[1], z-l[2]])
  124. elif isinstance(transformation, Str):
  125. trname = transformation.name
  126. lambda_transformation = CoordSys3D._get_transformation_lambdas(trname)
  127. if parent is not None:
  128. if parent.lame_coefficients() != (S.One, S.One, S.One):
  129. raise ValueError('Parent for pre-defined coordinate '
  130. 'system should be Cartesian.')
  131. lambda_lame = CoordSys3D._get_lame_coeff(trname)
  132. lambda_inverse = CoordSys3D._set_inv_trans_equations(trname)
  133. elif isinstance(transformation, Lambda):
  134. if not CoordSys3D._check_orthogonality(transformation):
  135. raise ValueError("The transformation equation does not "
  136. "create orthogonal coordinate system")
  137. lambda_transformation = transformation
  138. lambda_lame = CoordSys3D._calculate_lame_coeff(lambda_transformation)
  139. lambda_inverse = None
  140. else:
  141. lambda_transformation = lambda x, y, z: transformation(x, y, z)
  142. lambda_lame = CoordSys3D._get_lame_coeff(transformation)
  143. lambda_inverse = None
  144. if variable_names is None:
  145. if isinstance(transformation, Lambda):
  146. variable_names = ["x1", "x2", "x3"]
  147. elif isinstance(transformation, Str):
  148. if transformation.name == 'spherical':
  149. variable_names = ["r", "theta", "phi"]
  150. elif transformation.name == 'cylindrical':
  151. variable_names = ["r", "theta", "z"]
  152. else:
  153. variable_names = ["x", "y", "z"]
  154. else:
  155. variable_names = ["x", "y", "z"]
  156. if vector_names is None:
  157. vector_names = ["i", "j", "k"]
  158. # All systems that are defined as 'roots' are unequal, unless
  159. # they have the same name.
  160. # Systems defined at same orientation/position wrt the same
  161. # 'parent' are equal, irrespective of the name.
  162. # This is true even if the same orientation is provided via
  163. # different methods like Axis/Body/Space/Quaternion.
  164. # However, coincident systems may be seen as unequal if
  165. # positioned/oriented wrt different parents, even though
  166. # they may actually be 'coincident' wrt the root system.
  167. if parent is not None:
  168. obj = super().__new__(
  169. cls, Str(name), transformation, parent)
  170. else:
  171. obj = super().__new__(
  172. cls, Str(name), transformation)
  173. obj._name = name
  174. # Initialize the base vectors
  175. _check_strings('vector_names', vector_names)
  176. vector_names = list(vector_names)
  177. latex_vects = [(r'\mathbf{\hat{%s}_{%s}}' % (x, name)) for
  178. x in vector_names]
  179. pretty_vects = ['%s_%s' % (x, name) for x in vector_names]
  180. obj._vector_names = vector_names
  181. v1 = BaseVector(0, obj, pretty_vects[0], latex_vects[0])
  182. v2 = BaseVector(1, obj, pretty_vects[1], latex_vects[1])
  183. v3 = BaseVector(2, obj, pretty_vects[2], latex_vects[2])
  184. obj._base_vectors = (v1, v2, v3)
  185. # Initialize the base scalars
  186. _check_strings('variable_names', vector_names)
  187. variable_names = list(variable_names)
  188. latex_scalars = [(r"\mathbf{{%s}_{%s}}" % (x, name)) for
  189. x in variable_names]
  190. pretty_scalars = ['%s_%s' % (x, name) for x in variable_names]
  191. obj._variable_names = variable_names
  192. obj._vector_names = vector_names
  193. x1 = BaseScalar(0, obj, pretty_scalars[0], latex_scalars[0])
  194. x2 = BaseScalar(1, obj, pretty_scalars[1], latex_scalars[1])
  195. x3 = BaseScalar(2, obj, pretty_scalars[2], latex_scalars[2])
  196. obj._base_scalars = (x1, x2, x3)
  197. obj._transformation = transformation
  198. obj._transformation_lambda = lambda_transformation
  199. obj._lame_coefficients = lambda_lame(x1, x2, x3)
  200. obj._transformation_from_parent_lambda = lambda_inverse
  201. setattr(obj, variable_names[0], x1)
  202. setattr(obj, variable_names[1], x2)
  203. setattr(obj, variable_names[2], x3)
  204. setattr(obj, vector_names[0], v1)
  205. setattr(obj, vector_names[1], v2)
  206. setattr(obj, vector_names[2], v3)
  207. # Assign params
  208. obj._parent = parent
  209. if obj._parent is not None:
  210. obj._root = obj._parent._root
  211. else:
  212. obj._root = obj
  213. obj._parent_rotation_matrix = rotation_matrix
  214. obj._origin = origin
  215. # Return the instance
  216. return obj
  217. def _sympystr(self, printer):
  218. return self._name
  219. def __iter__(self):
  220. return iter(self.base_vectors())
  221. @staticmethod
  222. def _check_orthogonality(equations):
  223. """
  224. Helper method for _connect_to_cartesian. It checks if
  225. set of transformation equations create orthogonal curvilinear
  226. coordinate system
  227. Parameters
  228. ==========
  229. equations : Lambda
  230. Lambda of transformation equations
  231. """
  232. x1, x2, x3 = symbols("x1, x2, x3", cls=Dummy)
  233. equations = equations(x1, x2, x3)
  234. v1 = Matrix([diff(equations[0], x1),
  235. diff(equations[1], x1), diff(equations[2], x1)])
  236. v2 = Matrix([diff(equations[0], x2),
  237. diff(equations[1], x2), diff(equations[2], x2)])
  238. v3 = Matrix([diff(equations[0], x3),
  239. diff(equations[1], x3), diff(equations[2], x3)])
  240. if any(simplify(i[0] + i[1] + i[2]) == 0 for i in (v1, v2, v3)):
  241. return False
  242. else:
  243. if simplify(v1.dot(v2)) == 0 and simplify(v2.dot(v3)) == 0 \
  244. and simplify(v3.dot(v1)) == 0:
  245. return True
  246. else:
  247. return False
  248. @staticmethod
  249. def _set_inv_trans_equations(curv_coord_name):
  250. """
  251. Store information about inverse transformation equations for
  252. pre-defined coordinate systems.
  253. Parameters
  254. ==========
  255. curv_coord_name : str
  256. Name of coordinate system
  257. """
  258. if curv_coord_name == 'cartesian':
  259. return lambda x, y, z: (x, y, z)
  260. if curv_coord_name == 'spherical':
  261. return lambda x, y, z: (
  262. sqrt(x**2 + y**2 + z**2),
  263. acos(z/sqrt(x**2 + y**2 + z**2)),
  264. atan2(y, x)
  265. )
  266. if curv_coord_name == 'cylindrical':
  267. return lambda x, y, z: (
  268. sqrt(x**2 + y**2),
  269. atan2(y, x),
  270. z
  271. )
  272. raise ValueError('Wrong set of parameters.'
  273. 'Type of coordinate system is defined')
  274. def _calculate_inv_trans_equations(self):
  275. """
  276. Helper method for set_coordinate_type. It calculates inverse
  277. transformation equations for given transformations equations.
  278. """
  279. x1, x2, x3 = symbols("x1, x2, x3", cls=Dummy, reals=True)
  280. x, y, z = symbols("x, y, z", cls=Dummy)
  281. equations = self._transformation(x1, x2, x3)
  282. solved = solve([equations[0] - x,
  283. equations[1] - y,
  284. equations[2] - z], (x1, x2, x3), dict=True)[0]
  285. solved = solved[x1], solved[x2], solved[x3]
  286. self._transformation_from_parent_lambda = \
  287. lambda x1, x2, x3: tuple(i.subs(list(zip((x, y, z), (x1, x2, x3)))) for i in solved)
  288. @staticmethod
  289. def _get_lame_coeff(curv_coord_name):
  290. """
  291. Store information about Lame coefficients for pre-defined
  292. coordinate systems.
  293. Parameters
  294. ==========
  295. curv_coord_name : str
  296. Name of coordinate system
  297. """
  298. if isinstance(curv_coord_name, str):
  299. if curv_coord_name == 'cartesian':
  300. return lambda x, y, z: (S.One, S.One, S.One)
  301. if curv_coord_name == 'spherical':
  302. return lambda r, theta, phi: (S.One, r, r*sin(theta))
  303. if curv_coord_name == 'cylindrical':
  304. return lambda r, theta, h: (S.One, r, S.One)
  305. raise ValueError('Wrong set of parameters.'
  306. ' Type of coordinate system is not defined')
  307. return CoordSys3D._calculate_lame_coefficients(curv_coord_name)
  308. @staticmethod
  309. def _calculate_lame_coeff(equations):
  310. """
  311. It calculates Lame coefficients
  312. for given transformations equations.
  313. Parameters
  314. ==========
  315. equations : Lambda
  316. Lambda of transformation equations.
  317. """
  318. return lambda x1, x2, x3: (
  319. sqrt(diff(equations(x1, x2, x3)[0], x1)**2 +
  320. diff(equations(x1, x2, x3)[1], x1)**2 +
  321. diff(equations(x1, x2, x3)[2], x1)**2),
  322. sqrt(diff(equations(x1, x2, x3)[0], x2)**2 +
  323. diff(equations(x1, x2, x3)[1], x2)**2 +
  324. diff(equations(x1, x2, x3)[2], x2)**2),
  325. sqrt(diff(equations(x1, x2, x3)[0], x3)**2 +
  326. diff(equations(x1, x2, x3)[1], x3)**2 +
  327. diff(equations(x1, x2, x3)[2], x3)**2)
  328. )
  329. def _inverse_rotation_matrix(self):
  330. """
  331. Returns inverse rotation matrix.
  332. """
  333. return simplify(self._parent_rotation_matrix**-1)
  334. @staticmethod
  335. def _get_transformation_lambdas(curv_coord_name):
  336. """
  337. Store information about transformation equations for pre-defined
  338. coordinate systems.
  339. Parameters
  340. ==========
  341. curv_coord_name : str
  342. Name of coordinate system
  343. """
  344. if isinstance(curv_coord_name, str):
  345. if curv_coord_name == 'cartesian':
  346. return lambda x, y, z: (x, y, z)
  347. if curv_coord_name == 'spherical':
  348. return lambda r, theta, phi: (
  349. r*sin(theta)*cos(phi),
  350. r*sin(theta)*sin(phi),
  351. r*cos(theta)
  352. )
  353. if curv_coord_name == 'cylindrical':
  354. return lambda r, theta, h: (
  355. r*cos(theta),
  356. r*sin(theta),
  357. h
  358. )
  359. raise ValueError('Wrong set of parameters.'
  360. 'Type of coordinate system is defined')
  361. @classmethod
  362. def _rotation_trans_equations(cls, matrix, equations):
  363. """
  364. Returns the transformation equations obtained from rotation matrix.
  365. Parameters
  366. ==========
  367. matrix : Matrix
  368. Rotation matrix
  369. equations : tuple
  370. Transformation equations
  371. """
  372. return tuple(matrix * Matrix(equations))
  373. @property
  374. def origin(self):
  375. return self._origin
  376. def base_vectors(self):
  377. return self._base_vectors
  378. def base_scalars(self):
  379. return self._base_scalars
  380. def lame_coefficients(self):
  381. return self._lame_coefficients
  382. def transformation_to_parent(self):
  383. return self._transformation_lambda(*self.base_scalars())
  384. def transformation_from_parent(self):
  385. if self._parent is None:
  386. raise ValueError("no parent coordinate system, use "
  387. "`transformation_from_parent_function()`")
  388. return self._transformation_from_parent_lambda(
  389. *self._parent.base_scalars())
  390. def transformation_from_parent_function(self):
  391. return self._transformation_from_parent_lambda
  392. def rotation_matrix(self, other):
  393. """
  394. Returns the direction cosine matrix(DCM), also known as the
  395. 'rotation matrix' of this coordinate system with respect to
  396. another system.
  397. If v_a is a vector defined in system 'A' (in matrix format)
  398. and v_b is the same vector defined in system 'B', then
  399. v_a = A.rotation_matrix(B) * v_b.
  400. A SymPy Matrix is returned.
  401. Parameters
  402. ==========
  403. other : CoordSys3D
  404. The system which the DCM is generated to.
  405. Examples
  406. ========
  407. >>> from sympy.vector import CoordSys3D
  408. >>> from sympy import symbols
  409. >>> q1 = symbols('q1')
  410. >>> N = CoordSys3D('N')
  411. >>> A = N.orient_new_axis('A', q1, N.i)
  412. >>> N.rotation_matrix(A)
  413. Matrix([
  414. [1, 0, 0],
  415. [0, cos(q1), -sin(q1)],
  416. [0, sin(q1), cos(q1)]])
  417. """
  418. from sympy.vector.functions import _path
  419. if not isinstance(other, CoordSys3D):
  420. raise TypeError(str(other) +
  421. " is not a CoordSys3D")
  422. # Handle special cases
  423. if other == self:
  424. return eye(3)
  425. elif other == self._parent:
  426. return self._parent_rotation_matrix
  427. elif other._parent == self:
  428. return other._parent_rotation_matrix.T
  429. # Else, use tree to calculate position
  430. rootindex, path = _path(self, other)
  431. result = eye(3)
  432. for i in range(rootindex):
  433. result *= path[i]._parent_rotation_matrix
  434. for i in range(rootindex + 1, len(path)):
  435. result *= path[i]._parent_rotation_matrix.T
  436. return result
  437. @cacheit
  438. def position_wrt(self, other):
  439. """
  440. Returns the position vector of the origin of this coordinate
  441. system with respect to another Point/CoordSys3D.
  442. Parameters
  443. ==========
  444. other : Point/CoordSys3D
  445. If other is a Point, the position of this system's origin
  446. wrt it is returned. If its an instance of CoordSyRect,
  447. the position wrt its origin is returned.
  448. Examples
  449. ========
  450. >>> from sympy.vector import CoordSys3D
  451. >>> N = CoordSys3D('N')
  452. >>> N1 = N.locate_new('N1', 10 * N.i)
  453. >>> N.position_wrt(N1)
  454. (-10)*N.i
  455. """
  456. return self.origin.position_wrt(other)
  457. def scalar_map(self, other):
  458. """
  459. Returns a dictionary which expresses the coordinate variables
  460. (base scalars) of this frame in terms of the variables of
  461. otherframe.
  462. Parameters
  463. ==========
  464. otherframe : CoordSys3D
  465. The other system to map the variables to.
  466. Examples
  467. ========
  468. >>> from sympy.vector import CoordSys3D
  469. >>> from sympy import Symbol
  470. >>> A = CoordSys3D('A')
  471. >>> q = Symbol('q')
  472. >>> B = A.orient_new_axis('B', q, A.k)
  473. >>> A.scalar_map(B)
  474. {A.x: B.x*cos(q) - B.y*sin(q), A.y: B.x*sin(q) + B.y*cos(q), A.z: B.z}
  475. """
  476. origin_coords = tuple(self.position_wrt(other).to_matrix(other))
  477. relocated_scalars = [x - origin_coords[i]
  478. for i, x in enumerate(other.base_scalars())]
  479. vars_matrix = (self.rotation_matrix(other) *
  480. Matrix(relocated_scalars))
  481. return {x: trigsimp(vars_matrix[i])
  482. for i, x in enumerate(self.base_scalars())}
  483. def locate_new(self, name, position, vector_names=None,
  484. variable_names=None):
  485. """
  486. Returns a CoordSys3D with its origin located at the given
  487. position wrt this coordinate system's origin.
  488. Parameters
  489. ==========
  490. name : str
  491. The name of the new CoordSys3D instance.
  492. position : Vector
  493. The position vector of the new system's origin wrt this
  494. one.
  495. vector_names, variable_names : iterable(optional)
  496. Iterables of 3 strings each, with custom names for base
  497. vectors and base scalars of the new system respectively.
  498. Used for simple str printing.
  499. Examples
  500. ========
  501. >>> from sympy.vector import CoordSys3D
  502. >>> A = CoordSys3D('A')
  503. >>> B = A.locate_new('B', 10 * A.i)
  504. >>> B.origin.position_wrt(A.origin)
  505. 10*A.i
  506. """
  507. if variable_names is None:
  508. variable_names = self._variable_names
  509. if vector_names is None:
  510. vector_names = self._vector_names
  511. return CoordSys3D(name, location=position,
  512. vector_names=vector_names,
  513. variable_names=variable_names,
  514. parent=self)
  515. def orient_new(self, name, orienters, location=None,
  516. vector_names=None, variable_names=None):
  517. """
  518. Creates a new CoordSys3D oriented in the user-specified way
  519. with respect to this system.
  520. Please refer to the documentation of the orienter classes
  521. for more information about the orientation procedure.
  522. Parameters
  523. ==========
  524. name : str
  525. The name of the new CoordSys3D instance.
  526. orienters : iterable/Orienter
  527. An Orienter or an iterable of Orienters for orienting the
  528. new coordinate system.
  529. If an Orienter is provided, it is applied to get the new
  530. system.
  531. If an iterable is provided, the orienters will be applied
  532. in the order in which they appear in the iterable.
  533. location : Vector(optional)
  534. The location of the new coordinate system's origin wrt this
  535. system's origin. If not specified, the origins are taken to
  536. be coincident.
  537. vector_names, variable_names : iterable(optional)
  538. Iterables of 3 strings each, with custom names for base
  539. vectors and base scalars of the new system respectively.
  540. Used for simple str printing.
  541. Examples
  542. ========
  543. >>> from sympy.vector import CoordSys3D
  544. >>> from sympy import symbols
  545. >>> q0, q1, q2, q3 = symbols('q0 q1 q2 q3')
  546. >>> N = CoordSys3D('N')
  547. Using an AxisOrienter
  548. >>> from sympy.vector import AxisOrienter
  549. >>> axis_orienter = AxisOrienter(q1, N.i + 2 * N.j)
  550. >>> A = N.orient_new('A', (axis_orienter, ))
  551. Using a BodyOrienter
  552. >>> from sympy.vector import BodyOrienter
  553. >>> body_orienter = BodyOrienter(q1, q2, q3, '123')
  554. >>> B = N.orient_new('B', (body_orienter, ))
  555. Using a SpaceOrienter
  556. >>> from sympy.vector import SpaceOrienter
  557. >>> space_orienter = SpaceOrienter(q1, q2, q3, '312')
  558. >>> C = N.orient_new('C', (space_orienter, ))
  559. Using a QuaternionOrienter
  560. >>> from sympy.vector import QuaternionOrienter
  561. >>> q_orienter = QuaternionOrienter(q0, q1, q2, q3)
  562. >>> D = N.orient_new('D', (q_orienter, ))
  563. """
  564. if variable_names is None:
  565. variable_names = self._variable_names
  566. if vector_names is None:
  567. vector_names = self._vector_names
  568. if isinstance(orienters, Orienter):
  569. if isinstance(orienters, AxisOrienter):
  570. final_matrix = orienters.rotation_matrix(self)
  571. else:
  572. final_matrix = orienters.rotation_matrix()
  573. # TODO: trigsimp is needed here so that the matrix becomes
  574. # canonical (scalar_map also calls trigsimp; without this, you can
  575. # end up with the same CoordinateSystem that compares differently
  576. # due to a differently formatted matrix). However, this is
  577. # probably not so good for performance.
  578. final_matrix = trigsimp(final_matrix)
  579. else:
  580. final_matrix = Matrix(eye(3))
  581. for orienter in orienters:
  582. if isinstance(orienter, AxisOrienter):
  583. final_matrix *= orienter.rotation_matrix(self)
  584. else:
  585. final_matrix *= orienter.rotation_matrix()
  586. return CoordSys3D(name, rotation_matrix=final_matrix,
  587. vector_names=vector_names,
  588. variable_names=variable_names,
  589. location=location,
  590. parent=self)
  591. def orient_new_axis(self, name, angle, axis, location=None,
  592. vector_names=None, variable_names=None):
  593. """
  594. Axis rotation is a rotation about an arbitrary axis by
  595. some angle. The angle is supplied as a SymPy expr scalar, and
  596. the axis is supplied as a Vector.
  597. Parameters
  598. ==========
  599. name : string
  600. The name of the new coordinate system
  601. angle : Expr
  602. The angle by which the new system is to be rotated
  603. axis : Vector
  604. The axis around which the rotation has to be performed
  605. location : Vector(optional)
  606. The location of the new coordinate system's origin wrt this
  607. system's origin. If not specified, the origins are taken to
  608. be coincident.
  609. vector_names, variable_names : iterable(optional)
  610. Iterables of 3 strings each, with custom names for base
  611. vectors and base scalars of the new system respectively.
  612. Used for simple str printing.
  613. Examples
  614. ========
  615. >>> from sympy.vector import CoordSys3D
  616. >>> from sympy import symbols
  617. >>> q1 = symbols('q1')
  618. >>> N = CoordSys3D('N')
  619. >>> B = N.orient_new_axis('B', q1, N.i + 2 * N.j)
  620. """
  621. if variable_names is None:
  622. variable_names = self._variable_names
  623. if vector_names is None:
  624. vector_names = self._vector_names
  625. orienter = AxisOrienter(angle, axis)
  626. return self.orient_new(name, orienter,
  627. location=location,
  628. vector_names=vector_names,
  629. variable_names=variable_names)
  630. def orient_new_body(self, name, angle1, angle2, angle3,
  631. rotation_order, location=None,
  632. vector_names=None, variable_names=None):
  633. """
  634. Body orientation takes this coordinate system through three
  635. successive simple rotations.
  636. Body fixed rotations include both Euler Angles and
  637. Tait-Bryan Angles, see https://en.wikipedia.org/wiki/Euler_angles.
  638. Parameters
  639. ==========
  640. name : string
  641. The name of the new coordinate system
  642. angle1, angle2, angle3 : Expr
  643. Three successive angles to rotate the coordinate system by
  644. rotation_order : string
  645. String defining the order of axes for rotation
  646. location : Vector(optional)
  647. The location of the new coordinate system's origin wrt this
  648. system's origin. If not specified, the origins are taken to
  649. be coincident.
  650. vector_names, variable_names : iterable(optional)
  651. Iterables of 3 strings each, with custom names for base
  652. vectors and base scalars of the new system respectively.
  653. Used for simple str printing.
  654. Examples
  655. ========
  656. >>> from sympy.vector import CoordSys3D
  657. >>> from sympy import symbols
  658. >>> q1, q2, q3 = symbols('q1 q2 q3')
  659. >>> N = CoordSys3D('N')
  660. A 'Body' fixed rotation is described by three angles and
  661. three body-fixed rotation axes. To orient a coordinate system D
  662. with respect to N, each sequential rotation is always about
  663. the orthogonal unit vectors fixed to D. For example, a '123'
  664. rotation will specify rotations about N.i, then D.j, then
  665. D.k. (Initially, D.i is same as N.i)
  666. Therefore,
  667. >>> D = N.orient_new_body('D', q1, q2, q3, '123')
  668. is same as
  669. >>> D = N.orient_new_axis('D', q1, N.i)
  670. >>> D = D.orient_new_axis('D', q2, D.j)
  671. >>> D = D.orient_new_axis('D', q3, D.k)
  672. Acceptable rotation orders are of length 3, expressed in XYZ or
  673. 123, and cannot have a rotation about about an axis twice in a row.
  674. >>> B = N.orient_new_body('B', q1, q2, q3, '123')
  675. >>> B = N.orient_new_body('B', q1, q2, 0, 'ZXZ')
  676. >>> B = N.orient_new_body('B', 0, 0, 0, 'XYX')
  677. """
  678. orienter = BodyOrienter(angle1, angle2, angle3, rotation_order)
  679. return self.orient_new(name, orienter,
  680. location=location,
  681. vector_names=vector_names,
  682. variable_names=variable_names)
  683. def orient_new_space(self, name, angle1, angle2, angle3,
  684. rotation_order, location=None,
  685. vector_names=None, variable_names=None):
  686. """
  687. Space rotation is similar to Body rotation, but the rotations
  688. are applied in the opposite order.
  689. Parameters
  690. ==========
  691. name : string
  692. The name of the new coordinate system
  693. angle1, angle2, angle3 : Expr
  694. Three successive angles to rotate the coordinate system by
  695. rotation_order : string
  696. String defining the order of axes for rotation
  697. location : Vector(optional)
  698. The location of the new coordinate system's origin wrt this
  699. system's origin. If not specified, the origins are taken to
  700. be coincident.
  701. vector_names, variable_names : iterable(optional)
  702. Iterables of 3 strings each, with custom names for base
  703. vectors and base scalars of the new system respectively.
  704. Used for simple str printing.
  705. See Also
  706. ========
  707. CoordSys3D.orient_new_body : method to orient via Euler
  708. angles
  709. Examples
  710. ========
  711. >>> from sympy.vector import CoordSys3D
  712. >>> from sympy import symbols
  713. >>> q1, q2, q3 = symbols('q1 q2 q3')
  714. >>> N = CoordSys3D('N')
  715. To orient a coordinate system D with respect to N, each
  716. sequential rotation is always about N's orthogonal unit vectors.
  717. For example, a '123' rotation will specify rotations about
  718. N.i, then N.j, then N.k.
  719. Therefore,
  720. >>> D = N.orient_new_space('D', q1, q2, q3, '312')
  721. is same as
  722. >>> B = N.orient_new_axis('B', q1, N.i)
  723. >>> C = B.orient_new_axis('C', q2, N.j)
  724. >>> D = C.orient_new_axis('D', q3, N.k)
  725. """
  726. orienter = SpaceOrienter(angle1, angle2, angle3, rotation_order)
  727. return self.orient_new(name, orienter,
  728. location=location,
  729. vector_names=vector_names,
  730. variable_names=variable_names)
  731. def orient_new_quaternion(self, name, q0, q1, q2, q3, location=None,
  732. vector_names=None, variable_names=None):
  733. """
  734. Quaternion orientation orients the new CoordSys3D with
  735. Quaternions, defined as a finite rotation about lambda, a unit
  736. vector, by some amount theta.
  737. This orientation is described by four parameters:
  738. q0 = cos(theta/2)
  739. q1 = lambda_x sin(theta/2)
  740. q2 = lambda_y sin(theta/2)
  741. q3 = lambda_z sin(theta/2)
  742. Quaternion does not take in a rotation order.
  743. Parameters
  744. ==========
  745. name : string
  746. The name of the new coordinate system
  747. q0, q1, q2, q3 : Expr
  748. The quaternions to rotate the coordinate system by
  749. location : Vector(optional)
  750. The location of the new coordinate system's origin wrt this
  751. system's origin. If not specified, the origins are taken to
  752. be coincident.
  753. vector_names, variable_names : iterable(optional)
  754. Iterables of 3 strings each, with custom names for base
  755. vectors and base scalars of the new system respectively.
  756. Used for simple str printing.
  757. Examples
  758. ========
  759. >>> from sympy.vector import CoordSys3D
  760. >>> from sympy import symbols
  761. >>> q0, q1, q2, q3 = symbols('q0 q1 q2 q3')
  762. >>> N = CoordSys3D('N')
  763. >>> B = N.orient_new_quaternion('B', q0, q1, q2, q3)
  764. """
  765. orienter = QuaternionOrienter(q0, q1, q2, q3)
  766. return self.orient_new(name, orienter,
  767. location=location,
  768. vector_names=vector_names,
  769. variable_names=variable_names)
  770. def create_new(self, name, transformation, variable_names=None, vector_names=None):
  771. """
  772. Returns a CoordSys3D which is connected to self by transformation.
  773. Parameters
  774. ==========
  775. name : str
  776. The name of the new CoordSys3D instance.
  777. transformation : Lambda, Tuple, str
  778. Transformation defined by transformation equations or chosen
  779. from predefined ones.
  780. vector_names, variable_names : iterable(optional)
  781. Iterables of 3 strings each, with custom names for base
  782. vectors and base scalars of the new system respectively.
  783. Used for simple str printing.
  784. Examples
  785. ========
  786. >>> from sympy.vector import CoordSys3D
  787. >>> a = CoordSys3D('a')
  788. >>> b = a.create_new('b', transformation='spherical')
  789. >>> b.transformation_to_parent()
  790. (b.r*sin(b.theta)*cos(b.phi), b.r*sin(b.phi)*sin(b.theta), b.r*cos(b.theta))
  791. >>> b.transformation_from_parent()
  792. (sqrt(a.x**2 + a.y**2 + a.z**2), acos(a.z/sqrt(a.x**2 + a.y**2 + a.z**2)), atan2(a.y, a.x))
  793. """
  794. return CoordSys3D(name, parent=self, transformation=transformation,
  795. variable_names=variable_names, vector_names=vector_names)
  796. def __init__(self, name, location=None, rotation_matrix=None,
  797. parent=None, vector_names=None, variable_names=None,
  798. latex_vects=None, pretty_vects=None, latex_scalars=None,
  799. pretty_scalars=None, transformation=None):
  800. # Dummy initializer for setting docstring
  801. pass
  802. __init__.__doc__ = __new__.__doc__
  803. @staticmethod
  804. def _compose_rotation_and_translation(rot, translation, parent):
  805. r = lambda x, y, z: CoordSys3D._rotation_trans_equations(rot, (x, y, z))
  806. if parent is None:
  807. return r
  808. dx, dy, dz = [translation.dot(i) for i in parent.base_vectors()]
  809. t = lambda x, y, z: (
  810. x + dx,
  811. y + dy,
  812. z + dz,
  813. )
  814. return lambda x, y, z: t(*r(x, y, z))
  815. def _check_strings(arg_name, arg):
  816. errorstr = arg_name + " must be an iterable of 3 string-types"
  817. if len(arg) != 3:
  818. raise ValueError(errorstr)
  819. for s in arg:
  820. if not isinstance(s, str):
  821. raise TypeError(errorstr)
  822. # Delayed import to avoid cyclic import problems:
  823. from sympy.vector.vector import BaseVector