r5059 - /1.3/maths_fns/rotation_matrix.py -- February 20, 2008 - 18:53

Author: bugman
Date: Wed Feb 20 18:53:03 2008
New Revision: 5059

URL: http://svn.gna.org/viewcvs/relax?rev=5059&view=rev
Log:
Wrote a new function R_2vect() to generate a rotation matrix from 2 vectors.

This rotation matrix is not unique because of axially symmetry allowing the 
rotation axis to be any
vector within the symmetry plane, so the cross product with the axis-angle 
convention is used to
define this matrix.


Modified:
    1.3/maths_fns/rotation_matrix.py

Modified: 1.3/maths_fns/rotation_matrix.py
URL: 
http://svn.gna.org/viewcvs/relax/1.3/maths_fns/rotation_matrix.py?rev=5059&r1=5058&r2=5059&view=diff
==============================================================================
--- 1.3/maths_fns/rotation_matrix.py (original)
+++ 1.3/maths_fns/rotation_matrix.py Wed Feb 20 18:53:03 2008
@@ -21,8 +21,61 @@
 
###############################################################################
 
 # Python module imports.
-from numpy import dot
-from math import cos, sin
+from numpy import cross, dot
+from math import acos, cos, sin
+
+
+def R_2vect(R, vector_orig, vector_fin):
+    """Calculate the rotation matrix required to rotate from one vector to 
another.
+
+    For the rotation of one vector to another, there are an infinit series 
of rotation matrices
+    possible.  Due to axially symmetry, the rotation axis can be any vector 
lying in the symmetry
+    plane between the two vectors.  Hence the axis-angle convention will be 
used to construct the
+    matrix with the rotation axis defined as the cross product of the two 
vectors.  The rotation
+    angle is the arccosine of the dot product of the two unit vectors.
+
+    Given a unit vector parallel to the rotation axis, w = [x, y, z] and the 
rotation angle a,
+    the rotation matrix R is:
+
+              |  1 + (1-cos(a))*(x*x-1)   -z*sin(a)+(1-cos(a))*x*y   
y*sin(a)+(1-cos(a))*x*z |
+        R  =  |  z*sin(a)+(1-cos(a))*x*y   1 + (1-cos(a))*(y*y-1)   
-x*sin(a)+(1-cos(a))*y*z |
+              | -y*sin(a)+(1-cos(a))*x*z   x*sin(a)+(1-cos(a))*y*z   1 + 
(1-cos(a))*(z*z-1)  |
+
+    @param R:           The 3x3 rotation matrix to update.
+    @type R:            3x3 numpy array
+    @param vector_orig: The unrotated vector defined in the reference frame.
+    @type vector_orig:  numpy array, len 3
+    @param vector_fin:  The rotated vector defined in the reference frame.
+    @type vector_fin:   numpy array, len 3
+    """
+
+    # Convert the vectors to unit vectors.
+    vector_orig = vector_orig / norm(vector_orig)
+    vector_fin = vector_fin / norm(vector_fin)
+
+    # The rotation axis.
+    axis = cross(vector_orig, vector_fin)
+    x = axis[0]
+    y = axis[1]
+    z = axis[2]
+
+    # The rotation angle.
+    angle = acos(dot(vector_orig, vector_fin))
+
+    # Trig functions (only need to do this maths once!).
+    ca = cos(angle)
+    sa = sin(angle)
+
+    # Calculate the rotation matrix elements.
+    R[0,0] = 1.0 + (1.0 - ca)*(x**2 - 1.0)
+    R[0,1] = -z*sa + (1.0 - ca)*x*y
+    R[0,2] = y*sa + (1.0 - ca)*x*z
+    R[1,0] = z*sa+(1.0 - ca)*x*y
+    R[1,1] = 1.0 + (1.0 - ca)*(y**2 - 1.0)
+    R[1,2] = -x*sa+(1.0 - ca)*y*z
+    R[2,0] = -y*sa+(1.0 - ca)*x*z
+    R[2,1] = x*sa+(1.0 - ca)*y*z
+    R[2,2] = 1.0 + (1.0 - ca)*(z**2 - 1.0)
 
 
 def R_euler_zyz(matrix, alpha, beta, gamma):