# sympy/sympy

Merge pull request #1346 from angadhn/dev2

`Linear momentum function and tests.`
certik committed Jul 3, 2012
2 parents bdeee5b + f9d974c commit 48f42f8aabf988284f624453ef7c3fc8861325d3
 @@ -9,13 +9,17 @@ 'mlatex', 'kinematic_equations', 'inertia_of_point_mass', - 'partial_velocity'] + 'partial_velocity', + 'linearmomentum', + 'angularmomentum'] 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 import sympify, diff, sin, cos, Matrix def cross(vec1, vec2): @@ -443,3 +447,113 @@ def partial_velocity(vel_list, u_ind_list): list_of_pvlists = list_of_pvlists + [pvlist] i = i + 1 return list_of_pvlists + +def linearmomentum(bodylist, frame): + """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 + + bodylist : list + A list of all RigidBody's and Particle's in the system. + + frame : ReferenceFrame + The frame in which linear momentum is desired. + + Examples + ======== + + >>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame + >>> from sympy.physics.mechanics import RigidBody, outer, linearmomentum + >>> 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)) + >>> BL = [Pa, A] + >>> print linearmomentum(BL, N) + 10*N.x + 500*N.y + + """ + + if not isinstance(bodylist, list): + raise TypeError('System elements must be supplied as a list') + else: + linearmomentum_sys = 0 + for e in bodylist: + if isinstance(e, (RigidBody, Particle)): + lm_individual = e.linearmomentum(frame) + linearmomentum_sys = linearmomentum_sys + lm_individual + else: + raise TypeError('Body list must have only Particle' + 'or RigidBody') + return linearmomentum_sys + +def angularmomentum(bodylist, point, frame): + """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 + ========== + + bodylist : list + A list of all RigidBody's and Particle's in the system. + + point : Point + The point about which angular momentum of the system is desired. + + frame : ReferenceFrame + The frame in which angular momentum is desired. + + Examples + ======== + + >>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame + >>> from sympy.physics.mechanics import RigidBody, outer, angularmomentum + >>> 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)) + >>> BL = [Pa, A] + >>> print angularmomentum(BL, O, N) + 10*N.z + + """ + + if not isinstance(bodylist, list): + raise TypeError('System elements must be supplied as a list') + else: + angularmomentum_sys = 0 + for e in bodylist: + if isinstance(e, (RigidBody, Particle)): + angmom_individual = e.angularmomentum(point, frame) + angularmomentum_sys = angularmomentum_sys + angmom_individual + else: + raise TypeError('Body list must have only Particle' + 'or RigidBody') + return angularmomentum_sys
 @@ -66,3 +66,76 @@ def set_point(self, p): self._point = p point = property(get_point, set_point) + + def linearmomentum(self, frame): + """Linear momentum of the particle. + + The linear momentum L, of a particle P, with respect to frame N is + given by + + L = m * v + + where m is the mass of the particle, and v is the velocity of the + particle in the frame N. + + Parameters + ========== + + frame : ReferenceFrame + The frame in which linear momentum is desired. + + Examples + ======== + + >>> from sympy.physics.mechanics import Particle, Point, ReferenceFrame + >>> from sympy.physics.mechanics import dynamicsymbols + >>> m, v = dynamicsymbols('m v') + >>> N = ReferenceFrame('N') + >>> P = Point('P') + >>> A = Particle('A', P, m) + >>> P.set_vel(N, v * N.x) + >>> print A.linearmomentum(N) + m*v*N.x + + """ + + return self.mass * self.point.vel(frame) + + def angularmomentum(self, point, frame): + """Angular momentum of the particle about the point. + + The angular momentum H, about some point O of a particle, P, is given + by: + + H = r x m * v + + where r is the position vector from point O to the particle P, m is + the mass of the particle, and v is the velocity of the particle in + the inertial frame, N. + + Parameters + ========== + + point : Point + The point about which angular momentum of the particle is desired. + + frame : ReferenceFrame + The frame in which angular momentum is desired. + + Examples + ======== + + >>> from sympy.physics.mechanics import Particle, Point, ReferenceFrame + >>> from sympy.physics.mechanics import dynamicsymbols + >>> m, v, r = dynamicsymbols('m v r') + >>> N = ReferenceFrame('N') + >>> O = Point('O') + >>> A = O.locatenew('A', r * N.x) + >>> P = Particle('P', A, m) + >>> P.point.set_vel(N, v * N.y) + >>> print P.angularmomentum(O, N) + m*r*v*N.z + + """ + + return self.point.pos_from(point) ^ (self.mass * self.point.vel(frame))
 @@ -99,3 +99,84 @@ def set_inertia(self, I): self._inertia_point = I[1] inertia = property(get_inertia, set_inertia) + + def linearmomentum(self, frame): + """ Linear momentum of the rigid body. + + The linear momentum L, of a rigid body B, with respect to frame N is + given by + + L = M * v* + + where M is the mass of the rigid body and v* is the velocity of + the mass center of B in the frame, N. + + Parameters + ========== + + frame : ReferenceFrame + The frame in which linear momentum is desired. + + Examples + ======== + + >>> from sympy.physics.mechanics import Point, ReferenceFrame, outer + >>> from sympy.physics.mechanics import RigidBody, dynamicsymbols + >>> M, v = dynamicsymbols('M v') + >>> N = ReferenceFrame('N') + >>> P = Point('P') + >>> P.set_vel(N, v * N.x) + >>> I = outer (N.x, N.x) + >>> Inertia_tuple = (I, P) + >>> B = RigidBody('B', P, N, M, Inertia_tuple) + >>> B.linearmomentum(N) + M*v*N.x + + """ + + return self.mass * self.masscenter.vel(frame) + + def angularmomentum(self, point, frame): + """ Angular momentum of the rigid body. + + The angular momentum H, about some point O, of a rigid body B, in a + frame N is given by + + H = I* . omega + r* x (M * v) + + where I* is the central inertia dyadic of B, omega is the angular + velocity of body B in the frame, N, r* is the position vector from + point O to the mass center of B, and v is the velocity of point O in + the frame, N. + + Parameters + ========== + + point : Point + The point about which angular momentum is desired. + + frame : ReferenceFrame + The frame in which angular momentum is desired. + + Examples + ======== + + >>> from sympy.physics.mechanics import Point, ReferenceFrame, outer + >>> from sympy.physics.mechanics import RigidBody, dynamicsymbols + >>> M, v, r, omega = dynamicsymbols('M v r omega') + >>> N = ReferenceFrame('N') + >>> b = ReferenceFrame('b') + >>> b.set_ang_vel(N, omega * b.x) + >>> P = Point('P') + >>> P.set_vel(N, 1 * N.x) + >>> I = outer (b.x, b.x) + >>> Inertia_tuple = (I, P) + >>> B = RigidBody('B', P, b, M, Inertia_tuple) + >>> B.angularmomentum(P, N) + omega*b.x + + """ + + return ((self.inertia[0] & self.frame.ang_vel_in(frame)) + + (point.vel(frame) ^ -self.masscenter.pos_from(point)) * + self.mass)
 @@ -2,7 +2,10 @@ from sympy.physics.mechanics import (cross, dot, dynamicsymbols, express, ReferenceFrame, inertia, Point, kinematic_equations, Vector, - inertia_of_point_mass, partial_velocity) + inertia_of_point_mass, partial_velocity, + outer, Particle, + Point, RigidBody, angularmomentum, + linearmomentum) Vector.simp = True q1, q2, q3, q4, q5 = symbols('q1 q2 q3 q4 q5') @@ -300,3 +303,34 @@ def test_partial_velocity(): assert (partial_velocity(vel_list, u_list) == [[- r*L.y, 0, L.x], [r*L.x, 0, L.y], [0, 0, L.z], [L.x, L.x, 0], [cos(q2)*L.y - sin(q2)*L.z, cos(q2)*L.y - sin(q2)*L.z, 0]]) + +def test_linearmomentum(): + N = ReferenceFrame('N') + Ac = Point('Ac') + Ac.set_vel(N, 25 * N.y) + I = outer(N.x, N.x) + A = RigidBody('A', Ac, N, 20, (I, Ac)) + P = Point('P') + Pa = Particle('Pa', P, 1) + Pa.point.set_vel(N, 10 * N.x) + BL = [A, Pa] + assert linearmomentum(BL, N) == 10 * N.x + 500 * N.y + +def test_angularmomentum_and_linearmomentum(): + m, M, l1 = symbols('m M l1') + q1d = dynamicsymbols('q1d') + N = ReferenceFrame('N') + O = Point('O') + O.set_vel(N, 0 * N.x) + Ac = O.locatenew('Ac', l1 * N.x) + P = Ac.locatenew('P', l1 * N.x) + a = ReferenceFrame('a') + a.set_ang_vel(N, q1d * N.z) + Ac.v2pt_theory(O, N, a) + P.v2pt_theory(O, N, a) + Pa = Particle('Pa', P, m) + I = outer(N.z, N.z) + A = RigidBody('A', Ac, a, M, (I, Ac)) + BL =[Pa, A] + assert linearmomentum(BL, N) == 2 * m * q1d* l1 * N.y + M * l1 * q1d * N.y + assert angularmomentum(BL, O, N) == 4 * m * q1d * l1**2 * N.z + q1d * N.z
 @@ -1,8 +1,8 @@ from sympy import symbols -from sympy.physics.mechanics import Point, Particle +from sympy.physics.mechanics import Point, Particle, ReferenceFrame def test_particle(): - m, m2 = symbols('m m2') + m, m2, v1, v2, v3, r = symbols('m m2 v1 v2 v3 r') P = Point('P') P2 = Point('P2') p = Particle('pa', P, m) @@ -14,3 +14,19 @@ def test_particle(): # Test the point setter p.point = P2 assert p.point == P2 + # Test the linear momentum function + N = ReferenceFrame('N') + O = Point('O') + P2.set_pos(O, r * N.y) + P2.set_vel(N, v1 * N.x) + assert p.linearmomentum(N) == m2 * v1 * N.x + assert p.angularmomentum(O, N) == -m2 * r *v1 * N.z + P2.set_vel(N, v2 * N.y) + assert p.linearmomentum(N) == m2 * v2 * N.y + assert p.angularmomentum(O, N) == 0 + P2.set_vel(N, v3 * N.z) + assert p.linearmomentum(N) == m2 * v3 * N.z + assert p.angularmomentum(O, N) == m2 * r * v3 * N.x + P2.set_vel(N, v1 * N.x + v2 * N.y + v3 * N.z) + assert p.linearmomentum(N) == m2 * (v1 * N.x + v2 * N.y + v3 * N.z) + assert p.angularmomentum(O, N) == m2 * r * (v3 * N.x - v1 * N.z)