Changeset 2704


Ignore:
Timestamp:
Apr 12, 2006, 5:52:34 AM (18 years ago)
Author:
ole
Message:

Generalised angle computation and tested it

Location:
inundation
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/util.py

    r2701 r2704  
    399399
    400400
    401 def angle(v):
     401def angle(v1, v2):
    402402    """Temporary Interface to new location"""
    403403
     
    408408    warn(msg, DeprecationWarning)
    409409
    410     return NT.angle(v)
     410    return NT.angle(v1, v2)
    411411   
    412412def anglediff(v0, v1):
  • inundation/utilities/numerical_tools.py

    r2633 r2704  
    88#(this should move to somewhere central)
    99try:
    10     from scipy import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate   
     10    from scipy import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float   
    1111except:
    1212    #print 'Could not find scipy - using Numeric'
    13     from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate   
    14 
    15 # Getting an infinate number to use when using Numeric
     13    from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float   
     14
     15# Getting an infinite number to use when using Numeric
    1616INF = (array([1])/0.)[0]
    1717
    18 def angle(v):
    19     """Compute angle between e1 (the unit vector in the x-direction)
    20     and the specified vector
    21     """
    22 
     18
     19
     20def angle(v1, v2=None):
     21    """Compute angle between 2D vectors v1 and v2.
     22   
     23    If v2 is not specified it will default
     24    to e1 (the unit vector in the x-direction)
     25
     26    The angle is measured as a number in [0, 2pi] from v2 to v1.
     27    """
    2328    from math import acos, pi, sqrt
    24 
    25     l = sqrt( sum (array(v)**2))
    26     v1 = v[0]/l
    27     v2 = v[1]/l
    28 
    29 
    30     theta = acos(v1)
    31    
    32     #try:
    33     #   theta = acos(v1)
    34     #except ValueError, e:
    35     #    print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e)
    36     #
    37     #    #FIXME, hack to avoid acos(1.0) Value error
    38     #    # why is it happening?
    39     #    # is it catching something we should avoid?
    40     #    #FIXME (Ole): When did this happen? We need a unit test to
    41     #    #reveal this condition or otherwise remove the hack
    42     #   
    43     #    s = 1e-6
    44     #    if (v1+s >  1.0) and (v1-s < 1.0) :
    45     #        theta = 0.0
    46     #    elif (v1+s >  -1.0) and (v1-s < -1.0):
    47     #        theta = 3.1415926535897931
    48     #    print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
    49     #          %(v1, theta)
    50 
    51     if v2 < 0:
    52         #Quadrant 3 or 4
    53         theta = 2*pi-theta
    54 
     29 
     30    # Prepare two Numeric vectors
     31    if v2 is None:
     32        v2 = [1.0, 0.0] # Unit vector along the x-axis
     33       
     34    v1 = ensure_numeric(v1, Float)
     35    v2 = ensure_numeric(v2, Float)   
     36   
     37    # Normalise
     38    v1 = v1/sqrt(sum(v1**2))
     39    v2 = v2/sqrt(sum(v2**2))
     40   
     41    # Compute angle
     42    p = innerproduct(v1, v2)
     43    c = crossproduct_length(v1, v2)
     44   
     45    theta = acos(p)
     46   
     47    # Correct if v1 is in quadrant 3 or 4 with respect to v2 (as the x-axis)
     48    # If v2 was the unit vector [1,0] this would correspond to the test
     49    # if v1[1] < 0: theta = 2*pi-theta   
     50    # In general we use the sign of the cross product length
     51    if c > 0:
     52       #Quadrant 3 or 4
     53       theta = 2*pi-theta       
     54       
    5555    return theta
    5656
    57 
     57   
    5858def anglediff(v0, v1):
    5959    """Compute difference between angle of vector x0, y0 and x1, y1.
     
    7575    return a0-a1
    7676
    77 
     77def normal_vector(v):
     78    """Normal vector to v
     79    """
     80   
     81    return array([-v[1], v[0]], Float)
     82
     83   
     84def crossproduct_length(v1, v2):
     85    return v1[0]*v2[1]-v2[0]*v1[1]
     86   
     87       
    7888def mean(x):
    7989    """Mean value of a vector
     
    255265    pass
    256266
     267
     268   
     269def angle_obsolete(v):
     270    """Compute angle between e1 (the unit vector in the x-direction)
     271    and the specified vector v.
     272   
     273    Return a number in [0, 2pi]
     274    """
     275    from math import acos, pi, sqrt
     276 
     277    # Normalise v
     278    v = ensure_numeric(v, Float)
     279    v = v/sqrt(sum(v**2))
     280   
     281    # Compute angle
     282    theta = acos(v[0])
     283     
     284    if v[1] < 0:
     285       #Quadrant 3 or 4
     286        theta = 2*pi-theta
     287   
     288    return theta
     289   
  • inundation/utilities/test_numerical_tools.py

    r2534 r2704  
    2020
    2121
    22     def test_angle(self):
     22    def test_angle1(self):
     23        """Test angles between one vector and the x-axis
     24        """
     25        assert allclose(angle([1.0, 0.0])/pi*180, 0.0)     
    2326        assert allclose(angle([1.0, 1.0])/pi*180, 45.0)
     27        assert allclose(angle([0.0, 1.0])/pi*180, 90.0)         
     28        assert allclose(angle([-1.0, 1.0])/pi*180, 135.0)               
     29        assert allclose(angle([-1.0, 0.0])/pi*180, 180.0)
     30        assert allclose(angle([-1.0, -1.0])/pi*180, 225.0)
     31        assert allclose(angle([0.0, -1.0])/pi*180, 270.0)
     32        assert allclose(angle([1.0, -1.0])/pi*180, 315.0)               
     33               
     34                                                         
     35    def test_angle2(self):
     36        """Test angles between two arbitrary vectors
     37        """   
     38       
     39        assert allclose(angle([1.0, 0.0], [1.0, 1.0])/pi*180, 315.0)
     40        assert allclose(angle([1.0, 1.0], [1.0, 0.0])/pi*180, 45.0)
     41               
     42        assert allclose(angle([-1.0, -1.0], [1.0, 1.0])/pi*180, 180)   
     43        assert allclose(angle([-1.0, -1.0], [-1.0, 1.0])/pi*180, 90.0) 
     44       
     45        assert allclose(angle([-1.0, 0.0], [1.0, 1.0])/pi*180, 135.0)
     46        assert allclose(angle([0.0, -1.0], [1.0, 1.0])/pi*180, 225.0)   
     47       
     48        assert allclose(angle([1.0, -1.0], [1.0, 1.0])/pi*180, 270.0)   
     49        assert allclose(angle([1.0, 0.0], [0.0, 1.0])/pi*180, 270.0)   
     50               
     51       
     52                               
     53               
     54       
     55               
    2456
    2557
Note: See TracChangeset for help on using the changeset viewer.