Changeset 2704

Timestamp:

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

Author:

ole

Message:

Generalised angle computation and tested it

Location:

inundation

Files:

• pyvolution/util.py (modified) (2 diffs)
• utilities/numerical_tools.py (modified) (3 diffs)
• utilities/test_numerical_tools.py (modified) (1 diff)

inundation/pyvolution/util.py

@@ -399 +399 @@
 
 
-def angle(v):
+def angle(v1, v2):
     """Temporary Interface to new location"""
 
@@ -408 +408 @@
     warn(msg, DeprecationWarning) 
 
-    return NT.angle(v)
+    return NT.angle(v1, v2)
     
 def anglediff(v0, v1):

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
+    

inundation/utilities/test_numerical_tools.py

@@ -20 +20 @@
 
 
-    def test_angle(self):
+    def test_angle1(self):
+        """Test angles between one vector and the x-axis
+        """
+        assert allclose(angle([1.0, 0.0])/pi*180, 0.0)      
         assert allclose(angle([1.0, 1.0])/pi*180, 45.0)
+        assert allclose(angle([0.0, 1.0])/pi*180, 90.0)         
+        assert allclose(angle([-1.0, 1.0])/pi*180, 135.0)               
+        assert allclose(angle([-1.0, 0.0])/pi*180, 180.0)
+        assert allclose(angle([-1.0, -1.0])/pi*180, 225.0)
+        assert allclose(angle([0.0, -1.0])/pi*180, 270.0)
+        assert allclose(angle([1.0, -1.0])/pi*180, 315.0)               
+                
+                                                          
+    def test_angle2(self):
+        """Test angles between two arbitrary vectors
+        """    
+        
+        assert allclose(angle([1.0, 0.0], [1.0, 1.0])/pi*180, 315.0)
+        assert allclose(angle([1.0, 1.0], [1.0, 0.0])/pi*180, 45.0)
+                
+        assert allclose(angle([-1.0, -1.0], [1.0, 1.0])/pi*180, 180)    
+        assert allclose(angle([-1.0, -1.0], [-1.0, 1.0])/pi*180, 90.0)  
+        
+        assert allclose(angle([-1.0, 0.0], [1.0, 1.0])/pi*180, 135.0)
+        assert allclose(angle([0.0, -1.0], [1.0, 1.0])/pi*180, 225.0)   
+        
+        assert allclose(angle([1.0, -1.0], [1.0, 1.0])/pi*180, 270.0)   
+        assert allclose(angle([1.0, 0.0], [0.0, 1.0])/pi*180, 270.0)    
+                
+        
+                                
+                
+        
+