coordsysrect.py 36 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034
  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.matrices 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. i = -1
  433. for i in range(rootindex):
  434. result *= path[i]._parent_rotation_matrix
  435. i += 2
  436. while i < len(path):
  437. result *= path[i]._parent_rotation_matrix.T
  438. i += 1
  439. return result
  440. @cacheit
  441. def position_wrt(self, other):
  442. """
  443. Returns the position vector of the origin of this coordinate
  444. system with respect to another Point/CoordSys3D.
  445. Parameters
  446. ==========
  447. other : Point/CoordSys3D
  448. If other is a Point, the position of this system's origin
  449. wrt it is returned. If its an instance of CoordSyRect,
  450. the position wrt its origin is returned.
  451. Examples
  452. ========
  453. >>> from sympy.vector import CoordSys3D
  454. >>> N = CoordSys3D('N')
  455. >>> N1 = N.locate_new('N1', 10 * N.i)
  456. >>> N.position_wrt(N1)
  457. (-10)*N.i
  458. """
  459. return self.origin.position_wrt(other)
  460. def scalar_map(self, other):
  461. """
  462. Returns a dictionary which expresses the coordinate variables
  463. (base scalars) of this frame in terms of the variables of
  464. otherframe.
  465. Parameters
  466. ==========
  467. otherframe : CoordSys3D
  468. The other system to map the variables to.
  469. Examples
  470. ========
  471. >>> from sympy.vector import CoordSys3D
  472. >>> from sympy import Symbol
  473. >>> A = CoordSys3D('A')
  474. >>> q = Symbol('q')
  475. >>> B = A.orient_new_axis('B', q, A.k)
  476. >>> A.scalar_map(B)
  477. {A.x: B.x*cos(q) - B.y*sin(q), A.y: B.x*sin(q) + B.y*cos(q), A.z: B.z}
  478. """
  479. origin_coords = tuple(self.position_wrt(other).to_matrix(other))
  480. relocated_scalars = [x - origin_coords[i]
  481. for i, x in enumerate(other.base_scalars())]
  482. vars_matrix = (self.rotation_matrix(other) *
  483. Matrix(relocated_scalars))
  484. return {x: trigsimp(vars_matrix[i])
  485. for i, x in enumerate(self.base_scalars())}
  486. def locate_new(self, name, position, vector_names=None,
  487. variable_names=None):
  488. """
  489. Returns a CoordSys3D with its origin located at the given
  490. position wrt this coordinate system's origin.
  491. Parameters
  492. ==========
  493. name : str
  494. The name of the new CoordSys3D instance.
  495. position : Vector
  496. The position vector of the new system's origin wrt this
  497. one.
  498. vector_names, variable_names : iterable(optional)
  499. Iterables of 3 strings each, with custom names for base
  500. vectors and base scalars of the new system respectively.
  501. Used for simple str printing.
  502. Examples
  503. ========
  504. >>> from sympy.vector import CoordSys3D
  505. >>> A = CoordSys3D('A')
  506. >>> B = A.locate_new('B', 10 * A.i)
  507. >>> B.origin.position_wrt(A.origin)
  508. 10*A.i
  509. """
  510. if variable_names is None:
  511. variable_names = self._variable_names
  512. if vector_names is None:
  513. vector_names = self._vector_names
  514. return CoordSys3D(name, location=position,
  515. vector_names=vector_names,
  516. variable_names=variable_names,
  517. parent=self)
  518. def orient_new(self, name, orienters, location=None,
  519. vector_names=None, variable_names=None):
  520. """
  521. Creates a new CoordSys3D oriented in the user-specified way
  522. with respect to this system.
  523. Please refer to the documentation of the orienter classes
  524. for more information about the orientation procedure.
  525. Parameters
  526. ==========
  527. name : str
  528. The name of the new CoordSys3D instance.
  529. orienters : iterable/Orienter
  530. An Orienter or an iterable of Orienters for orienting the
  531. new coordinate system.
  532. If an Orienter is provided, it is applied to get the new
  533. system.
  534. If an iterable is provided, the orienters will be applied
  535. in the order in which they appear in the iterable.
  536. location : Vector(optional)
  537. The location of the new coordinate system's origin wrt this
  538. system's origin. If not specified, the origins are taken to
  539. be coincident.
  540. vector_names, variable_names : iterable(optional)
  541. Iterables of 3 strings each, with custom names for base
  542. vectors and base scalars of the new system respectively.
  543. Used for simple str printing.
  544. Examples
  545. ========
  546. >>> from sympy.vector import CoordSys3D
  547. >>> from sympy import symbols
  548. >>> q0, q1, q2, q3 = symbols('q0 q1 q2 q3')
  549. >>> N = CoordSys3D('N')
  550. Using an AxisOrienter
  551. >>> from sympy.vector import AxisOrienter
  552. >>> axis_orienter = AxisOrienter(q1, N.i + 2 * N.j)
  553. >>> A = N.orient_new('A', (axis_orienter, ))
  554. Using a BodyOrienter
  555. >>> from sympy.vector import BodyOrienter
  556. >>> body_orienter = BodyOrienter(q1, q2, q3, '123')
  557. >>> B = N.orient_new('B', (body_orienter, ))
  558. Using a SpaceOrienter
  559. >>> from sympy.vector import SpaceOrienter
  560. >>> space_orienter = SpaceOrienter(q1, q2, q3, '312')
  561. >>> C = N.orient_new('C', (space_orienter, ))
  562. Using a QuaternionOrienter
  563. >>> from sympy.vector import QuaternionOrienter
  564. >>> q_orienter = QuaternionOrienter(q0, q1, q2, q3)
  565. >>> D = N.orient_new('D', (q_orienter, ))
  566. """
  567. if variable_names is None:
  568. variable_names = self._variable_names
  569. if vector_names is None:
  570. vector_names = self._vector_names
  571. if isinstance(orienters, Orienter):
  572. if isinstance(orienters, AxisOrienter):
  573. final_matrix = orienters.rotation_matrix(self)
  574. else:
  575. final_matrix = orienters.rotation_matrix()
  576. # TODO: trigsimp is needed here so that the matrix becomes
  577. # canonical (scalar_map also calls trigsimp; without this, you can
  578. # end up with the same CoordinateSystem that compares differently
  579. # due to a differently formatted matrix). However, this is
  580. # probably not so good for performance.
  581. final_matrix = trigsimp(final_matrix)
  582. else:
  583. final_matrix = Matrix(eye(3))
  584. for orienter in orienters:
  585. if isinstance(orienter, AxisOrienter):
  586. final_matrix *= orienter.rotation_matrix(self)
  587. else:
  588. final_matrix *= orienter.rotation_matrix()
  589. return CoordSys3D(name, rotation_matrix=final_matrix,
  590. vector_names=vector_names,
  591. variable_names=variable_names,
  592. location=location,
  593. parent=self)
  594. def orient_new_axis(self, name, angle, axis, location=None,
  595. vector_names=None, variable_names=None):
  596. """
  597. Axis rotation is a rotation about an arbitrary axis by
  598. some angle. The angle is supplied as a SymPy expr scalar, and
  599. the axis is supplied as a Vector.
  600. Parameters
  601. ==========
  602. name : string
  603. The name of the new coordinate system
  604. angle : Expr
  605. The angle by which the new system is to be rotated
  606. axis : Vector
  607. The axis around which the rotation has to be performed
  608. location : Vector(optional)
  609. The location of the new coordinate system's origin wrt this
  610. system's origin. If not specified, the origins are taken to
  611. be coincident.
  612. vector_names, variable_names : iterable(optional)
  613. Iterables of 3 strings each, with custom names for base
  614. vectors and base scalars of the new system respectively.
  615. Used for simple str printing.
  616. Examples
  617. ========
  618. >>> from sympy.vector import CoordSys3D
  619. >>> from sympy import symbols
  620. >>> q1 = symbols('q1')
  621. >>> N = CoordSys3D('N')
  622. >>> B = N.orient_new_axis('B', q1, N.i + 2 * N.j)
  623. """
  624. if variable_names is None:
  625. variable_names = self._variable_names
  626. if vector_names is None:
  627. vector_names = self._vector_names
  628. orienter = AxisOrienter(angle, axis)
  629. return self.orient_new(name, orienter,
  630. location=location,
  631. vector_names=vector_names,
  632. variable_names=variable_names)
  633. def orient_new_body(self, name, angle1, angle2, angle3,
  634. rotation_order, location=None,
  635. vector_names=None, variable_names=None):
  636. """
  637. Body orientation takes this coordinate system through three
  638. successive simple rotations.
  639. Body fixed rotations include both Euler Angles and
  640. Tait-Bryan Angles, see https://en.wikipedia.org/wiki/Euler_angles.
  641. Parameters
  642. ==========
  643. name : string
  644. The name of the new coordinate system
  645. angle1, angle2, angle3 : Expr
  646. Three successive angles to rotate the coordinate system by
  647. rotation_order : string
  648. String defining the order of axes for rotation
  649. location : Vector(optional)
  650. The location of the new coordinate system's origin wrt this
  651. system's origin. If not specified, the origins are taken to
  652. be coincident.
  653. vector_names, variable_names : iterable(optional)
  654. Iterables of 3 strings each, with custom names for base
  655. vectors and base scalars of the new system respectively.
  656. Used for simple str printing.
  657. Examples
  658. ========
  659. >>> from sympy.vector import CoordSys3D
  660. >>> from sympy import symbols
  661. >>> q1, q2, q3 = symbols('q1 q2 q3')
  662. >>> N = CoordSys3D('N')
  663. A 'Body' fixed rotation is described by three angles and
  664. three body-fixed rotation axes. To orient a coordinate system D
  665. with respect to N, each sequential rotation is always about
  666. the orthogonal unit vectors fixed to D. For example, a '123'
  667. rotation will specify rotations about N.i, then D.j, then
  668. D.k. (Initially, D.i is same as N.i)
  669. Therefore,
  670. >>> D = N.orient_new_body('D', q1, q2, q3, '123')
  671. is same as
  672. >>> D = N.orient_new_axis('D', q1, N.i)
  673. >>> D = D.orient_new_axis('D', q2, D.j)
  674. >>> D = D.orient_new_axis('D', q3, D.k)
  675. Acceptable rotation orders are of length 3, expressed in XYZ or
  676. 123, and cannot have a rotation about about an axis twice in a row.
  677. >>> B = N.orient_new_body('B', q1, q2, q3, '123')
  678. >>> B = N.orient_new_body('B', q1, q2, 0, 'ZXZ')
  679. >>> B = N.orient_new_body('B', 0, 0, 0, 'XYX')
  680. """
  681. orienter = BodyOrienter(angle1, angle2, angle3, rotation_order)
  682. return self.orient_new(name, orienter,
  683. location=location,
  684. vector_names=vector_names,
  685. variable_names=variable_names)
  686. def orient_new_space(self, name, angle1, angle2, angle3,
  687. rotation_order, location=None,
  688. vector_names=None, variable_names=None):
  689. """
  690. Space rotation is similar to Body rotation, but the rotations
  691. are applied in the opposite order.
  692. Parameters
  693. ==========
  694. name : string
  695. The name of the new coordinate system
  696. angle1, angle2, angle3 : Expr
  697. Three successive angles to rotate the coordinate system by
  698. rotation_order : string
  699. String defining the order of axes for rotation
  700. location : Vector(optional)
  701. The location of the new coordinate system's origin wrt this
  702. system's origin. If not specified, the origins are taken to
  703. be coincident.
  704. vector_names, variable_names : iterable(optional)
  705. Iterables of 3 strings each, with custom names for base
  706. vectors and base scalars of the new system respectively.
  707. Used for simple str printing.
  708. See Also
  709. ========
  710. CoordSys3D.orient_new_body : method to orient via Euler
  711. angles
  712. Examples
  713. ========
  714. >>> from sympy.vector import CoordSys3D
  715. >>> from sympy import symbols
  716. >>> q1, q2, q3 = symbols('q1 q2 q3')
  717. >>> N = CoordSys3D('N')
  718. To orient a coordinate system D with respect to N, each
  719. sequential rotation is always about N's orthogonal unit vectors.
  720. For example, a '123' rotation will specify rotations about
  721. N.i, then N.j, then N.k.
  722. Therefore,
  723. >>> D = N.orient_new_space('D', q1, q2, q3, '312')
  724. is same as
  725. >>> B = N.orient_new_axis('B', q1, N.i)
  726. >>> C = B.orient_new_axis('C', q2, N.j)
  727. >>> D = C.orient_new_axis('D', q3, N.k)
  728. """
  729. orienter = SpaceOrienter(angle1, angle2, angle3, rotation_order)
  730. return self.orient_new(name, orienter,
  731. location=location,
  732. vector_names=vector_names,
  733. variable_names=variable_names)
  734. def orient_new_quaternion(self, name, q0, q1, q2, q3, location=None,
  735. vector_names=None, variable_names=None):
  736. """
  737. Quaternion orientation orients the new CoordSys3D with
  738. Quaternions, defined as a finite rotation about lambda, a unit
  739. vector, by some amount theta.
  740. This orientation is described by four parameters:
  741. q0 = cos(theta/2)
  742. q1 = lambda_x sin(theta/2)
  743. q2 = lambda_y sin(theta/2)
  744. q3 = lambda_z sin(theta/2)
  745. Quaternion does not take in a rotation order.
  746. Parameters
  747. ==========
  748. name : string
  749. The name of the new coordinate system
  750. q0, q1, q2, q3 : Expr
  751. The quaternions to rotate the coordinate system by
  752. location : Vector(optional)
  753. The location of the new coordinate system's origin wrt this
  754. system's origin. If not specified, the origins are taken to
  755. be coincident.
  756. vector_names, variable_names : iterable(optional)
  757. Iterables of 3 strings each, with custom names for base
  758. vectors and base scalars of the new system respectively.
  759. Used for simple str printing.
  760. Examples
  761. ========
  762. >>> from sympy.vector import CoordSys3D
  763. >>> from sympy import symbols
  764. >>> q0, q1, q2, q3 = symbols('q0 q1 q2 q3')
  765. >>> N = CoordSys3D('N')
  766. >>> B = N.orient_new_quaternion('B', q0, q1, q2, q3)
  767. """
  768. orienter = QuaternionOrienter(q0, q1, q2, q3)
  769. return self.orient_new(name, orienter,
  770. location=location,
  771. vector_names=vector_names,
  772. variable_names=variable_names)
  773. def create_new(self, name, transformation, variable_names=None, vector_names=None):
  774. """
  775. Returns a CoordSys3D which is connected to self by transformation.
  776. Parameters
  777. ==========
  778. name : str
  779. The name of the new CoordSys3D instance.
  780. transformation : Lambda, Tuple, str
  781. Transformation defined by transformation equations or chosen
  782. from predefined ones.
  783. vector_names, variable_names : iterable(optional)
  784. Iterables of 3 strings each, with custom names for base
  785. vectors and base scalars of the new system respectively.
  786. Used for simple str printing.
  787. Examples
  788. ========
  789. >>> from sympy.vector import CoordSys3D
  790. >>> a = CoordSys3D('a')
  791. >>> b = a.create_new('b', transformation='spherical')
  792. >>> b.transformation_to_parent()
  793. (b.r*sin(b.theta)*cos(b.phi), b.r*sin(b.phi)*sin(b.theta), b.r*cos(b.theta))
  794. >>> b.transformation_from_parent()
  795. (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))
  796. """
  797. return CoordSys3D(name, parent=self, transformation=transformation,
  798. variable_names=variable_names, vector_names=vector_names)
  799. def __init__(self, name, location=None, rotation_matrix=None,
  800. parent=None, vector_names=None, variable_names=None,
  801. latex_vects=None, pretty_vects=None, latex_scalars=None,
  802. pretty_scalars=None, transformation=None):
  803. # Dummy initializer for setting docstring
  804. pass
  805. __init__.__doc__ = __new__.__doc__
  806. @staticmethod
  807. def _compose_rotation_and_translation(rot, translation, parent):
  808. r = lambda x, y, z: CoordSys3D._rotation_trans_equations(rot, (x, y, z))
  809. if parent is None:
  810. return r
  811. dx, dy, dz = [translation.dot(i) for i in parent.base_vectors()]
  812. t = lambda x, y, z: (
  813. x + dx,
  814. y + dy,
  815. z + dz,
  816. )
  817. return lambda x, y, z: t(*r(x, y, z))
  818. def _check_strings(arg_name, arg):
  819. errorstr = arg_name + " must be an iterable of 3 string-types"
  820. if len(arg) != 3:
  821. raise ValueError(errorstr)
  822. for s in arg:
  823. if not isinstance(s, str):
  824. raise TypeError(errorstr)
  825. # Delayed import to avoid cyclic import problems:
  826. from sympy.vector.vector import BaseVector