Changeset 2704 for inundation/utilities/numerical_tools.py

Timestamp:

Apr 12, 2006, 5:52:34 AM (19 years ago)

Author:

ole

Message:

Generalised angle computation and tested it

File:

• inundation/utilities/numerical_tools.py (modified) (3 diffs)

inundation/utilities/numerical_tools.py

@@ -8 +8 @@
 #(this should move to somewhere central)
 try:
-    from scipy import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate    
+    from scipy import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float    
 except:
     #print 'Could not find scipy - using Numeric'
-    from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate    
-
-# Getting an infinate number to use when using Numeric
+    from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float    
+
+# Getting an infinite number to use when using Numeric
 INF = (array([1])/0.)[0]
 
-def angle(v):
-    """Compute angle between e1 (the unit vector in the x-direction)
-    and the specified vector
-    """
-
+
+
+def angle(v1, v2=None):
+    """Compute angle between 2D vectors v1 and v2.
+    
+    If v2 is not specified it will default 
+    to e1 (the unit vector in the x-direction)
+
+    The angle is measured as a number in [0, 2pi] from v2 to v1.
+    """
     from math import acos, pi, sqrt
-
-    l = sqrt( sum (array(v)**2))
-    v1 = v[0]/l
-    v2 = v[1]/l
-
-
-    theta = acos(v1)
-    
-    #try:
-    #   theta = acos(v1)
-    #except ValueError, e:
-    #    print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e)
-    #
-    #    #FIXME, hack to avoid acos(1.0) Value error
-    #    # why is it happening?
-    #    # is it catching something we should avoid?
-    #    #FIXME (Ole): When did this happen? We need a unit test to
-    #    #reveal this condition or otherwise remove the hack
-    #    
-    #    s = 1e-6
-    #    if (v1+s >  1.0) and (v1-s < 1.0) :
-    #        theta = 0.0
-    #    elif (v1+s >  -1.0) and (v1-s < -1.0):
-    #        theta = 3.1415926535897931
-    #    print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
-    #          %(v1, theta)
-
-    if v2 < 0:
-        #Quadrant 3 or 4
-        theta = 2*pi-theta
-
+  
+    # Prepare two Numeric vectors
+    if v2 is None:
+        v2 = [1.0, 0.0] # Unit vector along the x-axis
+        
+    v1 = ensure_numeric(v1, Float)
+    v2 = ensure_numeric(v2, Float)    
+    
+    # Normalise
+    v1 = v1/sqrt(sum(v1**2))
+    v2 = v2/sqrt(sum(v2**2))
+   
+    # Compute angle
+    p = innerproduct(v1, v2)
+    c = crossproduct_length(v1, v2)
+    
+    theta = acos(p)
+    
+    # Correct if v1 is in quadrant 3 or 4 with respect to v2 (as the x-axis) 
+    # If v2 was the unit vector [1,0] this would correspond to the test
+    # if v1[1] < 0: theta = 2*pi-theta    
+    # In general we use the sign of the cross product length
+    if c > 0: 
+       #Quadrant 3 or 4
+       theta = 2*pi-theta        
+       
     return theta
 
-
+    
 def anglediff(v0, v1):
     """Compute difference between angle of vector x0, y0 and x1, y1.
@@ -75 +75 @@
     return a0-a1
 
-
+def normal_vector(v):
+    """Normal vector to v
+    """
+    
+    return array([-v[1], v[0]], Float)
+
+    
+def crossproduct_length(v1, v2):
+    return v1[0]*v2[1]-v2[0]*v1[1]
+   
+       
 def mean(x):
     """Mean value of a vector
@@ -255 +265 @@
     pass
 
+
+    
+def angle_obsolete(v):
+    """Compute angle between e1 (the unit vector in the x-direction)
+    and the specified vector v.
+    
+    Return a number in [0, 2pi]
+    """
+    from math import acos, pi, sqrt
+  
+    # Normalise v
+    v = ensure_numeric(v, Float)
+    v = v/sqrt(sum(v**2))
+   
+    # Compute angle
+    theta = acos(v[0])
+     
+    if v[1] < 0: 
+       #Quadrant 3 or 4
+        theta = 2*pi-theta
+    
+    return theta
+