r26624 - /trunk/lib/geometry/vectors.py -- November 19, 2014 - 16:20

Author: bugman
Date: Wed Nov 19 16:20:53 2014
New Revision: 26624

URL: http://svn.gna.org/viewcvs/relax?rev=26624&view=rev
Log:
Created functions in the relax library for calculating the inter-vector angle 
for complex vectors.

This is in the lib.geometry.vectors module.  The function 
vector_angle_complex_conjugate() has been
created to calculate the angle between two complex vectors.  This uses the 
new auxiliary function
complex_inner_product() to calculate <v1|v2>.


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=26624&r1=26623&r2=26624&view=diff
==============================================================================
--- trunk/lib/geometry/vectors.py       (original)
+++ trunk/lib/geometry/vectors.py       Wed Nov 19 16:20:53 2014
@@ -24,9 +24,34 @@
 
 # Python module imports.
 from math import acos, atan2, cos, pi, sin
-from numpy import array, cross, dot, float64
+from numpy import array, cross, dot, float64, sqrt
 from numpy.linalg import norm
 from random import uniform
+
+
+def complex_inner_product(v1=None, v2_conj=None):
+    """Calculate the inner product <v1|v2> for the two complex vectors v1 
and v2.
+
+    This is calculated as::
+                   ___
+                   \
+        <v1|v2> =   >   v1i . v2i* ,
+                   /__
+                    i
+
+    where * is the complex conjugate.
+
+
+    @keyword v1:        The first vector.
+    @type v1:           numpy rank-1 complex array
+    @keyword v2_conj:   The conjugate of the second vector.  This is already 
in conjugate form to allow for non-standard definitions of the conjugate (for 
example Sm* = (-1)^m S-m).
+    @type v2_conj:      numpy rank-1 complex array
+    @return:            The value of the inner product <v1|v2>.
+    @rtype:             float
+    """
+
+    # Return the inner product.
+    return dot(v1, v2_conj)
 
 
 def random_unit_vector(vector):
@@ -114,6 +139,51 @@
 
     # Calculate and return the angle.
     return atan2(norm(cross(vector1, vector2)), dot(vector1, vector2))
+
+
+def vector_angle_complex_conjugate(v1=None, v2=None, v1_conj=None, 
v2_conj=None):
+    """Calculate the inter-vector angle between two complex vectors using 
the arccos formula.
+
+    The formula is::
+
+        theta = arccos(Re(<v1|v2>) / (|v1|.|v2|)) ,
+
+    where::
+                   ___
+                   \ 
+        <v1|v2> =   >   v1i . v2i* ,
+                   /__
+                    i
+
+    and::
+
+        |v1| = Re(<v1|v1>) .
+
+
+    @keyword v1:        The first vector.
+    @type v1:           numpy rank-1 complex array
+    @keyword v2:        The second vector.
+    @type v2:           numpy rank-1 complex array
+    @keyword v1_conj:   The conjugate of the first vector.  This is already 
in conjugate form to allow for non-standard definitions of the conjugate (for 
example Sm* = (-1)^m S-m).
+    @type v1_conj:      numpy rank-1 complex array
+    @keyword v2_conj:   The conjugate of the second vector.  This is already 
in conjugate form to allow for non-standard definitions of the conjugate (for 
example Sm* = (-1)^m S-m).
+    @type v2_conj:      numpy rank-1 complex array
+    @return:            The angle between 0 and pi.
+    @rtype:             float
+    """
+
+    # The inner products.
+    inner_v1v2 = complex_inner_product(v1=v1, v2_conj=v2_conj)
+    inner_v1v1 = complex_inner_product(v1=v1, v2_conj=v1_conj)
+    inner_v2v2 = complex_inner_product(v1=v2, v2_conj=v2_conj)
+
+    # The normalised inner product, correcting for truncation errors 
creating ratios slightly above 1.0.
+    ratio = inner_v1v2.real / (sqrt(inner_v1v1).real * sqrt(inner_v2v2).real)
+    if ratio > 1.0:
+        ratio = 1.0
+
+    # Calculate and return the angle.
+    return acos(ratio)
 
 
 def vector_angle_normal(vector1, vector2, normal):