[docs]classDyadic(object):"""A Dyadic object. See: http://en.wikipedia.org/wiki/Dyadic_tensor Kane, T., Levinson, D. Dynamics Theory and Applications. 1985 McGraw-Hill A more powerful way to represent a rigid body's inertia. While it is more complex, by choosing Dyadic components to be in body fixed basis vectors, the resulting matrix is equivalent to the inertia tensor. """def__init__(self,inlist):""" Just like Vector's init, you shouldn't call this unless creating a zero dyadic. zd = Dyadic(0) Stores a Dyadic as a list of lists; the inner list has the measure number and the two unit vectors; the outerlist holds each unique unit vector pair. """self.args=[]ifinlist==0:inlist=[]whilelen(inlist)!=0:added=0fori,vinenumerate(self.args):if((str(inlist[0][1])==str(self.args[i][1]))and(str(inlist[0][2])==str(self.args[i][2]))):self.args[i]=(self.args[i][0]+inlist[0][0],inlist[0][1],inlist[0][2])inlist.remove(inlist[0])added=1breakifadded!=1:self.args.append(inlist[0])inlist.remove(inlist[0])i=0# This code is to remove empty parts from the listwhilei<len(self.args):if((self.args[i][0]==0)|(self.args[i][1]==0)|(self.args[i][2]==0)):self.args.remove(self.args[i])i-=1i+=1def__add__(self,other):"""The add operator for Dyadic. """other=_check_dyadic(other)returnDyadic(self.args+other.args)def__and__(self,other):"""The inner product operator for a Dyadic and a Dyadic or Vector. Parameters ========== other : Dyadic or Vector The other Dyadic or Vector to take the inner product with Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, outer >>> N = ReferenceFrame('N') >>> D1 = outer(N.x, N.y) >>> D2 = outer(N.y, N.y) >>> D1.dot(D2) (N.x|N.y) >>> D1.dot(N.y) N.x """ifisinstance(other,Dyadic):other=_check_dyadic(other)ol=Dyadic(0)fori,vinenumerate(self.args):fori2,v2inenumerate(other.args):ol+=v[0]*v2[0]*(v[2]&v2[1])*(v[1]|v2[2])else:other=_check_vector(other)ol=Vector(0)fori,vinenumerate(self.args):ol+=v[0]*v[1]*(v[2]&other)returnoldef__div__(self,other):"""Divides the Dyadic by a sympifyable expression. """returnself.__mul__(1/other)__truediv__=__div__def__eq__(self,other):"""Tests for equality. Is currently weak; needs stronger comparison testing """ifother==0:other=Dyadic(0)other=_check_dyadic(other)if(self.args==[])and(other.args==[]):returnTrueelif(self.args==[])or(other.args==[]):returnFalsereturnset(self.args)==set(other.args)def__mul__(self,other):"""Multiplies the Dyadic by a sympifyable expression. Parameters ========== other : Sympafiable The scalar to multiply this Dyadic with Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, outer >>> N = ReferenceFrame('N') >>> d = outer(N.x, N.x) >>> 5 * d 5*(N.x|N.x) """newlist=[vforvinself.args]fori,vinenumerate(newlist):newlist[i]=(sympify(other)*newlist[i][0],newlist[i][1],newlist[i][2])returnDyadic(newlist)def__ne__(self,other):returnnotself.__eq__(other)def__neg__(self):returnself*-1def_latex(self,printer=None):ar=self.args# just to shorten thingsiflen(ar)==0:returnstr(0)ol=[]# output list, to be concatenated to a stringmlp=MechanicsLatexPrinter()fori,vinenumerate(ar):# if the coef of the dyadic is 1, we skip the 1ifar[i][0]==1:ol.append(' + '+mlp.doprint(ar[i][1])+r"\otimes "+mlp.doprint(ar[i][2]))# if the coef of the dyadic is -1, we skip the 1elifar[i][0]==-1:ol.append(' - '+mlp.doprint(ar[i][1])+r"\otimes "+mlp.doprint(ar[i][2]))# If the coefficient of the dyadic is not 1 or -1,# we might wrap it in parentheses, for readability.elifar[i][0]!=0:arg_str=mlp.doprint(ar[i][0])ifisinstance(ar[i][0],Add):arg_str='(%s)'%arg_strifarg_str.startswith('-'):arg_str=arg_str[1:]str_start=' - 'else:str_start=' + 'ol.append(str_start+arg_str+r" "+mlp.doprint(ar[i][1])+r"\otimes "+mlp.doprint(ar[i][2]))outstr=''.join(ol)ifoutstr.startswith(' + '):outstr=outstr[3:]elifoutstr.startswith(' '):outstr=outstr[1:]returnoutstrdef_pretty(self,printer=None):e=selfclassFake(object):baseline=0defrender(self,*args,**kwargs):self=ear=self.args# just to shorten thingsmpp=MechanicsPrettyPrinter()iflen(ar)==0:returnunicode(0)ol=[]# output list, to be concatenated to a stringfori,vinenumerate(ar):# if the coef of the dyadic is 1, we skip the 1ifar[i][0]==1:ol.append(u(" + ")+mpp.doprint(ar[i][1])+u("\u2a02 ")+mpp.doprint(ar[i][2]))# if the coef of the dyadic is -1, we skip the 1elifar[i][0]==-1:ol.append(u(" - ")+mpp.doprint(ar[i][1])+u("\u2a02 ")+mpp.doprint(ar[i][2]))# If the coefficient of the dyadic is not 1 or -1,# we might wrap it in parentheses, for readability.elifar[i][0]!=0:arg_str=mpp.doprint(ar[i][0])ifisinstance(ar[i][0],Add):arg_str=u("(%s)")%arg_strifarg_str.startswith(u("-")):arg_str=arg_str[1:]str_start=u(" - ")else:str_start=u(" + ")ol.append(str_start+arg_str+u(" ")+mpp.doprint(ar[i][1])+u("\u2a02 ")+mpp.doprint(ar[i][2]))outstr=u("").join(ol)ifoutstr.startswith(u(" + ")):outstr=outstr[3:]elifoutstr.startswith(" "):outstr=outstr[1:]returnoutstrreturnFake()def__rand__(self,other):"""The inner product operator for a Vector or Dyadic, and a Dyadic This is for: Vector dot Dyadic Parameters ========== other : Vector The vector we are dotting with Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, dot, outer >>> N = ReferenceFrame('N') >>> d = outer(N.x, N.x) >>> dot(N.x, d) N.x """other=_check_vector(other)ol=Vector(0)fori,vinenumerate(self.args):ol+=v[0]*v[2]*(v[1]&other)returnoldef__rsub__(self,other):return(-1*self)+otherdef__rxor__(self,other):"""For a cross product in the form: Vector x Dyadic Parameters ========== other : Vector The Vector that we are crossing this Dyadic with Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, outer, cross >>> N = ReferenceFrame('N') >>> d = outer(N.x, N.x) >>> cross(N.y, d) - (N.z|N.x) """other=_check_vector(other)ol=Dyadic(0)fori,vinenumerate(self.args):ol+=v[0]*((other^v[1])|v[2])returnoldef__str__(self,printer=None):"""Printing method. """ar=self.args# just to shorten thingsiflen(ar)==0:returnstr(0)ol=[]# output list, to be concatenated to a stringfori,vinenumerate(ar):# if the coef of the dyadic is 1, we skip the 1ifar[i][0]==1:ol.append(' + ('+str(ar[i][1])+'|'+str(ar[i][2])+')')# if the coef of the dyadic is -1, we skip the 1elifar[i][0]==-1:ol.append(' - ('+str(ar[i][1])+'|'+str(ar[i][2])+')')# If the coefficient of the dyadic is not 1 or -1,# we might wrap it in parentheses, for readability.elifar[i][0]!=0:arg_str=MechanicsStrPrinter().doprint(ar[i][0])ifisinstance(ar[i][0],Add):arg_str="(%s)"%arg_strifarg_str[0]=='-':arg_str=arg_str[1:]str_start=' - 'else:str_start=' + 'ol.append(str_start+arg_str+'*('+str(ar[i][1])+'|'+str(ar[i][2])+')')outstr=''.join(ol)ifoutstr.startswith(' + '):outstr=outstr[3:]elifoutstr.startswith(' '):outstr=outstr[1:]returnoutstrdef__sub__(self,other):"""The subtraction operator. """returnself.__add__(other*-1)def__xor__(self,other):"""For a cross product in the form: Dyadic x Vector. Parameters ========== other : Vector The Vector that we are crossing this Dyadic with Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, outer, cross >>> N = ReferenceFrame('N') >>> d = outer(N.x, N.x) >>> cross(d, N.y) (N.x|N.z) """other=_check_vector(other)ol=Dyadic(0)fori,vinenumerate(self.args):ol+=v[0]*(v[1]|(v[2]^other))returnol_sympystr=__str___sympyrepr=_sympystr__repr__=__str____radd__=__add____rmul__=__mul__

[docs]defexpress(self,frame1,frame2=None):"""Expresses this Dyadic in alternate frame(s) The first frame is the list side expression, the second frame is the right side; if Dyadic is in form A.x|B.y, you can express it in two different frames. If no second frame is given, the Dyadic is expressed in only one frame. Calls the global express function Parameters ========== frame1 : ReferenceFrame The frame to express the left side of the Dyadic in frame2 : ReferenceFrame If provided, the frame to express the right side of the Dyadic in Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, outer, dynamicsymbols >>> N = ReferenceFrame('N') >>> q = dynamicsymbols('q') >>> B = N.orientnew('B', 'Axis', [q, N.z]) >>> d = outer(N.x, N.x) >>> d.express(B, N) cos(q)*(B.x|N.x) - sin(q)*(B.y|N.x) """fromsympy.physics.mechanicsimportexpressreturnexpress(self,frame1,frame2)

[docs]defdoit(self,**hints):"""Calls .doit() on each term in the Dyadic"""returnsum([Dyadic([(v[0].doit(**hints),v[1],v[2])])forvinself.args],Dyadic(0))

[docs]classCoordinateSym(Symbol):""" A coordinate symbol/base scalar associated wrt a Reference Frame. Ideally, users should not instantiate this class. Instances of this class must only be accessed through the corresponding frame as 'frame[index]'. CoordinateSyms having the same frame and index parameters are equal (even though they may be instantiated separately). Parameters ========== name : string The display name of the CoordinateSym frame : ReferenceFrame The reference frame this base scalar belongs to index : 0, 1 or 2 The index of the dimension denoted by this coordinate variable Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, CoordinateSym >>> A = ReferenceFrame('A') >>> A[1] A_y >>> type(A[0]) <class 'sympy.physics.mechanics.essential.CoordinateSym'> >>> a_y = CoordinateSym('a_y', A, 1) >>> a_y == A[1] True """def__new__(cls,name,frame,index):obj=super(CoordinateSym,cls).__new__(cls,name)_check_frame(frame)ifindexnotinrange(0,3):raiseValueError("Invalid index specified")obj._id=(frame,index)returnobj@propertydefframe(self):returnself._id[0]def__eq__(self,other):#Check if the other object is a CoordinateSym of the same frame#and same indexifisinstance(other,CoordinateSym):ifother._id==self._id:returnTruereturnFalsedef__ne__(self,other):returnnotself.__eq__(other)def__hash__(self):returntuple((self._id[0].__hash__(),self._id[1])).__hash__()

[docs]classReferenceFrame(object):"""A reference frame in classical mechanics. ReferenceFrame is a class used to represent a reference frame in classical mechanics. It has a standard basis of three unit vectors in the frame's x, y, and z directions. It also can have a rotation relative to a parent frame; this rotation is defined by a direction cosine matrix relating this frame's basis vectors to the parent frame's basis vectors. It can also have an angular velocity vector, defined in another frame. """def__init__(self,name,indices=None,latexs=None,variables=None):"""ReferenceFrame initialization method. A ReferenceFrame has a set of orthonormal basis vectors, along with orientations relative to other ReferenceFrames and angular velocities relative to other ReferenceFrames. Parameters ========== indices : list (of strings) If custom indices are desired for console, pretty, and LaTeX printing, supply three as a list. The basis vectors can then be accessed with the get_item method. latexs : list (of strings) If custom names are desired for LaTeX printing of each basis vector, supply the names here in a list. Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, mlatex >>> N = ReferenceFrame('N') >>> N.x N.x >>> O = ReferenceFrame('O', indices=('1', '2', '3')) >>> O.x O['1'] >>> O['1'] O['1'] >>> P = ReferenceFrame('P', latexs=('A1', 'A2', 'A3')) >>> mlatex(P.x) 'A1' """ifnotisinstance(name,string_types):raiseTypeError('Need to supply a valid name')# The if statements below are for custom printing of basis-vectors for# each frame.# First case, when custom indices are suppliedifindicesisnotNone:ifnotisinstance(indices,(tuple,list)):raiseTypeError('Supply the indices as a list')iflen(indices)!=3:raiseValueError('Supply 3 indices')foriinindices:ifnotisinstance(i,string_types):raiseTypeError('Indices must be strings')self.str_vecs=[(name+'[\''+indices[0]+'\']'),(name+'[\''+indices[1]+'\']'),(name+'[\''+indices[2]+'\']')]self.pretty_vecs=[(u("\033[94m\033[1m")+name.lower()+u("_")+indices[0]+u("\033[0;0m\x1b[0;0m")),(u("\033[94m\033[1m")+name.lower()+u("_")+indices[1]+u("\033[0;0m\x1b[0;0m")),(u("\033[94m\033[1m")+name.lower()+u("_")+indices[2]+u("\033[0;0m\x1b[0;0m"))]self.latex_vecs=[(r"\mathbf{\hat{%s}_{%s}}"%(name.lower(),indices[0])),(r"\mathbf{\hat{%s}_{%s}}"%(name.lower(),indices[1])),(r"\mathbf{\hat{%s}_{%s}}"%(name.lower(),indices[2]))]self.indices=indices# Second case, when no custom indices are suppliedelse:self.str_vecs=[(name+'.x'),(name+'.y'),(name+'.z')]self.pretty_vecs=[(u("\033[94m\033[1m")+name.lower()+u("_x\033[0;0m\x1b[0;0m")),(u("\033[94m\033[1m")+name.lower()+u("_y\033[0;0m\x1b[0;0m")),(u("\033[94m\033[1m")+name.lower()+u("_z\033[0;0m\x1b[0;0m"))]self.latex_vecs=[(r"\mathbf{\hat{%s}_x}"%name.lower()),(r"\mathbf{\hat{%s}_y}"%name.lower()),(r"\mathbf{\hat{%s}_z}"%name.lower())]self.indices=['x','y','z']# Different step, for custom latex basis vectorsiflatexsisnotNone:ifnotisinstance(latexs,(tuple,list)):raiseTypeError('Supply the indices as a list')iflen(latexs)!=3:raiseValueError('Supply 3 indices')foriinlatexs:ifnotisinstance(i,string_types):raiseTypeError('Latex entries must be strings')self.latex_vecs=latexsself.name=nameself._var_dict={}#The _dcm_dict dictionary will only store the dcms of parent-child#relationships. The _dcm_cache dictionary will work as the dcm#cache.self._dcm_dict={}self._dcm_cache={}self._ang_vel_dict={}self._ang_acc_dict={}self._dlist=[self._dcm_dict,self._ang_vel_dict,self._ang_acc_dict]self._cur=0self._x=Vector([(Matrix([1,0,0]),self)])self._y=Vector([(Matrix([0,1,0]),self)])self._z=Vector([(Matrix([0,0,1]),self)])#Associate coordinate symbols wrt this frameifvariablesisnotNone:ifnotisinstance(variables,(tuple,list)):raiseTypeError('Supply the variable names as a list/tuple')iflen(variables)!=3:raiseValueError('Supply 3 variable names')foriinvariables:ifnotisinstance(i,string_types):raiseTypeError('Variable names must be strings')else:variables=[name+'_x',name+'_y',name+'_z']self.varlist=(CoordinateSym(variables[0],self,0), \
CoordinateSym(variables[1],self,1), \
CoordinateSym(variables[2],self,2))def__getitem__(self,ind):""" Returns basis vector for the provided index, if the index is a string. If the index is a number, returns the coordinate variable correspon- -ding to that index. """ifnotisinstance(ind,str):ifind<3:returnself.varlist[ind]else:raiseValueError("Invalid index provided")ifself.indices[0]==ind:returnself.xifself.indices[1]==ind:returnself.yifself.indices[2]==ind:returnself.zelse:raiseValueError('Not a defined index')def__iter__(self):returniter([self.x,self.y,self.z])def__str__(self):"""Returns the name of the frame. """returnself.name__repr__=__str__def_dict_list(self,other,num):"""Creates a list from self to other using _dcm_dict. """outlist=[[self]]oldlist=[[]]whileoutlist!=oldlist:oldlist=outlist[:]fori,vinenumerate(outlist):templist=v[-1]._dlist[num].keys()fori2,v2inenumerate(templist):ifnotv.__contains__(v2):littletemplist=v+[v2]ifnotoutlist.__contains__(littletemplist):outlist.append(littletemplist)fori,vinenumerate(oldlist):ifv[-1]!=other:outlist.remove(v)outlist.sort(key=len)iflen(outlist)!=0:returnoutlist[0]raiseValueError('No Connecting Path found between '+self.name+' and '+other.name)def_w_diff_dcm(self,otherframe):"""Angular velocity from time differentiating the DCM. """dcm2diff=self.dcm(otherframe)diffed=dcm2diff.diff(dynamicsymbols._t)angvelmat=diffed*dcm2diff.Tw1=trigsimp(expand(angvelmat[7]),recursive=True)w2=trigsimp(expand(angvelmat[2]),recursive=True)w3=trigsimp(expand(angvelmat[3]),recursive=True)return-Vector([(Matrix([w1,w2,w3]),self)])

[docs]defvariable_map(self,otherframe):""" Returns a dictionary which expresses the coordinate variables of this frame in terms of the variables of otherframe. If Vector.simp is True, returns a simplified version of the mapped values. Else, returns them without simplification. Simplification of the expressions may take time. Parameters ========== otherframe : ReferenceFrame The other frame to map the variables to Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, dynamicsymbols >>> A = ReferenceFrame('A') >>> q = dynamicsymbols('q') >>> B = A.orientnew('B', 'Axis', [q, A.z]) >>> A.variable_map(B) {A_x: B_x*cos(q(t)) - B_y*sin(q(t)), A_y: B_x*sin(q(t)) + B_y*cos(q(t)), A_z: B_z} """_check_frame(otherframe)if(otherframe,Vector.simp)inself._var_dict:returnself._var_dict[(otherframe,Vector.simp)]else:vars_matrix=self.dcm(otherframe)*Matrix(otherframe.varlist)mapping={}fori,xinenumerate(self):ifVector.simp:mapping[self.varlist[i]]=trigsimp(vars_matrix[i],method='fu')else:mapping[self.varlist[i]]=vars_matrix[i]self._var_dict[(otherframe,Vector.simp)]=mappingreturnmapping

[docs]deforient(self,parent,rot_type,amounts,rot_order=''):"""Defines the orientation of this frame relative to a parent frame. Parameters ========== parent : ReferenceFrame The frame that this ReferenceFrame will have its orientation matrix defined in relation to. rot_type : str The type of orientation matrix that is being created. Supported types are 'Body', 'Space', 'Quaternion', and 'Axis'. See examples for correct usage. amounts : list OR value The quantities that the orientation matrix will be defined by. rot_order : str If applicable, the order of a series of rotations. Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, Vector >>> from sympy import symbols >>> q0, q1, q2, q3, q4 = symbols('q0 q1 q2 q3 q4') >>> N = ReferenceFrame('N') >>> B = ReferenceFrame('B') Now we have a choice of how to implement the orientation. First is Body. Body orientation takes this reference frame through three successive simple rotations. Acceptable rotation orders are of length 3, expressed in XYZ or 123, and cannot have a rotation about about an axis twice in a row. >>> B.orient(N, 'Body', [q1, q2, q3], '123') >>> B.orient(N, 'Body', [q1, q2, 0], 'ZXZ') >>> B.orient(N, 'Body', [0, 0, 0], 'XYX') Next is Space. Space is like Body, but the rotations are applied in the opposite order. >>> B.orient(N, 'Space', [q1, q2, q3], '312') Next is Quaternion. This orients the new ReferenceFrame with Quaternions, defined as a finite rotation about lambda, a unit vector, by some amount theta. This orientation is described by four parameters: q0 = cos(theta/2) q1 = lambda_x sin(theta/2) q2 = lambda_y sin(theta/2) q3 = lambda_z sin(theta/2) Quaternion does not take in a rotation order. >>> B.orient(N, 'Quaternion', [q0, q1, q2, q3]) Last is Axis. This is a rotation about an arbitrary, non-time-varying axis by some angle. The axis is supplied as a Vector. This is how simple rotations are defined. >>> B.orient(N, 'Axis', [q1, N.x + 2 * N.y]) """_check_frame(parent)amounts=list(amounts)fori,vinenumerate(amounts):ifnotisinstance(v,Vector):amounts[i]=sympify(v)def_rot(axis,angle):"""DCM for simple axis 1,2,or 3 rotations. """ifaxis==1:returnMatrix([[1,0,0],[0,cos(angle),-sin(angle)],[0,sin(angle),cos(angle)]])elifaxis==2:returnMatrix([[cos(angle),0,sin(angle)],[0,1,0],[-sin(angle),0,cos(angle)]])elifaxis==3:returnMatrix([[cos(angle),-sin(angle),0],[sin(angle),cos(angle),0],[0,0,1]])approved_orders=('123','231','312','132','213','321','121','131','212','232','313','323','')rot_order=str(rot_order).upper()# Now we need to make sure XYZ = 123rot_type=rot_type.upper()rot_order=[i.replace('X','1')foriinrot_order]rot_order=[i.replace('Y','2')foriinrot_order]rot_order=[i.replace('Z','3')foriinrot_order]rot_order=''.join(rot_order)ifnotrot_orderinapproved_orders:raiseTypeError('The supplied order is not an approved type')parent_orient=[]ifrot_type=='AXIS':ifnotrot_order=='':raiseTypeError('Axis orientation takes no rotation order')ifnot(isinstance(amounts,(list,tuple))&(len(amounts)==2)):raiseTypeError('Amounts are a list or tuple of length 2')theta=amounts[0]axis=amounts[1]axis=_check_vector(axis)ifnotaxis.dt(parent)==0:raiseValueError('Axis cannot be time-varying')axis=axis.express(parent).normalize()axis=axis.args[0][0]parent_orient=((eye(3)-axis*axis.T)*cos(theta)+Matrix([[0,-axis[2],axis[1]],[axis[2],0,-axis[0]],[-axis[1],axis[0],0]])*sin(theta)+axis*axis.T)elifrot_type=='QUATERNION':ifnotrot_order=='':raiseTypeError('Quaternion orientation takes no rotation order')ifnot(isinstance(amounts,(list,tuple))&(len(amounts)==4)):raiseTypeError('Amounts are a list or tuple of length 4')q0,q1,q2,q3=amountsparent_orient=(Matrix([[q0**2+q1**2-q2**2-q3**2,2*(q1*q2-q0*q3),2*(q0*q2+q1*q3)],[2*(q1*q2+q0*q3),q0**2-q1**2+q2**2-q3**2,2*(q2*q3-q0*q1)],[2*(q1*q3-q0*q2),2*(q0*q1+q2*q3),q0**2-q1**2-q2**2+q3**2]]))elifrot_type=='BODY':ifnot(len(amounts)==3&len(rot_order)==3):raiseTypeError('Body orientation takes 3 values & 3 orders')a1=int(rot_order[0])a2=int(rot_order[1])a3=int(rot_order[2])parent_orient=(_rot(a1,amounts[0])*_rot(a2,amounts[1])*_rot(a3,amounts[2]))elifrot_type=='SPACE':ifnot(len(amounts)==3&len(rot_order)==3):raiseTypeError('Space orientation takes 3 values & 3 orders')a1=int(rot_order[0])a2=int(rot_order[1])a3=int(rot_order[2])parent_orient=(_rot(a3,amounts[2])*_rot(a2,amounts[1])*_rot(a1,amounts[0]))else:raiseNotImplementedError('That is not an implemented rotation')#Reset the _dcm_cache of this frame, and remove it from the _dcm_caches#of the frames it is linked to. Also remove it from the _dcm_dict of#its parentframes=self._dcm_cache.keys()forframeinframes:ifframeinself._dcm_dict:delframe._dcm_dict[self]delframe._dcm_cache[self]#Add the dcm relationship to _dcm_dictself._dcm_dict=self._dlist[0]={}self._dcm_dict.update({parent:parent_orient.T})parent._dcm_dict.update({self:parent_orient})#Also update the dcm cache after resetting itself._dcm_cache={}self._dcm_cache.update({parent:parent_orient.T})parent._dcm_cache.update({self:parent_orient})ifrot_type=='QUATERNION':t=dynamicsymbols._tq0,q1,q2,q3=amountsq0d=diff(q0,t)q1d=diff(q1,t)q2d=diff(q2,t)q3d=diff(q3,t)w1=2*(q1d*q0+q2d*q3-q3d*q2-q0d*q1)w2=2*(q2d*q0+q3d*q1-q1d*q3-q0d*q2)w3=2*(q3d*q0+q1d*q2-q2d*q1-q0d*q3)wvec=Vector([(Matrix([w1,w2,w3]),self)])elifrot_type=='AXIS':thetad=(amounts[0]).diff(dynamicsymbols._t)wvec=thetad*amounts[1].express(parent).normalize()else:try:fromsympy.polys.polyerrorsimportCoercionFailedfromsympy.physics.mechanics.functionsimportkinematic_equationsq1,q2,q3=amountsu1,u2,u3=dynamicsymbols('u1, u2, u3')templist=kinematic_equations([u1,u2,u3],[q1,q2,q3],rot_type,rot_order)templist=[expand(i)foriintemplist]td=solve(templist,[u1,u2,u3])u1=expand(td[u1])u2=expand(td[u2])u3=expand(td[u3])wvec=u1*self.x+u2*self.y+u3*self.zexcept(CoercionFailed,AssertionError):wvec=self._w_diff_dcm(parent)self._ang_vel_dict.update({parent:wvec})parent._ang_vel_dict.update({self:-wvec})self._var_dict={}

[docs]deforientnew(self,newname,rot_type,amounts,rot_order='',variables=None,indices=None,latexs=None):"""Creates a new ReferenceFrame oriented with respect to this Frame. See ReferenceFrame.orient() for acceptable rotation types, amounts, and orders. Parent is going to be self. Parameters ========== newname : str The name for the new ReferenceFrame rot_type : str The type of orientation matrix that is being created. amounts : list OR value The quantities that the orientation matrix will be defined by. rot_order : str If applicable, the order of a series of rotations. Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, Vector >>> from sympy import symbols >>> q1 = symbols('q1') >>> N = ReferenceFrame('N') >>> A = N.orientnew('A', 'Axis', [q1, N.x]) .orient() documentation:\n ======================== """newframe=ReferenceFrame(newname,variables,indices,latexs)newframe.orient(self,rot_type,amounts,rot_order)returnnewframe

[docs]defx(self):"""The basis Vector for the ReferenceFrame, in the x direction. """returnself._x

@property

[docs]defy(self):"""The basis Vector for the ReferenceFrame, in the y direction. """returnself._y

@property

[docs]defz(self):"""The basis Vector for the ReferenceFrame, in the z direction. """returnself._z

[docs]classVector(object):"""The class used to define vectors. It along with ReferenceFrame are the building blocks of describing a classical mechanics system in PyDy. Attributes ========== simp : Boolean Let certain methods use trigsimp on their outputs """simp=Falsedef__init__(self,inlist):"""This is the constructor for the Vector class. You shouldn't be calling this, it should only be used by other functions. You should be treating Vectors like you would with if you were doing the math by hand, and getting the first 3 from the standard basis vectors from a ReferenceFrame. The only exception is to create a zero vector: zv = Vector(0) """self.args=[]ifinlist==0:inlist=[]whilelen(inlist)!=0:added=0fori,vinenumerate(self.args):ifinlist[0][1]==self.args[i][1]:self.args[i]=(self.args[i][0]+inlist[0][0],inlist[0][1])inlist.remove(inlist[0])added=1breakifadded!=1:self.args.append(inlist[0])inlist.remove(inlist[0])i=0# This code is to remove empty frames from the listwhilei<len(self.args):ifself.args[i][0]==Matrix([0,0,0]):self.args.remove(self.args[i])i-=1i+=1def__hash__(self):returnhash(tuple(self.args))def__add__(self,other):"""The add operator for Vector. """other=_check_vector(other)returnVector(self.args+other.args)def__and__(self,other):"""Dot product of two vectors. Returns a scalar, the dot product of the two Vectors Parameters ========== other : Vector The Vector which we are dotting with Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, Vector, dot >>> from sympy import symbols >>> q1 = symbols('q1') >>> N = ReferenceFrame('N') >>> dot(N.x, N.x) 1 >>> dot(N.x, N.y) 0 >>> A = N.orientnew('A', 'Axis', [q1, N.x]) >>> dot(N.y, A.y) cos(q1) """ifisinstance(other,Dyadic):returnNotImplementedother=_check_vector(other)out=S(0)fori,v1inenumerate(self.args):forj,v2inenumerate(other.args):out+=((v2[0].T)*(v2[1].dcm(v1[1]))*(v1[0]))[0]ifVector.simp:returntrigsimp(sympify(out),recursive=True)else:returnsympify(out)def__div__(self,other):"""This uses mul and inputs self and 1 divided by other. """returnself.__mul__(1/other)__truediv__=__div__def__eq__(self,other):"""Tests for equality. It is very import to note that this is only as good as the SymPy equality test; False does not always mean they are not equivalent Vectors. If other is 0, and self is empty, returns True. If other is 0 and self is not empty, returns False. If none of the above, only accepts other as a Vector. """ifother==0:other=Vector(0)other=_check_vector(other)if(self.args==[])and(other.args==[]):returnTrueelif(self.args==[])or(other.args==[]):returnFalseframe=self.args[0][1]forvinframe:ifexpand((self-other)&v)!=0:returnFalsereturnTruedef__mul__(self,other):"""Multiplies the Vector by a sympifyable expression. Parameters ========== other : Sympifyable The scalar to multiply this Vector with Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, Vector >>> from sympy import Symbol >>> N = ReferenceFrame('N') >>> b = Symbol('b') >>> V = 10 * b * N.x >>> print(V) 10*b*N.x """newlist=[vforvinself.args]fori,vinenumerate(newlist):newlist[i]=(sympify(other)*newlist[i][0],newlist[i][1])returnVector(newlist)def__ne__(self,other):returnnotself.__eq__(other)def__neg__(self):returnself*-1def__or__(self,other):"""Outer product between two Vectors. A rank increasing operation, which returns a Dyadic from two Vectors Parameters ========== other : Vector The Vector to take the outer product with Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, outer >>> N = ReferenceFrame('N') >>> outer(N.x, N.x) (N.x|N.x) """other=_check_vector(other)ol=Dyadic(0)fori,vinenumerate(self.args):fori2,v2inenumerate(other.args):# it looks this way because if we are in the same frame and# use the enumerate function on the same frame in a nested# fashion, then bad things happenol+=Dyadic([(v[0][0]*v2[0][0],v[1].x,v2[1].x)])ol+=Dyadic([(v[0][0]*v2[0][1],v[1].x,v2[1].y)])ol+=Dyadic([(v[0][0]*v2[0][2],v[1].x,v2[1].z)])ol+=Dyadic([(v[0][1]*v2[0][0],v[1].y,v2[1].x)])ol+=Dyadic([(v[0][1]*v2[0][1],v[1].y,v2[1].y)])ol+=Dyadic([(v[0][1]*v2[0][2],v[1].y,v2[1].z)])ol+=Dyadic([(v[0][2]*v2[0][0],v[1].z,v2[1].x)])ol+=Dyadic([(v[0][2]*v2[0][1],v[1].z,v2[1].y)])ol+=Dyadic([(v[0][2]*v2[0][2],v[1].z,v2[1].z)])returnoldef_latex(self,printer=None):"""Latex Printing method. """ar=self.args# just to shorten thingsiflen(ar)==0:returnstr(0)ol=[]# output list, to be concatenated to a stringfori,vinenumerate(ar):forjin0,1,2:# if the coef of the basis vector is 1, we skip the 1ifar[i][0][j]==1:ol.append(' + '+ar[i][1].latex_vecs[j])# if the coef of the basis vector is -1, we skip the 1elifar[i][0][j]==-1:ol.append(' - '+ar[i][1].latex_vecs[j])elifar[i][0][j]!=0:# If the coefficient of the basis vector is not 1 or -1;# also, we might wrap it in parentheses, for readability.arg_str=MechanicsStrPrinter().doprint(ar[i][0][j])ifisinstance(ar[i][0][j],Add):arg_str="(%s)"%arg_strifarg_str[0]=='-':arg_str=arg_str[1:]str_start=' - 'else:str_start=' + 'ol.append(str_start+arg_str+'*'+ar[i][1].latex_vecs[j])outstr=''.join(ol)ifoutstr.startswith(' + '):outstr=outstr[3:]elifoutstr.startswith(' '):outstr=outstr[1:]returnoutstrdef_pretty(self,printer=None):"""Pretty Printing method. """e=selfclassFake(object):baseline=0defrender(self,*args,**kwargs):self=ear=self.args# just to shorten thingsiflen(ar)==0:returnunicode(0)ol=[]# output list, to be concatenated to a stringfori,vinenumerate(ar):forjin0,1,2:# if the coef of the basis vector is 1, we skip the 1ifar[i][0][j]==1:ol.append(u(" + ")+ar[i][1].pretty_vecs[j])# if the coef of the basis vector is -1, we skip the 1elifar[i][0][j]==-1:ol.append(u(" - ")+ar[i][1].pretty_vecs[j])elifar[i][0][j]!=0:# If the basis vector coeff is not 1 or -1,# we might wrap it in parentheses, for readability.arg_str=(MechanicsPrettyPrinter().doprint(ar[i][0][j]))ifisinstance(ar[i][0][j],Add):arg_str=u("(%s)")%arg_strifarg_str[0]==u("-"):arg_str=arg_str[1:]str_start=u(" - ")else:str_start=u(" + ")ol.append(str_start+arg_str+'*'+ar[i][1].pretty_vecs[j])outstr=u("").join(ol)ifoutstr.startswith(u(" + ")):outstr=outstr[3:]elifoutstr.startswith(" "):outstr=outstr[1:]returnoutstrreturnFake()def__ror__(self,other):"""Outer product between two Vectors. A rank increasing operation, which returns a Dyadic from two Vectors Parameters ========== other : Vector The Vector to take the outer product with Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, outer >>> N = ReferenceFrame('N') >>> outer(N.x, N.x) (N.x|N.x) """other=_check_vector(other)ol=Dyadic(0)fori,vinenumerate(other.args):fori2,v2inenumerate(self.args):# it looks this way because if we are in the same frame and# use the enumerate function on the same frame in a nested# fashion, then bad things happenol+=Dyadic([(v[0][0]*v2[0][0],v[1].x,v2[1].x)])ol+=Dyadic([(v[0][0]*v2[0][1],v[1].x,v2[1].y)])ol+=Dyadic([(v[0][0]*v2[0][2],v[1].x,v2[1].z)])ol+=Dyadic([(v[0][1]*v2[0][0],v[1].y,v2[1].x)])ol+=Dyadic([(v[0][1]*v2[0][1],v[1].y,v2[1].y)])ol+=Dyadic([(v[0][1]*v2[0][2],v[1].y,v2[1].z)])ol+=Dyadic([(v[0][2]*v2[0][0],v[1].z,v2[1].x)])ol+=Dyadic([(v[0][2]*v2[0][1],v[1].z,v2[1].y)])ol+=Dyadic([(v[0][2]*v2[0][2],v[1].z,v2[1].z)])returnoldef__rsub__(self,other):return(-1*self)+otherdef__str__(self,printer=None):"""Printing method. """ar=self.args# just to shorten thingsiflen(ar)==0:returnstr(0)ol=[]# output list, to be concatenated to a stringfori,vinenumerate(ar):forjin0,1,2:# if the coef of the basis vector is 1, we skip the 1ifar[i][0][j]==1:ol.append(' + '+ar[i][1].str_vecs[j])# if the coef of the basis vector is -1, we skip the 1elifar[i][0][j]==-1:ol.append(' - '+ar[i][1].str_vecs[j])elifar[i][0][j]!=0:# If the coefficient of the basis vector is not 1 or -1;# also, we might wrap it in parentheses, for readability.arg_str=MechanicsStrPrinter().doprint(ar[i][0][j])ifisinstance(ar[i][0][j],Add):arg_str="(%s)"%arg_strifarg_str[0]=='-':arg_str=arg_str[1:]str_start=' - 'else:str_start=' + 'ol.append(str_start+arg_str+'*'+ar[i][1].str_vecs[j])outstr=''.join(ol)ifoutstr.startswith(' + '):outstr=outstr[3:]elifoutstr.startswith(' '):outstr=outstr[1:]returnoutstrdef__sub__(self,other):"""The subraction operator. """returnself.__add__(other*-1)def__xor__(self,other):"""The cross product operator for two Vectors. Returns a Vector, expressed in the same ReferenceFrames as self. Parameters ========== other : Vector The Vector which we are crossing with Examples ======== >>> from sympy.physics.mechanics import ReferenceFrame, Vector >>> from sympy import symbols >>> q1 = symbols('q1') >>> N = ReferenceFrame('N') >>> N.x ^ N.y N.z >>> A = N.orientnew('A', 'Axis', [q1, N.x]) >>> A.x ^ N.y N.z >>> N.y ^ A.x - sin(q1)*A.y - cos(q1)*A.z """ifisinstance(other,Dyadic):returnNotImplementedother=_check_vector(other)ifother.args==[]:returnVector(0)def_det(mat):"""This is needed as a little method for to find the determinant of a list in python; needs to work for a 3x3 list. SymPy's Matrix won't take in Vector, so need a custom function. You shouldn't be calling this. """return(mat[0][0]*(mat[1][1]*mat[2][2]-mat[1][2]*mat[2][1])+mat[0][1]*(mat[1][2]*mat[2][0]-mat[1][0]*mat[2][2])+mat[0][2]*(mat[1][0]*mat[2][1]-mat[1][1]*mat[2][0]))outvec=Vector(0)ar=other.args# For brevityfori,vinenumerate(ar):tempx=v[1].xtempy=v[1].ytempz=v[1].ztempm=([[tempx,tempy,tempz],[self&tempx,self&tempy,self&tempz],[Vector([ar[i]])&tempx,Vector([ar[i]])&tempy,Vector([ar[i]])&tempz]])outvec+=_det(tempm)returnoutvec_sympystr=__str___sympyrepr=_sympystr__repr__=__str____radd__=__add____rand__=__and____rmul__=__mul__

[docs]defdoit(self,**hints):"""Calls .doit() on each term in the Vector"""ov=Vector(0)fori,vinenumerate(self.args):ov+=Vector([(v[0].applyfunc(lambdax:x.doit(**hints)),v[1])])returnov

[docs]defdt(self,otherframe):""" Returns a Vector which is the time derivative of the self Vector, taken in frame otherframe. Calls the global time_derivative method Parameters ========== otherframe : ReferenceFrame The frame to calculate the time derivative in """fromsympy.physics.mechanicsimporttime_derivativereturntime_derivative(self,otherframe)

[docs]defnormalize(self):"""Returns a Vector of magnitude 1, codirectional with self."""returnVector(self.args+[])/self.magnitude()

classMechanicsStrPrinter(StrPrinter):"""String Printer for mechanics. """def_print_Derivative(self,e):t=dynamicsymbols._tif(bool(sum([i==tforiine.variables]))&isinstance(type(e.args[0]),UndefinedFunction)):ol=str(e.args[0].func)fori,vinenumerate(e.variables):ol+=dynamicsymbols._strreturnolelse:returnStrPrinter().doprint(e)def_print_Function(self,e):t=dynamicsymbols._tifisinstance(type(e),UndefinedFunction):returnStrPrinter().doprint(e).replace("(%s)"%t,'')returne.func.__name__+"(%s)"%self.stringify(e.args,", ")classMechanicsLatexPrinter(LatexPrinter):"""Latex Printer for mechanics. """def_print_Function(self,expr,exp=None):func=expr.func.__name__t=dynamicsymbols._tifhasattr(self,'_print_'+func):returngetattr(self,'_print_'+func)(expr,exp)elifisinstance(type(expr),UndefinedFunction)and(expr.args==(t,)):name,sup,sub=split_super_sub(func)iflen(sup)!=0:sup=r"^{%s}"%"".join(sup)else:sup=r""iflen(sub)!=0:sub=r"_{%s}"%"".join(sub)else:sub=r""ifexp:sup+=r"^{%s}"%self._print(exp)returnr"%s"%(name+sup+sub)else:args=[str(self._print(arg))forarginexpr.args]# How inverse trig functions should be displayed, formats are:# abbreviated: asin, full: arcsin, power: sin^-1inv_trig_style=self._settings['inv_trig_style']# If we are dealing with a power-style inverse trig functioninv_trig_power_case=False# If it is applicable to fold the argument bracketscan_fold_brackets=self._settings['fold_func_brackets']and \
len(args)==1and \
notself._needs_function_brackets(expr.args[0])inv_trig_table=["asin","acos","atan","acot"]# If the function is an inverse trig function, handle the styleiffuncininv_trig_table:ifinv_trig_style=="abbreviated":func=funcelifinv_trig_style=="full":func="arc"+func[1:]elifinv_trig_style=="power":func=func[1:]inv_trig_power_case=True# Can never fold brackets if we're raised to a powerifexpisnotNone:can_fold_brackets=Falseifinv_trig_power_case:name=r"\operatorname{%s}^{-1}"%funcelifexpisnotNone:name=r"\operatorname{%s}^{%s}"%(func,exp)else:name=r"\operatorname{%s}"%funcifcan_fold_brackets:name+=r"%s"else:name+=r"\left(%s\right)"ifinv_trig_power_caseandexpisnotNone:name+=r"^{%s}"%expreturnname%",".join(args)def_print_Derivative(self,der_expr):# make sure it is an the right formder_expr=der_expr.doit()ifnotisinstance(der_expr,Derivative):returnself.doprint(der_expr)# check if expr is a dynamicsymbolfromsympy.core.functionimportAppliedUndeft=dynamicsymbols._texpr=der_expr.exprred=expr.atoms(AppliedUndef)syms=der_expr.variablestest1=notall([Trueforiinredifi.atoms()==set([t])])test2=notall([(t==i)foriinsyms])iftest1ortest2:returnLatexPrinter().doprint(der_expr)# done checkingdots=len(syms)base=self._print_Function(expr)base_split=base.split('_',1)base=base_split[0]ifdots==1:base=r"\dot{%s}"%baseelifdots==2:base=r"\ddot{%s}"%baseelifdots==3:base=r"\dddot{%s}"%baseiflen(base_split)isnot1:base+='_'+base_split[1]returnbaseclassMechanicsPrettyPrinter(PrettyPrinter):"""Pretty Printer for mechanics. """def_print_Derivative(self,deriv):# XXX use U('PARTIAL DIFFERENTIAL') here ?t=dynamicsymbols._tdots=0can_break=Truesyms=list(reversed(deriv.variables))x=Nonewhilelen(syms)>0:ifsyms[-1]==t:syms.pop()dots+=1else:breakf=prettyForm(binding=prettyForm.FUNC,*self._print(deriv.expr))ifnot(isinstance(type(deriv.expr),UndefinedFunction)and(deriv.expr.args==(t,))):dots=0can_break=Falsef=prettyForm(binding=prettyForm.FUNC,*self._print(deriv.expr).parens())ifdots==0:dots=u("")elifdots==1:dots=u("\u0307")elifdots==2:dots=u("\u0308")elifdots==3:dots=u("\u20db")elifdots==4:dots=u("\u20dc")uni_subs=[u("\u2080"),u("\u2081"),u("\u2082"),u("\u2083"),u("\u2084"),u("\u2085"),u("\u2086"),u("\u2087"),u("\u2088"),u("\u2089"),u("\u208a"),u("\u208b"),u("\u208c"),u("\u208d"),u("\u208e"),u("\u208f"),u("\u2090"),u("\u2091"),u("\u2092"),u("\u2093"),u("\u2094"),u("\u2095"),u("\u2096"),u("\u2097"),u("\u2098"),u("\u2099"),u("\u209a"),u("\u209b"),u("\u209c"),u("\u209d"),u("\u209e"),u("\u209f")]fpic=f.__dict__['picture']funi=f.__dict__['unicode']ind=len(funi)val=""foriinuni_subs:cur_ind=funi.find(i)if(cur_ind!=-1)and(cur_ind<ind):ind=cur_indval=iifind==len(funi):funi+=dotselse:funi=funi.replace(val,dots+val)iff.__dict__['picture']==[f.__dict__['unicode']]:fpic=[funi]f.__dict__['picture']=fpicf.__dict__['unicode']=funiif(len(syms))==0andcan_break:returnfforsym,numingroup(syms,multiple=False):s=self._print(sym)ds=prettyForm(*s.left('d'))ifnum>1:ds=ds**prettyForm(str(num))ifxisNone:x=dselse:x=prettyForm(*x.right(' '))x=prettyForm(*x.right(ds))pform=prettyForm('d')iflen(syms)>1:pform=pform**prettyForm(str(len(syms)))pform=prettyForm(*pform.below(stringPict.LINE,x))pform.baseline=pform.baseline+1pform=prettyForm(*stringPict.next(pform,f))returnpformdef_print_Function(self,e):t=dynamicsymbols._t# XXX works only for applied functionsfunc=e.funcargs=e.argsfunc_name=func.__name__prettyFunc=self._print(C.Symbol(func_name))prettyArgs=prettyForm(*self._print_seq(args).parens())# If this function is an Undefined function of t, it is probably a# dynamic symbol, so we'll skip the (t). The rest of the code is# identical to the normal PrettyPrinter codeifisinstance(func,UndefinedFunction)and(args==(t,)):pform=prettyForm(binding=prettyForm.FUNC,*stringPict.next(prettyFunc))else:pform=prettyForm(binding=prettyForm.FUNC,*stringPict.next(prettyFunc,prettyArgs))# store pform parts so it can be reassembled e.g. when poweredpform.prettyFunc=prettyFuncpform.prettyArgs=prettyArgsreturnpformclassMechanicsTypeError(TypeError):def__init__(self,other,type_str):super(MechanicsTypeError,self).__init__("Expected an instance of %s, ""instead received an object '%s' of type %s."%(type_str,other,type(other)))def_check_dyadic(other):ifnotisinstance(other,Dyadic):raiseTypeError('A Dyadic must be supplied')returnotherdef_check_frame(other):ifnotisinstance(other,ReferenceFrame):raiseMechanicsTypeError(other,"ReferenceFrame")def_check_vector(other):ifnotisinstance(other,Vector):raiseTypeError('A Vector must be supplied')returnother

[docs]defdynamicsymbols(names,level=0):"""Uses symbols and Function for functions of time. Creates a SymPy UndefinedFunction, which is then initialized as a function of a variable, the default being Symbol('t'). Parameters ========== names : str Names of the dynamic symbols you want to create; works the same way as inputs to symbols level : int Level of differentiation of the returned function; d/dt once of t, twice of t, etc. Examples ======== >>> from sympy.physics.mechanics import dynamicsymbols >>> from sympy import diff, Symbol >>> q1 = dynamicsymbols('q1') >>> q1 q1(t) >>> diff(q1, Symbol('t')) Derivative(q1(t), t) """esses=symbols(names,cls=Function)t=dynamicsymbols._tifhasattr(esses,'__iter__'):esses=[reduce(diff,[t]*level,e(t))foreinesses]returnesseselse:returnreduce(diff,[t]*level,esses(t))