New matrices support in sympy.physics.matrices

sympy/matrices/__init__.py

@@ -4,4 +4,5 @@
 """
 from matrices import Matrix, SMatrix, zero, zeronm, zeros, one, ones, eye, \
      hessian, randMatrix, GramSchmidt, wronskian, casoratian, \
-     list2numpy, matrix2numpy, DeferredVector, block_diag, symarray
+     list2numpy, matrix2numpy, DeferredVector, block_diag, symarray, \
+     rot_axis1, rot_axis2, rot_axis3

sympy/matrices/matrices.py

@@ -9,6 +9,7 @@
 from sympy.simplify import simplify
 from sympy.utilities import any, all
 from sympy.printing import sstr
+from sympy.functions.elementary.trigonometric import cos, sin
 
 
 import random
@@ -2203,3 +2204,94 @@ def symarray(prefix, shape):
         arr[index] = Symbol('%s_%s' % (prefix, '_'.join(map(str, index))))
     return arr
 
+def rot_axis3(theta):
+    """Returns a rotation matrix for a rotation of theta (in radians) about
+    the 3-axis.
+
+    Examples
+    --------
+
+    >>> from sympy import pi
+    >>> from sympy.matrices import rot_axis3
+
+    A rotation of pi/3 (60 degrees):
+    >>> theta = pi/3
+    >>> rot_axis3(theta)
+    [        1/2, 3**(1/2)/2, 0]
+    [-3**(1/2)/2,        1/2, 0]
+    [          0,          0, 1]
+
+    If we rotate by pi/2 (90 degrees):
+    >> rot_axis3(pi/2)
+    [ 0, 1, 0]
+    [-1, 0, 0]
+    [ 0, 0, 1]
+    """
+    ct = cos(theta)
+    st = sin(theta)
+    mat = ((ct,st,0),
+           (-st,ct,0),
+           (0,0,1))
+    return Matrix(mat)
+
+def rot_axis2(theta):
+    """Returns a rotation matrix for a rotation of theta (in radians) about
+    the 2-axis.
+
+    Examples
+    --------
+
+    >>> from sympy import pi
+    >>> from sympy.matrices import rot_axis2
+
+    A rotation of pi/3 (60 degrees):
+    >>> theta = pi/3
+    >>> rot_axis2(theta)
+    [       1/2, 0, -3**(1/2)/2]
+    [         0, 1,           0]
+    [3**(1/2)/2, 0,         1/2]
+
+    If we rotate by pi/2 (90 degrees):
+    >>> rot_axis2(pi/2)
+    [0, 0, -1]
+    [0, 1,  0]
+    [1, 0,  0]
+    """
+    ct = cos(theta)
+    st = sin(theta)
+    mat = ((ct,0,-st),
+           (0,1,0),
+           (st,0,ct))
+    return Matrix(mat)
+
+def rot_axis1(theta):
+    """Returns a rotation matrix for a rotation of theta (in radians) about
+    the 1-axis.
+
+    Examples
+    --------
+
+    >>> from sympy import pi
+    >>> from sympy.matrices import rot_axis1
+
+    A rotation of pi/3 (60 degrees):
+    >>> theta = pi/3
+    >>> rot_axis1(theta)
+    [1,           0,          0]
+    [0,         1/2, 3**(1/2)/2]
+    [0, -3**(1/2)/2,        1/2]
+
+    If we rotate by pi/2 (90 degrees):
+    >>> rot_axis1(pi/2)
+    [1,  0, 0]
+    [0,  0, 1]
+    [0, -1, 0]
+    """
+    ct = cos(theta)
+    st = sin(theta)
+    mat = ((1,0,0),
+           (0,ct,st),
+           (0,-st,ct))
+    return Matrix(mat)
+
+

sympy/matrices/tests/test_matrices.py

@@ -4,6 +4,7 @@
 from sympy.matrices.matrices import (ShapeError, MatrixError,
     matrix_multiply_elementwise)
 from sympy.utilities.pytest import raises
+from sympy.matrices import rot_axis1, rot_axis2, rot_axis3
 
 def test_division():
     x, y, z = symbols('x','y','z')
@@ -1202,3 +1203,26 @@ def test_creation_args():
     assert eye(3.) == eye(3)
     assert ones((3L, Integer(4))) == ones((3, 4))
     raises(TypeError, 'Matrix(1, 2)')
+
+def test_rotation_matrices():
+    #This tests the rotation matrices by rotating about an axis and back.
+    theta = pi/3
+    r3_plus = rot_axis3(theta)
+    r3_minus = rot_axis3(-theta)
+    r2_plus = rot_axis2(theta)
+    r2_minus = rot_axis2(-theta)
+    r1_plus = rot_axis1(theta)
+    r1_minus = rot_axis1(-theta)
+    assert r3_minus*r3_plus*eye(3)== eye(3)
+    assert r2_minus*r2_plus*eye(3)== eye(3)
+    assert r1_minus*r1_plus*eye(3)== eye(3)
+
+    # Check that the rotation matrix trace is correct
+    assert r1_plus.trace() == 1 + 2*cos(theta)
+    assert r2_plus.trace() == 1 + 2*cos(theta)
+    assert r3_plus.trace() == 1 + 2*cos(theta)
+
+    # Chcek that a rotation of zero doesn't change things.
+    assert rot_axis1(0) == eye(3)
+    assert rot_axis2(0) == eye(3)
+    assert rot_axis3(0) == eye(3)

sympy/physics/matrices.py

@@ -26,6 +26,38 @@ def msigma(i):
         raise IndexError("Invalid Pauli index")
     return Matrix(mat)
 
+def pat_matrix(m, dx, dy, dz):
+    """Returns the Parallel Axis Theorem matrix to translate the inertia
+    matrix a distance of (dx, dy, dz) for a body of mass m.
+
+    Examples
+    --------
+    If the point we want the inertia about is a distance of 2 units of
+    length and 1 unit along the x-axis we get.
+    >>> from sympy.physics.matrices import pat_matrix
+    >>> pat_matrix(2,1,0,0)
+    [0, 0, 0]
+    [0, 2, 0]
+    [0, 0, 2]
+
+    and if we want to move the inertia along a vector of (1,1,1) for the
+    same mass,
+    >>> pat_matrix(2,1,1,1)
+    [ 4, -2, -2]
+    [-2,  4, -2]
+    [-2, -2,  4]
+    """
+    dxdy = -dx*dy
+    dydz = -dy*dz
+    dxdz = -dx*dz
+    dx2 = dx**2
+    dy2 = dy**2
+    dz2 = dz**2
+    mat = ((dy2 + dz2,dxdy,dxdz),
+           (dxdy,dx2+dz2,dydz),
+           (dxdz, dydz, dy2+dx2))
+    return m*Matrix(mat)
+
 def mgamma(mu,lower=False):
     """Returns a Dirac gamma matrix gamma^mu in the standard
     (Dirac) representation.

sympy/physics/tests/test_physics_matrices.py

@@ -1,7 +1,29 @@
-from sympy.physics.matrices import msigma, mgamma, minkowski_tensor
-from sympy import zeros, eye, I
+from sympy.physics.matrices import msigma, mgamma, minkowski_tensor, pat_matrix
+from sympy import zeros, eye, I, Matrix
 
+def test_parallel_axis_theorem():
+    # This tests the parallel axis theorem matrix by comparing to test
+    # matrices.
 
+    # First case, 1 in all directions.
+    mat1 = Matrix(((2,-1,-1),(-1,2,-1),(-1,-1,2)))
+    assert pat_matrix(1,1,1,1) == mat1
+    assert pat_matrix(2,1,1,1) == 2*mat1
+
+    # Second case, 1 in x, 0 in all others
+    mat2 = Matrix(((0,0,0),(0,1,0),(0,0,1)))
+    assert pat_matrix(1,1,0,0) == mat2
+    assert pat_matrix(2,1,0,0) == 2*mat2
+
+    # Third case, 1 in y, 0 in all others
+    mat3 = Matrix(((1,0,0),(0,0,0),(0,0,1)))
+    assert pat_matrix(1,0,1,0) == mat3
+    assert pat_matrix(2,0,1,0) == 2*mat3
+
+    # Fourth case, 1 in z, 0 in all others
+    mat4 = Matrix(((1,0,0),(0,1,0),(0,0,0)))
+    assert pat_matrix(1,0,0,1) == mat4
+    assert pat_matrix(2,0,0,1) == 2*mat4
 
 def test_Pauli():
     #this and the following test are testing both Pauli and Dirac matrices