r22520 - /trunk/lib/geometry/vectors.py -- March 21, 2014 - 17:29

Author: bugman
Date: Fri Mar 21 17:29:24 2014
New Revision: 22520

URL: http://svn.gna.org/viewcvs/relax?rev=22520&view=rev
Log:
Changes to the lib.geometry.vectors.vector_angle() function.

This now expects the normal of the plane in which the angle is defined.  The 
original logic was not
functional, therefore the angle is forced to be negative if the cross product 
between the two
vectors points in the opposite direction as the normal.


Modified:
    trunk/lib/geometry/vectors.py

Modified: trunk/lib/geometry/vectors.py
URL: 
http://svn.gna.org/viewcvs/relax/trunk/lib/geometry/vectors.py?rev=22520&r1=22519&r2=22520&view=diff
==============================================================================
--- trunk/lib/geometry/vectors.py       (original)
+++ trunk/lib/geometry/vectors.py       Fri Mar 21 17:29:24 2014
@@ -74,23 +74,30 @@
     return vect / norm(vect)
 
 
-def vector_angle(vector1, vector2):
+def vector_angle(vector1, vector2, normal):
     """Calculate the directional angle between two N-dimensional vectors.
-
-    The angle between vectors A and B is calculated using the formula::
-
-        theta = arctan(AxB / A.B),
-
-    where the arctan function used is atan2, AxB is the cross product 
between the two vectors, and A.B is the dot product.
-
 
     @param vector1:     The first vector.
     @type vector1:      numpy rank-1 array
     @param vector2:     The second vector.
     @type vector2:      numpy rank-1 array
+    @param normal:      The vector defining the plane, to determine the sign.
+    @type normal:       numpy rank-1 array
     @return:            The angle between -pi and pi.
     @rtype:             float
     """
 
-    # Calculate and return the value.
-    return atan2(cross(vector1, vector2), dot(vector1, vector2))
+    # Normalise the vectors (without changing the original vectors).
+    vector1 = vector1 / norm(vector1)
+    vector2 = vector2 / norm(vector2)
+
+    # The cross product.
+    cp = cross(vector1, vector2)
+
+    # The angle.
+    angle = acos(dot(vector1, vector2))
+    if dot(cp, normal) < 0.0:
+        angle = -angle
+
+    # Return the signed angle.
+    return angle