import math
import numpy

def normalise(a):
    'Normalises a vector'
    a_norm = numpy.linalg.norm(a)
    return numpy.array([float(each)/a_norm for each in a])

def rotation_axis_angle(axis, angle):
    '''
    Returns the matrix for rotation by an angle around an axis
    '''
    if abs(angle) < 1E-7:
        return numpy.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
    elif abs(math.pi - angle) < 1E-7:
        return numpy.array([[1, 0, 0], [0, -1, 0], [0, 0, -1]])
    else:
        sin = math.sin(angle)
        cos = math.cos(angle)
        comp = 1 - cos
        x, y, z = normalise(axis)
        mat = numpy.array([[(cos + x*x*comp), (x*y*comp - z*sin), (x*z*comp + y*sin)],
               [(y*x*comp + z*sin), (cos + y*y*comp), (y*z*comp - x*sin)],
               [(z*x*comp - y*sin), (z*y*comp + x*sin), (cos + z*z*comp)]])
        should_be_I = mat.dot(mat.transpose())
        I = numpy.ma.identity(3)
        # There should be a way of doing this in numpy
        numpy.testing.assert_array_almost_equal(I, should_be_I, 3)
        return mat


def rotation_from_axes(ax1, ax2):  # To test
    '''
    Given two 3-vectors as tuples, gives the matrix that will rotate the first
    to the second
    '''
    # Probably a more numpy-ish way of doing this
    if max(numpy.absolute(ax1 - ax2)) < 1E-7:
        return numpy.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
    elif max(numpy.absolute(ax1 + ax2)) < 1E-7:
        return numpy.array([[1, 0, 0], [0, -1, 0], [0, 0, -1]])
    else:
        ang = angle_between(ax1, ax2)
        rot_ax = numpy.cross(ax1, ax2)
        return rotation_axis_angle(rot_ax, ang)

def angle_between(vec1, vec2):
    return math.acos(float(vec1.dot(vec2))/
                     (numpy.linalg.norm(vec1)*numpy.linalg.norm(vec2)))

def reflection_plane(vec1, vec2):
    '''
    Returns the Householder reflection matrix for reflection through 
    the plane containning two vectors
    '''
    norm = numpy.cross(vec1, vec2)
    a, b, c = normalise(norm)
    return numpy.array([[1 - 2*a*a, -2*a*b, -2*a*c],
                        [-2*a*b, 1-2*b*b, -2*b*c],
                        [-2*a*c, -2*b*c, 1-2*c*c]])
@@ -278,212 +278,8 @@
 '''%0A
-    if abs(angle) %3C 1E-7:%0A        return numpy.array(%5B%5B1, 0, 0%5D, %5B0, 1, 0%5D, %5B0, 0, 1%5D%5D)%0A    elif abs(math.pi - angle) %3C 1E-7:%0A        return numpy.array(%5B%5B1, 0, 0%5D, %5B0, -1, 0%5D, %5B0, 0, -1%5D%5D)%0A    else:%0A    
     
@@ -296,28 +296,24 @@
 .sin(angle)%0A
-    
     cos = ma
@@ -326,20 +326,16 @@
 (angle)%0A
-    
     comp
@@ -345,20 +345,16 @@
 1 - cos%0A
-    
     x, y
@@ -375,20 +375,16 @@
 e(axis)%0A
-    
     mat 
@@ -473,16 +473,24 @@
         
+        
 %5B(y*x*co
@@ -556,16 +556,24 @@
         
+        
 %5B(z*x*co
@@ -625,20 +625,16 @@
 omp)%5D%5D)%0A
-    
     shou
@@ -672,20 +672,16 @@
 ())%0A    
-    
 I = nump
@@ -701,67 +701,8 @@
 (3)%0A
-        # There should be a way of doing this in numpy%0A    
     
@@ -760,20 +760,16 @@
 e_I, 3)%0A
-    
     retu
@@ -1141,38 +1141,160 @@
 %3C 1E-7:%0A        
-return
+ang = angle_between(ax1, ax2)%0A        z = math.sqrt(1/(1 + (ax1%5B2%5D/ax1%5B1%5D)**2))%0A        y = math.sqrt(1 - z**2)%0A        rot_ax =
  numpy.array(%5B%5B1
@@ -1295,42 +1295,64 @@
 ay(%5B
-%5B1, 0, 0%5D, %5B0, -1, 0%5D, %5B0, 0, -1%5D%5D
+0, y, z%5D)%0A        return rotation_axis_angle(rot_ax, ang
 )%0A  

Fix edge-case handling in rotation matrix generation