__all__ = ['cross',
'dot',
'express',
'outer',
'inertia',
'mechanics_printing',
'mprint',
'msprint',
'mpprint',
'mlatex',
'kinematic_equations',
'inertia_of_point_mass',
'partial_velocity',
'linear_momentum',
'angular_momentum',
'kinetic_energy',
'potential_energy',
'Lagrangian']
from sympy.physics.mechanics.essential import (Vector, Dyadic, ReferenceFrame,
MechanicsStrPrinter,
MechanicsPrettyPrinter,
MechanicsLatexPrinter,
dynamicsymbols)
from sympy.physics.mechanics.particle import Particle
from sympy.physics.mechanics.rigidbody import RigidBody
from sympy.physics.mechanics.point import Point
from sympy import sympify, diff, sin, cos, Matrix
from sympy.core.basic import S
[docs]def cross(vec1, vec2):
"""Cross product convenience wrapper for Vector.cross(): \n"""
if not isinstance(vec1, (Vector, Dyadic)):
raise TypeError('Cross product is between two vectors')
return vec1 ^ vec2
cross.__doc__ += Vector.cross.__doc__
[docs]def dot(vec1, vec2):
"""Dot product convenience wrapper for Vector.dot(): \n"""
if not isinstance(vec1, (Vector, Dyadic)):
raise TypeError('Dot product is between two vectors')
return vec1 & vec2
dot.__doc__ += Vector.dot.__doc__
[docs]def express(vec, frame, frame2=None):
"""Express convenience wrapper for Vector.express(): \n"""
if not isinstance(vec, (Vector, Dyadic)):
raise TypeError('Can only express Vectors')
if isinstance(vec, Vector):
return vec.express(frame)
else:
return vec.express(frame, frame2)
express.__doc__ += Vector.express.__doc__
[docs]def outer(vec1, vec2):
"""Outer product convenience wrapper for Vector.outer():\n"""
if not isinstance(vec1, Vector):
raise TypeError('Outer product is between two Vectors')
return vec1 | vec2
outer.__doc__ += Vector.express.__doc__
[docs]def inertia(frame, ixx, iyy, izz, ixy=0, iyz=0, izx=0):
"""Simple way to create inertia Dyadic object.
If you don't know what a Dyadic is, just treat this like the inertia
tensor. Then, do the easy thing and define it in a body-fixed frame.
Parameters
==========
frame : ReferenceFrame
The frame the inertia is defined in
ixx : Sympifyable
the xx element in the inertia dyadic
iyy : Sympifyable
the yy element in the inertia dyadic
izz : Sympifyable
the zz element in the inertia dyadic
ixy : Sympifyable
the xy element in the inertia dyadic
iyz : Sympifyable
the yz element in the inertia dyadic
izx : Sympifyable
the zx element in the inertia dyadic
Examples
========
>>> from sympy.physics.mechanics import ReferenceFrame, inertia
>>> N = ReferenceFrame('N')
>>> inertia(N, 1, 2, 3)
(N.x|N.x) + 2*(N.y|N.y) + 3*(N.z|N.z)
"""
if not isinstance(frame, ReferenceFrame):
raise TypeError('Need to define the inertia in a frame')
ol = sympify(ixx) * (frame.x | frame.x)
ol += sympify(ixy) * (frame.x | frame.y)
ol += sympify(izx) * (frame.x | frame.z)
ol += sympify(ixy) * (frame.y | frame.x)
ol += sympify(iyy) * (frame.y | frame.y)
ol += sympify(iyz) * (frame.y | frame.z)
ol += sympify(izx) * (frame.z | frame.x)
ol += sympify(iyz) * (frame.z | frame.y)
ol += sympify(izz) * (frame.z | frame.z)
return ol
[docs]def inertia_of_point_mass(mass, pos_vec, frame):
"""Inertia dyadic of a point mass realtive to point O.
Parameters
==========
mass : Sympifyable
Mass of the point mass
pos_vec : Vector
Position from point O to point mass
frame : ReferenceFrame
Reference frame to express the dyadic in
Examples
========
>>> from sympy import symbols
>>> from sympy.physics.mechanics import ReferenceFrame, inertia_of_point_mass
>>> N = ReferenceFrame('N')
>>> r, m = symbols('r m')
>>> px = r * N.x
>>> inertia_of_point_mass(m, px, N)
m*r**2*(N.y|N.y) + m*r**2*(N.z|N.z)
"""
return mass * (((frame.x | frame.x) + (frame.y | frame.y) +
(frame.z | frame.z)) * (pos_vec & pos_vec) -
(pos_vec | pos_vec))
[docs]def mechanics_printing():
"""Sets up interactive printing for mechanics' derivatives.
The main benefit of this is for printing of time derivatives;
instead of displaying as Derivative(f(t),t), it will display f'
This is only actually needed for when derivatives are present and are not
in a physics.mechanics object.
Examples
========
>>> # 2 lines below are for tests to function properly
>>> import sys
>>> sys.displayhook = sys.__displayhook__
>>> from sympy import Function, Symbol, diff
>>> from sympy.physics.mechanics import mechanics_printing
>>> f = Function('f')
>>> t = Symbol('t')
>>> x = Symbol('x')
>>> diff(f(t), t)
Derivative(f(t), t)
>>> mechanics_printing()
>>> diff(f(t), t)
f'
>>> diff(f(x), x)
Derivative(f(x), x)
>>> # 2 lines below are for tests to function properly
>>> import sys
>>> sys.displayhook = sys.__displayhook__
"""
import sys
sys.displayhook = mprint
[docs]def mprint(expr, **settings):
r"""Function for printing of expressions generated in mechanics.
Extends SymPy's StrPrinter; mprint is equivalent to:
print sstr()
mprint takes the same options as sstr.
Parameters
==========
expr : valid sympy object
SymPy expression to print
settings : args
Same as print for SymPy
Examples
========
>>> from sympy.physics.mechanics import mprint, dynamicsymbols
>>> u1 = dynamicsymbols('u1')
>>> print(u1)
u1(t)
>>> mprint(u1)
u1
"""
outstr = msprint(expr, **settings)
import __builtin__
if (outstr != 'None'):
__builtin__._ = outstr
print(outstr)
def msprint(expr, **settings):
r"""Function for displaying expressions generated in mechanics.
Returns the output of mprint() as a string.
Parameters
==========
expr : valid sympy object
SymPy expression to print
settings : args
Same as print for SymPy
Examples
========
>>> from sympy.physics.mechanics import msprint, dynamicsymbols
>>> u1, u2 = dynamicsymbols('u1 u2')
>>> u2d = dynamicsymbols('u2', level=1)
>>> print("%s = %s" % (u1, u2 + u2d))
u1(t) = u2(t) + Derivative(u2(t), t)
>>> print("%s = %s" % (msprint(u1), msprint(u2 + u2d)))
u1 = u2 + u2'
"""
pr = MechanicsStrPrinter(settings)
return pr.doprint(expr)
[docs]def mpprint(expr, **settings):
r"""Function for pretty printing of expressions generated in mechanics.
Mainly used for expressions not inside a vector; the output of running
scripts and generating equations of motion. Takes the same options as
SymPy's pretty_print(); see that function for more information.
Parameters
==========
expr : valid sympy object
SymPy expression to pretty print
settings : args
Same as pretty print
Examples
========
Use in the same way as pprint
"""
mp = MechanicsPrettyPrinter(settings)
print(mp.doprint(expr))
[docs]def mlatex(expr, **settings):
r"""Function for printing latex representation of mechanics objects.
For latex representation of Vectors, Dyadics, and dynamicsymbols. Takes the
same options as SymPy's latex(); see that function for more information;
Parameters
==========
expr : valid sympy object
SymPy expression to represent in LaTeX form
settings : args
Same as latex()
Examples
========
>>> from sympy.physics.mechanics import mlatex, ReferenceFrame, dynamicsymbols
>>> N = ReferenceFrame('N')
>>> q1, q2 = dynamicsymbols('q1 q2')
>>> q1d, q2d = dynamicsymbols('q1 q2', 1)
>>> q1dd, q2dd = dynamicsymbols('q1 q2', 2)
>>> mlatex(N.x + N.y)
'\\mathbf{\\hat{n}_x} + \\mathbf{\\hat{n}_y}'
>>> mlatex(q1 + q2)
'q_{1} + q_{2}'
>>> mlatex(q1d)
'\\dot{q}_{1}'
>>> mlatex(q1 * q2d)
'q_{1} \\dot{q}_{2}'
>>> mlatex(q1dd * q1 / q1d)
'\\frac{q_{1} \\ddot{q}_{1}}{\\dot{q}_{1}}'
"""
return MechanicsLatexPrinter(settings).doprint(expr)
[docs]def kinematic_equations(speeds, coords, rot_type, rot_order=''):
"""Gives equations relating the qdot's to u's for a rotation type.
Supply rotation type and order as in orient. Speeds are assumed to be
body-fixed; if we are defining the orientation of B in A using by rot_type,
the angular velocity of B in A is assumed to be in the form: speed[0]*B.x +
speed[1]*B.y + speed[2]*B.z
Parameters
==========
speeds : list of length 3
The body fixed angular velocity measure numbers.
coords : list of length 3 or 4
The coordinates used to define the orientation of the two frames.
rot_type : str
The type of rotation used to create the equations. Body, Space, or
Quaternion only
rot_order : str
If applicable, the order of a series of rotations.
Examples
========
>>> from sympy.physics.mechanics import dynamicsymbols
>>> from sympy.physics.mechanics import kinematic_equations, mprint
>>> u1, u2, u3 = dynamicsymbols('u1 u2 u3')
>>> q1, q2, q3 = dynamicsymbols('q1 q2 q3')
>>> mprint(kinematic_equations([u1,u2,u3], [q1,q2,q3], 'body', '313'),
... order=None)
[-(u1*sin(q3) + u2*cos(q3))/sin(q2) + q1', -u1*cos(q3) + u2*sin(q3) + q2', (u1*sin(q3) + u2*cos(q3))*cos(q2)/sin(q2) - u3 + q3']
"""
# Code below is checking and sanitizing input
approved_orders = ('123', '231', '312', '132', '213', '321', '121', '131',
'212', '232', '313', '323', '1', '2', '3', '')
rot_order = str(rot_order).upper() # Now we need to make sure XYZ = 123
rot_type = rot_type.upper()
rot_order = [i.replace('X', '1') for i in rot_order]
rot_order = [i.replace('Y', '2') for i in rot_order]
rot_order = [i.replace('Z', '3') for i in rot_order]
rot_order = ''.join(rot_order)
if not isinstance(speeds, (list, tuple)):
raise TypeError('Need to supply speeds in a list')
if len(speeds) != 3:
raise TypeError('Need to supply 3 body-fixed speeds')
if not isinstance(coords, (list, tuple)):
raise TypeError('Need to supply coordinates in a list')
if rot_type.lower() in ['body', 'space']:
if rot_order not in approved_orders:
raise ValueError('Not an acceptable rotation order')
if len(coords) != 3:
raise ValueError('Need 3 coordinates for body or space')
# Actual hard-coded kinematic differential equations
q1, q2, q3 = coords
q1d, q2d, q3d = [diff(i, dynamicsymbols._t) for i in coords]
w1, w2, w3 = speeds
s1, s2, s3 = [sin(q1), sin(q2), sin(q3)]
c1, c2, c3 = [cos(q1), cos(q2), cos(q3)]
if rot_type.lower() == 'body':
if rot_order == '123':
return [q1d - (w1 * c3 - w2 * s3) / c2, q2d - w1 * s3 - w2 *
c3, q3d - (-w1 * c3 + w2 * s3) * s2 / c2 - w3]
if rot_order == '231':
return [q1d - (w2 * c3 - w3 * s3) / c2, q2d - w2 * s3 - w3 *
c3, q3d - w1 - (- w2 * c3 + w3 * s3) * s2 / c2]
if rot_order == '312':
return [q1d - (-w1 * s3 + w3 * c3) / c2, q2d - w1 * c3 - w3 *
s3, q3d - (w1 * s3 - w3 * c3) * s2 / c2 - w2]
if rot_order == '132':
return [q1d - (w1 * c3 + w3 * s3) / c2, q2d + w1 * s3 - w3 *
c3, q3d - (w1 * c3 + w3 * s3) * s2 / c2 - w2]
if rot_order == '213':
return [q1d - (w1 * s3 + w2 * c3) / c2, q2d - w1 * c3 + w2 *
s3, q3d - (w1 * s3 + w2 * c3) * s2 / c2 - w3]
if rot_order == '321':
return [q1d - (w2 * s3 + w3 * c3) / c2, q2d - w2 * c3 + w3 *
s3, q3d - w1 - (w2 * s3 + w3 * c3) * s2 / c2]
if rot_order == '121':
return [q1d - (w2 * s3 + w3 * c3) / s2, q2d - w2 * c3 + w3 *
s3, q3d - w1 + (w2 * s3 + w3 * c3) * c2 / s2]
if rot_order == '131':
return [q1d - (-w2 * c3 + w3 * s3) / s2, q2d - w2 * s3 - w3 *
c3, q3d - w1 - (w2 * c3 - w3 * s3) * c2 / s2]
if rot_order == '212':
return [q1d - (w1 * s3 - w3 * c3) / s2, q2d - w1 * c3 - w3 *
s3, q3d - (-w1 * s3 + w3 * c3) * c2 / s2 - w2]
if rot_order == '232':
return [q1d - (w1 * c3 + w3 * s3) / s2, q2d + w1 * s3 - w3 *
c3, q3d + (w1 * c3 + w3 * s3) * c2 / s2 - w2]
if rot_order == '313':
return [q1d - (w1 * s3 + w2 * c3) / s2, q2d - w1 * c3 + w2 *
s3, q3d + (w1 * s3 + w2 * c3) * c2 / s2 - w3]
if rot_order == '323':
return [q1d - (-w1 * c3 + w2 * s3) / s2, q2d - w1 * s3 - w2 *
c3, q3d - (w1 * c3 - w2 * s3) * c2 / s2 - w3]
if rot_type.lower() == 'space':
if rot_order == '123':
return [q1d - w1 - (w2 * s1 + w3 * c1) * s2 / c2, q2d - w2 *
c1 + w3 * s1, q3d - (w2 * s1 + w3 * c1) / c2]
if rot_order == '231':
return [q1d - (w1 * c1 + w3 * s1) * s2 / c2 - w2, q2d + w1 *
s1 - w3 * c1, q3d - (w1 * c1 + w3 * s1) / c2]
if rot_order == '312':
return [q1d - (w1 * s1 + w2 * c1) * s2 / c2 - w3, q2d - w1 *
c1 + w2 * s1, q3d - (w1 * s1 + w2 * c1) / c2]
if rot_order == '132':
return [q1d - w1 - (-w2 * c1 + w3 * s1) * s2 / c2, q2d - w2 *
s1 - w3 * c1, q3d - (w2 * c1 - w3 * s1) / c2]
if rot_order == '213':
return [q1d - (w1 * s1 - w3 * c1) * s2 / c2 - w2, q2d - w1 *
c1 - w3 * s1, q3d - (-w1 * s1 + w3 * c1) / c2]
if rot_order == '321':
return [q1d - (-w1 * c1 + w2 * s1) * s2 / c2 - w3, q2d - w1 *
s1 - w2 * c1, q3d - (w1 * c1 - w2 * s1) / c2]
if rot_order == '121':
return [q1d - w1 + (w2 * s1 + w3 * c1) * c2 / s2, q2d - w2 *
c1 + w3 * s1, q3d - (w2 * s1 + w3 * c1) / s2]
if rot_order == '131':
return [q1d - w1 - (w2 * c1 - w3 * s1) * c2 / s2, q2d - w2 *
s1 - w3 * c1, q3d - (-w2 * c1 + w3 * s1) / s2]
if rot_order == '212':
return [q1d - (-w1 * s1 + w3 * c1) * c2 / s2 - w2, q2d - w1 *
c1 - w3 * s1, q3d - (w1 * s1 - w3 * c1) / s2]
if rot_order == '232':
return [q1d + (w1 * c1 + w3 * s1) * c2 / s2 - w2, q2d + w1 *
s1 - w3 * c1, q3d - (w1 * c1 + w3 * s1) / s2]
if rot_order == '313':
return [q1d + (w1 * s1 + w2 * c1) * c2 / s2 - w3, q2d - w1 *
c1 + w2 * s1, q3d - (w1 * s1 + w2 * c1) / s2]
if rot_order == '323':
return [q1d - (w1 * c1 - w2 * s1) * c2 / s2 - w3, q2d - w1 *
s1 - w2 * c1, q3d - (-w1 * c1 + w2 * s1) / s2]
elif rot_type.lower() == 'quaternion':
if rot_order != '':
raise ValueError('Cannot have rotation order for quaternion')
if len(coords) != 4:
raise ValueError('Need 4 coordinates for quaternion')
# Actual hard-coded kinematic differential equations
e0, e1, e2, e3 = coords
w = Matrix(speeds + [0])
E = Matrix([[e0, -e3, e2, e1], [e3, e0, -e1, e2], [-e2, e1, e0, e3],
[-e1, -e2, -e3, e0]])
edots = Matrix([diff(i, dynamicsymbols._t) for i in [e1, e2, e3, e0]])
return list(edots.T - 0.5 * w.T * E.T)
else:
raise ValueError('Not an approved rotation type for this function')
[docs]def partial_velocity(vel_list, u_list, frame):
"""Returns a list of partial velocities.
For a list of velocity or angular velocity vectors the partial derivatives
with respect to the supplied generalized speeds are computed, in the
specified ReferenceFrame.
The output is a list of lists. The outer list has a number of elements
equal to the number of supplied velocity vectors. The inner lists are, for
each velocity vector, the partial derivatives of that velocity vector with
respect to the generalized speeds supplied.
Parameters
==========
vel_list : list
List of velocities of Point's and angular velocities of ReferenceFrame's
u_list : list
List of independent generalized speeds.
frame : ReferenceFrame
The ReferenceFrame the partial derivatives are going to be taken in.
Examples
========
>>> from sympy.physics.mechanics import Point, ReferenceFrame
>>> from sympy.physics.mechanics import dynamicsymbols
>>> from sympy.physics.mechanics import partial_velocity
>>> u = dynamicsymbols('u')
>>> N = ReferenceFrame('N')
>>> P = Point('P')
>>> P.set_vel(N, u * N.x)
>>> vel_list = [P.vel(N)]
>>> u_list = [u]
>>> partial_velocity(vel_list, u_list, N)
[[N.x]]
"""
if not hasattr(vel_list, '__iter__'):
raise TypeError('Provide velocities in an iterable')
if not hasattr(u_list, '__iter__'):
raise TypeError('Provide speeds in an iterable')
list_of_pvlists = []
for i in vel_list:
pvlist = []
for j in u_list:
vel = i.diff(j, frame)
pvlist += [vel]
list_of_pvlists += [pvlist]
return list_of_pvlists
[docs]def linear_momentum(frame, *body):
"""Linear momentum of the system.
This function returns the linear momentum of a system of Particle's and/or
RigidBody's. The linear momentum of a system is equal to the vector sum of
the linear momentum of its constituents. Consider a system, S, comprised of
a rigid body, A, and a particle, P. The linear momentum of the system, L,
is equal to the vector sum of the linear momentum of the particle, L1, and
the linear momentum of the rigid body, L2, i.e-
L = L1 + L2
Parameters
==========
frame : ReferenceFrame
The frame in which linear momentum is desired.
body1, body2, body3... : Particle and/or RigidBody
The body (or bodies) whose kinetic energy is required.
Examples
========
>>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
>>> from sympy.physics.mechanics import RigidBody, outer, linear_momentum
>>> N = ReferenceFrame('N')
>>> P = Point('P')
>>> P.set_vel(N, 10 * N.x)
>>> Pa = Particle('Pa', P, 1)
>>> Ac = Point('Ac')
>>> Ac.set_vel(N, 25 * N.y)
>>> I = outer(N.x, N.x)
>>> A = RigidBody('A', Ac, N, 20, (I, Ac))
>>> linear_momentum(N, A, Pa)
10*N.x + 500*N.y
"""
if not isinstance(frame, ReferenceFrame):
raise TypeError('Please specify a valid ReferenceFrame')
else:
linear_momentum_sys = S(0)
for e in body:
if isinstance(e, (RigidBody, Particle)):
linear_momentum_sys += e.linear_momentum(frame)
else:
raise TypeError('*body must have only Particle or RigidBody')
return linear_momentum_sys
[docs]def angular_momentum(point, frame, *body):
"""Angular momentum of a system
This function returns the angular momentum of a system of Particle's and/or
RigidBody's. The angular momentum of such a system is equal to the vector
sum of the angular momentum of its constituents. Consider a system, S,
comprised of a rigid body, A, and a particle, P. The angular momentum of
the system, H, is equal to the vector sum of the linear momentum of the
particle, H1, and the linear momentum of the rigid body, H2, i.e-
H = H1 + H2
Parameters
==========
point : Point
The point about which angular momentum of the system is desired.
frame : ReferenceFrame
The frame in which angular momentum is desired.
body1, body2, body3... : Particle and/or RigidBody
The body (or bodies) whose kinetic energy is required.
Examples
========
>>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
>>> from sympy.physics.mechanics import RigidBody, outer, angular_momentum
>>> N = ReferenceFrame('N')
>>> O = Point('O')
>>> O.set_vel(N, 0 * N.x)
>>> P = O.locatenew('P', 1 * N.x)
>>> P.set_vel(N, 10 * N.x)
>>> Pa = Particle('Pa', P, 1)
>>> Ac = O.locatenew('Ac', 2 * N.y)
>>> Ac.set_vel(N, 5 * N.y)
>>> a = ReferenceFrame('a')
>>> a.set_ang_vel(N, 10 * N.z)
>>> I = outer(N.z, N.z)
>>> A = RigidBody('A', Ac, a, 20, (I, Ac))
>>> angular_momentum(O, N, Pa, A)
10*N.z
"""
if not isinstance(frame, ReferenceFrame):
raise TypeError('Please enter a valid ReferenceFrame')
if not isinstance(point, Point):
raise TypeError('Please specify a valid Point')
else:
angular_momentum_sys = S(0)
for e in body:
if isinstance(e, (RigidBody, Particle)):
angular_momentum_sys += e.angular_momentum(point, frame)
else:
raise TypeError('*body must have only Particle or RigidBody')
return angular_momentum_sys
[docs]def kinetic_energy(frame, *body):
"""Kinetic energy of a multibody system.
This function returns the kinetic energy of a system of Particle's and/or
RigidBody's. The kinetic energy of such a system is equal to the sum of
the kinetic energies of its constituents. Consider a system, S, comprising
a rigid body, A, and a particle, P. The kinetic energy of the system, T,
is equal to the vector sum of the kinetic energy of the particle, T1, and
the kinetic energy of the rigid body, T2, i.e.
T = T1 + T2
Kinetic energy is a scalar.
Parameters
==========
frame : ReferenceFrame
The frame in which the velocity or angular velocity of the body is
defined.
body1, body2, body3... : Particle and/or RigidBody
The body (or bodies) whose kinetic energy is required.
Examples
========
>>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
>>> from sympy.physics.mechanics import RigidBody, outer, kinetic_energy
>>> N = ReferenceFrame('N')
>>> O = Point('O')
>>> O.set_vel(N, 0 * N.x)
>>> P = O.locatenew('P', 1 * N.x)
>>> P.set_vel(N, 10 * N.x)
>>> Pa = Particle('Pa', P, 1)
>>> Ac = O.locatenew('Ac', 2 * N.y)
>>> Ac.set_vel(N, 5 * N.y)
>>> a = ReferenceFrame('a')
>>> a.set_ang_vel(N, 10 * N.z)
>>> I = outer(N.z, N.z)
>>> A = RigidBody('A', Ac, a, 20, (I, Ac))
>>> kinetic_energy(N, Pa, A)
350
"""
if not isinstance(frame, ReferenceFrame):
raise TypeError('Please enter a valid ReferenceFrame')
ke_sys = S(0)
for e in body:
if isinstance(e, (RigidBody, Particle)):
ke_sys += e.kinetic_energy(frame)
else:
raise TypeError('*body must have only Particle or RigidBody')
return ke_sys
[docs]def potential_energy(*body):
"""Potential energy of a multibody system.
This function returns the potential energy of a system of Particle's and/or
RigidBody's. The potential energy of such a system is equal to the sum of
the potential energy of its constituents. Consider a system, S, comprising
a rigid body, A, and a particle, P. The potential energy of the system, V,
is equal to the vector sum of the potential energy of the particle, V1, and
the potential energy of the rigid body, V2, i.e.
V = V1 + V2
Potential energy is a scalar.
Parameters
==========
body1, body2, body3... : Particle and/or RigidBody
The body (or bodies) whose potential energy is required.
Examples
========
>>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
>>> from sympy.physics.mechanics import RigidBody, outer, potential_energy
>>> from sympy import symbols
>>> M, m, g, h = symbols('M m g h')
>>> N = ReferenceFrame('N')
>>> O = Point('O')
>>> O.set_vel(N, 0 * N.x)
>>> P = O.locatenew('P', 1 * N.x)
>>> Pa = Particle('Pa', P, m)
>>> Ac = O.locatenew('Ac', 2 * N.y)
>>> a = ReferenceFrame('a')
>>> I = outer(N.z, N.z)
>>> A = RigidBody('A', Ac, a, M, (I, Ac))
>>> Pa.set_potential_energy(m * g * h)
>>> A.set_potential_energy(M * g * h)
>>> potential_energy(Pa, A)
M*g*h + g*h*m
"""
pe_sys = S(0)
for e in body:
if isinstance(e, (RigidBody, Particle)):
pe_sys += e.potential_energy
else:
raise TypeError('*body must have only Particle or RigidBody')
return pe_sys
[docs]def Lagrangian(frame, *body):
"""Lagrangian of a multibody system.
This function returns the Lagrangian of a system of Particle's and/or
RigidBody's. The Lagrangian of such a system is equal to the difference
between the kinetic energies and potential energies of its constituents. If
T and V are the kinetic and potential energies of a system then it's
Lagrangian, L, is defined as
L = T - V
The Lagrangian is a scalar.
Parameters
==========
frame : ReferenceFrame
The frame in which the velocity or angular velocity of the body is
defined to determine the kinetic energy.
body1, body2, body3... : Particle and/or RigidBody
The body (or bodies) whose kinetic energy is required.
Examples
========
>>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
>>> from sympy.physics.mechanics import RigidBody, outer, Lagrangian
>>> from sympy import symbols
>>> M, m, g, h = symbols('M m g h')
>>> N = ReferenceFrame('N')
>>> O = Point('O')
>>> O.set_vel(N, 0 * N.x)
>>> P = O.locatenew('P', 1 * N.x)
>>> P.set_vel(N, 10 * N.x)
>>> Pa = Particle('Pa', P, 1)
>>> Ac = O.locatenew('Ac', 2 * N.y)
>>> Ac.set_vel(N, 5 * N.y)
>>> a = ReferenceFrame('a')
>>> a.set_ang_vel(N, 10 * N.z)
>>> I = outer(N.z, N.z)
>>> A = RigidBody('A', Ac, a, 20, (I, Ac))
>>> Pa.set_potential_energy(m * g * h)
>>> A.set_potential_energy(M * g * h)
>>> Lagrangian(N, Pa, A)
-M*g*h - g*h*m + 350
"""
if not isinstance(frame, ReferenceFrame):
raise TypeError('Please supply a valid ReferenceFrame')
for e in body:
if not isinstance(e, (RigidBody, Particle)):
raise TypeError('*body must have only Particle or RigidBody')
return kinetic_energy(frame, *body) - potential_energy(*body)